Multiphase processes in the EC-Earth model and their relevance to the atmospheric oxalate, sulfate, and iron cycles

Understanding how multiphase processes affect the iron-containing aerosol cycle is key to predicting ocean biogeochemistry changes and hence the feedback effects on climate. For this work, the EC-Earth Earth system model in its climate–chemistry configuration is used to simulate the global atmospheric oxalate (OXL), sulfate (SO${}_{\mathrm{4}}^{\mathrm{2}-}$), and iron (Fe) cycles after incorporating a comprehensive representation of the multiphase chemistry in cloud droplets and aerosol water. The model considers a detailed gas-phase chemistry scheme, all major aerosol components, and the partitioning of gases in aerosol and atmospheric water phases. The dissolution of Fe-containing aerosols accounts kinetically for the solution's acidity, oxalic acid, and irradiation. Aerosol acidity is explicitly calculated in the model, both for accumulation and coarse modes, accounting for thermodynamic processes involving inorganic and crustal species from sea salt and dust.

Simulations for present-day conditions (2000–2014) have been carried out with both EC-Earth and the atmospheric composition component of the model in standalone mode driven by meteorological fields from ECMWF's ERA-Interim reanalysis. The calculated global budgets are presented and the links between the (1) aqueous-phase processes, (2) aerosol dissolution, and (3) atmospheric composition are demonstrated and quantified. The model results are supported by comparison to available observations. We obtain an average global OXL net chemical production of 12.615 ± 0.064 Tg yr⁻¹ in EC-Earth, with glyoxal being by far the most important precursor of oxalic acid. In comparison to the ERA-Interim simulation, differences in atmospheric dynamics and the simulated weaker oxidizing capacity in EC-Earth overall result in a ∼ 30 % lower OXL source. On the other hand, the more explicit representation of the aqueous-phase chemistry in EC-Earth compared to the previous versions of the model leads to an overall ∼ 20 % higher sulfate production, but this is still well correlated with atmospheric observations.

The total Fe dissolution rate in EC-Earth is calculated at 0.806 ± 0.014 Tg yr⁻¹ and is added to the primary dissolved Fe (DFe) sources from dust and combustion aerosols in the model (0.072 ± 0.001 Tg yr⁻¹). The simulated DFe concentrations show a satisfactory comparison with available observations, indicating an atmospheric burden of ∼0.007 Tg, resulting in an overall atmospheric deposition flux into the global ocean of 0.376 ± 0.005 Tg yr⁻¹, which is well within the range reported in the literature. All in all, this work is a first step towards the development of EC-Earth into an Earth system model with fully interactive bioavailable atmospheric Fe inputs to the marine biogeochemistry component of the model.

1 Introduction

Clouds, fog, and deliquescent aerosols host chemical reactions involving inorganic and organic polar atmospheric compounds (Calvert et al., 1985; Chameides and Davis, 1983; Collett et al., 1999; Donaldson and Valsaraj, 2010; Jacob, 1986; Lelieveld and Crutzen, 1991). These reactions result in the production of species that can neither be formed via gas-phase processes directly, nor explained solely by primary sources. These compounds participate in chemical transformations across the gas, aqueous, and solid phases. Such multiphase processes have a significant impact on the atmospheric cycles of important inorganic species like sulfur (e.g., Hoyle et al., 2016; Seinfeld and Pandis, 2006; Tsai et al., 2010) and act as a complementary pathway for the formation of organic particulate matter (e.g., Lin et al., 2014; Liu et al., 2012; Myriokefalitakis et al., 2011). The produced inorganic and organic aerosols serve as cloud condensation nuclei and thus affect the Earth's energy balance (IPCC, 2013).

Multiphase processes may also impact the global carbon balance indirectly by altering the atmospheric cycles of species that act as nutrients for the marine biota (Hamilton et al., 2022; Kanakidou et al., 2018; Mahowald et al., 2017; Myriokefalitakis et al., 2020a). Nutrient availability in marine ecosystems is key for the primary production that modulates both the surface oceanic concentrations and the uptake of atmospheric CO₂ (e.g., Le Quéré et al., 2007, 2013; Gruber et al., 2019). A large portion of the global ocean is found, however, to be limited in iron (Krishnamurthy et al., 2009, 2010); therefore, the importance of iron (Fe) to oceanic productivity is well established (Hamilton et al., 2020; Kanakidou et al., 2020; Meskhidze et al., 2019; Tagliabue et al., 2016). Besides rivers and sea ice, in addition to sediment dissolution and hydrothermal vents, which are the main sources of bioavailable Fe in the ocean, the atmospheric deposition of nutrients is the most effective external pathway that provides Fe in the open ocean. Fe is a critical micronutrient for marine biota that is mainly utilized in its dissolved form (e.g., aqueous, colloidal, or nanoparticulate). Thus, the atmospheric processing of Fe-containing minerals, i.e., the conversion from insoluble to soluble that is readily available Fe for marine organisms, is a central step in the atmospheric and marine Fe cycles and directly connected to atmospheric multiphase processes.

Fe is mainly present in the atmosphere in crystalline lattices of aluminosilicates or as iron oxides in dust aerosols (∼ 95 %; Mahowald et al., 2009) and tends to be rather insoluble when emitted (up to ∼ 1 % solubility; Journet et al., 2008). In fact, observed high Fe solubility downwind of dust source regions can be only explained via the atmospheric processing of dust aerosols (Baker and Jickells, 2017; Oakes et al., 2012). Enhanced Fe solubility is observed for biomass burning aerosols (e.g., ranging 2 %–46 %; Bowie et al., 2009; Guieu et al., 2005​​​​​​​; Mahowald et al., 2018; Oakes et al., 2012; Paris et al., 2010), depending strongly on the source region and/or the type of burned wood. Significantly higher Fe solubilities are found, however, for anthropogenic combustion-related Fe-containing aerosols, especially for Fe in oil fly ash from industries and shipping, which is mainly in the form of ferric sulfates (Chen et al., 2012; Ito, 2013; Rathod et al., 2020; Schroth et al., 2009). The uncertainty in Fe-containing combustion aerosol solubility (e.g., Rathod et al., 2020) is nevertheless also reflected in modeling studies, with some models assuming relatively high solubility at emission (e.g., Hamilton et al., 2019; Myriokefalitakis et al., 2011) depending on the aerosol size, and others assuming an almost completely insoluble emitted Fe whose solubility is then enhanced during transport via atmospheric processing (Ito, 2015; Ito et al., 2021). Recent multimodel studies estimate an overall global dissolved Fe (DFe) production rate due to atmospheric processing of dust and combustion aerosols of 0.56 ± 0.29 Tg yr⁻¹ (Ito et al., 2019; Myriokefalitakis et al., 2018), indicating that a large uncertainty still remains in the impact of atmospheric processing on the mineral Fe solubilization processes.

During atmospheric transport, inorganic strong acids and organic ligands may coat mineral aerosols and eventually convert part of the contained insoluble Fe forms (e.g., hematite) to bioavailable forms of Fe for marine biota in the euphotic zone (e.g., free ferrous forms, inorganic soluble Fe, and organic Fe complexes). Mineral dissolution rates depend on the solution's acidity levels, the mineral surface concentration of organic ligands, sunlight, and ambient temperature (e.g., Hamer et al., 2003; Lanzl et al., 2012; Lasaga et al., 1994; Zhu et al., 1993). Although sulfate (SO${}_{\mathrm{4}}^{\mathrm{2}-}$) is the dominant aerosol species that controls the aerosol liquid water content and acidity, oxalate ((COO⁻)₂; hereafter OXL) acts as an organic ligand for the Fe-containing aerosol dissolution processes (e.g., Paris et al., 2011; Paris and Desboeufs, 2013) that can effectively break the Fe–O bonds at the mineral's surface via the formation of ligand-containing surface structures (Yoon et al., 2004). Despite the dominant role of acidity in the mineral Fe dissolution processes, modeling estimates (Ito, 2015; Johnson and Meskhidze, 2013; Myriokefalitakis et al., 2015) show the importance of OXL to atmospheric DFe concentrations (e.g., including the formation of Fe(II/III) oxalate complexes). The dissolution of Fe by OXL may further contribute to the organic-bounded pool of nutrients deposited into the ocean, and it thus affects the marine primary production, especially in oligotrophic subtropical gyres (e.g., up to 20 %; Myriokefalitakis et al., 2020a).

Notwithstanding their different roles and efficiencies in Fe solubilization processes, atmospheric observations demonstrate a strong correlation between SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and OXL concentrations (Yu et al., 2005), especially above clouds (Sorooshian et al., 2006), indicating common chemical production pathways despite the differences in their precursors and primary sources. SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and OXL are the most common species formed via aqueous-phase reactions of inorganic and organic origin, respectively, with modeling studies supporting the conclusion that more than 60 % of the sulfates (e.g., Liao et al., 2003) and about 90 % of oxalates (Lin et al., 2012; Liu et al., 2012; Myriokefalitakis et al., 2011) are produced in clouds. OXL is the dominant dicarboxylic acid (DCA) in the troposphere (e.g., Kawamura and Ikushima, 1993; Kawamura and Sakaguchi, 1999; Norton et al., 1983) and is formed primarily through cloud processing of glyoxal and other water-soluble products of alkenes and aromatics of anthropogenic, biogenic, and marine origin (Carlton et al., 2007; Warneck, 2003). OXL is mostly present in the troposphere in particulate form (Yang and Yu, 2008), with aerosol concentrations roughly 4 times larger than in the gas phase (Martinelango et al., 2007; Yao et al., 2002). OXL can be present in urban environments (Yang et al., 2009) and in remote regions (Sempére and Kawamura, 1994) and is produced during the photochemical aging of organic aerosols (Eliason et al., 2003). The observed correlation of OXL with ammonium (NH${}_{\mathrm{4}}^{+}$) (Martinelango et al., 2007) indicates that OXL is mostly present as a salt (i.e., ammonium oxalate; (NH₄)₂C₂O₄) in the atmosphere (Paciga et al., 2014). Ortiz-Montalvo et al. (2014) found that in the presence of NH${}_{\mathrm{4}}^{+}$ under cloud-relevant conditions, the OXL produced by the aqueous-phase glyoxal oxidation is efficiently converted to ammonium oxalate, with its vapor pressure being several orders of magnitude lower than that of oxalic acid. However, in the presence of metals, such as calcium (Ca²⁺) and magnesium (Mg²⁺) from dust and sea salt aerosols, most of the oxalic acid is found to be present in the form of metal complexes (Furukawa and Takahashi, 2011). Nevertheless, due to their different solubility, the stability of oxalate complexes can be rather diverse, while calcium and magnesium oxalates precipitate from the solution, other salts, such as sodium or ammonium oxalates, remain in a deliquescent form (Furukawa and Takahashi, 2011).

Laboratory and modeling studies support the conclusion that OXL is directly produced in atmospheric water via glyoxylic acid (GLX; HC(O)COOH) oxidation by hydroxyl (OH) and nitrate (NO₃) radicals. The estimated net global OXL production rate in atmospheric water ranges between 13 and 30 Tg yr⁻¹ (Lin et al., 2014; Liu et al., 2012; Myriokefalitakis et al., 2011). However, modeling studies where the OXL production is only based on the GLX aqueous-phase oxidation tend to underestimate its observed atmospheric concentrations (e.g., Lin et al., 2014; Myriokefalitakis et al., 2011). Based on laboratory experiments, Carlton et al. (2007) proposed that predictions of oxalic acid concentrations could be significantly improved when larger multifunctional compounds are allowed to be produced under elevated glyoxal concentrations in typical cloud conditions. These larger multifunctional products can act as precursors for the glyoxylic and oxalic acids via their rapid oxidation by OH radicals Carlton et al., 2007). When such reactions are included, models tend to predict a higher oxalate atmospheric load and thus better match the observations (e.g., Myriokefalitakis et al., 2011). Note that although small carbonyl compounds, such as glyoxal and methylglyoxal, can undergo oligomerization under concentrated acidic conditions (Ervens and Volkamer, 2010; Lim et al., 2010, 2013), the mechanism behind the production of larger multifunctional products in dilute solutions may be rather complex, e.g., for products with alcohol functional groups, covalently bonded oligomers, larger carboxylic acids, and other humic-like substance (HULIS) components (Altieri et al., 2006; Blando and Turpin, 2000; Cappiello et al., 2003; Carlton et al., 2007).

The involvement of Fe chemistry in the aqueous phase decreases the global OXL net production rates overall (by ∼ 57 %), despite the increase in dissolved OH radical sources and thus the oxidation of OXL precursors (Lin et al., 2014). Besides the dissolved H₂O₂ photolysis that drastically enhances the OH production in the solution during the daytime, the presence of transition metal ions (TMIs) may play a central role in aqueous-phase oxidizing capacity, especially under dark conditions (Tilgner et al., 2013; Tilgner and Herrmann, 2018). Among other metals, Fe is the most efficient for the aqueous-phase oxidizing capacity, since on one hand it contributes to the OH reactivity via the Fenton reaction and the direct Fe photolysis, and on the other hand its dissolved concentrations are high due to the mineral dust contribution. The metal oxalate complexes formed in the presence of Fe in the solution (Zuo and Deng, 1997), however, can also undergo Fenton reaction and further increase the dissolved OH source, particularly for air masses of continental origin (Bianco et al., 2020) where elevated concentrations of OXL precursors and Fe-containing aerosols from both lithogenic and pyrogenic sources can exist. The photolysis of Fe oxalate complex [Fe(C₂O₄)₂]⁻ eventually transforms C₂O${}_{\mathrm{4}}^{\mathrm{2}-}$ into CO₂ in the aqueous phase (Ervens et al., 2003). Overall, it is clear that the impact of the Fe redox chemistry on the OXL production (and vice versa) is a rather complex issue that can also affect the ligand-promoted dissolution process of the Fe-containing minerals under ambient atmospheric conditions.

For this work, we incorporate a comprehensive aqueous-phase chemistry scheme into a state-of-the-art global climate–chemistry model to simulate the atmospheric multiphase processes with respect to iron-containing aerosol dissolution. Section 2 provides an overview of the model, focusing mostly on the new implementations. In particular, we describe the multiphase chemistry scheme used to simulate the atmospheric OXL, SO${}_{\mathrm{4}}^{\mathrm{2}-}$, and Fe cycles, along with the respective developments for the primary soil and combustion sources applied in the model. In Sect. 3, we present the model-derived OXL-, SO${}_{\mathrm{4}}^{\mathrm{2}-}$-, and Fe-containing aerosol atmospheric concentrations and their evaluation with available observations, and in Sect. 4 we discuss the impact of the simulated aqueous-phase processes on the DFe deposition fluxes to the global ocean. Finally, in Sect. 5, we summarize the global implications of explicitly resolving multiphase chemistry in a climate–chemistry model for the atmospheric Fe cycle, along with the plans for future model development.

2 Model description

2.1 The EC-Earth3 Earth system Model

Our tropospheric multiphase chemistry developments have been implemented in the global Earth system model (ESM) EC-Earth3 (Döscher et al., 2021). EC-Earth3 took part in the Coupled Model Intercomparison Project phase 6 (CMIP6; Eyring et al., 2016). The atmospheric general circulation model (GCM) of EC-Earth3 is based on cycle 36r4 of the Integrated Forecast System (IFS) from the European Centre for Medium-Range Weather Forecasts (ECMWF), which includes the land surface model H-TESSEL (Balsamo et al., 2009). The ocean model is the Nucleus for European Modeling of the Ocean (NEMO) release 3.6 (Rousset et al., 2015), with sea ice processes represented by the Louvain-la-Neuve sea ice model (LIM) (Rousset et al., 2015; Vancoppenolle et al., 2009). The ESM presents the following two configurations: (1) the carbon cycle configuration that represents the marine biogeochemistry processes through PISCES (Aumont et al., 2015), the dynamic terrestrial vegetation through LPJ-Guess (Smith et al., 2001, 2014), and the atmospheric cycle of CO₂ through the Tracer Model version 5 release 3.0 (TM5-MP 3.0) and (2) the EC-Earth3-AerChem configuration (van Noije et al., 2021) that represents the atmospheric chemistry and transport of aerosols and reactive species (also through the TM5-MP 3.0). Most of the information exchange and interpolation between modules is handled through the Ocean Atmosphere Sea Ice Soil version 3 (OASIS3) coupler (Craig et al., 2017). For this work we rely on the EC-Earth3-AerChem branch specifically (van Noije et al., 2021).

EC-Earth3-AerChem includes TM5-MP to simulate tropospheric aerosols and the reactive greenhouse gases methane (CH₄) and ozone (O₃) and allows the coupling of those species to relevant processes in the atmospheric module IFS (e.g., radiation and clouds). The model can be executed in an atmospheric mode only, i.e., using prescribed sea surface temperature and sea ice concentration, or coupled to the NEMO-LIM ocean and sea ice model. In addition, TM5-MP can run as a standalone (offline) atmospheric chemistry and transport model (CTM) driven by meteorological and surface fields (Krol et al., 2005). The present work is structured around a recently released version of TM5-MP that incorporates a rather detailed gas-phase tropospheric chemistry scheme, the MOGUNTIA (Myriokefalitakis et al., 2020b). MOGUNTIA explicitly simulates the organic polar species that partition in the atmospheric aqueous phase and allows for a sophisticated parameterization of the multiphase processes needed for this study.

All major aerosol components such as sulfate, black carbon, organic aerosols, sea salt, and mineral dust aerosols are included in TM5-MP and are distributed (depending on the aerosol type) in seven lognormal modes, i.e., four soluble modes (i.e., nucleation, Aitken, accumulation, and coarse) and three insoluble modes (i.e., Aitken, accumulation, and coarse). The aerosol microphysics in the model is calculated by the modal aerosol scheme M7 (Aan de Brugh et al., 2011; Vignati et al., 2004), which represents both the evolution of the total particle number and mass of the different species in each mode. Ammonium, nitrate, and aerosol water are determined based on gas–particle partitioning. M7 uses seven lognormal size distributions with predefined geometric standard deviations, with four water-soluble modes (nucleation, Aitken, accumulation, and coarse) and three insoluble modes (Aitken, accumulation, and coarse). Note that the new developments of this work are added to the model on top of the aerosols already represented by M7 and that the new aerosol components are introduced using the existing modes. Primary emissions of anthropogenic, biogenic, and biomass burning processes are defined through a variety of datasets; the most updated being those produced for the CMIP6 project. Natural emissions of mineral dust, sea salt, marine dimethyl sulfide (DMS), and nitrogen oxides from lighting are calculated online, while other natural emissions are prescribed. Details on the various parameterizations used for the definition of the gas and aerosol emissions in the model can be found in van Noije et al. (2021).

2.2 The EC-Earth3-Iron model

EC-Earth3-Iron is the new version of the model developed and used for this work that builds on EC-Earth3-AerChem. The new features required to determine the global aqueous-phase OXL formation, the atmospheric acidity, and the Fe cycle in the atmosphere can be summarized as follows:

1. treatment of mineral dust emission that considers soil mineralogical composition variations to account for the emission of Fe-containing minerals (and calcite), along with a detailed speciation of anthropogenic combustion and biomass burning emissions to explicitly account for Fe both in soluble and insoluble forms;

2. acidity calculations for water contained in fine and coarse aerosols, as well as for cloud droplets;

3. a comprehensive aqueous phase chemistry scheme in cloud droplets and aerosol water;

4. an explicit description of the Fe-containing aerosol dissolution processes of mineral dust, anthropogenic combustion, and biomass burning aerosols.

2.2.1 Speciated emissions

EC-Earth3-Iron includes a characterization of the dust mineralogical composition at emission and explicitly traces the Fe and calcium-containing species. The relative amounts of eight different minerals, namely illite, kaolinite, montmorillonite, calcite, feldspars, quartz, gypsum, and hematite, are derived from the soil mineralogy atlas of Claquin et al. (1999), including the updates proposed in Nickovic et al. (2012). The atlas provides the soil mineralogical composition in arid and semi-arid regions of the world, distinguishing between two soil size classes (i.e., the clay size fraction up to 2 µm and the silt size fraction from 2 to 50 µm diameter). The mineral fractions emitted in the accumulation and coarse insoluble modes of TM5-MP are estimated from the soil mineralogy atlas based on the brittle fragmentation theory (BFT) from Kok (2011). BFT posits that the emitted particle size distribution is independent of wind and soil conditions and additionally allows for estimating the size-resolved mineral fractions (Pérez García-Pando et al., 2016; Perlwitz et al., 2015a, b). The resulting mineral mass fractions are then applied to the dust emission fluxes, as calculated online in the model, yielding the corresponding accumulation- and coarse-mode emission of each mineral. We note that although we derive the mineral dust fractions in each mode using BFT, we maintain the dependence of the ratio between the accumulation- and coarse-mode dust mass at emissions upon wind and soil conditions of the original dust emission scheme (Tegen et al., 2002).

In EC-Earth3-Iron, the different Fe-containing minerals are not prognostic variables (tracers). Instead, we trace the mineral dust Fe according to three dissolution classes, namely fast, intermediate, and slow Fe pools (Ito and Shi, 2016). No relationship of Fe dissolution with other elements is observed, however, for clays and feldspars, where the total Fe content of the minerals is very low (<0.54 %), and the Fe is in the form of impurities (Journet et al., 2008). For this, 0.1 % Fe content in total Fe-containing minerals is here assumed directly soluble as amorphous free iron impurities regardless of mineralogy (Ito and Shi, 2016). The emitted amounts of calcium (i.e., in calcite) and Fe (i.e., in illite, kaolinite, montmorillonite, feldspars, and hematite) are derived either from the average elemental compositions of minerals or based on experimental analyses (Journet et al., 2008; Nickovic et al., 2013). The respective average fractions applied to mineral dust sources of this work are listed in Table S1.

Fe is also emitted in the model from anthropogenic activities (including fossil and biomass fuels) and biomass burning (excluding biofuel combustion) following Ito et al. (2018). The Fe-containing fossil fuel and biofuel combustion emissions are estimated here by applying specific factors (i.e., per emission sector and per particle size) to the total particulate emissions (i.e., the sum of organic carbon, black carbon, and inorganic matter), as derived for this work based on estimates from Ito et al. (2018), for the Fe content in the sub-micrometer and super-micrometer combustion aerosols. The historical anthropogenic emissions are taken here from the Community Emissions Data System (Hoesly et al., 2018) and the historical fire emissions from the BB4CMIP6 dataset (van Marle et al., 2017). We note, however, that the estimate of Fe emission from metal smelting remains highly uncertain and that further work is needed (Rathod et al., 2020). As for the biomass burning, the iron fractions in the fine particles are related to the combustion stages of flaming (0.46 ± 0.51 %) and smoldering (0.06 ± 0.03 %) fires, while the averaged iron fraction is used for coarse particles (3.4 %) (Ito, 2011). The global mean ratio of 0.04 gFe gBC⁻¹ for biomass burning in fine particles is consistent with that of 0.032 in the review paper by Hamilton et al. (2022). Fe-containing aerosol combustion emissions are considered to be insoluble (Ito, 2015), except for ship oil combustion, which is assumed to be mostly soluble, i.e., ∼ 79 % on average for the years 2000–2014. We note that the value of 79 % represents the high solubility of iron emissions in oil fly ash (Ito et al., 2021). Rathod et al. (2020) proposed a lower solubility in emissions (i.e., 47.5 % for iron sulfates), with an upper value, however, at ∼ 90 %. The year-to-year variation in anthropogenic combustion Fe-emission fractions follows Ito et al. (2018). On the contrary, for biomass burning Fe-emission fractions no such variation is provided. The average Fe fractions (per sector) for the years 2000 to 2014 applied to the total particulate carbonaceous emissions are also listed in Table S1.

EC-Earth3-Iron also includes OXL primary emissions from natural and anthropogenic wood-burning processes that mainly account for its rapid formation in the sub-grid plumes not represented in the model. Indeed, OXL is well correlated with elemental carbon and levoglucosan (Cao et al., 2017; Cong et al., 2015), which are observed at significant levels during biomass burning episodes in the Amazon (Kundu et al., 2010), suggesting that oxalic acid could be either directly emitted or formed rapidly via combustion processes. During biomass burning episodes, enhanced emissions of ionic species have been generally measured, indicating an average OXL mass concentration measured in plumes of ∼ 0.04 %–0.07 % w/w​​​​​​​ (Yamasoe et al., 2000). Furthermore, domestic wood combustion is a potential OXL source (Schmidl et al., 2008) since measurements indicate an OXL contribution to the total particulate concentrations of ∼ 0.09 %–0.28 % w/w. Gasoline engines may also contribute to total dicarboxylic acid mass emitted to the atmosphere (Kawamura and Kaplan, 1987), although their direct contribution to ambient OXL concentrations is generally found to be low (Huang and Yu, 2007) and is therefore neglected here. All in all, primary OXL sources are quite uncertain and, given the current estimates, may only have a limited impact on the calculation of its atmospheric concentrations (e.g., Myriokefalitakis et al., 2011).

2.2.2 Thermodynamic equilibrium and atmospheric acidity calculations

The gas and particle equilibrium calculations of NH${}_{\mathrm{3}}/$NH${}_{\mathrm{4}}^{+}$ and HNO${}_{\mathrm{3}}/$NO${}_{\mathrm{3}}^{-}$ have been substantially revised in EC-Earth3-Iron. In EC-Earth3-AerChem, EQSAM (Metzger et al., 2002) is used to determine the partitioning of NH${}_{\mathrm{3}}/$NH${}_{\mathrm{4}}^{+}$ and HNO${}_{\mathrm{3}}/$NO${}_{\mathrm{3}}^{-}$. In EC-Earth3-Iron, the ISORROPIA II thermodynamic equilibrium model (Fountoukis and Nenes, 2007) replaces EQSAM to determine the equilibrium between the inorganic gas and the aerosol phases. ISORROPIA-II calculates the gas–liquid–solid equilibrium partitioning of the K⁺-Ca²⁺-Mg²⁺-NH${}_{\mathrm{4}}^{+}$-Na⁺-SO${}_{\mathrm{4}}^{\mathrm{2}-}$-NO${}_{\mathrm{3}}^{-}$-Cl⁻-H₂O aerosol system and is used in the forward mode, assuming that all aerosols are in a metastable (liquid) state. The inclusion of sea salt and dust aerosols in the aerosol thermodynamic calculations has been shown to nevertheless substantially affect the ion balance and thus the partitioning of HNO${}_{\mathrm{3}}/$NO${}_{\mathrm{3}}^{-}$and NH${}_{\mathrm{3}}/$NH${}_{\mathrm{4}}^{+}$ species, especially in areas with abundant mineral dust and/or sea spray aerosols (Athanasopoulou et al., 2008, 2016; Karydis et al., 2016). In EC-Earth3-Iron nitrate aerosols are calculated for both the accumulation and coarse modes, in contrast to the bulk aerosol approximation used in the EC-Earth3-AerChem. For this, kinetic limitations by mass transfer and transport between the gas and the particulate phases in accumulation and coarse modes (Pringle et al., 2010) are considered, with ISORROPIA-II then re-distributing the respective masses between the gas and the aerosol phases. We note that Ca²⁺ from calcite is simulated prognostically in the model based on mineralogy maps (Sect. 2.2.1), in contrast to other crustal elements in soils that are calculated by assuming constant mass ratios to dust concentrations of 1.2 %, 1.5 %, and 0.9 % for Na⁺, K⁺, and Mg²⁺, respectively (Karydis et al., 2016; Sposito, 1989). For sea spray aerosols, mean mass fractions of 55.0 % Cl⁻, 30.6 % Na⁺, 7.7 % SO${}_{\mathrm{4}}^{\mathrm{2}-}$, 3.7 % Mg²⁺, 1.2 % Ca²⁺, and 1.1% K⁺ (Seinfeld and Pandis, 2006) are also applied.

The acidity levels of deliquescent aerosols are calculated in the model based on thermodynamic processes for accumulation and coarse particles. Aerosol acidity impacts the scavenging efficiency and the dry deposition of inorganic reactive nitrogen species due to changes in the partitioning of total nitrate and ammonium between the gas and aerosol phases and between the various aerosol sizes (Pye et al., 2020). Acidity levels also play a fundamental role in the aqueous-phase chemistry by controlling the dissociation reactions and thus the reactivity of the chemical mechanism. Indeed, aqueous-phase species, such as organic and inorganic acids, are oxidized with higher rates when they are dissociated. Nevertheless, in the case of the forward and reverse reactions, they typically occur fast and thus the concentrations of the reactants and the products are generally assumed to be in equilibrium in the global model due to its relatively long time step and large model grid. Note, however, that recent modeling studies showed that the metastable assumption produces pH values that are different from the stable assumption (e.g., regionally up to 2 pH units in the presence of crustal elements over dust sources, and roughly 0.5 pH units globally; Karydis et al., 2021). However, work to date, such as in Bougiatioti et al. (2016), Guo et al. (2019, 2015) and others identified in the review of Pye et al. (2020), has shown that the metastable solution tends to provide semi-volatile partitioning of pH-sensitive species (e.g., NH${}_{\mathrm{3}}/$NH₄ and HNO${}_{\mathrm{3}}/$NO₃) and aerosol liquid water content that is closer to observations – at least for when the relative humidity is above 40 %. For this reason, we assume that the most plausible estimates of acidity are to be obtained with the metastable assumption, and we base our simulations on that.

Under ambient atmospheric conditions, the water vapor uptake on aerosols depends on both the inorganic and organic components, along with the meteorological conditions (e.g., the temperature and the relative humidity conditions). ISORROPIA II does not, however, include water associated with organic aerosols, possibly leading to an underestimation of the aerosol hygroscopicity, especially within the boundary layer where the contribution of water-soluble organics to total aerosol mass can be substantial. For this, we account here for a contribution of aerosol water from organic particles in the acidity calculations,using a hygroscopicity parameter κ_org=0.15 (Bougiatioti et al., 2016). In more detail, the particulate water due to the organics (W_org) that is added to the aerosol water associated with the inorganic aerosol as calculated from ISORROPIA-II (W_inorg) is determined in the model as follows:

$$\begin{array}{}\text{(1)}& W_{\mathrm{org}}=m_{\mathrm{s}}\cdot {\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{s}}}}\cdot {\displaystyle \frac{{\mathit{\kappa }}_{\mathrm{org}}}{(\frac{\mathrm{1}}{\mathrm{RH}}-\mathrm{1})}},\end{array}$$

where m_s is the soluble organic mass concentration (µg m⁻³) as simulated by the TM5-MP chemistry scheme, ρ_w is the water density (1 kg m⁻³), ρ_s is the organic aerosol density (1.4 kg m⁻³), and RH (0–1) is the relative humidity.

Cloud acidity is also an important factor for simulating the multiphase processes in the atmosphere. The in-cloud proton concentration is initially determined by the electro-neutrality of strong acids and bases (i.e., H₂SO₄, SO${}_{\mathrm{4}}^{\mathrm{2}-}$, methanesulfonate (MS⁻), HNO₃, NO${}_{\mathrm{3}}^{-}$, and NH${}_{\mathrm{4}}^{+}$), and then the subsequent dissociations of CO₂, SO₂, and NH₃ (Jeuken et al., 2001) are solved iteratively in the model. For the cloud acidity calculations, the liquid water content, and the respective cloud cover fraction (i.e., 0–1) are obtained from meteorology. Note, however, that the effect of mineral dust (especially calcium) on cloud proton concentrations is neglected here. This assumption may result in some overestimation of cloud acidity, although the overall impact should be small, particularly in dusty areas with a low presence of clouds. Another limitation is the omission of light gaseous organic acids (such as formic and acetic acids) in the cloud pH calculations, possibly leading to some underestimation in cloud acidity where their concentration is important.

2.2.3 The aqueous-phase chemistry scheme

The aqueous-phase chemistry scheme used in this work is based to a large extent on the Chemical Aqueous Phase Radical Mechanism (CAPRAM) (e.g., Deguillaume et al., 2004; Ervens et al., 2003; Herrmann et al., 2000, 2015). However, CAPRAM includes more than 70 aqueous-phase species, 34 equilibria for compounds that are present both in the gas and the aqueous phases, along with numerous photolytic and aqueous-phase reactions, also covering a large series of acid–base and metal–complex equilibria. Note that various updates may further extend the mechanism by including, among other processes, the oxidation of aromatic hydrocarbons (Hoffmann et al., 2018), the multiphase oxidation of DMS (Hoffmann et al., 2016), and the tropospheric multiphase halogen chemistry (Bräuer et al., 2013). For this, some reactions are considered here in a more simplified way based on various assumptions published in the literature. Indeed, the level of chemical complexity of such a detailed mechanism is beyond the computational resources available for three-dimensional global climate–chemistry simulations, and thus simplifications that preserve however the essential features of the aqueous mechanism are needed.

Aqueous-phase chemical transformations are considered at the interface and in the bulk, initiated mainly by free radicals and oxidants produced both via photochemical reactions and in dark conditions (Bianco et al., 2020). The sources of OH radicals in the aqueous phase, however, strongly differ from those in the gas phase, primarily because of the presence of ionic species and TMIs in the solution. OH radicals are the main oxidant in the aqueous phase, either produced directly in the aqueous medium or diffused from the gas phase (i.e., via a gas-to-liquid transfer). However, aqueous-phase oxidation can also be induced by non-radical species, such as ozone (O₃) and hydrogen peroxide (H₂O₂). A characteristic example is the formation of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ in cloud droplets, via the oxidation of dissolved sulfur dioxide (SO₂) by O₃ and H₂O₂, with H₂O₂ nevertheless being the most effective oxidant (Seinfeld and Pandis, 2006), especially when the solution becomes acidic. Upon the absorption of SO₂ in cloud droplets, the establishment of the equilibrium between the dissolved sulfur species in oxidation state four, i.e., SO₂•H₂O, HSO${}_{\mathrm{3}}^{-}$ (pKa₁=1.9), and SO${}_{\mathrm{3}}^{\mathrm{2}-}$ (pKa₂=7.2) (hereafter also referred to as S(IV)) is calculated in the model. Thus, depending on the availability of oxidants and the solution's acidity, the different S(IV) species can participate in the formation of S(VI) (i.e., dissolved sulfur in oxidation state six).

In EC-Earth3-Iron, the aqueous-phase sulfur scheme is applied both in cloud droplets and aerosol water, replacing the S(VI) production through the dissolved S(IV) oxidation in cloud droplets previously included in the EC-Earth-AerChem (van Noije et al., 2014, 2021). In more detail, besides the two classic reactions of bisulfite and sulfite with hydrogen peroxide and ozone included in EC-Earth3-AerChem, additional reactions of S(IV) oxidation via methyl hydroperoxide (CH₃O₂H), peroxyacetic acid, and with the hydroperoxyl radical (HO₂)/superoxide radical anion (O${}_{\mathrm{2}}^{-}$) are considered. Nevertheless, in acidic solutions, the oxidation by peroxides, and especially H₂O₂, is significantly more important than other oxidants (Herrmann, 2003; Jacob, 1986). H₂O₂ is produced in the gas phase and can be rapidly dissolved in the liquid phase due to its high solubility. The dissolved H₂O₂ (as well as the organic peroxides, such as CH₃OOH) can react rapidly with the HSO${}_{\mathrm{3}}^{-}$. However, the pH-independent reaction of HSO${}_{\mathrm{3}}^{-}$ with CH₃OOH (or other organic peroxides) is expected to be less important than H₂O₂ under typical cloud conditions due to the much lower solubility of CH₃OOH. Note that the dissociation of H₂O₂ is neglected here since it is not expected to significantly influence the total H₂O₂ concentrations under typical tropospheric conditions (Herrmann, 2003; Jacob, 1986). In contrast, at a higher pH, the S(IV) oxidation by ozone tends to dominate the S(IV) oxidation (Seinfeld and Pandis, 2006). O₃ oxidizes rapidly all three S(IV) forms in the aqueous phase, becoming significant at pH higher than 4 (Seinfeld and Pandis, 2006), even in the absence of light. S(IV) oxidation by O₃ is also predicted to dominate S(VI) formation during winter in arctic regions due to the lack of photochemical production of OH and H₂O₂ at high latitudes, as well as the high anthropogenic SO₂ emissions in the Northern Hemisphere (Alexander et al., 2009). Laboratory studies indicate that S(IV) compounds may be also oxidized in the aqueous phase via other pathways. For example, the aqueous S(VI) production can be enhanced by TMIs (Harris et al., 2013), such as the Mn(II) catalyzed oxidation of S(IV) by dissolved O₂. In a global modeling study, Alexander et al. (2009) attributed 9 %–17 % of the total S(VI) production to the latter mechanism. However, such reactions would require several oxysulfur radicals as intermediates (e.g., Deguillaume et al., 2004; Herrmann et al., 2005), like a free radical chain mechanism initiated by reactions of HSO${}_{\mathrm{3}}^{-}$, SO${}_{\mathrm{3}}^{\mathrm{2}-}$ with radicals and radical anions, or TMIs catalyzed via oxidation of several S(IV) compounds, which is not considered in our model. Thus, in the case of the sulfate radical anion (SO${}_{\mathrm{4}}^{-}$) production via the Fe(III) sulfate complex [Fe(SO₄)]⁺ photolysis (Table S2), the sulfate radical anion is simply added to the S(VI) pool.

Gas-phase organics can be also oxidized in the interstitial cloud space, form water-soluble compounds like aldehydes, and rapidly partition into the droplets. In the presence of oxidants such as OH and NO₃ radicals in the solution, the dissolved organics undergo chemical conversions and form low-volatility organics that remain, at least partly, in the particulate phase upon droplet evaporation (Blando and Turpin, 2000). The dissolved OH radicals react with organic compounds in the aqueous phase by hydrogen abstraction or electron transfer, forming alkyl radicals (R), which in the presence of dissolved oxygen further form peroxyl radicals (RO₂). The OH oxidation of organic compounds in the aqueous phase can lead to either fragmentation or the formation of oxidized organic species, resulting overall in CO₂. However, the recombination of organic radicals can also be a favorable pathway when the water evaporates, and thus the aqueous solution becomes more concentrated. Box model simulations have shown that the cloud processing of polar products from isoprene oxidation can be an important contributor to secondary organic aerosol (SOA) production (Lim et al., 2005). Indeed, laboratory measurements show that the aqueous-phase photooxidation of C2 and C3​​​​​​​ carbonyl compounds (Perri et al., 2009, 2010), such as glyoxal (Carlton et al., 2007, 2009), methylglyoxal (Altieri et al., 2008), glycolaldehyde, pyruvic acid (Carlton et al., 2006), and acetic acid (Tan et al., 2012) leads to the production of low-volatility DCAs, which are commonly found in atmospheric aerosols and clouds (Sorooshian et al., 2006).

In EC-Earth3-Iron, gas-phase species can be reversibly transferred to the aqueous phase and oxidized by radicals and radical anions. The partitioning of 15 organic species that exist in both phases are considered in the aqueous-phase mechanism, namely methyl-peroxy radical (CH₃O₂), methyl hydroperoxide (CH₃O₂H), formaldehyde (HCHO), methanol (CH₃OH), formic acid (HCOOH), acetaldehyde (CH₃CHO), glycolaldehyde (GLYAL; HOCH₂CHO), glyoxal (GLY; CH(O)CH(O)), ethanol (CH₃CH₂OH), acetic acid (CH₃COOH), methylglyoxal (MGLY; CH₃C(O)CHO), hydroxyacetone (HYAC; CH₃C(O)CH₂OH), pyruvic acid (PRV; CH₃C(O)COOH), GLX, and oxalic acid (H₂C₂O₄). The aqueous-phase oxidation is taking place by the OH and NO₃ radicals, as well as the CO${}_{\mathrm{3}}^{-}$ radical anion. OH is either produced by photolytic reactions of dissolved compounds or via a direct transfer from the gas phase into the solution, as well as by Fenton reaction (Deguillaume et al., 2010). NO₃ radicals are transferred from the gas phase, while the CO${}_{\mathrm{3}}^{-}$ radical anion is produced mainly via the oxidation of hydrated CO₂. In general, the aqueous-phase oxidation largely proceeds via OH radicals, followed by NO₃ radicals under dark conditions, while the CO${}_{\mathrm{3}}^{-}$ radical has an overall small impact on the oxidizing capacity of the solution.

Upon their transfer to the solution, aldehydes are considered to be in equilibrium with the corresponding diols. The hydrated aldehydes are oxidized via H-atom abstraction with radicals (OH, NO₃) or radical anions (CO${}_{\mathrm{3}}^{-}$), followed by the elimination of HO₂ in reaction with O₂, leading overall to the formation of organic acids. Alcohols, such as CH₃OH and C₂H₅OH, are also oxidized via an H-atom abstraction; the resulting α-hydroxy-alkyl radicals, however, are not explicitly resolved, but the direct formation of aldehydes (e.g., formaldehyde and acetaldehyde) is considered via the respective peroxyl radical reactions with molecular oxygen to yield HO₂. Moreover, the glycolic acid (HOCH₂COOH) production via glycolaldehyde oxidation is not also explicitly described in the aqueous-phase scheme, and only the direct production of GLX is considered (Lin et al., 2012; Myriokefalitakis et al., 2011). This assumption is expected to have a negligible impact on the overall chemical mechanism since the glycolic acid is rapidly oxidized into glyoxylic acid with its net in-cloud production being rather small (Liu et al., 2012).

After cloud evaporation, OXL and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ are considered to reside entirely in the particulate phase of the model. This approximation may nevertheless result in an overestimate of OXL (pKa₁=1.23; (COO⁻)₂, pKa₂=4.19) concentrations, since low levels of gas-phase oxalic acid have been also observed in the atmosphere under favorable conditions (e.g., Baboukas et al., 2000; Martinelango et al., 2007). Note that other products, such as pyruvate, glyoxylate, and the oligomers from GLY and MGLY, are also considered to reside in the particulate phase upon cloud evaporation (Lim et al., 2005; Lin et al., 2012; Liu et al., 2012) and are thus added directly to the SOA pool of the model. However, in contrast to OXL and the low-volatility oligomers, the pyruvic and glyoxylic acids are allowed to be partially transferred back to the gas phase of the model when the cloud droplets evaporate.

For the present work, the aqueous reaction rate coefficients are taken (where available) from the available literature of the CAPRAM schemes and supplemented with reaction rates from laboratory and modeling studies (i.e., Carlton et al., 2007; Deguillaume et al., 2009; Lim et al., 2005; Sedlak and Hoigné, 1993). For the sulfur chemistry, the aqueous reaction rates are taken from Seinfeld and Pandis (2006). In the case of missing experimental data for temperature dependencies, the rate constants for T=298 K are only applied in chemistry calculations. O₃, H₂O₂, NO₃, HONO$/$NO${}_{\mathrm{2}}^{-}$, HNO${}_{\mathrm{3}}/$NO${}_{\mathrm{3}}^{-}$, CH₃O₂H, Fe³⁺, [Fe(SO₄)]⁺, and [Fe(OXL)₂]⁻ are photolyzed in the aqueous phase. Aqueous photolysis frequencies (where available) are taken from the gas-phase chemistry. For Fe species (e.g., Fe³⁺, [Fe(SO₄)]⁺, [Fe(OXL)₂]⁻), their maximum (i.e., noon at 51^∘ N) photolysis frequencies as proposed by Ervens et al. (2003) are scaled based on the gas-phase H₂O₂ photolysis rates. A list of all aqueous and photochemical reactions included in the chemical scheme of this study is presented in Table S2, with the respective equilibrium reactions shown in Table S3.

2.2.4 The iron solubilization scheme

A three-stage kinetic approach (Shi et al., 2011) is applied to describe the solubilization of the Fe-containing dust mineral pools (Ito and Shi, 2016), representing: (1) a rapid dissolution of ferrihydrite on the surface of minerals (i.e., fast pool), (2) an intermediate stage dissolution of nano-sized Fe oxides from the surface of minerals (i.e., intermediate pool), and (3) the Fe release from heterogeneous inclusion of nano-Fe grains in the internal mixture of various Fe-containing minerals, such as aluminosilicates, hematite, and goethite (i.e., slow pool). A separate Fe pool for combustion aerosols (Ito, 2015) is also considered in the model.

The dissolved Fe in the model is produced via dissolution processes in aerosol water and cloud droplets depending on the acidity levels of the solution (i.e., proton-promoted dissolution scheme), the OXL concentration (i.e., ligand-promoted dissolution scheme), and irradiation (photo-reductive dissolution scheme), following Ito (2015) and Ito and Shi (2016). The Fe release from different types of minerals thus depends on the solution acidity (pH) and the temperature (T), as well as on the degree of solution saturation. In more detail, the dissolution rates for each of the three dissolution processes considered can be empirically described (e.g., Ito, 2015; Ito and Shi, 2016; Lasaga et al., 1994) as follows:

$$\begin{array}{}\text{(2)}& {\mathrm{RFe}}_i=K_i\left(\mathrm{pH},T\right)\cdot \mathit{\alpha }({\mathrm{H}}^{+}{)}^{m_i}\cdot f_i\cdot g_i\end{array}$$

where K_i (mol Fe $g_i^{-\mathrm{1}}$ s⁻¹) is the Fe release rate due to the dissolution process i, α(H⁺) is the H⁺ activity of the solution, and m_i is the empirical reaction order for protons derived from experimental data. The functions f_i and g_i represent the suppression of the different dissolution rates due to the solution saturation state as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }{f}_i=\mathrm{1}-(a_{{\mathrm{Fe}}^{\mathrm{3}+}}\cdot a_{{\mathrm{H}}^{+}}^{-n_i})/K_{{\mathrm{eq}}_i},\text{(4)}& {\displaystyle }{g}_i=\mathrm{0.17}\cdot \ln \left({\displaystyle \frac{a_{\mathrm{OXL}}}{a_{{\mathrm{Fe}}^{\mathrm{3}+}}}}\right)+\mathrm{0.63},\end{array}$$

where ${\mathit{\alpha }}_{{\mathrm{H}}^{+}}$, ${\mathit{\alpha }}_{{\mathrm{Fe}}^{\mathrm{3}+}}$​​​​​​​, and α_OXL stand for the solution's activities of protons, ferric cations, and OXL, respectively, as calculated each time step in the model, and K_eqi (mol² kg⁻²) is the equilibrium constant. The activation energy that accounts for the temperature dependence is derived as a function of acidity based on soil measurements (Bibi et al., 2014; Ito and Shi, 2016), i.e.

$$\begin{array}{}\text{(5)}& E_{\mathrm{pH}}=-\mathrm{1.56}\times {\mathrm{10}}^{\mathrm{3}}\cdot \mathrm{pH}+\mathrm{1.08}\times {\mathrm{10}}^{\mathrm{4}}\end{array}$$

Overall, the net Fe dissolution rate results from the sum of the three rates. All parameters used for the calculation of dissolution rates for this work are presented in Table S4.

2.3 The chemistry solver

All concentrations of gas, aqueous, and aerosol species evolve dynamically in the model. The ordinary differential equations that govern the production and destruction terms due to chemical reaction and interphase mass transfer in the model are as follows:

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{d}G}{\mathrm{d}t}}=R_{\mathrm{G}}-\mathrm{LWC}{k}_{\mathrm{mt}}G+{\displaystyle \frac{k_{\mathrm{mt}}}{\mathrm{HRT}}}A,\text{(7)}& {\displaystyle \frac{\mathrm{d}A}{\mathrm{d}t}}=R_{\mathrm{A}}+\mathrm{LWC}{k}_{\mathrm{mt}}G-{\displaystyle \frac{k_{\mathrm{mt}}}{\mathrm{HRT}}}A,\end{array}$$

where G indicates gas-phase concentrations (molec. cm⁻³ of air), A indicates aqueous-phase concentrations (molec. cm⁻³ of air), R_G indicates gas-phase reaction terms (molec. cm⁻³ of air per second), R_A indicates aqueous-phase reaction terms (molec. cm⁻³ of air per second), LWC stands for liquid water content (cubic centimeter of water per cubic centimeter of air), k_mt indicates the mass transfer coefficient (s⁻¹), H indicates the Henry's Law coefficient (mol L⁻¹ atm⁻¹), R indicates the ideal gas constant (L atm mol⁻¹ K⁻¹), and T is temperature (K)

The mass transfer between the gas and aqueous phases (Lelieveld and Crutzen, 1991; Schwartz, 1986) is applied only for those species that exist in both phases and is represented in the mechanism by two separate reactions, i.e., one reaction for transfer from the gas to the aqueous phase and one for the transfer from the aqueous to the gas phase. All Henry's law solubility constants (H) used in this work are taken from Sander (2015) and are presented in Table S5.

The mass transfer coefficient (k_mt) for a species is calculated as follows:

$$\begin{array}{}\text{(8)}& k_{\mathrm{mt}}={\left({\displaystyle \frac{r^{\mathrm{2}}}{\mathrm{3}{D}_{\mathrm{g}}}}+{\displaystyle \frac{\mathrm{4}r}{\mathrm{3}\mathit{\upsilon }\mathit{\alpha }}}\right)}^{-\mathrm{1}},\end{array}$$

where r is the effective droplet or aqueous aerosol radius (m), D_g is the gas-phase diffusion coefficient (m² s⁻¹), υ the mean molecular speed (m s⁻¹), and α the mass accommodation coefficient (dimensionless). The cloud droplet effective radius may vary between ∼ 3.6 and 16.5 µm for remote clouds, 1 and 15 µm for continental clouds, and ∼ 1 and 25 µm for polluted clouds (Herrmann, 2003). For this work, the effective radius of cloud droplets (ranging between 4 and 30 µm in the model) is calculated online based on the cloud liquid water content and the cloud droplet number concentration (van Noije et al., 2021). The effective radii (i.e., the ratio of the third to the second wet aerosol moments) for the accumulation and coarse deliquescence particles are based on the respective M7 calculations. According to Eq. (8), the gas transfer to small droplets is faster, owing to the larger surface-to-volume ratio of smaller droplets. However, sensitivity model simulations using different droplet radii showed that varying droplet sizes result only in small changes in the chemical production of aqueous-phase species (Lelieveld and Crutzen, 1991; Liu et al., 2012; Myriokefalitakis et al., 2011).

The mean molecular speed of a gaseous species is calculated as follows:

$$\begin{array}{}\text{(9)}& \mathit{\upsilon }=\sqrt{\left({\displaystyle \frac{\mathrm{8}{R}_{\mathrm{g}}T}{\mathit{\pi }\mathrm{MW}}}\right)},\end{array}$$

where MW​​​​​​​ is the respective molecular weight (kg mol⁻¹) and R_g is the ideal gas constant (J mol⁻¹ K⁻¹) (Herrmann et al., 2000). The D_g and α used for this study are also presented in Table S5.

KPP version 2.2.3 (Damian et al., 2002; Sandu and Sander, 2006) was used to generate the Fortran 90 code for the numerical integration of the aqueous-phase chemical mechanism. For this, a separate model driver was developed to arrange the respective couplings to the TM5-MP I/O requirements (e.g., species that partition in the aqueous phase, the reaction and dissolution rates, and the photolysis coefficients). The Rosenbrock solver is used in this work as the numerical integrator since it is found to be rather robust and capable of integrating very stiff sets of equations (Sander et al., 2019). However, as for the case of the gas-phase mechanism's coupling (Myriokefalitakis et al., 2020b), minor changes needed to be applied to the original KPP code. For instance, the aqueous and photolysis reactions are not calculated inside KPP but directly provided through calculations in the aqueous chemistry driver. In contrast, for the Fe dissolution scheme, the suppressions of the mineral dissolution rates due to the solution saturation are calculated online by KPP (see Eqs. 3 and 4).

2.4 Simulations

We performed a range of present-day simulations, including experiments using EC-Earth3-Iron atmosphere-only runs (hereafter referred to as EC-Earth) and TM5-MP standalone driven by ERA-Interim (Dee et al., 2011) reanalysis fields (hereafter referred to as ERA-Interim), covering the period 2000–2014. For the EC-Earth simulation, TM5-MP is coupled to the IFS atmospheric dynamics. We used prescribed sea surface temperature and sea ice concentration fields from a set of input files through the AMIP interface (Taylor et al., 2000). Thus, for the atmosphere and chemistry modules, our setup follows the EC-Earth3-AerChem standard configuration in CMIP6 experiments. The IFS horizontal resolution is T255 (i.e., a spacing of roughly 80 km), 91 layers are used in the vertical direction up to 0.01 hPa, and a time step of 45 min is applied. Respectively, TM5-MP (both for the online and offline configurations) has a horizontal resolution of 3^∘ in longitude by 2^∘ in latitude and 34 layers in the vertical direction up to 0.1 hPa (∼ 60 km).

The ERA-Interim setup allows for constraining the model with the assimilated observed atmospheric circulation data and is therefore used for budget analysis and comparison with other estimates from the literature. ERA-Interim is further used to explore uncertainties regarding the aqueous-phase chemistry scheme. Specifically, an additional simulation is performed to identify the potential importance of glyoxal-derived oligomers and high molecular weight species in the aqueous phase (Carlton et al., 2007) on the OXL production rates and the respective ambient concentrations. In this sensitivity simulation (hereafter referred to as ERA-Interim(sens)), the OXL formation via formation of species of high molecular weight from glyoxal oxidation is neglected. Comparisons between the corresponding 15-year climatologies from the EC-Earth and ERA-Interim simulations are used to identify uncertainties in the aqueous-phase production terms of OXL, the iron-dissolution rates, and finally the atmospheric concentrations and deposition rates of Fe-containing aerosols due to the applied meteorology (i.e., online vs. offline). Note that the same emission datasets are used both in the ERA-Interim-driven and the EC-Earth experiments, thus only natural primary sources depending on meteorology may differ (see Sect. 2.1). A summary of the simulations is listed in Table 1.

Table 1 Overview of the simulations performed for this study.

2.5 Observations

A general evaluation of the modeled aerosol optical depth (AOD) at 550 nm allows for characterizing EC-Earth3-Iron's ability to reproduce the aerosol fields. The Aerosol Robotic Network (AERONET) version 3 (Giles et al., 2019) level 2.0 direct sun retrievals at a monthly basis are used to calculate annual mean AOD values for the 2000–2014 period. However, the model's coarse horizontal resolution hinders the representation of high-altitude locations; thus, following Huneeus et al. (2011), we exclude sites above 1000 m a.s.l., leaving 738 locations with information available during the simulated period. In addition, we perform a specific evaluation of mineral dust, which constitutes a key modulator of the outcome of our new developments as a source of Fe and Ca. To that end, we apply two additional filters to the AERONET data mentioned above, also following Huneeus et al. (2011), to identify dust-dominated sites. First, we exclude those sites where the monthly mean Ångström exponent is above 0.4 more than 2 months in the selected period. To further discriminate dust from sea salt, a minimum threshold of 0.2 for AOD at 550 nm is considered (i.e., if more than half of the retrieved AOD is above that threshold, the site is considered as dominated by dust). This filtering allows identifying a subset of stations potentially dominated by dust aerosols; however, it cannot ensure that there is no influence of other aerosol types in the monthly retrievals. Therefore, the evaluation of AOD at 550 nm at those sites is taken as a proxy for the dust optical depth, acknowledging that other aerosols may also be present.

Pure dust measurements of surface concentration and deposition complement our evaluation of the model. The modeled annual mean surface dust concentration for the years 2000–2014 is compared to climatological observations from the Rosenstiel School of Marine and Atmospheric Science (RSMAS) of the University of Miami (Arimoto et al., 1995; Prospero, 1996, 1999; Prospero et al., 1989) and the African Aerosol Multidisciplinary Analysis (AMMA) international program (Marticorena et al., 2010) observations. The 23 available sites cover locations close to sources (e.g., the AMMA stations over the Sahelian dust transect), in transport regions (e.g., stations from RSMAS in the Atlantic), and remote regions (e.g., RSMAS sites close to Antarctica). The modeled dust deposition fluxes are compared to the compilation of observations for the modern climate in Albani et al. (2014), including measurements at 110 locations, and the mass fraction for particles with a diameter lower than 10 µm is used to keep the observed mass fluxes within the range of the modeled sizes.

The simulated OXL and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations are compared against measurements for representative sites, such as the eastern Mediterranean (Finokalia, Greece; Koulouri et al., 2008), central Europe (Puy de Dome, France; Legrand et al., 2007), and the northern Atlantic Ocean (Azores, Portugal; Legrand et al., 2007). Simulated monthly mean surface concentrations of OXL are also compared against a range of observations (n=143) from remote sites around the world, as compiled in Myriokefalitakis et al. (2011). Moreover, SO${}_{\mathrm{4}}^{\mathrm{2}-}$ monthly mean surface concentrations over Europe and the USA are also compared against observations (n=3828) obtained from the European Monitoring and Evaluation Programme (EMEP; http://www.emep.int, last access 11 June 2021) and the Interagency Monitoring of Protected Visual Environments (IMPROVE; http://vista.cira.colostate.edu/improve/, last access 11 June 2021​​​​​​​), respectively, as compiled in Daskalakis et al. (2016). The simulated Fe-containing aerosol concentrations are evaluated against cruise measurements covering a period from late 1999 up to early 2015, as compiled by Myriokefalitakis et al. (2018) and Ito et al. (2019​​​​​​​), and include daily observations for fine, coarse, and total suspended particles.

Statistical parameters are here used to demonstrate the model's ability to represent atmospheric observations. These are the correlation coefficient (R) that reflects the strength of the linear relationship between model results and observations (i.e., the ability of the model to simulate the observed variability), the normalized mean bias (nMB), and the normalized root-mean-square error (nRMSE) as a measure of the mean deviation of the model from the observations due to random and systematic errors. The equations used for the statistical analysis of model results are provided in the Supplement (Eqs. S1–S3), and the locations (and regions) of the various observations used for evaluating the model for this work are presented in Fig. 1.

Figure 1 Site location map of observations (a) for AOD (AERONET, red dots; AERONET-DUST, blue squares), dust surface concentration (RSMAS, purple squares; AMMA, orange diamonds), and dust deposition rates (several sources compiled in Albani et al.,2014, green triangles) and (b) for surface oxalate (OXL, blue triangles), surface sulfate (green diamonds), and cruise aerosol Fe concentrations (red circles).

2.6 Model performance

The coupling of the aqueous-phase chemistry scheme along with the description of the atmospheric iron cycle for this work increases the model runtime. Here EC-Earth3-Iron uses 109 transported and 33 non-transported tracers, which are significantly larger numbers than in the EC-Earth3-AerChem configuration (i.e., 69 transported and 21 non-transported tracers). Note, however, that the EC-Earth3-Iron model used for this work employs the MOGUNTIA gas-phase chemistry scheme configuration, in contrast to the modified Carbon Bond Mechanism 2005 (mCB05) configuration (Huijnen et al., 2010; Williams et al., 2013, 2017) used in EC-Earth3-AerChem, which is overall found to be ∼ 27 % more expensive computationally (Myriokefalitakis et al., 2020b). In the Marenostrum4 supercomputer architecture (two Intel Xeon Platinum 8160 24C at 2.1 GHz), the EC-Earth3-AerChem configuration (van Noije et al., 2021) simulates 1.85 years per day of simulation time (SYPD) with 187 CPUs, while reaching a comparable performance (i.e., 1.41 SYPD) with the EC-Earth3-Iron configurations requires 432 CPUs. This means that the EC-Earth3-Iron corresponds to 7353 computation hours per year (CHPY) overall, which is roughly 3 times larger than the standard EC-Earth3-AerChem.

3 Results

3.1 Budget calculations

The chemical production and destruction terms of OXL and its precursors, along with the Fe-containing aerosols' dissolution rates from combustion (FeC) and mineral dust (FeD), their emissions, and their removal terms from the atmosphere, are presented for EC-Earth and ERA-Interim model configurations in this section. Additionally, we discuss differences compared to sensitivity simulations. Due to the common formation pathways of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ and OXL in the atmosphere, the SO${}_{\mathrm{4}}^{\mathrm{2}-}$ budget calculations are also presented and discussed. All calculations are presented as a mean (± standard error) for the years 2000–2014.

3.1.1 Oxalate

The annual net chemistry production of OXL (Table 2a) in EC-Earth is 12.615 ± 0.064 Tg yr⁻¹, which is lower than in ERA-Interim (18.116 ± 0.071 Tg yr⁻¹). The difference is explained by a higher oxidizing capacity in ERA-Interim than in EC-Earth. ERA-Interim calculates higher OH concentrations in the tropical and subtropical troposphere (Fig. S1b). In contrast, zonal mean OH levels in EC-Earth are slightly higher in the extratropics, causing a more efficient oxidation of the OXL precursors such as GLY (Fig. S1d), GLYAL (Fig. S1f), MGLY (Fig. S1h), and CH₃COOH (Fig. S1j) at higher latitudes, especially in the Southern Hemisphere (SH). Note that van Noije et al. (2014) also showed that the simulated oxidizing capacity in the previous version of EC-Earth (EC-Earth v2.4) was lower compared to a respective ERA-Interim configuration in large parts of the troposphere, due to the simulated lower temperatures (cold biases) and specific humidities. However, since sea surface temperatures (SSTs) and sea ice concentrations are prescribed in our EC-Earth atmosphere-only simulations, the long-term means of tropospheric temperatures and water vapor are not expected to differ significantly from ERA-Interim close to the surface levels, as also indicated by the low differences in the OH levels of the two simulations at low altitudes (Fig. S1b).

Table 2 Global budgets, atmospheric burdens, and lifetimes, averaged for the period 2000–2014, of (a) oxalate (OXL), (b) sulfate, and dissolved Fe-containing aerosols from (c) combustion processes (FeC) and (d) mineral dust (FeD) for EC-Earth, ERA-Interim, and ERA-Interim(sens) simulations.

The production term of OH via the H₂O pathway in EC-Earth is ∼ 5 % lower than in ERA-Interim due to a lower amount of water vapor being available to react with O(¹D). In addition, a ∼ 6 % lower OH production through the H₂O₂ photolysis is simulated in EC-Earth. Note that H₂O₂ is an important driver of the aqueous-phase oxidizing capacity in the model, with about 80 % of the OH radicals in the liquid phase being produced by photolysis of the dissolved H₂O₂. In more detail, the lower atmospheric abundance of the gas-phase H₂O₂ in EC-Earth (∼ 11 %) leads to smaller H₂O₂ uptake in the aqueous phase (∼ 13 %) and thus to a slower oxidation of OXL precursors due to the respective lower dissolved OH radical production (∼ 19 %). Overall, the total OH production is ∼ 7 % lower in EC-Earth, which corresponds to a ∼ 18 % lower aqueous-phase OH production, resulting in a ∼ 30 % lower OXL net chemistry production compared to ERA-Interim.

The total OXL production is 15.5 Tg yr⁻¹ in Lin et al. (2014) and 14.5 Tg yr⁻¹ in Liu et al. (2012), both of these values are lower than our ERA-Interim estimates (20.8 Tg yr⁻¹; Table 2a) but close to EC-Earth (15.0 Tg yr⁻¹; Table 2a). The main reason for the lower chemistry production of other published estimates compared to our results is the contribution of the aqueous-phase glyoxal oxidation scheme proposed by Carlton et al. (2007) that is applied in our simulations. The oxidation of the glyoxal-derived high-molecular-weight products formed mainly in the cloud droplets is calculated to contribute significantly to the global OXL production in our model (Table 2a). This result is in line with Carlton et al. (2007), who indicated that the GLX pathway may not be the primary pathway for oxalic acid formation, but this is instead attributed to the rapid oxidation of GLY multifunctional products via the OH radicals (i.e., ​​​​​​​3.1 × 10¹⁰ L mol⁻¹ s⁻¹; Table S2). However, for ERA-Interim(sens), where no such reactions are considered, the total OXL chemical production is calculated on average 11.5 Tg yr⁻¹ (Table 2); i.e., closer to the estimates of Lin et al. (2014) and Liu et al. (2012). On the other hand, our ERA-Interim net chemistry production calculations are close to the estimates of Myriokefalitakis et al. (2011) (i.e., 21.2 Tg yr⁻¹) when no potential effects of the ionic strength (e.g., Herrmann, 2003) on OXL precursors are considered, although this is still lower since no Fe chemistry was considered in that latter study. Indeed, the enhanced aqueous-phase oxidation capacity due to the Fenton reaction increases both the production and the destruction terms of OXL in our model, leading to ∼ 7 % lower net OXL production and a lower (∼ 8 %) atmospheric abundance, respectively. Nonetheless, our calculations indicate that Fe chemistry impacts on OXL net production drastically, increasing the destruction of the dissolved oxalic acid by at least ∼ 50 %. The potential primary sources (0.373 ± 0.005 Tg yr⁻¹) accounted for in the model (Table 2a) do not, however, significantly contribute to the simulated OXL atmospheric levels, and only a small fraction of OXL is calculated to be formed in aerosol water (∼ 6 %) for all simulations in this work.

Focusing further on the atmospheric sinks of OXL, roughly 13 % in ERA-Interim and 16 % in EC-Earth of the produced oxalic acid is oxidized into CO₂ in the aqueous phase, mainly via the photolysis of the [Fe(OXL)₂]⁻ complex (∼ 55 %) and via OH radicals (∼ 45 %). The fraction of the total produced OXL that is destroyed in the aqueous-phase is higher than in Liu et al. (2012) by ∼ 7 %, where no Fe chemistry was considered, but lower compared to Lin et al. (2014) and Myriokefalitakis et al. (2011), where roughly 30 % of the produced OXL is oxidized into CO₂ in the aqueous phase. Finally, a total average deposition rate of 18.5 Tg yr⁻¹ is calculated in ERA-Interim, primarily due to wet scavenging (∼ 99 %), resulting in a global atmospheric lifetime of 5.7 d, which is close to Liu et al. (2012) and Lin et al. (2014) but higher compared to Myriokefalitakis et al. (2011) (∼ 3 d); this is probably because of the more intense OXL production at higher altitudes in our model.

The major pathways of global OXL production, both in ERA-Interim and EC-Earth, are the oxidation of glyoxal (∼ 74 %), followed by glycolaldehyde (∼ 11 %), methylglyoxal (∼ 8 %), and acetic acid (∼ 7 %). Glyoxylic acid is nevertheless an important intermediate species because it is directly converted to OXL in the aqueous phase upon oxidation. Other important findings concerning the chemical budgets are summarized below.

1. Glyoxal. About 70 Tg yr⁻¹ GLY is produced in the gas-phase in ERA-Interim, similar to Lin et al. (2014), while in EC-Earth it is calculated 3 % lower. The global gas-phase production of the present work is higher than other global model estimates, e.g., about 56 Tg yr⁻¹ (Myriokefalitakis et al., 2008), 40 Tg yr⁻¹ (Fu et al., 2009, 2008), and 21 Tg yr⁻¹ (Liu et al., 2012). This difference can be explained by the more comprehensive isoprene chemistry of the gas-phase scheme used here (Myriokefalitakis et al., 2020b). Indeed, isoprene secondary oxidation products (e.g., epoxides) are significant precursors of GLY in the atmosphere (Knote et al., 2014) and the contribution of isoprene epoxides (IEPOX) from the gas-phase isoprene oxidation is here considered as a pathway of GLY formation. Note that the oxidation of other biogenic hydrocarbons, like terpenes and other reactive organics, may also result in GLY formation, since their chemistry is lumped on the first-generation peroxy radicals of isoprene in the model (Myriokefalitakis et al., 2020b). Besides the biogenic hydrocarbon oxidation, the model considers GLY formation due to the oxidation of other organic species (e.g., Warneck, 2003), such as acetylene (4.8 Tg yr⁻¹) and aromatics (18.8 Tg yr⁻¹). In the gas phase, other hydrocarbons, like ethene, further contribute to the atmospheric production of GLY via their oxidation products, mainly glycolaldehyde (5.4 Tg yr⁻¹). However, as in many modeling studies, additional primary and/or secondary glyoxal sources might be still missing in our model. Indeed, the elevated glyoxal concentrations over oceans that have been observed from space (e.g., Wittrock et al., 2006) would require at least 20 Tg yr⁻¹ of extra marine sources to reconcile model simulations with satellite retrievals (Myriokefalitakis et al., 2008). Great uncertainties, however, still exist on these oceanic sources (Alvarado et al., 2020; Sinreich et al., 2010), and therefore the only glyoxal primary sources accounted for in the model are from biofuel combustion and biomass burning processes (e.g., Christian et al., 2003; Fu et al., 2008; Hays et al., 2002), overall resulting in about 7 Tg yr⁻¹ on average in the model (Myriokefalitakis et al., 2020b). Glyoxal is rapidly destroyed in the atmosphere via photolysis (∼ 70 %), followed by its oxidation in the gas phase (∼ 15 %) and the aqueous phase (∼ 15 %). Roughly 5.4 Tg yr⁻¹ of glyoxal is produced in the aqueous phase via the dissolved GLYAL oxidation in ERA-Interim, close to the Liu et al. (2012) calculations but somehow higher compared to EC-Earth. Overall, the net cloud uptake of glyoxal in ERA-Interim is 6.3 Tg yr⁻¹, which is higher than the estimates from Liu et al. (2012) (1.6 Tg yr⁻¹). As expected, this increase is due to the applied glyoxal oxidation scheme in the aqueous phase of our base simulations. Finally, 4.2 Tg yr⁻¹ of glyoxal is removed from the atmosphere via wet scavenging (∼ 73 %) and dry deposition (∼ 27 %).

2. Glycolaldehyde. GLYAL is also a significant species for OXL atmospheric abundance since its oxidation directly produces GLY both in the gas and the aqueous phase. In ERA-Interim, the gas-phase production is 92.5 Tg yr⁻¹ on a global scale, with the primary sources accounting for 5.4 Tg yr⁻¹ (Myriokefalitakis et al., 2020b) on average. In EC-Earth, the gas-phase production is ∼ 1 % lower. GLYAL is destroyed via gas-phase photolysis (∼ 55 %) and by OH radicals in the gas phase (∼ 35 %) and the aqueous phase (∼ 10 %). Ethene oxidation products contribute ∼ 39 % to GLYAL production, but isoprene chemistry dominates its chemical production in the model. The only source of GLYAL in the aqueous phase is nevertheless the transfer from the gas phase. The dissolved GLYAL is oxidized to produce GLY (∼ 60 %) and GLX (∼ 40 %), overall resulting in a net aqueous uptake of 8.3 Tg yr⁻¹ in ERA-Interim, close to the estimates of Liu et al. (2012), but almost 40 % higher than in Lin et al. (2014). This higher uptake of GLYAL in the aqueous phase is due to the respective higher (∼ 14 %) gas-phase production in our model. Note that in ERA-Interim the net aqueous uptake of GLYAL is calculated ∼ 24 % lower compared to ERA-Interim.

3. Methylglyoxal. The global annual mean gas-phase production of MGLY in ERA-Interim is 237 Tg yr⁻¹ on average, with the primary sources accounting for 4.6 Tg yr⁻¹. The gas-phase production is higher than the 160–169 Tg yr⁻¹ reported by other modeling studies (Fu et al., 2008; Lin et al., 2014; Liu et al., 2012) owing to the contribution of oxidation products considered in the gas-phase isoprene chemistry scheme (Myriokefalitakis et al., 2020b). Roughly 56 % of MGLY is produced via the gas-phase oxidation of HYAC with OH radicals, which is lower than the estimated ∼ 75 % in Fu et al. (2008). The remaining MGLY production is due to isoprene oxidation products, i.e., ∼ 10 % from IEPOX oxidation, and ∼ 7 % from methyl vinyl ketone (MVK) and methacrolein (MACR) oxidation. In the aqueous phase, MGLY is produced via the dissolved HYAC oxidation (13.0 Tg yr⁻¹) and then further oxidized by OH radicals (11.6 Tg yr⁻¹) into pyruvic acid (PRV), methylglyoxal oligomers (MGLYOLI), and to a lesser extent into GLX. Note that the calculated contribution of dissolved HYAC to the aqueous-phase production of MGLY is higher compared to the nearly negligible rates in Liu et al. (2012) because of the higher gas-phase production of HYAC in our model. MGLY is chemically destroyed in the model mainly by gas-phase photolysis (∼ 60 %), the OH radicals in the gas phase (∼ 35 %), and via oxidation in the aqueous phase (∼ 5 %).

4. Pyruvic and acetic acids. The chemical production of PRV is 14.7 Tg yr⁻¹ in ERA-Interim and 16.7 Tg yr⁻¹ in EC-Earth. PRV is mainly produced by terpene oxidation via O₃ (∼ 51 %) in the gas phase followed by methyl vinyl ketone (MVK) oxidation (∼ 5 %). In the aqueous phase, PRV is solely produced from MGLY oxidation (6.5 Tg yr⁻¹) and subsequently oxidized to CH₃COOH. PRV is mainly removed via photolysis in the gas phase and via oxidation by OH radicals in the aqueous phase (∼ 30 %). However, more than half of the produced PRV in the aqueous phase directly contributes to the SOA mass of the model upon cloud evaporation. The gas-phase production of acetic acid is 44.3 Tg yr⁻¹, with the primary sources accounting for approximately 23.9 Tg yr⁻¹. In the aqueous phase, roughly 3 Tg yr⁻¹ of CH₃COOH is produced via PRV oxidation. Note that the net uptake of CH₃COOH (0.7 Tg yr⁻¹) is calculated in the model similar to the Lin et al. (2014) estimates but smaller than the 6.7 Tg yr⁻¹ calculated by Liu et al. (2012).

5. Glyoxylic acid. The GLX production rate is 7.1 Tg yr⁻¹ in EC-Earth and is ∼ 30 % lower in ERA-Interim. About 55 % of the produced GLX is directly oxidized to oxalic acid in the aqueous phase and ∼ 25 % is added directly to the SOA pool. Upon cloud evaporation, part of the produced GLX is also transferred in the gas phase, where it is either oxidized by OH radicals (∼ 60 %), photolyzed (∼ 33 %), or deposited (∼ 7 %). Due to the destruction of GLX in the gas phase, its total production is lower (∼ 60 %) compared to the production estimates in Lin et al. (2014) and Liu et al. (2012). For the EC-Earth and ERA-Interim, most of the produced GLX in the aqueous phase is derived from the oxidation of GLYAL (∼ 48 %), followed by the oxidation of CH₃COOH (30 %), GLY and its oligomeric products (∼ 15 %), and MGLY. The relative contributions in our calculations differ from the estimates in Lin et al. (2014), where GLX is primarily produced by GLY oxidation (∼ 77 %) followed by GLYAL (∼ 14 %), MGLY (∼ 1 %), and acetic acid (∼ 8 %). These differences are also caused by the direct contribution of the GLY oxidation products to the OXL formation. On the other hand, in the ERA-Interim(sens) simulation the calculated fractions agree well with other published estimates, where GLY overall dominates (∼ 60 %) the GLX production in the aqueous phase.

Figure 2a presents the annual mean (average 2000–2014) net chemistry production rates of OXL in EC-Earth and the respective absolute differences compared to ERA-interim (Fig. 2c). The maximum OXL production rates are calculated around the tropics and in the Southern Hemisphere, where both biogenic emissions (mainly isoprene) and the liquid cloud water are substantially enhanced (Fig. S2a). The Amazon region appears as the largest source of OXL, along with central Africa and Southeast Asia. At higher latitudes (>45^∘ N) the lower cloud liquid water content and vegetation cover lead to a lower OXL production over Asia and North America. However, over highly populated regions in the Northern Hemisphere, such as in Europe, the US, and China, enhanced OXL production rates are calculated due to its anthropogenic precursors. Furthermore, a significant source of OXL is calculated downwind of land areas, such as the South Pacific and the tropical Atlantic Ocean due to the long-range transport of OXL precursors.

Figure 2 Annual mean net chemical production rates for (a) oxalate (mg m⁻² yr⁻¹) and (d) sulfate (mg S m⁻² yr⁻¹) as calculated for the EC-Earth simulation averaged for the period 2000–2014, and the respective absolute differences to the ERA-Interim simulation (c, d).

The illustrated differences in OXL production between ERA-interim and EC-Earth (Fig. 2c) are caused due to the adopted atmospheric dynamics (i.e., online calculated versus offline), as both simulations use identical prescribed anthropogenic and biogenic emissions (see van Noije et al., 2021). For instance, EC-Earth calculates higher cloud water concentrations at ∼ 800–600 hPa around the tropics (30^∘ S–30^∘ N) compared to ERA-Interim. In contrast, lower concentrations are derived aloft, with ERA-Interim presenting enhanced cloud water concentrations at ∼ 400 hPa (Fig. S2b). Moreover, due to the lower OH concentrations in the tropical and subtropical troposphere (Fig. S1b), EC-Earth gives lower OXL production rates, especially over intense biogenic emission areas. Overall, the difference in the oxidizing capacity of the atmosphere between the two configurations significantly impacts the aqueous-phase OXL production efficiency in the model.

3.1.2 Sulfate

Sulfate (SO${}_{\mathrm{4}}^{\mathrm{2}-}$) is the main inorganic aerosol species produced in the aqueous phase, and similar to OXL its production in the model mainly occurs in cloud droplets. In addition, these two species largely reside in the aerosol accumulation mode of the model (roughly 99 % for SO${}_{\mathrm{4}}^{-\mathrm{2}}$ and 97 % for OXL). SO${}_{\mathrm{4}}^{\mathrm{2}-}$ is a key species for determining atmospheric acidity, and therefore here we also present the sulfate budget in conjunction with that of OXL. Sulfate is produced both in cloud droplets and in aerosol water, with the production in aerosol water having a negligible contribution on a global scale. In contrast to OXL, for which no gas-phase production is considered, the gas-phase oxidation of SO₂ via OH radicals contributes to the total SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations with about 12.0 Tg S yr⁻¹ (Table 2b). Our global estimate of the gaseous sulfuric acid (H₂SO₄) production is higher than in EC-Earth v2.4 (7.8 Tg S yr⁻¹ averaged for the years 2000–2009; van Noije et al., 2014) but slightly lower than in the EC-Earth3-AerChem AMIP simulations (van Noije et al., 2021) used for the CMIP6 experiments (available in https://esg-dn1.nsc.liu.se/search/cmip6-liu/, last access: 11 June 2021) where 12.9 Tg S yr⁻¹ of H₂SO₄ are produced (averaged for the years 2000–2014). These differences can be directly attributed to the OH radical production rates in the gas phase between the new and the previous chemistry versions of the atmospheric model, as have been discussed in Myriokefalitakis et al. (2020b). Despite the generally lower gas-phase OH radical levels (Fig. S1b), the slightly higher (∼ 8 %) global H₂SO₄ gas-phase production rate in EC-Earth than in ERA-Interim (Table 2b) can be attributed to the higher (∼ 6%) DMS emissions in EC-Earth (Fig. S4b) that contribute to the atmospheric SO₂ levels over the ocean.

The aqueous-phase SO${}_{\mathrm{4}}^{\mathrm{2}-}$ chemistry production from the oxidation of dissolved SO₂ is 39.8 Tg S yr⁻¹ in EC-Earth (Table 2b), which is higher than in EC-Earth v2.4 (29.3 Tg S yr⁻¹; van Noije et al., 2014) and EC-Earth3-AerChem (32.5 Tg S yr⁻¹). The higher SO${}_{\mathrm{4}}^{\mathrm{2}-}$ chemical production is mainly due to the higher SO₂ aqueous-phase oxidation rates by H₂O₂. In more detail, our calculations show that ∼ 84 % of the global SO${}_{\mathrm{4}}^{\mathrm{2}-}$ production in EC-Earth is due to the dissolved SO₂ oxidation via H₂O₂; 33.3 Tg S yr⁻¹ is produced due to H₂O₂, which is higher compared to 23.9 Tg S yr⁻¹ in van Noije et al. (2014). The dissolved SO₂ oxidation via O₃ (6.4 Tg S yr⁻¹) is also higher than in EC-Earth v2.4 (5.4 Tg S yr⁻¹). However, the contribution of CH₃O₂H to the SO${}_{\mathrm{4}}^{\mathrm{2}-}$ aqueous-phase production is small (0.05 Tg S yr⁻¹) in the model, with the HO₂ contribution being practically negligible on the global scale (0.02 Tg S yr⁻¹) for all simulations performed in this study. A total annual mean deposition rate of 52.5 Tg S yr⁻¹ is simulated in EC-Earth, with wet scavenging dominating the total deposition rate (∼ 93 %). Note that 2.5 % of the sulfur in the SO₂ emissions (1.6 Tg S yr⁻¹) is assumed to be in the form of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ for all simulations, which accounts for its formation in the sub-grid plumes (Aan de Brugh et al., 2011; Huijnen et al., 2010). Overall, a global SO${}_{\mathrm{4}}^{\mathrm{2}-}$ lifetime over deposition of 4.8 d is calculated in EC-Earth, which is lower than in ERA-Interim (6.6 d) but similar to the EC-Earth v2.4 estimate (4.9 d).

Figure 2b also shows the annual mean SO${}_{\mathrm{4}}^{\mathrm{2}-}$ net chemistry production rates in EC-Earth. High SO${}_{\mathrm{4}}^{\mathrm{2}-}$ production rates are calculated downwind of major anthropogenic SO₂ emission hotspots, such as central Europe, the eastern US, India, Russia, and eastern Asia. Furthermore, relatively high production rates due to biomass burning and volcanic eruptions are calculated in South America, southern Africa, and Indonesia. Significant SO${}_{\mathrm{4}}^{\mathrm{2}-}$ production is calculated over almost all oceanic regions due to the SO₂ production via the gas-phase oxidation of marine DMS emissions (Fig. S4a). Compared to the ERA-Interim simulation, however, the SO${}_{\mathrm{4}}^{\mathrm{2}-}$ production rates in EC-Earth are on average slightly higher over land in the tropics and extratropics. This increase can be attributed to combined effects that result in differences in chemical production and deposition rates (Table 2b). Some differences over oceans are nevertheless expected due to the differences in DMS concentrations, since DMS emissions are calculated online in the model based on sea surface temperature and wind velocity (Fig. S4b).

3.1.3 Iron

In EC-Earth, the total Fe (TFe) soil emissions result in 59.33 ± 1.22 Tg yr⁻¹, while in ERA-Interim they are 48.96 ± 0.95 Tg yr⁻¹. This difference results from the differences in wind speed between EC-Earth and ERA-Interim. EC-Earth produces higher dust emissions over large parts of the Middle East and Asia compared to ERA-Interim (Fig. S4f), which explains the differences in TFe emissions (TFeC emissions do not differ). However, most of the dissolved Fe from mineral dust in the model originates from atmospheric dissolution processes. In EC-Earth, FeD is primarily dissolved due to aerosol acidity at 0.31 Tg yr⁻¹, followed by the ligand-promoted dissolution that additionally produces 0.17 Tg yr⁻¹, while the photoinduced processes have a small impact on the global dissolved Fe release from dust, with 0.05 Tg yr⁻¹ (Table 2d). Fe primarily resides (98.4 %) in the slow pool of Fe-containing dust aerosols in the model, in particular in the coarse mode, with about 1.0 % being emitted as nano-sized iron oxides (intermediate pool) and 0.5 % as ferrihydrite (fast pool). Thus, most of the dissolved Fe release originates from the heterogeneous inclusion of nano-Fe grains in the internal mixture of various Fe-containing minerals such as aluminosilicates, hematite, and goethite (∼ 66 %), followed by nano-sized iron oxides (∼ 24 %) and to a lesser extent by ferrihydrite (∼ 10 %). Note, however, that the Fe release from aluminosilicates, hematite, and goethite particles is a slower process compared to the other soil classes considered in the model, as dictated by the three-stage approach applied for this study (Table S4).

Fe emissions from combustion processes are estimated at 2.518 ± 0.105 Tg yr⁻¹ in both simulations, with 0.012 Tg yr⁻¹ being emitted as dissolved from the primary oil combustion processes. Roughly 0.274 ± 0.010 Tg yr⁻¹ are released through Fe dissolution from combustion aerosols in EC-Earth, in good agreement with ERA-Interim (0.285 ± 0.011 Tg yr⁻¹). The acid-promoted dissolution contributes ∼ 17 % and photo-reductive processes ∼ 16 % to the Fe release from combustion particles, thus most Fe release comes from the ligand-promoted dissolution. This result is in line with laboratory studies (e.g., Chen and Grassian, 2013), where the contribution of oxalate-promoted dissolution is several times larger than the proton-promoted pathway under highly acidic dark conditions. According to our calculations, the relative contribution of atmospheric processing to the combustion aerosol Fe solubilization (∼ 11 %) is significantly higher compared to that of crystalline dust minerals (∼ 1 %), in agreement with laboratory (e.g., Chen et al., 2012; Fu et al., 2012) and modeling (e.g., Ito, 2015; Ito and Shi, 2016) studies.

The annual mean dissolution rates of FeC and FeD in EC-Earth are presented in Fig. 3. For combustion aerosols, the maximum dissolution rates occur downwind of biomass burning sources and highly populated regions, such as South America and Central Africa, the Middle East, India, and China. High dissolution rates are more likely to coincide with high OXL concentrations (Ito, 2015). Indeed, the model calculates important dissolution rates near regions where the OXL production rates are enhanced (Fig. 2a), such as over the Amazon basin and central Africa, as well as downwind of these regions, as the combustion aerosols are transported to the open ocean, in agreement with observations (e.g., Sholkovitz et al., 2012). For the mineral dust aerosols, most of the FeD dissolution fluxes occur downwind of the major dust source regions (e.g., the Sahara and the Gobi Desert), where the atmospheric transport of anthropogenic pollutants, such as SO_x and NO_x, enhances atmospheric acidity; e.g., the Fe release from the dust minerals due to proton-promoted dissolution processes is enhanced over the Middle East. Significant dissolution rates are also simulated over the Atlantic Ocean at the outflow of the Sahara, as well as at the outflow of Asian desert regions to the Pacific Ocean. High rates due to the contribution of the organic ligand-promoted dissolution processes are calculated downwind of central Africa and the equatorial Atlantic Ocean, where the oxidation of biogenic hydrocarbons in the presence of cloudiness leads to enhanced OXL aqueous-phase formation rates. On the contrary, the efficiency of ligand-promoted dissolution is substantially suppressed near dust source regions due to the low OXL availability (Fig. 2a).

Figure 3 Annual mean dissolution rates (mg m⁻² yr⁻¹) of combustion (a) and mineral dust (b) aerosols, as calculated for the EC-Earth simulation averaged for the period 2000–2014, and the respective absolute differences to the ERA-Interim simulation (c, d).

Figure 3 also presents the absolute differences between the ERA-Interim and EC-Earth annual mean Fe dissolution rates. ERA-Interim has significantly lower dissolution rates in the tropics (e.g., central Africa) and around the Equator, both for FeC (Fig. 3c) and FeD (Fig. 3d). This decrease is attributed both to the differences in atmospheric dynamics between the two model configurations and the suppression of the organic dissolution processes with lower OXL production. Indeed, Fig. 2c shows that in ERA-Interim, OXL production rates increase in the tropics, impacting Fe dissolution rates. In contrast, FeC dissolution rates increase in ERA-Interim over the Arabian Peninsula, India, and eastern Asia, due to fluctuations in OXL production and aerosol acidity. EC-Earth also shows lower FeD dissolution rates over the northern Pacific in the outflow of Asia. These differences are due to a higher aerosol acidity (i.e., up to ∼ 1 pH unit; Fig. S3d, f) in ERA-Interim due to changes in the buffering capacity of dust promoted by the higher calcite emissions in EC-Earth (Fig. S4f). This is especially the case for coarse dust aerosols where the majority of the Fe resides. EC-Earth also shows differences (positive or negative) with ERA-Interim in the Fe dissolution rates over oceanic regions (Fig. 3e, d), likely due to differences in SO${}_{\mathrm{4}}^{\mathrm{2}-}$ production over oceans from marine DMS emissions (Fig. S4b) and the impact of sea salt emissions (Fig. S4d) upon the buffering capacity of the solution. All in all, the total DFe atmospheric source in EC-Earth, accounting for both primary emissions and atmospheric processing, is 0.806 ± 0.014 Tg yr⁻¹ for the present day, well within the range of estimates presented in the model intercomparison study (0.7 ± 0.3 Tg yr⁻¹) in Myriokefalitakis et al. (2018).

3.2 Evaluation of new model features against observations

All developments described in this work have been implemented over the EC-Earth3-AerChem model version, which has been proven to simulate the atmospheric aerosol cycles better than other global models and to reproduce satisfactory the optical properties (Gliß et al., 2021). Thus, we do not expect substantial changes in EC-Earth3-Iron ability to represent the aerosol cycle or their optical properties compared to EC-Earth3-AerChem. However, owing to the significant differences in the gas-phase and aqueous chemistry between versions, we provide an overall assessment of the aerosol optical depth (AOD). In addition, as one of the novelties of this work is to consider explicitly how dust composition affects the atmospheric iron burden and alters acidity (e.g., through calcite), a comparison of dust fields with in situ observations is also provided. Finally, simulations of specific species key to our developments, such as oxalate, sulfate, and total and soluble iron are also evaluated.

3.2.1 AOD, dust concentration, and deposition

The annual mean AOD at 550 nm modeled in EC-Earth for 2000–2014 compares favorably with AERONETv3 direct-sun level 2.0 data (Fig. 4a). Overall, the model presents an nMB of −9 % and an nRMSE of 46 %, considering information from 738 AERONET sites. The regional analysis suggests a slightly better behavior in northern hemisphere regions (e.g., North America, Europe, East Asia) dominated by anthropogenic aerosols (normalized errors and biases, below 45 % and ± 10 %, respectively). The largest deviations from observed AOD occur in the Southern Hemisphere (e.g., South Africa, Australia, and Oceania) or remote regions. Over dust-dominated regions (e.g., North Africa, West Asia, and the Middle East) the model also behaves well (with normalized errors and biases below 45 % and ± 10 %, respectively). Selecting specifically dust-dominated sites for the comparison (Fig. 4b) and following the criteria explained in Sect. 2.5, EC-Earth slightly overestimates the retrieved AOD at 550 nm over North Africa (nMB = 21 %) and shows underestimations over sources in West Asia and the Middle East, as well as in transport regions such as Central America. In general, EC-Earth's ability to reproduce the annual mean AOD at 550 nm holds for dusty sites, with a normalized mean bias of 3 %, and a normalized root-mean-square error of 37 % over 38 sites. Overall, the average optical depth for dust at 550 nm (annual mean over the 2000–2014) yields a value of 0.032 ± 0.005, which falls well in the range of observationally based estimates based on in situ measurements, satellites, and global models (i.e., 0.030 ± 0.005; Ridley et al., 2016).

Figure 4 Comparison of (a) the modeled annual mean AOD at 550 nm against AERONET retrievals for all available stations covering the 2000–2014 period, (b) the same but for selected dusty AERONET sites, (c) the same but for modeled annual mean dust surface concentration for 2000–2014 compared to climatological mean values from RSMAS sites and AMMA campaign, and (d) the same but for modeled annual dust deposition flux averaged for the period 2000–2014 against observations compiled in Albani et al. (2014) from several sources.

The comparison of model outputs with climatologies of dust surface concentration from the RSMAS and the AMMA campaign (Fig. 4c) yields slightly poorer results, with an nMB of 19 % and an nRMSE of 81 %, as an average of the 23 sites available. EC-Earth best reproduces dust surface concentrations over source regions, such as North Africa (nMB = 21.9 %, nRMSE = 37.7 %), shows underestimations in transport areas (e.g., Central America: nMB $=-\mathrm{37}$ %; nRMSE = 39 %) and poorly represents the surface concentration in remote regions (e.g., the South Pacific and Southern Ocean, with nMB up to −98 % and nRMSE up to 113 %). The evaluation of the dust deposition field (Fig. 4d) shows both positive and negative biases over source and transport regions (see Table A1), with the deposited mass being generally underestimated, except for the Southern Ocean where the model tends to overestimate the observations. EC-Earth3-Iron may thus share the difficulties of many global models in representing the long-range transport of dust, particularly coarse particles downwind of dust sources (e.g., Adebiyi and Kok, 2020). As minerals in dust constitute the primary source of TFe to the atmosphere, the aforementioned discrepancies with respect to observations (e.g., higher concentrations over dust source areas and an underestimation of the dust deposition rates) are expected to also affect the representation of dust-related DFe in the model.

3.2.2 Oxalate

The averaged OXL surface concentrations in EC-Earth for the boreal winter (December, January, and February, i.e., DJF) and summer (June, July, and August, i.e., JJA) are presented in Fig. 5. OXL surface concentrations are distributed roughly between 60^∘ S and 60^∘ N, mainly in regions where intensive volatile organic compound (VOC) emissions from anthropogenic and biogenic sources coexist with cloud water. The highest OXL concentrations are calculated over tropical Africa, the Amazon Basin, eastern Asia, the eastern United States, and Europe, clearly showing the strong impact of OXL precursors (e.g., glyoxal) and the availability of cloud water. In the Northern Hemisphere, OXL concentrations are generally calculated higher in summer and lower in winter, indicating a strong impact of temperature and photochemistry on the production rate of oxalic acid in the aqueous phase. During DJF, the model calculates lower OXL concentrations over midlatitude and high-latitude regions, such as East Asia, central Europe, and the northern US. Over these highly populated regions, the aerosol water content is enhanced, following the increased SO${}_{\mathrm{4}}^{-\mathrm{2}}$ production due to anthropogenic activities, and the aqueous-phase OXL production in deliquesce particles also contributes to OXL atmospheric concentrations. Furthermore, high OXL concentrations are calculated in the tropics for both seasons due to the photochemical activity and the intense sources of biogenic VOCs in these regions.

Figure 5 Oxalate (OXL) surface concentrations (µg m⁻³) for the boreal winter (DJF; a) and boreal summer (JJA; b), as simulated for the EC-Earth simulation averaged for the period 2000–2014, and the respective absolute differences to the ERA-Interim simulation for surface (c, d) and zonal mean (e, f).

Figure 5 further presents the differences of OXL concentrations between EC-Earth and ERA-Interim. In general, OH levels in ERA-Interim are higher, which causes a more efficient oxidation of OXL precursors for both seasons. Moreover, ERA-Interim shows higher concentrations around the intertropical convergence zone (ITCZ) due to differences in meteorology between the two simulations, as discussed above. During boreal winter, some differences are observed over the subtropics of the Northern Hemisphere. Although OXL concentrations are very low over these latitudes, the relatively strong increase in liquid water that serves as a medium for OXL production, both for clouds (Fig. S2b) in higher altitudes (∼ 400 hPa) and at the surface for deliquescent particles (Fig. S2d, e) in the ERA-Interim simulation. In the vertical, OXL concentrations are distributed in the model from the surface to ∼ 400 hPa with a maximum at around 900 hPa. The zonal mean differences, however, indicate strong increases in the Southern Hemisphere (30^∘ S–0^∘​​​​​​​) during boreal winter (Fig. 5e). Compared to EC-Earth, ERA-Interim calculates higher OXL concentrations in the upper troposphere for both seasons (Fig. 5e, f), mostly due to more efficient transport of OXL precursors by deep convection into the tropical and extratropical upper troposphere. In the lower and middle troposphere, higher concentrations are calculated in ERA-Interim depending on the location and the season. The concentrations in EC-Earth are also lower than in ERA-Interim in the Northern Hemisphere (NH) extratropics during boreal summer due to a lower chemical production (Fig. 5f).

Figure 6 presents the comparison of the different model simulations performed for this work with OXL surface observations. OXL concentrations show a strong seasonal dependence, with maxima during the warm season due to the intense photochemical activity combined with the higher precursor abundance. Over the Mediterranean, and specifically the eastern part which is characterized by the long-range transport of air pollution and from surrounding urban centers (Kanakidou et al., 2011), the model underestimates the observed concentrations during winter at the Finokalia station in all simulations (Fig. 6a), either due to missing OXL primary and secondary sources or a too strong removal. During summer, ERA-Interim satisfactorily simulates the observed OXL levels, also representing the observed trend, which indicates that the model reproduces the mixing and aging of the air masses in the region under favorable meteorological conditions and intense solar radiation. EC-Earth calculates lower OXL concentrations than ERA-Interim due to the lower oxidizing capacity, thus underestimating the observed concentrations for all seasons. On the other hand, ERA-Interim(sens) tends to underestimate the observations for all seasons more than the other simulations, further indicating the important role of the secondary sources to OXL atmospheric concentrations in the region. At the Puy de Dôme site (Fig. 6c), which is located at 1450 m a.s.l., ERA-Interim underestimates the observed OXL concentrations, although it simulates them more realistically compared to EC-Earth, especially during summer (Fig. 2c). The seasonal variation in the area can be explained by the stronger upward transport of air masses during summer (Legrand et al., 2007), thus increasing the OXL production in the region. However, the model fails to represent the observed OXL levels, possibly due to missing sources. The importance of other production pathways not related to the aqueous-phase GLX oxidation is demonstrated in the comparison of the observed OXL levels with the ERA-Interim(sens) simulation. Again, ERA-Interim(sens) deviates more strongly from the measured values than other simulations. Nevertheless, this indicates that other species may further contribute to OXL production, such as the decay of longer diacids (e.g., azelaic and malonic acids) (Legrand et al., 2007) that are currently not included in the model. Another reason may be the impact of the enhanced cloud LWC in the region, implying a more intense cloud processing compared to other surface sites and thus a faster oxidation of oxalic acid into CO₂ (Ervens et al., 2004) in the model. Finally, at the Azores (Fig. 6e), a site that is characterized by a marine environment, the model tends to underestimate the observed OXL concentrations most of the time, with ERA-Interim again presenting a better skill than other simulations. EC-Earth underestimates the observed concentrations more than ERA-Interim, especially during summertime, and ERA-Interim(sens) simulates the lowest OXL concentrations. The observed OXL levels in the region, however, can be explained either by the transport of pollutants from the continents or the photochemical production in the region. Thus, the illustrated differences against the observations between EC-Earth and ERA-Interim can be attributed to differences in the oxidizing capacity and in simulated transport of EC-Earth, such as the vertical mixing in the troposphere (e.g., van Noije et al., 2014) that has a further impact on OXL precursors like glyoxal. Furthermore, since the long-range transport is found to be relatively constant in summer and winter in the region (Legrand et al., 2007), other species of marine origin, such as the unsaturated fatty acids (e.g., linoleic and oleic acids) may also contribute as precursors to the OXL production, especially during summer, but are not included in the model. All in all, we acknowledge that other formation pathways of OXL, both primary or secondary, may exist in the atmosphere (e.g., Baboukas et al., 2000); for example, higher DCAs (such as malonic, succinic, glutaric, and adipic acids) may act as precursors for smaller dicarboxylic acids like OXL, both in the gas phase (e.g., Kawamura and Ikushima, 1993) and the aqueous phase (e.g., Ervens et al., 2004; Lim et al., 2005; Sorooshian et al., 2006) and could further contribute to atmospheric OXL concentrations.

Figure 6 Comparison of daily mean observations (black line) of OXL (ng m⁻³; a, c, e, g) and nss-SO${}_{\mathrm{4}}^{\mathrm{2}-}$ (µg m⁻³; b, d, f, h) with the EC-Earth (orange line and circles) and the ERA-Interim (light blue line and triangles) simulations for Finokalia (Greece) (a), (b) for the period July 2004–July 2006 (Koulouri et al., 2008); for Puy de Dome (France) (c), (d) and Azores (Portugal)(e), (f) for the period September 2002–September 2004 (Legrand et al., 2007); and scatterplot comparisons for observations around the globe (d), (f); the solid line represents the 1:1 correspondence, and the dashed lines show the 10:1 and 1:10 relationships, respectively. For completeness, the comparisons for the sensitivity simulation ERA-Interim(sens) for OXL and the EC-Earth(AerChem-AMIP) (green line and squares) for sulfate are also presented. Gray-shaded areas represent the standard deviation of the observations and the color-coded error bars represent the model's standard error of the multi-annual mean for the individual observational period.​​​​​​​

In Fig. 6g, OXL observations reported in the literature are compared with monthly mean simulations. Due to the relatively low resolution of the global model (i.e., 3^∘ × 2^∘ in longitude by latitude), the spatial variability of urban emissions cannot be well resolved. Therefore, urban stations are omitted, and the comparison is limited to locations representative of background concentrations. All simulations tend to underestimate OXL observations, with lower biases in ERA-Interim (i.e., nMB $=-\mathrm{46}$ %, nRSME = 110 %). As expected, for almost all sites the ERA-Interim simulation calculates the highest OXL concentrations and the ERA-Interim(sens) the lowest (i.e., nMB $=-\mathrm{74}$ %, nRSME = 125 %). The latter indicates that additional production pathways need to be considered in modeling studies to capture the observed OXL concentrations. EC-Earth underestimates the observed concentrations more than ERA-Interim (i.e., nMB $=-\mathrm{64}$ % and nRSME = 117 %) but less than ERA-Interim(sens), also highlighting the importance of the atmospheric oxidating capacity and the atmospheric dynamics in the OXL production. A summary of statistics for the evaluation of the simulated OXL concentrations for the different simulations is presented in Table A2. Overall, our analysis indicates that the model either misses OXL sources (primary and secondary) or overestimates OXL sinks, especially during winter. Thus, under relatively low temperatures and irradiation, the model may not be representative of the fast secondary OXL production in wood-burning plumes or the secondary production through species produced by the oxidation of emitted from vehicles and other anthropogenic activities, such as ethane and aromatic hydrocarbons.

3.2.3 Sulfate

The averaged SO${}_{\mathrm{4}}^{\mathrm{2}-}$ surface concentrations as calculated in EC-Earth for the boreal winter and summer are presented in Fig. 7. During both seasons, high SO${}_{\mathrm{4}}^{\mathrm{2}-}$concentrations are simulated near or downwind major anthropogenic emission hotspots, where the vast majority of the surface SO₂ emissions from anthropogenic origin occur (e.g., Tsai et al., 2010). Enhanced SO${}_{\mathrm{4}}^{\mathrm{2}-}$ surface concentrations are also calculated downwind of biomass burning and volcanic eruptions, showing the overall impact of SO₂ primary sources and the abundance of cloud water over these latitudes. Over the remote oceans, however, DMS oxidation may significantly contribute to the SO${}_{\mathrm{4}}^{\mathrm{2}-}$ surface concentrations, as is also the case for the sulfur emissions over major shipping routes. In Fig. 7, the differences between EC-Earth and EC-Earth-AerChem in the averaged SO${}_{\mathrm{4}}^{\mathrm{2}-}$ surface and zonal mean concentrations are also presented for both boreal winter (Fig. 7c, e) and summer (Fig. 7d, f). Considering, however, that the EC-Earth version developed for this work is based on the EC-Earth-AerChem model version, the illustrated differences are solely due to the applied chemistry schemes in the model. During boreal winter, the largest differences appear over eastern Asia, the Middle East, India, and central Europe. In EC-Earth-AerChem, the gas-phase SO₂ oxidation by OH radicals is roughly 8 % larger than in EC-Earth (see Sect. 3.1.2), leading to an overall lower conversion of SO₂ to sulfate in the aqueous phase. Moreover, in our simulations SO${}_{\mathrm{4}}^{\mathrm{2}-}$ is also produced in deliquesced particles, which partly contribute to the SO${}_{\mathrm{4}}^{\mathrm{2}-}$ atmospheric concentrations, especially over highly populated regions during the boreal winter when sulfur emissions due to anthropogenic activities are enhanced. On the contrary, during boreal summer some differences are illustrated (Fig. 7d) mostly over the Middle East where the EC-Earth-AerChem simulation results in lower SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations, and over eastern Asia, where some higher concentrations in EC-Earth-AerChem occur. These negative and positive differences are attributed to differences in oxidizing capacity between the two models in both the gas and the aqueous phase. Finally, the zonal mean differences indicate higher concentrations in EC-Earth in the lower troposphere of the Northern Hemisphere during boreal winter (Fig. 7e) than for the SO${}_{\mathrm{4}}^{\mathrm{2}-}$ surface concentrations. For boreal summer no important differences are presented here either (Fig. 7f).

Figure 7 Sulfate (SO₄) surface concentrations (µg m⁻³) for boreal winter (DJF; a) and boreal summer (JJA; b), as simulated for the EC-Earth simulation averaged for the period 2000–2014, and the respective absolute differences with the EC-Earth(AerChem-AMIP) simulation for surface (c, d) and zonal mean (e, f).

Figure 6 further presents the model comparison with SO${}_{\mathrm{4}}^{\mathrm{2}-}$ surface observations. Sulfate concentrations maximize under intense photochemical activity and high SO₂ atmospheric levels; generally also suggesting a faster formation rate compared to oxalic acid (Ervens et al., 2004; Legrand et al., 2007). At the Finokalia station in the Mediterranean (Fig. 6b), the model overestimates the observed SO${}_{\mathrm{4}}^{\mathrm{2}-}$concentrations during boreal winter and summer in ERA-Interim, probably due to too high SO₂ background concentrations because of a too strong long-range transport from surrounding regions. During winter, EC-Earth better reproduces the observations, probably due to the lower oxidizing capacity compared to ERA-Interim, but in late spring and early summer the active photochemistry in the region leads again to an overestimation of the observed concentrations. The EC-Earth-AerChem simulation leads to generally lower concentrations compared to our EC-Earth simulation, tending to somehow underestimate the observed concentrations except for autumn. At Puy de Dôme (Fig. 6d), ERA-Interim, which has a higher cloud liquid water content aloft and a more intense oxidizing capacity compared to EC-Earth, overestimates the observed SO${}_{\mathrm{4}}^{\mathrm{2}-}$concentrations in almost in all seasons. In contrast, EC-Earth better simulates the measured SO${}_{\mathrm{4}}^{\mathrm{2}-}$ concentrations, although it seems to underestimate the observations during summer. At that site, the EC-Earth-AerChem calculations agree well with EC-Earth, although concentrations are again slightly lower. At the Azores site (Fig. 6f), ERA-Interim also simulates the observed concentrations well, following the observed annual cycle. In contrast, both EC-Earth and EC-Earth-AerChem underestimate the SO${}_{\mathrm{4}}^{\mathrm{2}-}$ observations, especially during spring and summer. Note that the SO${}_{\mathrm{4}}^{\mathrm{2}-}$ production in this marine site is attributed to the SO₂ atmospheric levels both from air masses advected from industrialized regions and the local production due to the oxidation of marine DMS emissions. Thus, the differences between ERA-Interim and the other EC-Earth simulations presented in this study may indicate slower aging of the polluted air masses transported in the region. Finally, a comparison of the model's monthly mean predictions with a compilation of SO${}_{\mathrm{4}}^{\mathrm{2}-}$ observations (n=3828) around the globe (Daskalakis et al., 2016) is presented in Fig. 6h. EC-Earth tends to overestimate the available SO${}_{\mathrm{4}}^{\mathrm{2}-}$ observations (Table A2), presenting positive biases (i.e., nMB = 16 %, nRSME = 55 %) that are slightly lower than in ERA-Interim (i.e., nMB = 23 %, nRSME = 57 %). In contrast, EC-Earth-AerChem tends to slightly underestimate (Table A2) the observed concentrations (i.e., nMB $=-\mathrm{5}$ %, nRSME = 57 %), showing a slightly lower correlation coefficient (R=0.70) than EC-Earth (i.e., R=0.76) and ERA-Interim (i.e., R=0.75).

3.2.4 Dissolved iron

Figures 8 and 9 present the averaged dissolved FeC (DFeC) and FeD (DFeD) surface concentrations, respectively, for DJF and JJA. EC-Earth calculates an annual global DFeC atmospheric burden of 0.002 Tg, while ERA-Interim calculates slightly higher concentrations (0.003 Tg) due to the more intense ligand-promoted dissolution rates (Table 2c). Elevated DFeC concentrations during boreal winter (Fig. 8a) are calculated over central Africa, eastern Asia, and India, where significant DFe concentrations (∼ 0.01–0.1 µg m⁻³) are associated with biomass burning and anthropogenic combustion emissions. During boreal summer (Fig. 8b), the maximum DFeC concentrations are calculated in the northern latitudes, in particular over the Mediterranean Basin, the Middle East, the western US, and China. The increase in the surface dissolved Fe concentrations over these regions, ranging from ∼ 0.01 to 0.1 µg m⁻³, clearly highlights the anthropogenic contribution due to the enhanced solubilization of Fe when the combustion aerosols are mixed with acidic and organic pollutants during atmospheric transport (Fig. 3b). Due to the intense biomass burning in the Southern Hemisphere (i.e., South America, central Africa, and Indonesia), the enhanced OXL production rates over such regions (Fig. 2a) lead to higher dissolved Fe concentrations. Figure 8a and b demonstrate that the geographic pattern of the DFeC concentrations may change overall from boreal winter to summer, following the biomass burning activity and the atmospheric processing of Fe-containing combustion aerosols of anthropogenic origin.

Figure 8 Dissolved iron surface concentrations (µg m⁻³) from combustion aerosols (DFeC) for the boreal winter (DJF; a) and summer (JJA; b) seasons for the EC-Earth simulation, averaged for the period 2000–2014, and the respective absolute differences to the ERA-Interim simulation (c, d).

Figure 9 Dissolved iron surface concentrations (µg m⁻³) from mineral dust (DFeD) for the boreal winter (DJF; a) and summer (JJA; b) seasons for the EC-Earth simulation, averaged for the period 2000–2014, and the respective absolute differences to the ERA-Interim simulation (c, d).

Mineral dust emissions mostly occur in the midlatitudes of the Northern Hemisphere (Fig. 9a, b), relatively close to where the vast majority of the population exists, and the anthropogenic emissions of acidic compounds dominate. For both seasons, high DFe concentrations from dust occur over the midlatitudes of the Northern Hemisphere, where the major dust sources are located. However, the equatorial maximum during boreal winter tends to shift to the north during boreal summer following the migration of the ITCZ (Fig. 9b). DFe from mineral dust aerosols maximize over the major dust source regions, with surface concentrations of roughly 0.1–1 µg m⁻³ (Fig. 9a, b) overall dominating the Fe burden. The outflow from these regions transports DFe over the global oceans, where secondary maxima of ∼ 0.01–0.1 µg m⁻³ are calculated, mainly over the Northern Hemisphere in the tropical Atlantic Ocean. The dissolved Fe associated with Saharan dust is nevertheless attributed to the long-range transport and the atmospheric processing that converts the insoluble Fe minerals to soluble forms.

The differences between the EC-Earth and ERA-Interim are illustrated in the averaged DFeC (Fig. 8c, d) and DFeD (Fig. 9c, d) surface concentrations. The differences in DFeC between the EC-Earth and ERA-Interim are well correlated with those of OXL concentrations (Fig. 5c, d), indicating the strong impact of ligand-promoted dissolution on the DFeC atmospheric load (Table 2c). Note, however, that the most important relative differences between the two simulations are calculated over regions with low dissolved FeC concentrations, and consequently the total burden does not change significantly (Table 2c). The differences in the DFe concentrations associated with mineral dust aerosols follow the general anomaly pattern of the two model configurations due to differences in transport that lead overall to higher (∼ 15 % globally) dust emissions in EC-Earth (Fig. S4f). We note, however, that the annual mean dust emission in EC-Earth3-Iron amounts to 1257 ± 26 Tg yr⁻¹ (years 2000–2014 averaged), which falls in the lower range of the AEROCOM phase III models (Gliß et al., 2021) and is also at the low end of the range estimated by Kok et al. (2021) from inverse modeling (1200–2900 Tg yr⁻¹) for dust aerosol with a geometric diameter ≤10 µm.

Figure 10 presents a comparison of the different model simulations with cruise observations of dissolved Fe concentrations. The spatial distributions of the DFe observations for the accumulation, coarse, and total suspended aerosols are shown in Fig. 10a–c, respectively. The median (mean ± standard deviation) DFe concentration in the accumulation mode amounts to 0.96 (3.68 ± 6.44) ng m⁻³ in the observations, while in EC-Earth and ERA-Interim it is 3.67 (5.79 ± 5.17) ng m⁻³ and 3.12 (5.25 ± 4.54) ng m⁻³, respectively. The respective observed concentration of DFe in the coarse mode is 0.85 (5.03 ± 10.90) ng m⁻³, while in EC-Earth and ERA-Interim the calculated values are around 1.74 (6.18 ± 7.29) ng m⁻³ and 1.04 (5.33 ± 7.30 ng m⁻³). Finally, the concentration of DFe in total suspended particles (tsp) is 1.93 (7.83 ± 19.03) ng m⁻³ in the observations, and it is 4.37 (12.51 ± 19.09) ng m⁻³ and 4.39 (11.63 ± 17.12) ng m⁻³ in EC-Earth and ERA-Interim, respectively. The correlation coefficients between the median values of the DFe cruise observations and the model results for EC-Earth and ERA-Interim are calculated as 0.49 (nMB = 38 %, nRMSE = 164 %) and 0.58 (nMB = 23 %, nRMSE = 145 %) for the accumulation aerosols, while for the coarse mode these values are 0.46 (nMB $=-\mathrm{9}$ %, nRMSE = 193 %) and 0.59 (nMB $=-\mathrm{25}$ %, nRMSE = 177 %) for EC-Earth and ERA-Interim simulations, respectively (Fig. S5). A summary of statistics for the evaluation of the simulated DFe concentrations for the EC-Earth and ERA-Interim simulations is also presented in Table A2.

Figure 10 Observed (averaged for each model's grid-cell) dissolved iron (DFe) concentrations (ng m⁻³) of (a) accumulation aerosols, (b) coarse aerosols, and (c) total suspended particles (tsp). The respective absolute differences to the ERA-Interim simulation (d, e, f) and the comparison to observations (black line with crosses) in latitudinal order (g, e, f) with the EC-Earth (orange line with circles) and ERA-Interim (light blue line with triangles) simulations; the grey-shaded areas correspond to the standard deviation of the observations and the color-coded shaded areas/error bars correspond to the model's standard error of the multi-annual mean for the individual observational period.

The spatial distributions of the absolute differences between the simulated DFe concentrations in EC-Earth and the observations are also presented in Fig. 10d–e. The model shows a general overestimation of the observed DFe concentrations around the tropics (up to ∼ 10 ng m⁻³) but an underestimation at mid-to-high latitudes (up to ∼ 5 ng m⁻³). It is, however, unclear why the model is unable to capture the observed concentrations over such regions. Reasons could be the misrepresentation of coarse dust emissions, missing anthropogenic primary Fe emissions, weak secondary sources of the dissolved Fe (especially in the southern latitudes), or even systematic errors in the transport of coarse particles. Over the Pacific, EC-Earth better predicts the average DFe concentrations. For completeness, Fig. S5 also presents a comparison of the simulated and observed TFe concentrations, showing that the model better captures the TFe concentrations in the accumulation mode (Fig. S6g) than in coarse mode (Fig. S6h). Considering, however, that the TFe is mostly dominated by primary sources, the calculated differences to the observed concentrations (Fig. S6e–f) downwind continental sources should mainly depict errors in the emission parameterizations or a misrepresentation in the mineralogical composition of the larger Fe-containing soil particles.

The overestimation around the tropics and the northern latitudes is further illustrated by a comparison of the model predictions with observations as a function of latitude (binned at 2^∘) (Fig. 10g–i). Although the dissolved Fe in the accumulation mode (Fig. 10g) is simulated well over the Southern Ocean, the model strongly overestimates the observed concentrations in the tropics, especially around the Equator, as well as in the northern extratropics. In contrast, for the coarse and the total suspended aerosols, an underestimation of the DFe aerosol observations in the southern latitudes (i.e., around 30–60^∘ S) is clearly visible for both simulations (Fig. 10h and i), along with an overestimation around the tropics that is similar to that for the accumulation DFe aerosols. All in all, the differences between the two model configurations are relatively small (see Table A2).

4 Discussion

The ocean is a critical component of the Earth's climate system, and Fe plays a key role in the efficiency of the biological carbon pump. For this reason, accurate estimates of the bioavailable Fe inputs to the ocean are a prerequisite for climate simulations. Our work attempts to properly simulate the effects of atmospheric multiphase processes on the chemical sources and sinks of the Fe-containing aerosols in an Earth system model. Indeed, the atmospheric processing of Fe-containing aerosols is rather important for the geographical pattern of DFe deposition fluxes into the ocean, especially in remote regions away from land sources. In agreement with other studies (e.g., Hamilton et al., 2022; Ito et al., 2019; Mahowald et al., 2005; Myriokefalitakis et al., 2018; Scanza et al., 2018), we find that mineral dust is the principal source of atmospheric Fe in EC-Earth (95 ± 1 %), with most of the remaining sources attributable to biomass burning. Focusing on the bioavailable fraction for marine biota, OXL aqueous-phase production is shown to be an important driver of aerosol DFe release, contributing ∼ 44 % to the aerosol Fe dissolution, along with the atmospheric acidity that accounts for ∼ 45 % of total DFe secondary sources in the model. Therefore, the realistic representation of the atmospheric OXL concentrations is a prerequisite to properly simulate the atmospheric mineral Fe dissolution processes.

Figure 11 Annual mean Fe-containing aerosol solubility at deposition fluxes (%) as simulated for the EC-Earth simulation averaged for the period 2000–2014 for (a) mineral dust aerosols; (b) the sum of solid fuel combustion, liquid fuel combustion, and open biomass burning aerosols; and (c) the sum of all aerosol sources.​​​​​​​

Present-day simulations indicate that 61.816 ± 1.295 Tg yr⁻¹ of Fe in EC-Earth is deposited to the Earth's surface, which is towards the low end (40–140 Tg yr⁻¹) of the model intercomparison study by Myriokefalitakis et al. (2018). The amount of total Fe deposited to the global ocean is calculated to be 12.937 ± 0.308 Tg yr⁻¹ in EC-Earth, which is about 50 % lower than recent estimates by Hamilton et al. (2019), owing to the significantly larger (almost double) mineral dust emission flux in that study. However, the large variability in global models can be partly attributed to the different mineral dust size ranges considered in the models. Indeed, since most of the Fe mass is associated with coarse dust aerosols (∼ 91 % in this work), models that additionally account for super-coarse mineral dust emission sources (i.e., >10 µm in diameter) eventually calculate a higher TFe source and thus increased TFe global deposition rates. In our EC-Earth simulations, roughly 0.878 ± 0.015 Tg yr⁻¹ of DFe is calculated to be deposited globally (Table 2) in the range of estimates presented in Myriokefalitakis et al. (2018) (0.8 ± 0.2 Tg yr⁻¹). Focusing on the marine environment, about 40 % (0.376 ± 0.005 Tg yr⁻¹) of the simulated DFe is deposited into the global ocean, indicating that a large fraction of Fe atmospheric inputs to the global ocean results from the dissolution of atmospheric aerosols. Our results are close to the high-end of other global estimates (0.173–0.419 Tg yr⁻¹) as presented in the model intercomparison study of Myriokefalitakis et al. (2018), slightly higher than the respective DFe deposition fluxes in Ito et al. (2021) (0.271 Tg yr⁻¹) but somehow lower compared to Hamilton et al. (2019) estimates (roughly 0.5 Tg yr⁻¹), although the significant differences in the dust emission fluxes. Thus, even though we do not consider a super-coarse mode of dust in our simulations, the DFe deposition rates over the remote ocean are not severely impacted (Myriokefalitakis et al., 2018) by the size of the emitted minerals, but instead by the atmospheric processing during long-range transport. Nevertheless, reaching a firm conclusion in that respect will need further work.

The Fe-containing dust aerosols dominate (∼ 70 %) the total deposition fluxes over the ocean in the model, although combustion sources are calculated to have a significant impact on the Fe inputs to remote oceanic regions, such as the Pacific and the Southern Oceans, in agreement with other studies (e.g., Hamilton et al., 2020). The maximum DFe deposition fluxes occur in EC-Earth downwind of the main desert source regions, with high deposition rates being simulated in the outflow of tropical biomass burning regions (such as South America, Africa, and Indonesia), as well as over highly populated regions due to the Fe released from anthropogenic combustion processes in the presence of polluted air masses (such as in India and China). Overall, the average Fe solubility at the deposition of combustion aerosols is found here ∼ 19 % (Fig. 11a), much higher compared to the solubility of mineral dust aerosols (∼ 2 %; Fig. 11b), clearly indicating the importance of atmospheric processing on the potential bioavailable inputs to the global ocean. We further note that although a relatively high Fe solubility is applied here for oil fly ash emission (∼ 79 %), a sensitivity simulation (not shown) using a solubility for ship oil emissions of 47.5 %, as proposed by Rathod et al. (2020), leads to overall an only slight decrease (up to ∼ 2 %) in Fe solubility calculations, mainly in the northern Atlantic Ocean, and does not substantially affect our results. Indeed, our simulations show high Fe solubilities far from continental regions (Fig. S6g–i), such as the tropical Pacific and Atlantic Oceans (Fig. 11c), due to aerosol aging and lower Fe concentrations. Note that an evaluation of Fe solubility model calculation over oceanic regions (based on cruise measurements, where available) is also provided here (Fig. S5), showing overall a general overestimation of the observed values (Table A2) with nevertheless both positive and negative biases.

5 Summary and conclusions

This work documents the implementation of a detailed multiphase chemistry scheme in the EC-Earth3 Earth system model, aiming to provide consistent estimates of the atmospheric concentrations of the Fe-containing aerosols, along with the species that modulates its atmospheric processing, i.e., OXL and SO${}_{\mathrm{4}}^{\mathrm{2}-}$. For this, a comprehensive description of the atmospheric Fe cycle is included in the model, accounting for (1) an explicit soil mineralogy, (2) the contribution of combustion emissions, and (3) an atmospheric dissolution scheme that accounts for atmospheric acidity, ambient levels of OXL, and photoinduced processes. The multiphase chemistry scheme simulates the aqueous-phase processes of the troposphere for inorganic and organic compounds, along with the Fenton reaction. The KPP software is used in the model to integrate the aqueous phase and the dissolution equations, which adds flexibility to the code. Overall, simulations of tropospheric chemistry and aerosols for present-day conditions (2000–2014) have been realized, and budget calculations for OXL, SO${}_{\mathrm{4}}^{\mathrm{2}-}$, and the DFe-containing aerosols have been presented.

Model simulations have been performed both as a coupled system with IFS, as well as driven by offline meteorological fields from the ERA-Interim reanalysis. Budget analysis has shown that glyoxal is the main precursor of OXL in the atmosphere (∼ 74 %) and that the potential primary sources in the model (0.373 ± 0.005 Tg yr⁻¹) have a negligible impact on OXL concentrations. Non-traditional, laboratory-derived, aqueous-phase production pathways of OXL via glyoxal oxidation are also accounted for in our simulations. We have shown that when such pathways are omitted in the simulations, the calculated global OXL atmospheric concentrations are substantially lowered (∼ 43 %). For the ERA-Interim setup, the OXL net chemical production is calculated as 18.116 ± 0.071 Tg yr⁻¹, resulting in an atmospheric burden of 0.330 Tg on average. For the online-coupled system, however, the OXL net chemical production is ∼ 30 % lower, mainly attributed to a lower atmospheric OH abundance in EC-Earth due to biases that can be generally found in climate–chemistry models, such as those for temperature and humidity, especially at the higher altitudes. Overall, the simulated oxidizing capacity, along with the contribution of secondary sources other than the GLX oxidation, are shown to have a significant impact on the OXL atmospheric abundance. However, we acknowledge that other formation pathways of OXL, primary or secondary, may exist in the gas and aqueous phases of the atmosphere that could further contribute to OXL levels.

The dissolution of dust and combustion aerosols dominates on the dissolved Fe fraction in the model, calculated in EC-Earth at 0.806 ± 0.014 Tg yr⁻¹, in good agreement with the ERA-Interim simulation and well in the range of other model estimates. Furthermore, a broad evaluation of the EC-Earth proves the models' ability to represent AOD, particularly over regions that are dominated by anthropogenic aerosol but also over selected dusty sites, in line with previous EC-Earth evaluations. However, dust is underestimated over most remote areas, implying the EC-Earth shares the difficulties of many global models in representing the long-range transport of dust, especially the coarse particles. The model also underestimates the in situ OXL measurements, especially during winter, nevertheless indicating that additional sources (primary and/or secondary) are needed. Model comparisons with cruise measurements demonstrate a strong link between atmospheric DFe concentrations and atmospheric composition. The model seems to better capture the observations in the accumulation mode than in the coarse mode. This is attributed to differences either in the atmospheric processing between the accumulation and coarse particles (i.e., aerosol water content and acidity levels), the misrepresentation in the aerosol sizes (e.g., surface area/volume ratios), or to systematic errors in the mineralogical composition of the emitted Fe-containing soil particles. For this, several developments are planned by the EC-Earth consortium, aiming to improve the representation of the dust aerosols' size distribution, the description of the mineralogical composition, and the Fe-content in combustion sources, which are expected to reduce the existing uncertainties in the model for a more accurate simulation of the atmospheric Fe cycle.

Emphasizing on the biogeochemistry-related implications of this study, EC-Earth calculates a global annual present-day DFe deposition flux of 0.878 ± 0.015 Tg yr⁻¹, which is well within the range of estimates of other global modeling studies. About 40 % of the DFe deposition fluxes are calculated to occur over the global ocean (0.376 ± 0.005 Tg yr⁻¹), with a strong spatial and temporal variability. The highest annual mean DFe inputs to the global ocean are associated with aerosols of soil origin, especially downwind of the major dust source regions. In addition, Fe-containing combustion aerosols are calculated to have a significant contribution downwind of biomass burning source regions and highly populated areas in the Northern Hemisphere. It is further demonstrated that over the open ocean the Fe solubility at deposition for aerosols of combustion origin is about an order of magnitude higher than that of mineral dust origin, suggesting that the relative contributions of the primary sources can significantly affect bioavailable aerosol fraction and may thus play an important role in oceanic areas where the phytoplankton growth is limited by Fe supply, such as the Southern Ocean.

It is widely recognized that a combined approach considering both Fe atmospheric processing and deposition over oceans should ideally be used in Earth system models for the assessment of the impact of nutrient-containing aerosol deposition on marine productivity. A deeper understanding of the atmospheric Fe cycle is thus needed for a better description of the biogeochemistry implications in the presence of a changing climate. Such types of knowledge, however, should be obtained by extensive model evaluation with observations, especially over the remote regions of the world like the Southern Ocean, where currently the largest discrepancies between models and measurements exist. Therefore, a comprehensive calculation of the Fe physicochemical transformations is necessary to predict the strength of DFe inputs to the ocean, despite the complexity of the related atmospheric multiphase processes. The present study thus aims to complement the marine biogeochemistry component of EC-Earth using a fully coupled calculation scheme for atmospheric dissolved Fe fluxes into the global ocean. That new EC-Earth model version is expected to eventually allow for a better representation of the marine biogeochemistry perturbations in past and future climates and air quality.

Appendix A

Table A1 Summary of statistics for all points and per region (as depicted in Fig. 1) for the evaluation of (a) the modeled annual mean AOD at 550 nm against AERONET version 3 level 2.0 retrievals for all available stations covering the 2000–2014 period, (b) the same but for selected dust-dominated AERONET sites (characterized as described in Sect. 2.5), (c) the modeled annual mean dust surface concentration for 2000–2014 compared to climatological mean values from RSMAS sites and AMMA campaign, and (d) the modeled annual dust deposition flux averaged for the period 2000–2014 against observations as compiled in Albani et al. (2014) from several sources. The number of stations (n), the Pearson correlation coefficients (R) between the simulated and measured monthly mean concentrations, the normalized mean bias (nMB), and the normalized root-mean-square errors (nRMSEs) are indicated for the EC-Earth simulation.

Table A2 Summary of statistics for the evaluation of the simulated concentrations (as depicted in Fig. 6) of (a) oxalate and (b) sulfate. Statistics are also given for (c), (d), (e) the dissolved Fe-containing aerosols,; (f), (g), (h) the total Fe-containing aerosols; and (i), (j), (k) the derived aerosol solubility (SFe = %DFe/TFe) for (c), (f), (i) the accumulation mode; (d), (g), (j) the coarse mode; and (e), (h), (k) the total suspended particles (tsp), respectively. The number of stations (n), the Pearson correlation coefficients (R) between the simulated and the measured concentrations, the normalized mean bias (nMB) and the normalized root-mean-square errors (nRMSEs) are indicated for the EC-Earth and ERA-Interim simulation. The results for the sensitivity simulation EC-Earth(sens) for oxalate evaluation, as well as the EC-Earth(AerChem-AMIP) for the sulfate evaluation, are also shown.