The atmospheric composition component of the ICON modeling framework: ICON-ART version 2025.10

Accurate and efficient modeling of atmospheric composition, including aerosols and trace gases and their interactions with radiation, clouds, and dynamics is essential for improving predictions and understanding of air quality, weather, climate, and related health impacts. The Aerosols and Reactive Trace gases (ART) component extends the ICOsahedral Nonhydrostatic (ICON) modeling framework by enabling online, fully coupled simulations of atmospheric composition processes across scales. ICON-ART includes modules for emissions, transport, gas-phase chemistry, and aerosol microphysics in both the troposphere and stratosphere, allowing for the investigation of feedbacks between atmospheric composition and physical processes from the large-eddy to global scale.

This paper presents an updated overview of the ICON-ART framework as implemented in version 2025.10, highlighting recent developments in emission parameterizations, chemical mechanisms, aerosol processes, and coupling to the physical core of ICON via aerosol–radiation and aerosol–cloud interactions. We summarize the structure of the code infrastructure and demonstrate the model’s flexibility and scalability across a wide range of applications. ICON-ART provides a unified and modular platform for research and operational use in atmospheric composition, bridging the gap between regional air quality modeling and global Earth system simulations.

1 Introduction

Atmospheric composition focuses on the variations in and processes affecting trace gases and aerosols, which in turn influence air quality, weather, and climate. These components and their interactions with radiation and clouds are vital for processes such as cloud formation, radiative forcing, and precipitation patterns. Therefore, accurately simulating atmospheric composition is essential for improving predictions and understandings related to weather, renewable energy, climate change, air pollution, and associated health impacts.

The mutual feedback between the chemical and physical states of the atmosphere across different scales has driven the integration of atmospheric composition modeling into weather and climate models (Baklanov, 2010; Grell and Baklanov, 2011; Brasseur and Kumar, 2021). Beginning in the early 2010s, several global and regional-scale model systems have been developed to account for the interactions between atmospheric composition and the physical state of the atmosphere (Baklanov et al., 2014). This led to the development of online-coupled models ranging from global hydrostatic chemistry-climate models like ECHAM/HAMMOZ (Schultz et al., 2018), EMAC (Jöckel et al., 2006), and WACCM-X (Liu et al., 2018) to regional non-hydrostatic models, such as WRF-Chem (Grell et al., 2005) and COSMO-ART (Vogel et al., 2009). One particular focus of these developments has been on how aerosols, trace gases, and their interactions with radiation and clouds influence weather patterns, air quality, and climate. Recent advancements in modeling aerosol-radiation interaction (ARI) (Rieger et al., 2017; Hoshyaripour et al., 2019; Yang et al., 2020; Oh et al., 2024), aerosol-cloud interaction (ACI) (Glotfelty et al., 2019; Zhang et al., 2022; Seifert et al., 2023; Samanta et al., 2024), and the integration of chemistry into weather forecasting models (Kukkonen et al., 2012; Hodzic and Madronich, 2018; Deroubaix et al., 2024) have significantly improved our ability to predict these complex systems. However, only a few models can account for local, regional, and global weather and climate processes within a single modeling framework.

The ICON model has been developed and widely used for weather and climate prediction across scales. It solves the 3D non-hydrostatic and compressible Navier–Stokes equations on an icosahedral-triangular grid (Gassmann and Herzog, 2008), facilitating precise predictions across scales (Zängl et al., 2015; Heinze et al., 2017; Giorgetta et al., 2018). The ART module, integrated into the ICON framework, enables comprehensive modeling of atmospheric composition. It handles emissions, transport, and transformations of trace gases and aerosols, incorporating gas-phase chemistry and aerosol dynamics in the troposphere and stratosphere (Rieger et al., 2015; Weimer et al., 2017; Schröter et al., 2018). ICON-ART has been successfully used to investigate mutual feedbacks between the chemical and physical states of the atmosphere across different scales ranging from large-eddy simulations (Muth et al., 2025) to regional weather (Rieger et al., 2017; Seifert et al., 2023) and global climate (Weimer et al., 2021). While ICON can be coupled to external aerosol and chemistry models, ART represents the natively and fully integrated atmospheric composition component developed and maintained within the ICON modeling framework.

Previous works have described ICON-ART with respect to the basic equations, parameterizations and numerical methods (Rieger et al., 2015), tracer framework (Schröter et al., 2018) and chemistry processes (Weimer et al., 2017; Schröter et al., 2018). This paper provides an updated overview of the ICON-ART framework version 2025.10 (ICON Partnership, 2025), highlighting the recently developed components and features that enable its role in atmospheric composition modeling. We discuss the emission parameterizations in Sect. 2, followed by the description of chemistry and aerosol processes in Sects. 3 and 4, respectively. Then we explain interactions within the broader ICON modeling framework through aerosol-radiation and aerosol-cloud interactions. A summary of the code infrastructure is given in Sect. 6 followed by conclusions and outlook in Sect. 7. To facilitate readability, all acronyms employed in this paper are provided in the Appendix.

2 Emission processes

ART accounts for the emission of gases and particulate matter from both natural and anthropogenic sources, as listed in Table 1. The emission processes for wildfires, desert dust, sea salt, dimethyl sulfide (DMS), pollen, and the online emission module (OEM) are taken from the original implementations in the COSMO-ART framework. Other emission parameterizations have been developed and implemented within the ICON-ART framework. In both cases, significant modifications have been made for scientific or technical reasons, which are detailed below. For further information on the parameterizations, see the references in Table 1. Processes that are already well established in ICON-ART, such as desert dust and biogenic VOC emissions, remain unchanged in version 2025.10 and are therefore only briefly summarized in Table 1, with full details provided in the cited model description papers.

(Weimer et al., 2017)(Jähn et al., 2020)Weimer et al. (2017); Jähn et al. (2020)(Kaiser et al., 2012)(Freitas et al., 2007)Walter et al. (2016)(Folch et al., 2016)Bruckert et al. (2022)(Vogel et al., 2006)Rieger et al. (2017)(Monahan et al., 1986; Smith and Harrison, 1998; Mårtensson et al., 2003)(Grythe et al., 2014)Lundgren et al. (2013); Rieger et al. (2015)(Lana et al., 2011)(Guenther et al., 2012)Weimer et al. (2017)(Zink et al., 2013)Zink et al. (2013)Rieger et al. (2015)Rieger et al. (2015)

Table 1 Emissions in the ICON-ART model system including the main technical/scientific basis and references of the parameterization and the first published implementation in the ART framework. The new features presented in this work are shown in bold.

2.1 Online Emission module

The Online Emission Module (OEM) was first developed for COSMO-ART (Jähn et al., 2020) and then adapted to ICON-ART and refactored for improved computational performance. OEM enables efficient processing of emissions which can be represented by adding up individual source categories, with the emissions from each category being fixed in space but varying over time. This holds for most applications requiring input from anthropogenic emission inventories such as the Emission Database for Global Atmospheric Research (EDGAR) (Crippa et al., 2018) or the regional inventory CAMS-REG (Kuenen et al., 2022). These inventories provide gridded emissions divided into individual source categories. An example for such a source classification is the Gridded Nomenclature For Reporting (GNFR), which distinguishes 12 source categories from public power to agricultural emissions (Super et al., 2020). The main steps for integrating these inventories into the model are (i) re-mapping to the ICON grid, (ii) application of temporal scaling factors, (iii) vertical distribution of the emissions onto ICON's terrain-following hybrid sigma-z coordinates, and (iv) summing up over all source categories. In most model systems including WRF-Chem, CHIMERE or CAMx, (e.g., Menut et al., 2024; Woo et al., 2012), these steps are performed externally using a pre-processing software such as HERMES (Guevara et al., 2019), which generates a large number of files each containing the emission field representative of a given time interval (e.g. hourly), which are then read into the model during runtime. OEM, instead, allows reading the emissions together with all temporal and vertical scaling information only once during model initialization. The temporal and vertical scalings are then applied online during the simulation. This greatly reduces the need for data pre-processing and simplifies the setup of new simulations. As shown in Jähn et al. (2020), the additional time required for the online computations is fully compensated by the time saved through less frequent file access.

Inputs for OEM can be produced with the Python package emiproc (Constantin et al., 2025), which is able to process multiple emission inventories to prepare the inputs for a range of atmospheric transport models including ICON-ART. The tool maps the emissions in a mass-conserving way onto the ICON grid and creates all additional inputs including source and country-specific temporal and vertical profiles and the corresponding country masks. The tool is also able to merge multiple inventories, for example to embed a high-resolution national inventory into a coarser global inventory.

To compute the emission of a species X at time t in a given grid cell, OEM performs the following operation

$$\begin{array}{}\text{(1)}& E_X(z,t)=\sum _{s=\mathrm{0}}^{N_s}{E}_{X,s}\cdot w_{X,s}\left(t\right)\cdot v_{X,s}\left(z\right),\end{array}$$

with E_X,s the emission of species X from source category s in that grid cell from the inventory, w_X,s(t) the temporal scaling factor at time t for source category s, and v_X,s(z) the vertical scaling factor for source category s at vertical level z. N_s is the total number of categories. The temporal factor w_X,s(t) is computed as the product of three different scaling factors describing diurnal, day-of-week, and seasonal variability

$$\begin{array}{}\text{(2)}& w_{X,s}\left(t\right)=w_{X,s,h}\left(h\right(t\left)\right)\cdot w_{X,s,d}\left(d\right(t\left)\right)\cdot w_{X,s,m}\left(m\right(t\left)\right)\end{array}$$

with h(t) being the hour of the day, d(t) the day of the week, and m(t) the month of the year. Alternatively, a separate scaling factor can be defined for each hour of the year to represent, for example, heating emissions varying with outdoor temperatures. Furthermore, different temporal scaling factors can be provided for different countries (or regions) together with a country mask. Diurnal factors $w_{X,s,h}\left(h\right(t\left)\right)$ are computed with respect to local time.

In addition to temporal and vertical factors, it may be necessary to provide speciation factors, which describe the fractional contribution of individual model species X to the total emissions of a family of species $\stackrel{\mathrm{̃}}{X}$ in the inventory. Examples are non-methane volatile organic compounds (NMVOC) and nitrogen oxides (NO_x), for which inventories typically provide only the total emissions of the family but not of the individual compounds. Speciation factors cannot be supplied as input to OEM, but are dealt with by emiproc. Based on an emission field for the family $\stackrel{\mathrm{̃}}{X}$ and a set of speciation factors, emiproc generates and writes out emission fields for all model species X. XML tags and namelist settings for OEM are described in Appendix E.

2.2 Wildfires

The ICON-ART model system incorporates a one-dimensional, sub-grid-scale plume-rise model developed by Freitas et al. (2006, 2007, 2010), similar to its implementation in COSMO-ART by Walter et al. (2016). This plume-rise model is suited for applications with horizontal resolutions on the order of 10 to 100 km. For these applications, the model calculates plume height based on buoyancy, atmospheric stratification, and flow conditions, accounting for processes that occur on scales much smaller than the horizontal spacing of the ICON. The model uses an internal vertical grid spacing of 100 m with 200 vertical layers. Environmental conditions (pressure, humidity, temperature, wind speed) are provided by ICON and transferred to the plume-rise model as initial and environmental conditions for each active fire grid point to determine plume height.

Fire size and intensity, based on the ICON land use class, vegetation type and density, determine heat release and initial buoyancy. The lower boundary condition assumes a virtual buoyancy source below the surface, resulting in high vertical velocity. Final buoyancy is limited by turbulent and dynamic entrainment, with additional buoyancy from latent heat release during condensation. The plume top is defined where vertical velocity inside the plume drops below 1 m s⁻¹ after an equilibrium state between the surroundings and the heat source is achieved. Heat flux values, dependent on vegetation type, are taken from Freitas et al. (2006). They were corrected in comparison to the implementation in COSMO-ART (Walter et al., 2016) such that the plume-rise model exactly reproduces the results in Freitas et al. (2010).

In addition to the original numerical solver introduced by Freitas et al. (2007), a more efficient and stable first-order implicit solver for the same equations was developed and implemented. The solver relies on a Godunov type scheme to solve the inviscid Burgers equation for the vertical velocity inside the plume and on upwind schemes for the other modeled plume variables such as temperature, specific humidity, specific cloud water, specific rain water, specific ice, horizontal entrainment velocity, and plume radius. An internal grid is not needed for this solver since it copes with unevenly spaced grids like the vertical atmosphere discretization in numerical weather forecast models such as ICON while at the same time allowing much longer plume internal time-steps. The solver comes with an a priori condition that determines whether a height calculation is needed or if the plume stays at the minimal emission height. The plume heights strongly vary with the time of day since a diurnal cycle function d is applied to fire intensity and size. The function proposed by Kaiser et al. (2009); Andela et al. (2015) and applied by Walter et al. (2016) in COSMO-ART is also used in ICON-ART. The plume-rise model returns plume bottom and top heights, with a parabolic emission profile f describing the vertical distribution of emissions. The emission rate E of a species in $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ is calculated for each grid cell based on height z and time t.

$$\begin{array}{}\text{(3)}& E_X(z,t)=M_X\left(t\right)\cdot d\left(t\right)\cdot f(z,t)\end{array}$$

M_X is the daily mean emission flux of species X from CAMS GFAS (Kaiser et al., 2012), d is the diurnal cycle, and f is the parabolic emission profile between upper and lower injection heights.

CAMS GFAS relies on Fire Radiative Power (FRP) from the NASA MOD14 product, which includes thermal radiation observations from the MODIS instrument (Kaiser et al., 2012). To address data gaps due to, e.g. cloud cover, fire data is assimilated using a Kalman filter and statistics. FRP density is updated based on previous and current observations, with sampling limited to four times per day to represent the diurnal fire cycle (Kaiser et al., 2012). GFAS data is provided at a resolution of 0.1°.

2.3 Volcanic eruptions

The rise of volcanic plumes during eruptions depends on both volcanic and atmospheric conditions. Volcanic conditions are the exit temperature, exit velocity, exit volatile fraction, and the vent diameter. The exit velocity and the vent diameter control the mass eruption rate (MER). The height of the plume depends on volcanic conditions due to effects on the MER as well as on the plume density and atmospheric conditions. Due to the complexity of plume dynamics, simple relationships (e.g., Mastin et al., 2009), which only depend on plume height and are often used in dispersion models, can lead to large uncertainties in emissions (e.g., Marti et al., 2017; Bruckert et al., 2022).

Volcanic emissions in ICON-ART are calculated online using the 1-D volcanic plume model FPlume by Folch et al. (2016), which considers the volcanic conditions as well as processes during the plume rise such as ambient air entrainment, plume bending due to wind, particle wet aggregation, energy supply due to water phase changes, and particle fallout and re-entrainment. Bruckert et al. (2022) described the coupling of FPlume with ICON-ART in detail. In short, FPlume requires atmospheric profiles for temperature, pressure, density, zonal and meridional wind speed, and specific humidity at the volcanic vent. In addition to meteorological data, FPlume needs estimates of the exit temperature, exit velocity, and exit volatile fraction, which depend on the type and setting of the volcano. FPlume can either calculate the MER based on a given height or the height based on a given MER. During every time step when the volcano is active, FPlume first calculates the plume properties, i.e., the total MER in the case of a given plume height or plume height in the case of a given MER. Second, the fraction of very fine ash, which is relevant for long-range transport, is determined based on the plume height and the total MER by using the relationship of Gouhier et al. (2019). Third, very fine ash is emitted along a profile that has initially been defined by Suzuki (1983) and applied by Marti et al. (2017) for the coupling of FPlume to NMMB-MONARCH-ASH transport model (Nonhydrostatic Multiscale Model on the B-grid – Multiscale Online Nonhydrostatic AtmospheRe CHemistry model – ASH). The distribution of the very fine ash mass into the ICON-ART modes is prescribed in the FPlume input file.

SO₂ can be emitted alongside with ash using the same timing and profile, however, the MER of SO₂ is prescribed in the FPlume input file, as FPlume does not differentiate between gaseous compounds and water vapor. Bruckert et al. (2025) extended the coupling of ICON-ART and FPlume by water vapor emission. Here, the MER of water vapor is calculated from the exit water mass fraction and the total MER and is emitted through the same profile as ash. This simplification neglects the entrainment of water vapor and was so far only tested for the water-rich 2022 Hunga eruption (Bruckert et al., 2025).

The advantages of the coupling are (1) a more accurate MER and therefore better agreement with observed mass column loadings as shown for the 2019 Raikoke eruption, Kuril Islands (Bruckert et al., 2022) and (2) an easy consideration of eruption phases, which allows a comparison to observations in the near-field of complex volcanic emissions such as the 2021 La Soufrière eruption, St. Vincent (Bruckert et al., 2023).

2.4 Sea salt

For the emission of sea salt, ICON-ART offers two options. The first option is based on the parameterizations of Monahan et al. (1986), Mårtensson et al. (2003), and Smith and Harrison (1998) (hereafter referred to as MMS scheme) as described by Rieger et al. (2015). Another option is based on Grythe et al. (2014) (hereafter referred to as G14 scheme) which is described below.

Numerous studies have shown that sea surface temperature (SST) significantly influences the emission rates of sea salt by affecting the physical properties of foam and droplet formation at the ocean surface. To develop a globally applicable source function that can also represent realistic sea salt concentrations in tropical regions, incorporating this temperature dependence is essential (Grythe et al., 2014). To address this, we implemented an alternative scheme for sea salt emission in ICON-ART based on Grythe et al. (2014), which in addition to wind speed includes SST as a key parameter:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}F(D_{\mathrm{p}},U_{\mathrm{10}},T)}{\mathrm{d}{D}_{\mathrm{p}}}}=\\ & T_w\cdot (\mathrm{235}\cdot U_{\mathrm{10}}^{\mathrm{3.5}}\exp \left(-\mathrm{0.55}{\left(\ln \left({\displaystyle \frac{D_{\mathrm{p}}}{\mathrm{0.1}}}\right)\right)}^{\mathrm{2}}\right)\\ & +\mathrm{0.2}\cdot U_{\mathrm{10}}^{\mathrm{3.5}}\exp \left(-\mathrm{1.5}{\left(\ln \left({\displaystyle \frac{D_{\mathrm{p}}}{\mathrm{3}}}\right)\right)}^{\mathrm{2}}\right)\\ & +\mathrm{6.8}\times {\mathrm{10}}^{-\mathrm{3}}\cdot U_{\mathrm{10}}^{\mathrm{3}}\exp \left(-\mathrm{1}{\left(\ln \left({\displaystyle \frac{D_{\mathrm{p}}}{\mathrm{30}}}\right)\right)}^{\mathrm{2}}\right)).\end{array}\end{array}$$

In this equation, $F(D_{\mathrm{p}},U_{\mathrm{10}},T)$ is the sea salt emission flux in $\mathrm{particles}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, D_p is the dry particle diameter in µm, U₁₀ is the 10 m wind speed in m s⁻¹, T is the SST in K, and T_w is an empirical temperature correction factor (dimensionless), accounting for the effect of SST on sea salt emissions. This parameterization is applied to three sea salt modes with median diameters of 0.1, 3.0 and 30 µm, with standard deviations of 1.9, 2.0 and 1.7, respectively (Grythe et al., 2014). This enhancement enables a more accurate representation of the regional and seasonal variability of sea salt concentrations, especially in tropical and subtropical ocean regions where conventional parameterizations often underestimate emissions. The new parameterization thus represents an important step towards a more physically consistent modeling of marine aerosol sources on a global scale. We note that the factor 10⁻³ in the last term of the equation is mistakenly missing in the original publication by Grythe et al. (2014) and is corrected in Eq. (4).

MMS and G14 sea-salt emission schemes differ not only in their whitecap formulations but also in their treatment of particle-size distributions and SST-dependent scaling (Grythe et al., 2014; Barthel et al., 2019; Li et al., 2024). Barthel et al. (2019) demonstrated that SST corrections can substantially reduce coarse-mode concentrations and may even have a larger impact than switching between source functions. They also found the strongest divergences for particles larger than PM_2.5, with SST effects further amplifying these differences. These insights highlight that the structural contrasts between MMS and G14 schemes, particularly the inclusion of SST dependence and the size-resolved flux formulation, can significantly influence emitted mass. While a quantitative evaluation is beyond the scope of this study, this context helps to clarify the expected behavior of the new G14 implementation.

2.5 DMS

DMS has multiple sources with one large contribution from oceanic biogenic emissions (Lana et al., 2011). As such, these emissions are highly dependent on the exchange between ocean and atmosphere and hence on wind speed. Therefore, an online calculation is needed to account for DMS emissions from the ocean.

Lana et al. (2011) provided a 1° × 1° monthly climatology of DMS ocean surface concentrations (C_w) based on measured data. These climatological values are converted in ICON-ART to the DMS emission flux F_DMS using the following parameterization (Lundgren et al., 2013; Ullwer, 2018):

$$\begin{array}{}\text{(5)}& F_{\mathrm{DMS}}={\displaystyle \frac{C_{\mathrm{w}}\cdot M_{\mathrm{DMS}}}{\mathrm{3.6}\times {\mathrm{10}}^{\mathrm{2}}}}\left(\mathrm{0.222}{u}_{\mathrm{10}\mathrm{m}}^{\mathrm{2}}+\mathrm{0.333}{u}_{\mathrm{10}\mathrm{m}}\right).\end{array}$$

In this equation, M_DMS=6.21 g mol⁻¹ is the DMS molar weight, C_w has the unit mol L⁻¹ and u_10 m is the wind speed at the altitude of 10 m above sea level in m s⁻¹. The emission flux is added to the DMS tracer mixing ratio using the same procedure as for other emissions (Weimer et al., 2017).

2.6 Pollen

The pollen emission model is based on the EMPOL approach (Zink et al., 2013). The basic idea is that the emission process can be divided into two steps: (1) opening of the anthers and the accumulation of pollen in a reservoir; (2) release of pollen from the reservoir into the air, driven by turbulence.

The opening of the anthers is controlled by phenology, temperature and humidity. The phenology is modeled using a temperature sum approach, as described by Pauling et al. (2014). In general, warm and dry conditions favor opening of the anthers. To calculate the amount of pollen released into the reservoir, a plant distribution map is required. When turbulence is sufficiently strong, the pollen reservoir is emptied and the pollen becomes airborne. This process is parameterized using the Turbulent Kinetic Energy (TKE). Once airborne, pollen is treated as a passive tracer and can be removed by dry and wet deposition. Re-suspension of pollen is not considered by the model.

The pollen emission model was originally designed for COSMO-ART and later implemented in ICON-ART. Currently, implemented plant species include hazel, alder, birch, grasses and ambrosia. Each species is handled slightly differently within the model. Recent developments include an emission implementation for hazel and the integration of real-time pollen data. The approach is described by Adamov and Pauling (2023). The real-time pollen data is used in two ways. First, the start of the flowering season is set to the date when the first pollen are observed. Second, the emission flux is scaled so that the modeled concentrations match the observed values.

3 Chemistry processes

In ART different types of tracers and solvers for them can be chosen. Passive tracers are inert and are only changed by emission and transport processes. Region tracers are a specific subset of passive tracers that originate from predefined regions. Chemtracers are tracers that have a simplified solver scheme, e.g. Linoz for ozone. The most complex tracers in ICON are those participating in chemical reactions described by a coupled system of Ordinary Differential Equations (ODE). These tracers are called meccatracers, because they are solved by the atmospheric chemistry module MECCA as described in Sect. 3.2. An overview is given in Table 2.

Vogel et al. (2015)Rieger et al. (2015)Rieger et al. (2015)Rieger et al. (2015)Rieger et al. (2015)Rieger et al. (2015)Mclinden et al. (2000)Schröter et al. (2018)Hsu and Prather (2010)Ramezani Ziarani et al. (2025)Diekmann (2021)Diekmann (2021)Funke et al. (2016); Matthes et al. (2017)Ramezani Ziarani et al. (2025)Weimer et al. (2017)Weimer et al. (2017)Sander et al. (2019)Schröter et al. (2018)

Table 2 Types of chemistry in ART including reference to their main technical/scientific descriptions and the first published implementation in the ICON-ART framework.

3.1 Simplified chemistry options

3.1.1 Region tracers

Understanding the origin and subsequent pathways of air masses is fundamental to interpreting atmospheric composition and evaluating the transport characteristics of atmospheric models. To facilitate these investigations within the ICON-ART modeling system, a suite of passive “region tracers” has been incorporated.

These tracers serve as computationally efficient diagnostic tools, each representing a distinct geographical source region. The defined regions encompass a variety of scales and types, including continental regions (e.g., Europe, North America, East Asia), oceanic basins (e.g., Tropical Pacific, Tropical Atlantic) and broad hemispheric backgrounds. Besides region tracers described in Vogel et al. (2015), additional tracers used in the PHILEAS campaign and ASCCI campaign are implemented. Those tracers are set to 1 inside and 0 outside their source region at the lowest level within ICON-ART.

Conceptually, air originating within the lowest model layer of a specific source region is “tagged” with its corresponding tracer. Those tracers are then transported throughout the model domain solely by the simulated atmospheric dynamics. These tracers are inert; they do not undergo chemical transformation or deposition processes. As a result, the concentration of a specific region tracer at any location and time within the simulation directly quantifies the fractional contribution of air that originated from that source region. In contrast to Vogel et al. (2015) in ICON-ART a land-sea mask is added to the tracers.

Figure 1 Left: source regions for region tracers used during the ASCCI airplane measurement campaign; Right: Distribution of tracer originating from Europe in 300 hPa, as well as the charted flight path of Flight 09 of the ASCCI campaign.

An example usage of those tracers used during the ASCCI campaign is shown in Fig. 1. On the left, the source regions of the tracers are displayed. On the right, the distribution at 300 hPa of the region tracer originating from Europe is shown for 23 March 2025. The simulation was started on 1 December 2024 and reinitialized each day with the current meteorological conditions provided by DWD. In the figure a distinct narrow filament can be seen over the north-west of Scandinavia. The blue lines with red dots depict the corresponding flight path of the research aircraft HALO, where in-situ and remote sampling of the structure has been conducted.

3.1.2 Lifetime tracers

Lifetime tracers in ICON-ART should be used when a very fast chemistry, but not the highest accuracy is needed. Some tracers in ICON-ART get a special treatment when they have specific names and when ”lifetime“ is chosen as solver. Tracers with lifetimes changing vertically as a function of pressure or overlaying O₂ column are those presented in Fig. 2. The corresponding formula are summarized in Appendix G. For other tracers with “lifetime” chosen as solver, the globally constant lifetime given in an XML file is used.

Figure 2 Lifetimes from left to right for the tracers TRCH4, TRCO, TRH2O and TRCO2.

Combined effects on TRCO2 compared to a tracer named CO₂ of the lifetime approach and CO₂ deposition in the ocean can be found Fig. 7.

3.1.3 Linoz v2: simple atmospheric ozone scheme

One way to calculate ozone concentrations in ICON-ART is the implementation of the LINearized OZone scheme as described by Mclinden et al. (2000). It is a fast online calculation of ozone as a chemical tracer using a linear approximation for the ozone change depending on temperature, overhead ozone column, local mixing ratio, and optionally, as described in more detail below, polar depletion. An in-depth explanation of the implementation in ICON-ART can be found in Schröter et al. (2018). Since the publishing of the aforementioned paper, several bug fixes have been applied to LINOZ as well as an adjustment to include the near-surface relaxation to a globally constant ozone volume mixing ratio of 25 ppb, as proposed already by Mclinden et al. (2000). Comparing the same setup (40 km horizontal resolution, 90 height levels, 48 h lead time, 10 test runs) with and without LINOZ shows an increase in computational runtime for the entire model run of up to 7 %.

LINOZ provides the possibility to include an additional parametrization for polar ozone loss due to catalytic ozone depletion by chlorine and bromine. This parametrization is activated in model cells that are located in polar regions ($\left|\mathrm{lat}\right|>\mathrm{45}\mathit{°}$), have a temperature below 195 K and a solar zenith angle smaller than a fixed value of 85° (cold tracer parametrization) or 90° (lifetime parametrization). The cold tracer parametrization contains an additional tracer that emulates the prolonged activation time of these species after their initial activation. Figure 3 shows ozone columns when using LINOZ with polar loss parametrization and coldtracer activated. In addition to that, the chlorine and bromine loading used in the parametrization of LINOZ has been adapted according to Hossaini et al. (2019). Ozone is initalized with CAMS EAC4 data (Inness et al., 2019) and then runs according to the LINOZ calculation for four months with meteorological reinitialization every 24 h.

Figure 3 Comparison of Ozone Columns Calculated by LINOZ (top) and Satellite Data (bottom). The top row shows a model run that was initialized daily (start date 11 February 2022, 40 km horizontal resolution, 90 height levels) with meteorological reanalysis data by DWD, while ozone was initialized from CAMS EAC4 data (Inness et al., 2019) on 11 February 2022 and then passed on. The bottom row shows satellite data provided by NASA Ozone Watch (NASA, 2025).

3.1.4 Linoz v3 and UBCNOy: Solar forcing of stratospheric ozone

To include variable solar forcing by energetic electron precipitation (EEP) and spectral solar irradiance (SSI) via stratospheric ozone into ICON-ART, the following adaptations were made:

• UBCNOy. An upper boundary condition of NO_y (NO, NO₂, NO₃, 2 N₂O₅, HNO₃, HNO₄, ClNO₃) was implemented at three model levels below the model top boundary. UBCNOy is based on a semi-empirical model of the auroral and magnetospheric electron precipitation into the mesosphere and lower thermosphere using the geomagnetic Ap-Index as prognostic variable and constrained by MIPAS satellite observations (Funke et al., 2016) also part of the solar forcing recommendations for chemistry-climate models, e.g., for CMIP6 and CMIP7 (Matthes et al., 2017; Funke et al., 2024). The upper boundary NO_y was added to the stratospheric NO_y background derived from SimNOy (see Sect. 3.1.5). The SimNOy tables were extended to the mesopause by using output from the EMAC model to account for the mesospheric lifetime of NO_y.

• Linoz v3. To account for the impact of NO_y on stratospheric ozone, Linoz version 3 (Hsu and Prather, 2010) was implemented into ART, using the linearized terms of temperature, NO_y, and ozone column.

• Solar Spectral Irradiance variability was incorporated by generating two distinct sets of LINOZ coefficients corresponding to solar maximum and solar minimum conditions. The model interpolates between these coefficient sets based on the daily F10.7 index (solar radio flux at 10.7 cm).

A comparison of NO_y and ozone from model experiments with and without UBCNOy and with constant solar maximum respectively solar minimum SSI is provided in Fig. 4, highlighting the impact of EEP-NO_y on NO_y, and the combined impact of EPP-NO_y and variable SSI on ozone. A more detailed description and evaluation against observations are given in Ramezani Ziarani et al. (2025).

Figure 4 Left: NO_y averaged over 8 years of model time. Colors: model experiment with variable EPP-NO_y for 2002–2009 and constant solar maximum SSI. Lines: no EPP-NO_y and constant solar minimum SSI. Right: percentage difference of ozone for the same period and model experiments, highlighting enhanced ozone formation around the tropical stratopause for solar maximum conditions, as well as strong ozone loss in polar regions due to the EPP-NO_y. Model experiments following Ramezani Ziarani et al. (2025) using ICON-ART version 2025.10

3.1.5 Stratospheric N₂O-NO_y scheme

Tropospheric N₂O is the main source of reactive nitrogen compounds in the stratosphere, which are summarised as NO_y (Seinfeld and Pandis, 2016): NO_y = NO + NO₂ + NO₃ + HNO₃ + HNO₄ + 2 N₂O₅ + ClONO₂ +…

N₂O is chemically inert in the troposphere and very long-lived due to its estimated lifetime of 120 years (Prather et al., 2015), so that it can be transported into the stratosphere. 90 % of the stratospheric N₂O is destroyed via photolysis:

$$\begin{array}{}\text{(R1)}& {\mathrm{N}}_{\mathrm{2}}\mathrm{O}+hv\to {\mathrm{N}}_{\mathrm{2}}+{\mathrm{O}}^{\mathrm{1}}\left(\mathrm{D}\right)\end{array}$$

The remaining N₂O molecules react with excited oxygen atoms O¹(D), which originate from ozone photolysis:

$$\begin{array}{}\text{(R2)}& {\displaystyle }{\mathrm{N}}_{\mathrm{2}}\mathrm{O}+{\mathrm{O}}^{\mathrm{1}}\left(\mathrm{D}\right)\to {\mathrm{N}}_{\mathrm{2}}+{\mathrm{O}}_{\mathrm{2}}\text{(R3)}& {\displaystyle }{\mathrm{N}}_{\mathrm{2}}\mathrm{O}+{\mathrm{O}}^{\mathrm{1}}\left(\mathrm{D}\right)\to \mathrm{NO}+\mathrm{NO}\end{array}$$

Reaction (R3) leads to the production of two NO molecules, which initiate the NO_y cycle through photolysis and oxidation:

$$\begin{array}{}\text{(R4)}& {\displaystyle }\mathrm{NO}+hv\to \mathrm{N}+{\mathrm{O}}^{\mathrm{3}}\left(\mathrm{P}\right)\text{(R5)}& {\displaystyle }\mathrm{N}+{\mathrm{O}}_{\mathrm{2}}\to \mathrm{NO}+{\mathrm{O}}^{\mathrm{3}}\left(\mathrm{P}\right)\end{array}$$

NO_y is destroyed via

$$\begin{array}{}\text{(R6)}& {\displaystyle }\mathrm{N}+\mathrm{NO}\to {\mathrm{N}}_{\mathrm{2}}+{\mathrm{O}}^{\mathrm{3}}\left(\mathrm{P}\right)\text{(R7)}& {\displaystyle }\mathrm{N}+{\mathrm{NO}}_{\mathrm{2}}\to {\mathrm{N}}_{\mathrm{2}}\mathrm{O}+{\mathrm{O}}^{\mathrm{3}}\left(\mathrm{P}\right)\end{array}$$

Since the production and destruction of NO according to the Reactions (R3) and (R5)–(R6) are in a first approximation the only chemical sources and sinks for NO_y, the simulation of NO_y with N₂O as a source takes place via the Reactions (R1)–(R7) (Olsen et al., 2001).

The N₂O-NO_y scheme of Olsen et al. (2001) is based on five parameters, which were determined in a photochemical box model at 20 pressure altitudes between 14–52 km, 18 latitudes between 85° S–85° N and for 12 months and are available as parameter tables:

• C1: 24 h – mean of N₂O loss frequency

• C2: 24 h – mean of NO photolysis frequency per NO_y molecule

• C3: proportion of NO formation during N₂O degradation

• C4: proportion of N reacting with NO or NO₂

• C5: N₂O production rate

In order to interpolate the coefficients onto the ICON grid, the nearest neighbour method is used for horizontal interpolation and linear interpolation is used for vertical interpolation based on geometric height. As the coefficients are only defined between the heights of 14–52 km, the corresponding boundary values at 14 and 52 km are used for the layers below and above.

3.2 Detailed chemistry mechanisms

ICON-ART supports comprehensive and scalable gas-phase chemistry simulations through the integration of the atmospheric chemistry module MECCA (Module Efficiently Calculating the Chemistry of the Atmosphere), which is part of the CAABA (Chemistry As A Boxmodel Application) framework (Sander et al., 2019). MECCA provides a flexible and extensible platform for representing detailed chemical processes in the troposphere and stratosphere, including oxidation pathways, radical chemistry, and heterogeneous reactions.

The numerical integration of the chemical system is performed using the Kinetic Pre-Processor (KPP) (Sandu and Sander, 2006), which generates optimized code for solving systems of ordinary differential equations representing chemical kinetics. ICON-ART uses the Rosenbrock solver (Sandu et al., 1997) within KPP for efficient and stable time integration, particularly suitable for stiff chemical systems.

Photolysis rates, a critical component of atmospheric chemistry, are computed using CloudJ (Prather, 2015), an advanced photolysis scheme that accounts for the effects of clouds, aerosols, and molecular absorption in a multi-wavelength approach. CloudJ can be used in offline or online configurations, depending on the application and computational cost considerations.

A key feature of the ICON-ART tracer framework is the use of MECCA as an external preprocessor, which enables users to define and compile custom chemical mechanisms tailored to specific scientific questions or case studies. This modular approach allows for the integration of predefined standard mechanisms as well as user-defined chemical schemes. The integration of MECCA into ICON-ART has been previously described in detail by Schröter et al. (2018). In this section, we summarize the available chemistry mechanisms currently implemented in ICON-ART and provide guidance on their configuration and application domains.

3.2.1 MOZART-4 Chemistry

In CAABA/MECCA version 4.0 (available in the supplementary material of Sander et al., 2019), the MOZART-4 chemical mechanism (Model for Ozone and Related chemical Tracers, version 4) from Emmons et al. (2010) has been integrated. Figure 5 presents a schematic overview of the steps required to run an ICON-ART simulation with MOZART-4 chemistry. This includes processing the chemical reaction mechanism in MECCA and generating Fortran90 code via KPP for numerical integration. The implementation of MOZART-4 into ART utilizes additional scripts provided in the Supplement. A detailed user guide for these preprocessing steps is available in Appendix D.

Figure 5 Schematic of the processing steps for integrating a custom chemistry mechanism into ART.

3.2.2 MOZART-T1 Chemistry

In Emmons et al. (2020), a new MOZART tropospheric chemistry scheme (MOZART-T1) was introduced, incorporating several improvements over the previous version (MOZART-4), including improved oxidation of isoprene and terpenes, refined organic nitrate speciation, and more accurate aromatic speciation and oxidation, among others. As a result, MOZART-T1 provides an improved representation of ozone in the troposphere. For more details on how to implement the MOZART-T1 chemistry scheme in ART, the reader is referred to Appendix D.

Figure 6 presents ground-level ozone pollution from a full chemistry simulation using the MOZART-T1 mechanism. The figure shows the modeled mean afternoon ozone mixing ratios for winter and summer 2019, alongside observed values from the EMEP monitoring network (European Monitoring and Evaluation Programme, data retrieved from Tørseth et al., 2012). The ICON-ART simulations were conducted in limited-area mode over Europe using a R3B07 grid (Δx≈ 13.2 km). Initial and boundary conditions for trace gases and aerosols were derived from CAM-Chem output (Buchholz et al., 2024). Anthropogenic emissions were generated from the TNO CAMS-REG-v4 air pollutant inventory (Kuenen et al., 2022).

Figure 6 demonstrates that the model accurately reproduces seasonal variations in ozone levels and effectively captures spatial variability. A comparison of hourly afternoon ozone concentrations with data from EMEP monitoring stations shows a RMSE of 12.4 ppb and a mean bias of −0.27 ppb for winter, and 11.85 and 1.92 ppb, respectively, for summer. The corresponding correlation coefficients are r=0.57 (winter) and r=0.62 (summer). These results are in good agreement with findings from other studies employing comparable modeling setups, for instance Mar et al. (2016), who used the WRF-Chem model with the MOZART-4 chemical mechanism and compared simulated hourly O₃ concentrations with the same EMEP monitoring network. Similar to our findings, their study also reported a slight overestimation of O₃ concentrations during summer.

Figure 6 Mean afternoon ground-level O₃ mixing ratios from ICON-ART simulations for winter (JF, left) and summer (JJA, right) 2019. Filled circles indicate observations from EMEP monitoring stations. Elevated sites and stations with less than 75 % valid data were excluded.

A full-chemistry simulation with MOZART-T1 increases the total runtime by roughly a factor of 10 compared to an ICON simulation without ART (tested on an HPC system with AMD Rome nodes with two AMD Epyc 7742 64-core CPU sockets each), reflecting not only the computational cost of the chemical mechanism but also the additional overhead from tracer transport, emissions, deposition processes, and model output.

3.3 Scheme for polar stratospheric clouds

Polar stratospheric clouds (PSCs) are an essential part for polar ozone depletion (e.g., Solomon, 1999), and therefore have to be accounted for in atmopsheric chemistry models. In ICON-ART, the three major known types of PSCs are implemented. Nitric acid trihydrate (NAT) particles are calculated based on the scheme introduced by van den Broek et al. (2004) using a fixed measured size distribution where each size bin is transported as a passive tracer. Size changes are based on the saturation conditions for NAT particles, based on Hanson and Mauersberger (1988). The calculation of supercooled ternary solution (STS) particles uses the diagnostic scheme by Carslaw et al. (1995) with an adapted volume concentration (Hervig and Deshler, 1998). For ice particles, it has been found that lifting the altitude where the ICON microphysics are calculated to the lower stratosphere leads to realistic formation and transport of these particles (Weimer, 2019). Further details can be found in Weimer et al. (2021) where the scheme has also been applied to local grid refinements around the Antarctica.

3.4 Removal processes

3.4.1 Parameterization of CO₂ deposition in the ocean

CO₂ is the major driver of anthropogenic climate change (Solomon et al., 2009), so that it should be accounted for in atmospheric chemistry simulation. One large sink of CO₂ is the deposition in the ocean which is implemented in ICON-ART in a simplified way. Jacob (1999) conclude that 50 % of the emitted CO₂ remains in the atmosphere, which we incorporated in ICON-ART as a possibility for deposition in the ocean. We assume an average emission of $E_{{\mathrm{CO}}_{\mathrm{2}}}=\mathrm{2.77}\times {\mathrm{10}}^{-\mathrm{9}}$ kg m⁻² s⁻¹, divide this value by the fraction of ocean surface on Earth r_ocean=0.707 and multiply it by the deposition factor r_depo=0.5 based on Jacob (1999) to get the effective deposition $f_{\mathrm{eff},{\mathrm{CO}}_{\mathrm{2}}}$ of CO₂ in the ocean:

$$\begin{array}{}\text{(6)}& f_{\mathrm{eff},{\mathrm{CO}}_{\mathrm{2}}}={\displaystyle \frac{E_{{\mathrm{CO}}_{\mathrm{2}}}{r}_{\mathrm{depo}}}{r_{\mathrm{ocean}}}}\end{array}$$

This deposition is then converted to mixing ratio and subtracted from the atmospheric CO₂ mixing ratio in the lowest model layer.

The constant value of $E_{{\mathrm{CO}}_{\mathrm{2}}}$ is based on the external emission datasets of the year 2012. It is increasing and should be adapted in the future to be time-dependent based on the current emissions of CO₂. However, taking global averages needed for this case will decrease the model performance, because all parallel processes have to wait for each other, which is why a constant value is used as a first approximation. Another possible improvement of this could be to use the yearly global growth rates of CO₂. This deposition is automatically applied when the tracer is called TRCO2. Combined effect of this deposition and the non-globally-uniform lifetime is shown in Fig. 7.

Figure 7 Percentage change in lifetime tracer for column-averaged carbon dioxide when tracer is called TRCO2 instead of CO₂. Effects due to lifetime parameterization and deposition in the ocean after 2 years of simulation.

3.4.2 Dry deposition

The gaseous dry deposition parameterization follows the widely used resistance-in-series approach introduced by Wesely and Hicks (1977). This method expresses the deposition velocity as $v_{\mathrm{d}}=-F/C\left(z_{\text{ref}}\right)$, where F is the flux to the surface, and C(z_ref) is the concentration at a reference height z_ref. The deposition velocity is then modeled using resistances

$$\begin{array}{}\text{(7)}& v_{\mathrm{d}}={\displaystyle \frac{\mathrm{1}}{(r_{\mathrm{a}}+r_{\mathrm{b}}+r_{\mathrm{c}})}},\end{array}$$

where r_a, r_b and r_c are the aerodynamic, boundary layer and canopy resistances, respectively. The aerodynamic resistance (r_a) is the same for all gases and depends on surface properties, wind speed, and atmospheric stability. The boundary layer resistance (r_b) accounts for gas transport through the quasi-laminar layer to the surface and is influenced by the molecular diffusion coefficient (D_i) of the chemical species i in air (e.g. Seinfeld and Pandis, 2016).

The canopy resistance parameterization in ICON-ART follows Baer (1992), incorporating the influence of plant species composition, vegetation physiological state, and the chemical properties of the trace gases. Based on Baer (1992), ICON-ART distinguishes n_LU=7 deposition land use classes, each with tabulated plant specific constants and leaf area indices (LAI). The total canopy resistance is calculated as a combination of stomatal (r_st), mesophyll (r_mes), cuticular (r_cut), and soil (r_soil) resistances, as described in Eq. (10) below.

The uptake of trace gases through stomata occurs by diffusion and is therefore inversely proportional to their molecular diffusivities. In ICON-ART, stomatal opening is influenced by photosynthetically active radiation (I_PAR) and temperature, following the multiplicative model by Jarvis (1976):

$$\begin{array}{}\text{(8)}& r_{\text{st}}=r_{\text{st,min}}\left(\mathrm{1}+{\displaystyle \frac{b}{I_{\text{PAR}}}}\right){\displaystyle \frac{\mathrm{1}}{f_T{f}_D}}.\end{array}$$

Here, r_st,min is an empirically determined minimum resistance (Körner et al., 1979), and b an empirical constant, both tabulated for each deposition land use class. The correction factor f_T accounts for temperature effects, particularly stomatal closure at excessively low or high temperatures. The term

$$\begin{array}{}\text{(9)}& f_D={\displaystyle \frac{D_i}{D_{\mathrm{0}}}}\end{array}$$

represents the dependence on the diffusivity D_i of the trace gas, where D₀ is the molecular diffusivity of the reference species for which r_st,min was determined (usually H₂O or CO₂). Once inside the stomata, trace gases are either deposited onto the hydrated surface of mesophyll cells or undergo chemical destruction. These processes are represented by the mesophyll resistance r_mes, which thus depends on the solubility and reactivity of the gas. Total stomatal resistance is then given by the sum $r_{\mathrm{s}}=r_{\text{st}}+r_{\text{mes}}$. In ICON-ART, two additional uptake processes are considered: direct deposition onto the leaf cuticle and deposition to soils, represented by r_cut and r_soil, respectively. Both processes depend on the solubility and reactivity of the trace gases and are parametrized using empirical values for SO₂ and O₃. The total canopy resistance in each grid cell is computed as

$$\begin{array}{}\text{(10)}& r_{\mathrm{c}}=\sum _{j=\mathrm{1}}^{n_{\text{LU}}}{\displaystyle \frac{x_j}{\left(\frac{\mathrm{1}}{r_{\mathrm{s}}\left(j\right)}+\frac{\mathrm{1}}{r_{\text{cut}}\left(j\right)}\right)\text{LAI}\left(j\right)+\frac{\mathrm{1}}{r_{\text{soil}}\left(j\right)}}},\end{array}$$

where x_j represents the land use fraction. This formula applies to any trace gas species, except for H²SO₄, SO₂ and O₃, for which adapted parametrizations are implemented.

3.4.3 Wet deposition

The removal of trace gases from the atmosphere via wet deposition is typically divided into two processes: in-cloud and below-cloud scavenging. In ICON-ART, both processes are parameterized using scavenging ratios, following the formulation by Simpson et al. (2012). Specifically, the in-cloud scavenging of a soluble trace gas with mixing ratio χ is described by

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\mathrm{d}\mathit{\chi }}{\mathrm{d}t}}=-\mathit{\chi }{\displaystyle \frac{W_{\text{in}}P}{\mathrm{\Delta }z}}{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\text{W}}}},\end{array}$$

where W_in is the in-cloud scavenging ratio, P is the surface precipitation rate (kg m⁻² s⁻¹), Δz is the scavenging depth (assumed to be 1000 m), and ρ_W is the density of water. Below-cloud scavenging is treated analogously, with the in-cloud ratio W_in replaced by the below-cloud scavenging ratio W_out. Scavenging ratios for the main soluble trace gases are provided in the supplementary material of Simpson et al. (2012).

4 Aerosol processes

Aerosol processes in ART are represented using a flexible log-normal modal framework (Rieger et al., 2015). Each mode can belong to the Aitken, accumulation, coarse, or giant size ranges, with mean diameters spanning from below 0.01 µm to above 10 µm. ART allows both externally (single-component) and internally mixed (multi-component) modes. The transport equations in ICON-ART are Hesselberg-averaged (indicated by a hat) meaning a variable Ψ can be decomposed into a barycentric mean with respect to the air density ρ_a and its fluctuations:

$$\begin{array}{}\text{(12)}& \widehat{\mathrm{\Psi }}=\mathrm{\Psi }-{\mathrm{\Psi }}^{\prime \prime }={\displaystyle \frac{\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathrm{a}}\mathrm{\Psi }}}{\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathrm{a}}}}}\end{array}$$

The bar over a variable indicates Reynolds-averaging. The prognostic equations for number density $\left({\widehat{\mathrm{\Psi }}}_{\mathrm{0},l}\right)$ and mass mixing ratio $\left({\widehat{\mathrm{\Psi }}}_{\mathrm{3},l}\right)$ are solved at every fast physics time step and are given then by:

$$\begin{array}{}\text{(13)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial \left(\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathrm{a}}}{\widehat{\mathrm{\Psi }}}_{\mathrm{0},l}\right)}{\partial t}}& =-\mathrm{\nabla }\cdot \left(\widehat{v}\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathrm{a}}}{\widehat{\mathrm{\Psi }}}_{\mathrm{0},l}\right)-\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathrm{a}}{v}^{\prime \prime }{\mathrm{\Psi }}_{\mathrm{0},l}^{\prime \prime }}\right)\\ & -{\displaystyle \frac{\partial }{\partial z}}\left(v_{\text{sed},\mathrm{0},l}\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathrm{a}}}{\widehat{\mathrm{\Psi }}}_{\mathrm{0},l}\right)-{\mathrm{Wa}}_{\mathrm{0},l}-{\mathrm{Ca}}_{\mathrm{0},l}\\ & -{\mathrm{Nu}}_{\mathrm{0},l}-{\mathrm{Em}}_{\mathrm{0},l}\end{array}\end{array}$$

$$\begin{array}{}\text{(14)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial \left(\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathrm{a}}}{\widehat{\mathrm{\Psi }}}_{\mathrm{3},l}\right)}{\partial t}}& =-\mathrm{\nabla }\cdot \left(\widehat{v}\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathrm{a}}}{\widehat{\mathrm{\Psi }}}_{\mathrm{3},l}\right)-\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathrm{a}}{v}^{\prime \prime }{\mathrm{\Psi }}_{\mathrm{3},l}^{\prime \prime }}\right)\\ & -{\displaystyle \frac{\partial }{\partial z}}\left(v_{\text{sed},\mathrm{3},l}\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathrm{a}}}{\widehat{\mathrm{\Psi }}}_{\mathrm{3},l}\right)-{\mathrm{Wa}}_{\mathrm{3},l}-{\mathrm{Ca}}_{\mathrm{3},l}\\ & -{\mathrm{Nu}}_{\mathrm{3},l}-{\mathrm{Em}}_{\mathrm{3},l}-{\mathrm{Co}}_{\mathrm{3},l}-{\mathrm{Ch}}_{\mathrm{3},l}-{\mathrm{Eq}}_{\mathrm{3},l}\end{array}\end{array}$$

where ρ_a is the air density, $\mathrm{\nabla }\cdot \left(\widehat{v}\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathrm{a}}}{\widehat{\mathrm{\Psi }}}_{M,l}\right)$ and $\mathrm{\nabla }\cdot \left(\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathrm{a}}{v}^{\prime \prime }{\mathrm{\Psi }}_{M,l}^{\prime \prime }}\right)$ denote the changes of the Mth moment of mode l due to advection and turbulent fluxes, respectively, $\frac{\partial }{\partial z}\left(v_{\text{sed},M,l}\stackrel{\mathrm{‾}}{{\mathit{\rho }}_{\mathrm{a}}}{\widehat{\mathrm{\Psi }}}_{M,l}\right)$ describes the sedimentation flux with $v_{\text{sed},M,l}$ being the sedimentation velocity, Wa_M,l denotes washout due to wet scavenging below-cloud, Nu_M,l denotes the nucleation, and Em_M,l denotes the emission flux. The terms Co_3,l, Ch_3,l and Eq_3,l are relevant for the third moment only and denote condensation, chemical transformation and equilibrium gas-aerosol partitioning, respectively. The aerosol dynamics processes (nucleation, coagulation, condensation, chemical transformation and gas-aerosol partitioning) are relevant when considering internally mixed aerosols and are briefly explained in the following section. Dry deposition is considered as a lower boundary condition for tracer transport in the lowest model layer. The surface flux is determined using a one-dimensional turbulence scheme, where dry deposition is represented through a parameterized deposition velocity (Rieger et al., 2015).

The standard deviation of the modes is kept constant during the whole simulation. The median diameter of modes can change during atmospheric transport and be diagnosed from the zeroth and third moment (Rieger et al., 2015; Muser et al., 2020).

4.1 Aerosol types and properties

It is crucial for an aerosol model to account for key dimensions of aerosol variability, such as chemical composition, size distribution, and mixing state (Riemer, 2002). ART is capable of representing these aspects, allowing for a comprehensive description of aerosol interactions and impacts. It considers a diverse range of aerosol types, including carbonaceous particles, mineral dust, volcanic ash, sea salt, and water-soluble species such as sulfates, nitrates, and ammonium. Table 3 categorizes these aerosols based on their properties such as size modes, mixing states, and solubility. It distinguishes between externally and internally mixed aerosols, with some species exhibiting a mixed state depending on atmospheric conditions. The listed aerosols span a range of size modes from Aitken to giant, reflecting their diverse sources and transport characteristics. While sulfate, nitrate, ammonium, and water are fully or partially soluble and internally mixed, other aerosols, such as carbonaceous particles, dust, volcanic ash, and sea salt, tend to be less soluble and can exist in both external and internal mixtures. If the mixed model includes only soluble species, the composition is assumed to be volume-averaged. However, if insoluble species are present in the mixed mode, the model assumes a core-shell structure. This assumption is critical for the aerosol optical properties, which are addressed in Sect. 5.1.

Table 3 Aerosol components considered in ICON-ART, along with their standard mode sizes, mixing states, and solubility. Choices of components, their modes and mixing state are flexible and can be user-defined.

^* Note that giant here does not refer to giant dust particles with diameters >62.5 µm (Adebiyi et al., 2023).

4.2 Aerosol dynamics

The ICON-ART model considers aerosol dynamical processes using its AERODYN module (Muser et al., 2020), which simulates key processes such as nucleation, condensation, coagulation, and gas-aerosol partitioning. AERODYN allows to consider a flexible number of log-normal modes (up to 10) to represent various aerosol sizes (Aitken, accumulation, coarse) and three different mixing states (soluble, insoluble, mixed) as well as an insoluble giant mode, which is not involved in the aerosol dynamic processes. Table 3 summarizes all possible compounds in the different modes. As aerosol dynamics changes the sizes and mixing states of particles, two mechanisms shift particles to other modes: (1) a shift to larger modes, when a threshold diameter is exceeded, and (2) a shift to the corresponding mixed mode when the soluble components in the modes exceed a 5 % mass threshold (Muser et al., 2020).

Figure 8 Schematic of the modes and aerosol dynamics processes in ICON-ART. Modes and processes can be configured by the user as needed, allowing for flexible configuration based on modeling requirements.

Nucleation forms new soluble Aitken particles from gaseous H₂SO₄ increasing the zeroth and third moment of the Aitken mode. The parametrization used in AERODYN is based on Kerminen and Wexler (1995). It calculates a critical H₂SO₄ concentration which depends on temperature and relative humidity and above which it is assumed that H₂SO₄ nucleates. In addition to forming new particles through nucleation, H₂SO₄ can condense in all modes (except the insoluble giant mode) leading to an increase of the third moment of a mode and a particle growth. The condensation is parameterized based on Whitby et al. (1991) and was adapted in ART from Riemer (2002).

Intermodal and intramodal coagulation can be activated between all aerosol modes. Coagulation increases the median diameter and reduces the total number concentration within the aerosol size distribution. The parameterizations for the coagulation rates of the zeroth and third moments are based on Riemer (2002) and references therein, particularly Whitby et al. (1991).

Although the coagulation matrix is flexible and can be user-defined, several key aspects should be considered regarding the assignment of particles to modes after coagulation:

• In the case of intramodal coagulation, particles remain within the same mode.

• For intermodal coagulation, particles are assigned to the mode with the larger median diameter.

• If a mixed mode is involved, the resulting particles are assigned to the mixed mode.

• Coagulation between an insoluble and a soluble mode results in particles initially remaining in the insoluble mode. These particles may later be transferred to the mixed mode if the mass of soluble material exceeds a predefined threshold.

The equations for nucleation, condensation, and coagulation are given in Appendix B.

The ISORROPIA-2 model by Fountoukis and Nenes (2007) is coupled to ART to calculate gas-aerosol partitioning according to thermodynamics equilibrium. Its aerosol system considers the species potassium (K⁺), sodium (Na⁺), magnesium (Mg²⁺), calcium (Ca²⁺), ammonium (NH${}_{\mathrm{4}}^{+}$), sulfate (SO${}_{\mathrm{4}}^{\mathrm{2}-}$), nitrate (NO${}_{\mathrm{3}}^{-}$), chloride (Cl⁻), water (H₂O) and derives an equilibrium state for these species in the gas, liquid, and solid phase. The processes are called in the following order: coagulation, aerosol-gas partitioning, condensation, nucleation, and shifting.

4.3 Aerosol water uptake

The water uptake by aerosols in the atmosphere depends on the ambient relative humidity and the chemical composition of the aerosol. The chemical composition is determined by the molality, which is the amount of a component per kg of solvent. As molality or relative humidity increases, the liquid water content also increases. Conversely, crystallization can occur when the relative humidity falls below the efflorescence relative humidity (ERH). The relative humidity at which no further water uptake occurs depends on the aerosol components and their efflorescence properties. Therefore, determining the crystallization point for an internally mixed aerosol is significantly more complex compared to an externally mixed aerosol.

The thermodynamic model ISORROPIA-2 (Fountoukis and Nenes, 2007) is used to calculate the water uptake by aerosols. However, this comes with computational overhead, which is undesirable, especially when focusing solely on water uptake on aerosols like sea salt. Due to the high hygroscopicity of sea salt aerosol, which increases aerosol particle mass and consequently affects processes such as removal mechanisms and optical properties, an alternative method for water uptake has been implemented in ICON-ART. This alternative method, previously applied in the COSMO-ART, is described in Lundgren et al. (2013). In contrast to ISORROPIA-2, this approach bypasses full thermodynamic equilibrium calculations and instead uses an empirical hygroscopic growth formulation for sea salt. This greatly reduces computational cost while retaining the key impact of water uptake on particle mass and related aerosol processes.

4.4 Subpollen particles (SPPs)

As explained in Sect. 2.6, ART employs a parameterization for pollen emission (based on EMPOL-parameterization, Zink et al., 2013) and uses it to forecast alder, birch, grass and ragweed pollen at MeteoSwiss and DWD operationally. Due to their large size, pollen are not considered in the majority of processes in ICON-ART. Although they seem to be quite efficient in nucleating ice or as cloud condensation nuclei, they do not reach relevant altitudes in sufficiently large numbers to actually impact microphysical processes. In recent years, so-called subpollen particles (SPPs), released by those large pollen grains rupturing and bursting, have received increasing attention, since their significantly smaller size enables them to reach those altitudes in larger numbers. To reflect this increased interest, Werchner et al. (2022) implemented a pollen bursting parameterization into ICON-ART based on physical assumptions, observations and processes (according to Zhou, 2014), enabling it to emit SPPs.

The parameterization's driver is the turgor pressure that builds up inside each pollen grain. Once this turgor pressure reaches a pollen specific critical value that the pollen walls can no longer withstand the pollen bursts and releases SPPs into the atmosphere. The turgor pressure's temporal development is formulated to be:

$$\begin{array}{}\text{(15)}& p_{\text{T}}\left(t_{\text{m}}\right)=p_{\text{a}}+\mathrm{\Delta }\mathit{\pi }\left[\mathrm{1}-\exp \left(-{\displaystyle \frac{\mathrm{3}Ek}{r}}{\displaystyle \frac{{\mathit{\rho }}_{\text{w}}}{{\mathit{\rho }}_{\mathrm{0}}}}{t}_{\text{m}}\right)\right].\end{array}$$

In Eq. (15), p_T and p_a are the turgor pressure and the ambient pressure, respectively, in Pa, t_m is the model time step in s, Δπ is the difference in water potential in Pa, E is the compression module of water in N m⁻², k is the pollen grain's water permeability in m s⁻¹ Pa⁻¹, r is the pollen grain's radius in m and ρ_w and ρ₀ are water and pollen density, respectively, in kg m⁻³.

A complete description of this parameterization is given in Werchner et al. (2022). Werchner et al. (2022) reported that SPP concentrations (used only as INP, not as CCN) vary between 10²×10⁶ m⁻³, with a mean value of 4×10³ m⁻³ (particularly relevant in warmer levels suitable for biological ice nucleation) while mean pollen concentration amount to $\mathrm{3.4}\times {\mathrm{10}}^{-\mathrm{3}}$ m⁻³.

5 Interactions and feedback mechanisms

5.1 Aerosol-radiation interaction

ICON uses ecRad (Hogan and Bozzo, 2018) as the standard radiation scheme for numerical weather prediction (Rieger et al., 2019). To account for ARI, ART computes the local radiative transfer parameters based on the aerosol optical properties – mass extinction coefficient, single scattering albedo, and asymmetry parameter – as well as the prognostic aerosol mass concentrations at each grid point and model level. These parameters are then passed as input to the radiation scheme (Rieger et al., 2017).

In addition to the prognostic ARI calculations, ART includes forward operators to diagnose aerosol optical depth (AOD) and attenuated backscatter at different wavelengths. These diagnostics are derived by multiplying the prognosed aerosol mass concentrations with mass extinction and backscatter coefficients (Hoshyaripour et al., 2019).

ART provides two approaches for specifying aerosol optical properties. In the first, these properties are precomputed offline and stored in lookup tables for use in both prognostic and diagnostic calculations (Hoshyaripour et al., 2019; Muser et al., 2020). This method accounts for the variability in aerosol size distribution by applying polynomial fits to the optical properties (Gasch et al., 2017). However, it is limited in capturing variations in particle composition and mixing state (e.g., coating), particularly when aerosol dynamics lead to the formation of internally mixed particles during simulations.

To address these limitations, the second approach computes aerosol optical properties online using MieAI (Kumar et al., 2024), enabling a more flexible and physically consistent representation. Further details on both approaches are provided below.

5.1.1 Prescribed aerosol optics

In ICON-ART, the prescribed optical properties for aerosols were originally implemented as hard-coded tables within the model code (Gasch et al., 2017; Rieger et al., 2017). This static approach limited flexibility and made updates or extensions to the optical property datasets cumbersome. To overcome these limitations, the optical properties have been transitioned to external NetCDF files, allowing for a more modular and user-configurable setup. This change enables users to easily update or switch between different aerosol optical property datasets without modifying the source code, facilitates consistency across simulations, and supports the use of more detailed and physically representative optical data, including properties for varying size distributions, compositions, and wavelengths.

The NetCDF database consists of extinction coefficients, single scattering albedo values, asymmetry parameters, and backscattering coefficients in three different wavelength groups:

• 30 prognostic wavebands, as used by ecRad

• 9 diagnostic wavelengths from AERONET: 340.0, 380.0, 440.0, 500.0, 550.0 , 675.0, 870.0, 1020.0 and 1064.0 nm

• 3 diagnostic wavelengths for lidar (ceilometer or satellites): 355.0, 532.0 and 1064.0 nm

An overview of all available modes and their microphysical properties can be found in Appendix C.

5.1.2 Online aerosol optics with MieAI

Accurate and efficient estimation of aerosol-radiation interaction is paramount for precise weather and climate prediction. The interaction between aerosols and radiation is governed by their optical properties, which are determined by various aerosol attributes such as morphology, size distribution, and chemical composition. These properties undergo significant changes due to chemical reactions and microphysical processes as aerosol particles are transported through the atmosphere. Traditional methods for computing aerosol optical properties rely on precomputed look-up tables (LUTs), which are computationally inexpensive but prone to substantial errors.

To address this limitation, we have integrated MieAI, a recently developed neural network-based approach, to compute the optical properties of internally mixed aerosol particles with ICON-ART (Kumar et al., 2024). MieAI is a fully connected feed forward neural network with 4 hidden layers, each layer having 64 neurons and uses Gaussian Error Linear Unit (GELU) as the activation function. It emulates Mie theory, enabling on-the-fly optical properties calculations that account for variations in aerosol size distribution and chemical composition. MieAI takes aerosol mass concentration of aerosol components from ICON-ART simulations as input and outputs key optical properties, including the mass extinction coefficient, single scattering albedo, and asymmetry parameter. These estimated optical properties are then used in radiative transfer calculations to calculate the radiative effects of aerosols more efficiently and accurately. Fortran-Keras bridge (FKB) was used for online coupling of MieAI with ICON-ART (Ott et al., 2020). See Kumar et al. (2024) for detailed discussion about MieAI.

Figure 9 Comparison of net shortwave radiative flux estimated using MieAI against those estimated using Look-up table (LUT) approach for a case study involving the La Soufrière volcanic eruption (denoted by the black triangle) event simulated using ICON-ART. Here, panel (a) shows the net SW flux estimated using LUT, (b) shows the same estimated using MieAI and (c) shows the absolute difference between them. The volcanic plume is depicted in panel (d) whereas panel (e) shows the mixed mode aerosols within the plume. Panel (f) zooms panel (c) over the plume region.

An example usecase is shown in Fig. 9 for an ICON-ART simulation involving the La Soufrière volcanic eruption event in April 2021. Here, Fig. 9a shows the net shortwave flux at the surface estimated using traditional LUT approach whereas Fig. 9b shows the same using MieAI apporach and the difference between them are shown in Fig. 9c. As can be seen clearly from Fig. 9c and f, MieAI estimates lower SW flux (> 8 W m⁻²) at the surface over the region containing the volcanic plume depicted in Fig. 9d. A large part of this reduction in SW flux is due to the presence of mixed modes in the plume captured realistically by MieAI, causing enhanced extinction of solar radiation reaching the surface.

5.1.3 Multiple radiation call

Apart from including aerosol-radiation interactions in the modeling of atmospheric dynamics and composition, its diagnostic quantification is key to understanding aerosol impacts. To quantify the direct radiative effect (DRE), it is important to impose a technique that does not alter meteorological conditions, which would induce feedback processes and ultimately modify the aerosol effect. To diagnose the DRE, we therefore implemented a radiation multiple call scheme in ICON-ART, as it is also implemented in other models (e.g. HadCM (Woodward, 2001), LMDz-INCA (Balkanski et al., 2007), WRF-Chem (Zhao et al., 2013), MONARCH (Klose et al., 2021)). This approach executes calls of the radiation routine not just once (as needed to include aerosol-radiation interactions in the simulation), but multiple times for diagnostic purposes. In the first call, all aerosols are neglected and the radiation fluxes are stored. In optional intermediate calls, the effect of a single aerosol mode (and thus aerosol species) on the DRE can be quantified by either including or omitting the respective mode (see below). In the last call, all aerosol-radiation interactions are accounted for, just as in the conventional model run. Therefore, the model integration is not impacted by execution of the multiple call scheme. The aerosol DRE is then directly obtained as the difference of the radiation fluxes between the two (or multiple) calls.

In ICON-ART, three different types of the radiation multiple call are available, following the concept implemented in Klose et al. (2021):

Type I – Double call:

two calls, one without aerosols, one with all aerosols

Type II – Inclusive multiple call:

first call without aerosols, one call per mode considering only that mode, last call with all aerosols

Type III – Exclusive multiple call:

first call without aerosols, one call per mode considering all aerosols except that mode, last call with all aerosols.

The double call includes in total n=2 calls per radiation timestep. The inclusive and exclusive multiple calls entail $n=n_{\text{modes}}+\mathrm{2}$ calls, with n_modes being the number of aerosol modes, as for each mode there is one intermediate call. The DRE for each mode depends on the aerosol distribution of all modes in the column, such that the total DRE for all modes is not exactly equal to the sum of the individual contributions per mode, especially in regions with high aerosol load. Thus, call types II and III yield similar, but not equal results, especially in those regions. As call type III includes effects of other aerosol modes, we consider this type preferential over the conceptually simpler type II.

For the first two call types (double call and inclusive multiple call), the DRE for a specific mode m is evaluated for both shortwave (wb=sw) and longwave (wb=lw) radiation, and at the surface (lev=sfc) and the top of atmosphere (lev=toa), DRE${}_{\text{wb}}^{\text{m}}{|}_{\text{lev}}$. It is calculated from the radiative flux $F_{\text{wb}}^{\text{m}}{|}_{\text{lev}}$ as

$$\begin{array}{}\text{(16)}& {\mathrm{DRE}}_{\text{wb}}^{\text{m}}{|}_{\text{lev}}=F_{\text{wb}}^{\text{m}}{|}_{\text{lev}}-F_{\text{wb}}^{\text{none}}{|}_{\text{lev}},\end{array}$$

where $F_{\text{wb}}^{\text{none}}{|}_{\text{lev}}$ denotes the flux without any aerosol-radiation interactions, obtained from the first call. For call type III, the exclusive multiple call, the DRE for all aerosols is calculated in the same manner as

$$\begin{array}{}\text{(17)}& {\mathrm{DRE}}_{\text{wb}}^{\text{all}}{|}_{\text{lev}}=F_{\text{wb}}^{\text{all}}{|}_{\text{lev}}-F_{\text{wb}}^{\text{none}}{|}_{\text{lev}}.\end{array}$$

To obtain the DRE for individual modes, the situation is slightly different for type III. The flux obtained from the call for a specific mode m, $F_{\text{wb}}^{\text{all-m}}{|}_{\text{lev}}$, includes the effect of all modes except the desired mode. Thus, the reference flux for the DRE in this case is the one including the effect of all aerosols, yielding

$$\begin{array}{}\text{(18)}& {\mathrm{DRE}}_{\text{wb}}^{\text{all}}{|}_{\text{lev}}=F_{\text{wb}}^{\text{all}}{|}_{\text{lev}}-F_{\text{wb}}^{\text{all-m}}{|}_{\text{lev}}.\end{array}$$

All calculated DREs are available as output, yielding four diagnostic variables (2 wavebands and two levels) for the double call (type I) and another four per mode for the multiple calls, e.g. 16 for the case of a multiple call with the three default dust modes dust_acc, dust_coa and dust_gia (cf. Table 3). They are available both as instantaneous values and as accumulated quantities. For the latter, note that the accumulation occurs between times at which the DRE is evaluated, i.e. with the radiation time step.

An example of the short wave DRE for dust is presented in Fig. 10. Shown is the annual average DRE for all three default dust modes (Fig. 10a–c), and the total DRE for all modes (Fig. 10d), obtained using the exclusive multiple call. As expected, the main contribution to the short wave dust DRE stems from the accumulation dust mode, which yields a strong cooling effect of up to 8.4 W m⁻² near the horn of Africa. In contrast, the giant dust mode features a positive short wave DRE at the top of the atmosphere over land in areas of high dust loading due to the stronger absorption. Dust in the coarse mode contributes with an intermediate signal to the total DRE, which is overall negative in our example. Note that the size mode terminology used here follows the mode definitions given in Table 3 and differs from the dust size classification proposed in Adebiyi et al. (2023).

Figure 10 Direct radiative effect (DRE) for shortwave radiation at the top of atmosphere. Panels show the DRE for the different dust modes, i.e. (a) accumulation mode, (b) coarse mode, (c) giant mode (see Table 3), and (d) all aerosols. Note that the color scales differ between panels.

5.2 Aerosol-cloud interaction

To investigate aerosol–cloud interactions (ACI), ART is coupled with the two-moment cloud microphysics scheme of ICON (Seifert and Beheng, 2006), which provides a detailed and physically consistent representation of cloud processes. This scheme includes six hydrometeor categories – cloud droplets, ice crystals, rain, snow, graupel, and hail – and predicts both their mass and number concentrations through a set of prognostic budget equations. The evolution of liquid-phase hydrometeors is governed by processes such as condensational growth, autoconversion, accretion, self-collection, breakup, and freezing. In contrast, ice-phase hydrometeors are subject to diffusional growth, aggregation, self-collection, riming, secondary ice production (ice multiplication), and melting. The two-moment formulation enhances the model’s ability to capture the sensitivity of cloud development and precipitation formation to aerosol perturbations, making it well-suited for studies of ACI (Rieger et al., 2017; Gruber et al., 2019).

5.2.1 CCN and IN activation

A physically-based parameterization for cloud condensation nuclei (CCN) activation, following Abdul-Razzak et al. (1998) and Abdul-Razzak and Ghan (2000), is implemented within the two-moment cloud microphysics scheme of ICON-ART. This approach enables a more realistic coupling between aerosols and cloud droplet formation by using grid-scale vertical velocity and a parameterized subgrid contribution, available moisture, and soluble aerosol concentrations to compute the number concentration and mass mixing ratio of newly formed cloud droplets. The activation is computed only under appropriate atmospheric conditions, such as the presence of updrafts and supersaturation with respect to liquid water. With this implementation, microphysical computations are removed from the ART module; ART now supplies the activation parameterization and outputs the newly formed cloud droplets, which are then further processed by ICON’s native two-moment microphysics scheme Seifert and Beheng (2006). Initial testing and validation of the parameterization are performed successfully.

Figures 11 and 12 show preliminary results from idealized simulations of a warm bubble. The model setup follows the Weisman Klemp test case (Weisman and Klemp, 1982). A predefined sea salt concentration of 2×10⁷ # kg⁻¹ is uniformly distributed throughout the domain and equally distributed between the accumulation and coarse mode.

Figure 11 displays the number of activated cloud condensation nuclei (nccn) in # kg⁻¹, accumulated over 640 s from the start of the simulation. Sea salt aerosols are activated within the updraft region generated by the warm bubble. Figure 12 illustrate the ratio of activated particles to available sea salt aerosols as a function of vertical velocity, for (a) accumulation and (b) coarse mode. The results indicate that a substantial fraction of sea salt in coarse mode gets activated, whereas, only a small portion of sea salt in accumulation mode undergoes activation. This outcome aligns with Köhler theory, which predicts that larger particles are more likely activated due to their lower critical supersaturation.

Figure 11 Number of activated sea salt aerosols from an idealized warm bubble simulation based on the Weisman Klemp test case (Weisman and Klemp, 1982).

Figure 12 Ratio of activated particles to available sea salt aerosols as a function of vertical velocity, shown for (a) accumulation mode and (b) coarse mode.

For heterogeneous ice nucleation, the calculation of the total surface area of all dust aerosol modes is fully implemented on the ART side and serves as the basis for determining the number of ice nucleating particles according to the ice nucleation active site (INAS) density approach (Hoose and Möhler, 2012), following the empirical parameterizations for immersion freezing and deposition ice nucleation (e.g., Ullrich et al., 2017). This enables the prognostic treatment of heterogeneous ice nucleation in accordance with the evolving aerosol size distribution and composition within the model. Coupling of this INAS-based ice nucleation treatment with the two-moment microphysical scheme of ICON Seifert and Beheng (2006) is currently in progress. Once completed, this coupling will allow for a consistent representation of aerosol-cloud interactions, including the formation of ice crystals from mineral dust under mixed-phase and cirrus conditions.

5.2.2 Dusty cirrus parametrization

In addition to the ACI on the grid-scale, ICON-ART includes a sub-grid parametrization of dusty cirrus. This special parametrization is beneficial or even necessary because dusty cirrus clouds can form due to a small-scale mixing instability, which is difficult to represent in model configurations that are not eddy resolving (Seifert et al., 2023).

The dusty cirrus parametrization uses the mass concentration of mineral dust c_mode assuming three dust modes with mode∈{dustA, dustB, dustC}. Humidity is quantified by the ice saturation ratio $s_{\text{ice}}=p_{\mathrm{v}}/p_{\text{sat,ice}}$, where p_v is the vapor pressure and p_sat,ice is the saturation vapor pressure over ice. Atmospheric stability is characterized by the temperature lapse rate

$$\begin{array}{}\text{(19)}& {\mathit{\gamma }}_k={\left.{\displaystyle \frac{\partial T}{\partial z}}\right|}_k\approx {\displaystyle \frac{T_k-T_{k+\mathrm{1}}}{\mathrm{\Delta }z}}\end{array}$$

Note that ICON uses top-down indices, i.e., level k+1 is below level k. Dusty cirrus occurs in model level k if the following conditions are fulfilled:

$$\begin{array}{}\text{(20)}& {\displaystyle }{T}_k<\mathrm{240}\text{ K}\text{(21)}& {\displaystyle }{\widehat{c}}_{\text{dust},k}=\underset{j=k+\mathrm{1}}{\overset{k+N+\mathrm{1}}{max}}\left(c_{\text{dustB},j}+\mathrm{2}{c}_{\text{dustC},j}\right)>c_{\text{dust}}^{*}\text{(22)}& {\displaystyle }{\widehat{s}}_{\text{ice},k}=\underset{j=k}{\overset{k+N}{max}}{s}_{\text{ice},j}>s_{\text{ice}}^{*}\text{(23)}& {\displaystyle }{\widehat{\mathit{\gamma }}}_k=\underset{j=k}{\overset{k+\mathrm{1}}{min}}{\mathit{\gamma }}_j<{\mathit{\gamma }}^{*}\end{array}$$

with empirically determined thresholds $c_{\text{dust}}^{*}=\mathrm{50}$ µg kg⁻³, $s_{\text{ice}}^{*}=\mathrm{0.7}$, and ${\mathit{\gamma }}^{*}=-\mathrm{6.5}$ K km⁻¹ and with N=4 corresponding to a vertical depth of approximately 1500 m. That the scheme uses information from 4 layers below to predict dusty cirrus layers corresponds to the convective overturning and the mixing instability. Note that $c_{\text{dust}}^{*}$ does not include that smallest dust mode dustA, and the largest mode dustC has double the weight of dustB. The fact that the larger dust modes, dustB and dustC, are better predictors for the occurrence of a dusty cirrus than dustA is consistent with the increased ability of large mineral dust particles to act as INPs, whereas smaller particles are less relevant for the formation of ice clouds by heterogeneous nucleation (DeMott et al., 2010, 2015). For the sub-grid dusty cirrus a cloud fraction of one is assumed and the sub-grid ice water content is set to 80 mg m⁻³ with a linear tapering at the boundaries.

Figure 13 shows the results of a global ICON-ART simulation of a Saharan dust event over Europe leading to dusty cirrus cloud formation. The simulation uses the R03B06 icosahedral grid with a equivalent grid spacing of approximately 20 km. The satellite observations show an extended cirrus deck over central Europe with outgoing longwave radiation (OLR) of 200 W m⁻² and lower. Without the dusty cirrus parameterization underestimates the extent of the dusty cirrus cloud. With the sub-grid parameterization the simulation is improved but OLR is underestimated at the boundary of the cirrus cloud deck. This suggests that the simplistic threshold based parametrization needs to be further improved in the future.

Figure 13 Comparison of global ICON-ART simulation for 12:00 UTC of 6 May 2022 with CERES Level 2 satellite data of outgoing longwave radiation at the top of atmosphere. Observations (left), ICON-ART without dusty cirrus parametrization (center), and ICON-ART with dusty cirrus parametrization (right).

Dusty-cirrus parameterization performs robustly from mesoscale model resolutions (∼ 10–20 km) down to convection-resolving scales (∼ 1–3 km), as long as dust mass and number concentrations are physically reasonable. Coarser grids capture the large-scale cirrus response, while finer or convection-resolving scales benefit from enhanced representation of vertical motions and aerosol gradients.

6 Code infrastructure and periphery

6.1 Coupling ART with ICON

ART is currently designed as a submodel for ICON such that the development and organization of ART is independent from the base model in most aspects. Its coupling to ICON, however, is extensive, as Fig. 14 shows. Apart from several points within the time loop where ART provides additional processes, ART also influences some ICON processes directly. Beyond this process-based coupling, ART actively uses some structures introduced in the ICON code. The most prominent example for this is the tracer array, which ART extends to accommodate the additional tracers, whose set is defined by the user with ICON-ART's flexible tracer framework (Schröter et al., 2018). Accessing ICON's tracer array allows the ART tracers to be treated the same way as the base ICON tracers with regard to processes like advection or turbulent diffusion. Additionally, ART uses and extends structures that represent tracer metadata and requires grid information provided by ICON.For the bulk of the communication between ICON and ART, a set of interface modules are used. The interfaces enable ICON to call process and other routines in ART altering the model state, and ART to retrieve crucial information and structures to fulfill its tasks. The separation of the different modules is based on the general stages of the model (examples: mo_art_init_interface.f90, mo_art_general_interface.f90 and mo_art_diagnostics_interface.f90) and the processes provided (examples: mo_art_sedimentation_interface.f90, mo_art_aerodyn_interface.f90, mo_art_coagulation_interface.f90 and mo_art_reaction_interface.f90). While previously the interfacing between ICON and ART was managed in the host model, it moved to the ART codebase in recent developments. This reduces the additional infrastructural burden when developing ART, since the ICON codebase is only altered if actually required.

Figure 14 Schematic of the coupling of ICON–ART. The sequence in which processes of ICON are executed is illustrated by the blue boxes. Processes of ART are illustrated by the orange boxes. An orange frame around a blue box indicates, that the according code is part of the ICON tracer framework but ART tracers are treated inside this framework. The gray and black loops indicate the sequences of the time integration. Some examples of input and output fields are also show at the top and bottom of the loop.

6.2 Flexible tracer and mapping framework

ART showed great success using its flexible tracer framework introduced by Schröter et al. (2018) in such a way that further development was conducted to extend this flexibility to processes. One of these processes is coagulation, which requires a coagulation-matrix describing potential mode shifts when coagulation occurs. This matrix is to be provided to ART via a corresponding XML-file, while also adding a new XML-tag for the tracers indicating whether or not they should be considered for coagulation. An important part for the consideration of aerosols and tracers is their emission into the atmosphere. Emission parameterizations were initially assigned hardcoded by their and the corresponding tracers' name. This put a strain to the actual flexibility achieved, so a new method was introduced. By storing information (like number of emission fluxes or sizes of the emitted material) about a specific emission parameterization and assigning substances to it in a new XML-structure, a mapping routine can map the emission fluxes to the correct mode of the specific substance before calling the emission parameterization itself. Also the new OEM uses XMLs to control its behaviour in ART. The usability of the underlying XML-structure for providing the model with information also goes beyond the actual tracers, modes and processes, since diagnostic variables can also be defined with it.

Table 4 Overview of standard ICON-ART configurations with the relevant tracers, and processes.

6.3 GPU and exascale readiness

Since the development of ART has been largely driven by scientific motivations, particularly in the implementation of new processes and features, the resulting historically grown codebase contains numerous cross-dependencies between modules. This complexity makes it difficult to trace the code flow and hinders further development and maintenance. To address this issue, modularization efforts have been undertaken in recent years, streamlining the codebase, reducing inter-module dependencies, and providing a more intuitive interface. An additional benefit of these modularization efforts is the simplified parallelization and the potential to extend ART for use on new heterogeneous and future exascale HPC architectures.

Due to both its scientific design and historically evolved codebase, most processes, features, and modules have been developed exclusively for CPU usage. Although CPU parallelization is implemented and widely used across large parts of ICON, ART rarely exhibits comprehensive parallelization of its features. In recent years, the use of GPUs has gained increasing interest in scientific fields and is likely to continue doing so. Hence, more of the ART codebase needs to be ported to heterogeneous HPC architectures with GPUs to remain competitive and future-proof - due to its operational use at MeteoSwiss pollen-related ART-parts were already ported. As a first step, the online emissions module (OEM) and the vegetation photosynthesis and respiration model (VPRM) of ART have been refactored and ported using OpenACC compiler directives, following ICON's GPU porting approach in an attempt to accelerate the porting process while maintaining a unified codebase for both CPUs and GPUs. VPRM was introduced to enable the simulation of atmospheric carbon dioxide, which is not only affected by anthropogenic emissions but also by exchange with the biosphere. A first application of VPRM in ICON-ART was demonstrated by Ponomarev et al. (2026). Inverse emission estimation with ICON-ART as first demonstrated by Steiner et al. (2024) and Thanwerdas et al. (2025) is now possible with a GPU-accelerated version of ICON-ART and has already been successfully tested and validated on Switzerland’s new HPC system, ALPS, leveraging NVIDIA’s Grace-Hopper superchips. A systematic quantitative evaluation of CPU versus GPU performance for ICON-ART with OEM, including detailed speedup analyses, is currently underway and will be presented in a forthcoming dedicated study.

6.4 Standard configurations and simulation workflow

In ICON-ART, a set of standard experiment configurations has been established to ensure consistent testing and validation across model releases. Each configuration, see Table 4, combines specific model components, tracer species, and physical or chemical processes relevant to targeted research or operational applications. For example, “Dust” focuses on mineral dust and its radiative feedback, while “Volcano” emphasizes volcanic ash and aerosol-radiation interactions. “Nat. Aero.” incorporates a wide range of aerosol tracers and processes without radiative coupling, providing a comprehensive base for sensitivity studies. “Pollen” addresses both pollen and subpollen particles, representing bio-aerosols, and their relevant processes and the chemistry experiments (“Simp. Chem.” and “Comp. Chem.”) are tailored for atmospheric composition forecasts, differing in the complexity of the chemistry scheme being used. “OEM” is configured for regional applications with full chemistry and online emission treatment. These standard configurations are fully supported, regularly maintained, and systematically tested with every ICON-ART release to ensure model integrity, reproducibility, and performance across use cases.

All standard experiments are also integrated into the workflow management solution auto-icon (Baer et al., 2025a, b), which enables users to orchestrate a simulation workflow covering the retrieval of input data, compilation, and execution of the model as well as post-processing and visualization tasks. This integration enables a quick and smooth start for new users as well as easy-to-use and flexible testing and evaluation capabilities with the provided experiments during model development. In more complex model applications, auto-icon can significantly enhance throughput by executing all non-dependent tasks in parallel and wrapping individual tasks into larger allocations on the HPC system, which was shown to reduce the time-to-solution by 7 % or more (Marciani et al., 2025).

7 Conclusion and outlook

The current version, ICON-ART 2025.10, supports a wide range of applications, from air quality and weather forecasting to climate and Earth system simulations. It features a modular and flexible design, enabling the use of both simplified and detailed chemistry mechanisms, interactive aerosols, and online coupling with dynamics, radiation, and cloud microphysics.

Beside scientific research, ICON-ART is currently used in operational applications, including mineral dust forecasting at the Deutscher Wetterdienst (DWD), as well as pollen forecasting at both DWD and MeteoSwiss. These activities demonstrate the model’s robustness, flexibility, and ability to support real-time decision-making in environmental and health-related services.

Ongoing developments focus on further enhancing the coupling between atmospheric composition and the physical core of ICON. This includes the integration of ART with the seamless ICON configuration (Müller et al., 2025) for unified weather and climate simulations, which will allow for consistent aerosol and trace gas treatment across spatial and temporal scales. In addition, the implementation of organic aerosol processes is underway to improve the representation of both primary and secondary organic aerosol formation and their interactions with radiation and clouds.

Together, these efforts aim to strengthen ICON-ART as a versatile and comprehensive tool for studying the interactions between atmospheric composition and the Earth system in a fully coupled, scale-aware modeling framework.

Appendix A: Acronyms

The acronyms used in this study are summarized in Table A1.

Table A1 Summary of key acronyms used in this study, grouped by topic area.

Appendix B: Aerosol Dynamics

B1 Nucleation

Based on the parametrization by Kerminen and Wexler (1995), new sulfate particles form when a critical concentration of H₂SO₄ (given in µg m⁻³) is reached:

$$\begin{array}{}\text{(B1)}& c_{\mathrm{crit}}=\mathrm{0.16}\exp \left(\mathrm{0.1}T-\mathrm{3.5}{\displaystyle \frac{\mathrm{RH}}{\mathrm{100}}}-\mathrm{27.7}\right)\end{array}$$

with temperature T in K and relative humidity RH in %. The concentration above this threshold nucleates into the soluble Aitken mode with a mass rate of

$$\begin{array}{}\text{(B2)}& {\mathrm{Nu}}_{\mathrm{3},{\mathrm{sol}}_{\mathrm{Ait}}}={\displaystyle \frac{c_{{\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}}-\frac{c_{\mathrm{crit}}}{{\mathit{\rho }}_p}}{\mathrm{\Delta }t}}.\end{array}$$

$c_{{\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}}$ is the actual concentration of H₂SO₄. The number concentration (in kg⁻¹) follows from the assumed size distribution as

$$\begin{array}{}\text{(B3)}& {\mathrm{Nu}}_{\mathrm{0},{\mathrm{sol}}_{\mathrm{Ait}}}={\displaystyle \frac{\mathrm{6}}{\mathit{\pi }{\mathit{\rho }}_p}}\cdot {\displaystyle \frac{\exp \left(\mathrm{4.5}\cdot {\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_{{\mathrm{sol}}_{\mathrm{Ait}}}\right)\right)}{d_{\mathrm{0},{\mathrm{sol}}_{\mathrm{Ait}}}^{\mathrm{3}}}}\cdot {\mathrm{Nu}}_{\mathrm{3},{\mathrm{sol}}_{\mathrm{Ait}}}.\end{array}$$

B2 Condensation

The condensation rate of the third moment of a mode l is based on Whitby et al. (1991); Riemer (2002) and is given by

$$\begin{array}{}\text{(B4)}& \stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{o}}_{\mathrm{3},l}={\displaystyle \frac{\mathrm{6}}{\mathit{\pi }}}{\mathit{\chi }}_{\mathrm{T}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\mathit{\chi }\left(d_{\mathrm{3},\mathrm{0}}\right){\mathrm{\Psi }}_{\mathrm{3},l}\left(d\right)\text{d}{d}_{\mathrm{3},l}={\displaystyle \frac{\mathrm{6}}{\mathit{\pi }}}{\mathit{\chi }}_{\mathrm{T}}{I}_l.\end{array}$$

χ_T is independent of the particle size and depends on thermodynamic variables. It will be eliminated later. χ(d_3,l) differs for different size regimes. The expression for the near-continuum and the free-molecular regime are given by

$$\begin{array}{}\text{(B5)}& {\mathit{\chi }}_{\mathrm{nc}}\left(d_l\right)=\mathrm{2}\mathit{\pi }{D}_{\mathrm{g},{\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}}{d}_l\end{array}$$

and

$$\begin{array}{}\text{(B6)}& {\mathit{\chi }}_{\mathrm{fm}}\left(d_l\right)={\displaystyle \frac{\mathit{\pi }\mathit{\alpha }{\stackrel{\mathrm{‾}}{c}}_{{\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}}}{\mathrm{4}}}{d}_l^{\mathrm{2}}\end{array}$$

with $D_{\mathrm{g},{\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}}$ the diffusion coefficient of H₂SO₄, α the accommodation coefficient, and ${\stackrel{\mathrm{‾}}{c}}_{{\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}}$ the mean molecular velocity of H₂SO₄. Using these two formulations, the integral I_l in Eq. (B4) is evaluated separately as

$$\begin{array}{}\text{(B7)}& I_l^{\mathrm{nc}}=\mathrm{2}\mathit{\pi }{D}_{\mathrm{g},{\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}}{M}_{\mathrm{1},l}\end{array}$$

and

$$\begin{array}{}\text{(B8)}& I_l^{\mathrm{fm}}={\displaystyle \frac{\mathit{\pi }\mathit{\alpha }{\stackrel{\mathrm{‾}}{c}}_{{\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}}}{\mathrm{4}}}{M}_{\mathrm{2},l}.\end{array}$$

With the harmonic mean of the expressions for the near-continuum and free-molecular regime, Eq. (B4) can be written as

$$\begin{array}{}\text{(B9)}& \stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{o}}_{\mathrm{3},l}={\displaystyle \frac{\mathrm{6}}{\mathit{\pi }}}{\mathit{\chi }}_{\mathrm{T}}{\displaystyle \frac{I_{l,\mathrm{fm}}\cdot I_{l,\mathrm{nc}}}{I_{l,\mathrm{fm}}+I_{l,\mathrm{nc}}}}\end{array}$$

Assuming a much slower condensation of H₂SO₄ than its production leads to an equilibrium state in which the total condensation rate equals the production rate of H₂SO₄:

$$\begin{array}{}\text{(B10)}& \stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{o}}_{\mathrm{3}}=\sum _l\stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{o}}_{\mathrm{3},l}={\dot{M}}_{\mathrm{3}}\end{array}$$

with ${\dot{M}}_{\mathrm{3}}$ the third moment production rate of gaseous H₂SO₄. A third-moment condensation rate can be formulated as

$$\begin{array}{}\text{(B11)}& \stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{o}}_{\mathrm{3},l}=\sum _l\stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{o}}_{\mathrm{3},l}={\dot{M}}_{\mathrm{3}}{\mathrm{\Omega }}_l\end{array}$$

using a dimensionless coefficient defined as:

$$\begin{array}{}\text{(B12)}& {\mathrm{\Omega }}_l={\displaystyle \frac{\stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{o}}_{\mathrm{3},l}}{\stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{o}}_{\mathrm{3}}}}={\displaystyle \frac{\stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{o}}_{\mathrm{3},l}}{{\dot{M}}_{\mathrm{3}}}}={\displaystyle \frac{I_l}{{\sum }_l{I}_l}}\end{array}$$

Thus, the condensation rate is independent from χ_T now and can be written for the third moment as:

$$\begin{array}{}\text{(B13)}& {\mathrm{Co}}_{\mathrm{3},l}={\displaystyle \frac{\mathit{\pi }}{\mathrm{6}}}{\displaystyle \frac{{\mathit{\rho }}_{{\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}}}{{\mathit{\rho }}_{\mathrm{a}}}}\dot{M_{\mathrm{3}}}{\mathrm{\Omega }}_l\end{array}$$

In ICON-ART, we approximate $\dot{M_{\mathrm{3}}}$ by the H₂SO₄ mass mixing ratio and the model time step. Thus, for Eq. (B13) it follows:

$$\begin{array}{}\text{(B14)}& {\mathrm{Co}}_{\mathrm{3},l}={\displaystyle \frac{c_{{\mathrm{H}}_{\mathrm{2}}{\mathrm{SO}}_{\mathrm{4}}}}{\mathrm{\Delta }t}}{\mathrm{\Omega }}_l\end{array}$$

B3 Coagulation

The parametrizations for the coagulation follow Riemer (2002) and references therein (mainly Whitby et al., 1991). The rates for the zeroth and third moment are given by

$$\begin{array}{}\text{(B15)}& {\mathrm{Ca}}_{\mathrm{0},i}={\mathrm{Ca}}_{\mathrm{0},ii}+{\mathrm{Ca}}_{\mathrm{0},ij}\end{array}$$

and

$$\begin{array}{}\text{(B16)}& {\mathrm{Ca}}_{\mathrm{3},i}={\mathrm{Ca}}_{\mathrm{3},ij}\end{array}$$

with the Ca_k,ii and Ca_k,ij for the intra- and inter-modal coagulation rates of the kth moment, respectively, and i and j referring to two different modes. In ICON-ART, only coagulation due to Brownian motion is considered so far and reads for a system of two modes i and j as (Whitby et al., 1991)

$$\begin{array}{}\text{(B17)}& \stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{a}}_{\mathrm{0},ii}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\mathit{\beta }(d_{\mathrm{1}},d_{\mathrm{2}}){\mathit{\psi }}_{\mathrm{0},i}\left(d_{\mathrm{1}}\right){\mathit{\psi }}_{\mathrm{0},i}\left(d_{\mathrm{2}}\right)\text{d}{d}_{\mathrm{1}}\text{d}{d}_{\mathrm{2}},\end{array}$$

$$\begin{array}{}\text{(B18)}& \stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{a}}_{\mathrm{0},ij}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\mathit{\beta }(d_{\mathrm{1}},d_{\mathrm{2}}){\mathit{\psi }}_{\mathrm{0},i}\left(d_{\mathrm{1}}\right){\mathit{\psi }}_{\mathrm{0},j}\left(d_{\mathrm{2}}\right)\text{d}{d}_{\mathrm{1}}\text{d}{d}_{\mathrm{2}},\end{array}$$

and

$$\begin{array}{}\text{(B19)}& \stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{a}}_{\mathrm{3},ij}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{d}_{\mathrm{1}}^{\mathrm{3}}\mathit{\beta }(d_{\mathrm{1}},d_{\mathrm{2}}){\mathit{\psi }}_{\mathrm{0},i}\left(d_{\mathrm{1}}\right){\mathit{\psi }}_{\mathrm{0},j}\left(d_{\mathrm{2}}\right)\text{d}{d}_{\mathrm{1}}\text{d}{d}_{\mathrm{2}}.\end{array}$$

Usually only numerical solutions exist for β for all sizes, but analytical solutions for certain size regimes exist. The method by Whitby et al. (1991) uses analytic solutions for the integrals in Eqs. (B17) to (B19) for the near-continuum and free-molecular regime. The following line only derive the coagulation rates for inter-modal coagulation. The rates for intra-modal coagulation can be derived analogously.

The coagulation coefficient in the near-continuum is given by

$$\begin{array}{}\text{(B20)}& {\mathit{\beta }}_{\mathrm{nc}}(d_i,D_j)=\mathrm{2}\mathit{\pi }(D_i+D_j)(d_i+d_j).\end{array}$$

with D_l the diffusion coefficient of particles in mode l:

$$\begin{array}{}\text{(B21)}& D_l={\displaystyle \frac{k_{\mathrm{B}}T{C}_l}{\mathrm{3}\mathit{\pi }\mathit{\mu }{d}_l}}\end{array}$$

Here, k_B denotes the Boltzmann constant, T the temperature, μ the dynamic viscosity of air, C_l the Cunningham correction factor, and the particle diameter d_l. The Cunningham factor C_l is approximated by

$$\begin{array}{}\text{(B22)}& C_l\approx \mathrm{1}+\mathrm{1.246}K{n}_l\end{array}$$

with $K{n}_l=\frac{\mathrm{2}{\mathit{\lambda }}_{\mathrm{air}}}{d_l}$ the Knudsen number and λ_l the mean free path of air. The coagulation rate in the near continuum regime can now be written as

$$\begin{array}{}\text{(B23)}& \begin{array}{rlr}\stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{a}}_{\mathrm{0},ij}& =\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\mathit{\beta }}_{\mathrm{nc}}(d_{\mathrm{1}},d_{\mathrm{2}}){\mathit{\psi }}_{\mathrm{0},i}\left(d_{\mathrm{1}}\right){\mathit{\psi }}_{\mathrm{0},j}\left(d_{\mathrm{2}}\right)\text{d}{d}_{\mathrm{1}}\text{d}{d}_{\mathrm{2}}& \\ & ={\mathrm{\Psi }}_{\mathrm{0},i}{\mathrm{\Psi }}_{\mathrm{0},j}{\mathit{\rho }}_{\mathrm{a}}^{\mathrm{2}}{K}_{\mathrm{nc}}{d}_{\mathrm{0},i}^{\mathrm{3}}[\mathrm{2}{e}^{\frac{\mathrm{36}}{\mathrm{8}}l{n}^{\mathrm{2}}\left({\mathit{\sigma }}_i\right)}\\ & +a_{i}K{n}_{gi}\left(e^{\frac{\mathrm{16}}{\mathrm{8}}l{n}^{\mathrm{2}}\left({\mathit{\sigma }}_i\right)}+{\displaystyle \frac{d_{\mathrm{0},j}}{d_{\mathrm{0},i}}}{e}^{\frac{\mathrm{4}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_i\right)}{e}^{\frac{\mathrm{4}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_j\right)}\right)\\ & +a_{j}K{n}_{gj}\left(e^{\frac{\mathrm{36}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_j\right)}{e}^{\frac{\mathrm{4}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_j\right)}\right.\\ & \left.+{\displaystyle \frac{d_{\mathrm{0},i}}{d_{\mathrm{0},j}}}{e}^{\frac{\mathrm{64}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_j\right)}{e}^{\frac{\mathrm{16}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_i\right)}\right)+\left({\displaystyle \frac{d_{\mathrm{0},i}}{d_{\mathrm{0},j}}}+{\displaystyle \frac{d_{\mathrm{0},j}}{d_{\mathrm{0},i}}}\right)\\ & e^{\frac{\mathrm{4}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_j\right)}{e}^{\frac{\mathrm{4}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_i\right)}]\end{array}\end{array}$$

with the factor $K_{\mathrm{nc}}=\frac{\mathrm{2}{k}_{B}T}{\mathrm{3}\mathit{\mu }}$. Analogously, the intr-modal coagulation rates $\stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{a}}_{\mathrm{0},ii}^{\mathrm{nc}}$ and $\stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{a}}_{\mathrm{0},jj}^{\mathrm{nc}}$ are derived. The coagulation rate for the third moment in the near continuum regime reads as:

$$\begin{array}{}\text{(B24)}& \begin{array}{rlr}\stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{a}}_{\mathrm{3},ij}& =\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{d}_i^{\mathrm{3}}{\mathit{\beta }}_{\mathrm{nc}}(d_{\mathrm{1}},d_{\mathrm{2}}){\mathit{\psi }}_{\mathrm{0},i}\left(d_{\mathrm{1}}\right){\mathit{\psi }}_{\mathrm{0},j}\left(d_{\mathrm{2}}\right)\text{d}{d}_{\mathrm{1}}\text{d}{d}_{\mathrm{2}}& \\ & ={\mathrm{\Psi }}_{\mathrm{0},i}{\mathrm{\Psi }}_{\mathrm{0},j}{\mathit{\rho }}_{\mathrm{a}}^{\mathrm{2}}{K}_{\mathrm{nc}}[\mathrm{2}+a_{i}K{n}_{gi}\left(e^{\frac{\mathrm{4}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_i\right)}\right.\\ & \left.+{\displaystyle \frac{d_{\mathrm{0},j}}{d_{\mathrm{0},i}}}{e}^{\frac{\mathrm{16}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_i\right)}{e}^{\frac{\mathrm{4}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_j\right)}\right)\\ & +a_{j}K{n}_{gj}\left(e^{\frac{\mathrm{4}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_j\right)}+{\displaystyle \frac{d_{\mathrm{0},i}}{d_{\mathrm{0},j}}}{e}^{\frac{\mathrm{16}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_j\right)}{e}^{\frac{\mathrm{4}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_i\right)}\right)\\ & +{\displaystyle \frac{d_{\mathrm{0},j}}{d_{\mathrm{0},i}}}{e}^{\frac{\mathrm{16}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_j\right)}{e}^{\frac{\mathrm{4}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_i\right)}+{\displaystyle \frac{d_{\mathrm{0},i}}{d_{\mathrm{0},j}}}{e}^{\frac{\mathrm{64}}{\mathrm{8}}l{n}^{\mathrm{2}}\left({\mathit{\sigma }}_j\right)}{e}^{\frac{\mathrm{4}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_i\right)}]\end{array}\end{array}$$

For the free -molecular regime the approximation of the coagulation coefficient is given by

$$\begin{array}{}\text{(B25)}& \begin{array}{rl}& {\mathit{\beta }}_{\mathrm{fm}}(d_i,d_j)=\sqrt{{\displaystyle \frac{\mathrm{6}{k}_{B}T}{{\mathit{\rho }}_{p,i}+{\mathit{\rho }}_{p,j}}}}\\ & \left(\sqrt{d_i}+\mathrm{2}{\displaystyle \frac{d_j}{\sqrt{d_i}}}+{\displaystyle \frac{d_j^{\mathrm{2}}}{d_i^{\mathrm{3}/\mathrm{2}}}}+{\displaystyle \frac{d_i^{\mathrm{2}}}{d_j^{\mathrm{3}/\mathrm{2}}}}+\mathrm{2}{\displaystyle \frac{d_i}{\sqrt{d_j}}}+\sqrt{d_j}\right).\end{array}\end{array}$$

The resulting coagulation rate in the free-molecular regime becomes for the zeroth moment

$$\begin{array}{}\text{(B26)}& \begin{array}{rlr}\stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{a}}_{\mathrm{0},ij}^{\mathrm{fm}}& =\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\mathit{\beta }}_{\mathrm{fm}}(d_{\mathrm{1}},d_{\mathrm{2}}){\mathit{\psi }}_{\mathrm{0},i}\left(d_{\mathrm{1}}\right){\mathit{\psi }}_{\mathrm{0},j}\left(d_{\mathrm{2}}\right)\text{d}{d}_{\mathrm{1}}\text{d}{d}_{\mathrm{2}}& \\ & ={\mathrm{\Psi }}_{\mathrm{0},i}{\mathrm{\Psi }}_{\mathrm{0},j}{\mathit{\rho }}_{\mathrm{a}}^{\mathrm{2}}\sqrt{{\displaystyle \frac{\mathrm{6}{k}_{B}T}{{\mathit{\rho }}_{p,i}+{\mathit{\rho }}_{p,j}}}}{b}_{\mathrm{0}}\sqrt{d_{\mathrm{0},i}}\\ & [e^{\frac{\mathrm{1}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_i\right)}+\sqrt{{\displaystyle \frac{d_{\mathrm{0},j}}{d_{\mathrm{0},i}}}}{e}^{\frac{\mathrm{1}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_j\right)}\\ & +\mathrm{2}{\displaystyle \frac{d_{\mathrm{0},j}}{d_{\mathrm{0},i}}}{e}^{\frac{\mathrm{1}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_i\right)}{e}^{\frac{\mathrm{4}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_j\right)}\\ & +{\displaystyle \frac{d_{\mathrm{0},j}^{\mathrm{2}}}{d_{\mathrm{0},i}^{\mathrm{2}}}}{e}^{\frac{\mathrm{9}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_i\right)}{e}^{\frac{\mathrm{16}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_j\right)}\\ & +{\left(\sqrt{{\displaystyle \frac{d_{\mathrm{0},i}}{d_{\mathrm{0},j}}}}\right)}^{\mathrm{3}}{e}^{\frac{\mathrm{16}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_i\right)}{e}^{\frac{\mathrm{9}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_j\right)}\\ & +\mathrm{2}\sqrt{{\displaystyle \frac{d_{\mathrm{0},i}}{d_{\mathrm{0},j}}}}{e}^{\frac{\mathrm{4}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_i\right)}{e}^{\frac{\mathrm{1}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_j\right)}]\end{array}\end{array}$$

and for the third moment

$$\begin{array}{}\text{(B27)}& \begin{array}{rlr}\stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{a}}_{\mathrm{0},ij}^{\mathrm{fm}}& =\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{d}_i^{\mathrm{3}}{\mathit{\beta }}_{\mathrm{fm}}(d_{\mathrm{1}},d_{\mathrm{2}}){\mathit{\psi }}_{\mathrm{0},i}\left(d_{\mathrm{1}}\right){\mathit{\psi }}_{\mathrm{0},j}\left(d_{\mathrm{2}}\right)\text{d}{d}_{\mathrm{1}}\text{d}{d}_{\mathrm{2}}& \\ & ={\mathrm{\Psi }}_{\mathrm{0},i}{\mathrm{\Psi }}_{\mathrm{0},j}{\mathit{\rho }}_{\mathrm{a}}^{\mathrm{2}}\sqrt{{\displaystyle \frac{\mathrm{6}{k}_{B}T}{{\mathit{\rho }}_{p,i}+{\mathit{\rho }}_{p,j}}}}{b}_{\mathrm{3}}(d_{\mathrm{0},i}{)}^{\frac{\mathrm{7}}{\mathrm{2}}}\\ & [e^{\frac{\mathrm{49}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_i\right)}+\sqrt{{\displaystyle \frac{d_{\mathrm{0},j}}{d_{\mathrm{0},i}}}}{e}^{\frac{\mathrm{36}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_i\right)}{e}^{\frac{\mathrm{1}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_j\right)}\\ & +\mathrm{2}{\displaystyle \frac{d_{\mathrm{0},j}}{d_{\mathrm{0},i}}}{e}^{\frac{\mathrm{25}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_i\right)}{e}^{\frac{\mathrm{4}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_j\right)}\\ & +{\displaystyle \frac{d_{\mathrm{0},j}^{\mathrm{2}}}{d_{\mathrm{0},i}^{\mathrm{2}}}}{e}^{\frac{\mathrm{9}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_i\right)}{e}^{\frac{\mathrm{16}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_j\right)}\\ & +{\left(\sqrt{{\displaystyle \frac{d_{\mathrm{0},i}}{d_{\mathrm{0},j}}}}\right)}^{\mathrm{3}}{e}^{\frac{\mathrm{100}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_i\right)}{e}^{\frac{\mathrm{9}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_j\right)}\\ & +\mathrm{2}\sqrt{{\displaystyle \frac{d_{\mathrm{0},i}}{d_{\mathrm{0},j}}}}{e}^{\frac{\mathrm{64}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_i\right)}{e}^{\frac{\mathrm{1}}{\mathrm{8}}{\ln }^{\mathrm{2}}\left({\mathit{\sigma }}_j\right)}].\end{array}\end{array}$$

The values for b₀ and b₃ are given by Whitby et al. (1991) with b₀=0.9 and b₃=0.9 for inter-modal coagulation and b₀=0.8 for intramodal coagulation.

Finally a harmonic mean is applied, similar to the procedure for the condensation, to receive the coagulation rate for the full size range. For the conversion to coagulation rate of the number concentration and mass mixing ratio, the density of air and mode density are considered:

$$\begin{array}{}\text{(B28)}& {\mathrm{Ca}}_{\mathrm{0},l}={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{a}}}}{\displaystyle \frac{\stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{a}}_{\mathrm{0},ij}^{\mathrm{nc}}\cdot \stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{a}}_{\mathrm{0},ij}^{\mathrm{fm}}}{\stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{a}}_{\mathrm{0},ij}^{\mathrm{nc}}+\stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{a}}_{\mathrm{0},ij}^{\mathrm{fm}}}}\end{array}$$

and

$$\begin{array}{}\text{(B29)}& {\mathrm{Ca}}_{\mathrm{3},l}={\displaystyle \frac{{\mathit{\rho }}_{p,l}}{{\mathit{\rho }}_{\mathrm{a}}}}{\displaystyle \frac{\stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{a}}_{\mathrm{3},ij}^{\mathrm{nc}}\cdot \stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{a}}_{\mathrm{3},ij}^{\mathrm{fm}}}{\stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{a}}_{\mathrm{3},ij}^{\mathrm{nc}}+\stackrel{\mathrm{̃}}{\mathrm{C}}{\mathrm{a}}_{\mathrm{3},ij}^{\mathrm{fm}}}}\end{array}$$

Appendix C: Database of prescribed aerosol optics

To represent the optical properties of aerosols in the model, we use a set of predefined modes that account for typical particle size distributions and compositions. These modes are based on lognormal distributions, described by their count median diameter (cmd) and geometric standard deviation (σg). In addition to the size distribution, each mode specifies the chemical composition of the particle core and, if present, a coating layer. The coating fraction is defined as the ratio of the coating thickness to the core radius. The data are provided in the model repository and summarized in Table C1.

Table C1 List of available modes, the corresponding parameters of the particle size distribution and the particle composition, for whitch the optical properties are precalculated. The assumed particle size distribution is characterized by the count median diameter cmd and the geometric standard deviation σ_g. “core” and “coating” describe the constituents of the particles' core and its coating, if one is present. The coating fraction (coat. frac.) is the ratio of the thickness of the coating shell and the radius of the core. Refractive indices are taken from Gordon et al. (2022). Data for aerosols with non-sphericity approximation are taken from Hoshyaripour et al. (2019).

Appendix D: Preprocessing for MECCA chemistry simulations

Figure 5 illustrates a schematic of the preprocessing steps required to integrate the MOZART-4 chemistry mechanism into ART. The steps are as follows: Within MECCA (caaba_4.0/mecca), the eqn folder contains various chemical reaction mechanisms formatted for KPP, with the corresponding gas-phase species defined in gas.spc. To activate MOZART-4 chemistry, the xmecca script must be executed using the mozart.bat batch file, which utilizes the KPP equation file eqn/mozart/mozart.eqn. The batch file also allows users to select a subset of the full mechanism – for example, enabling all gas-phase reactions while excluding those involving halogens. Based on the selected mechanism (mecca.eqn) and species list (mecca.spc), KPP generates Fortran90 files (messy_mecca_kpp_*.f90) for the numerical integration of the ODEs derived from the reaction mechanism.

The MECCA preprocessing has been extended to couple KPP routines from the box model to ART. A shell script (create_icon_code.sh), along with an AWK-based routine (write_tracers_to_xml.awk), generates the necessary files and copies them to the corresponding ART directory. Additionally, the tracer XML file mecca_tracers.xml is created. This process uses the table caaba_4.0/mecca/process_gas.tbl to define advection properties: if advection is enabled, the advection template hadv52aero is selected; otherwise, the <transport> tag is set to off. Users can modify or extend this table to accommodate their specific mechanism and setup. The molecular weight of the tracers is computed based on the chemical formulas provided in gas.eqn. The scripts for these steps are provided in the supplementary material.

To integrate the MOZART-T1 mechanism into ART, the following steps should be followed. Using the reactions and rate coefficients provided in the supplementary material of Emmons et al. (2020), users must create an equation file, mozart-t1.eqn, containing the relevant chemical reactions. The syntax for this file is outlined in the CAABA/MECCA user manual available in the manual directory of caaba_4.0. Additionally, the species file gas.spc must be updated accordingly. For convenience, both mozart-t1.eqn and the modified gas.spc file for MOZART-T1 are provided in the supplementary material.

Next, the batch file mozart.bat must be adjusted to select mozart-t1.eqn as the active mechanism. The required KPP files can then be generated as described earlier. Finally, to define the advection properties of tracers in the XML file, users need to extend the table caaba_4.0/mecca/process_gas.tbl by including all MOZART-T1 species with their corresponding advection properties.

Appendix E: XML tags and namelist settings for the Online Emissions Module

The Online Emissions Module (OEM) has its own namelist section oemcntrl_nml, where the paths to the different input files (in netcdf format) need to be defined. Table E1 lists all possible entries. The last column describes in which case the file is required. The file hour_of_year_nc, for example, is only required if for one of the tracers in the tracer XML file the value of tscale is set to 2.

Table E1 List of all oemctrl_nml namelist parameters controlling the input for OEM.

In order to prepare a tracer for OEM, a number of tags need to be defined in the tracer XML file. An overview of all OEM tags is presented in Table E2. In order to be able to add emissions to a tracer using OEM, the tag oem_type needs to be set for this tracer to “emis”. The tag oem_cat lists all the emission categories, which contribute to emissions of this tracer. The tags oem_tp and oem_vp need to be lists of the same length as oem_cat and provide the names of the temporal and vertical profiles that should be used for each emission category. Temporal profiles can either be composed of the product of hour-of-day, day-of-week and month-of-year profiles (if tag oem_tscale=1) or of a single profile with scaling factors for each hour-of-year (oem_tscale=2).

Table E2 List of all XML tags available to control the OEM tracers.

Appendix F: XML tags for deposition of gases

F1 Dry deposition

The dry deposition scheme is available for both parameterized and full chemistry simulations via the tracer XML file by setting the following tags. To activate dry deposition for a specific tracer, one needs to set:

<idep type="int">1</idep>

Additionally, the reactivity of the tracer must be specified using the tag <reac>, which should be a value between 0 and 1. The solubility of the tracer is defined by the tag <heff>, which represents its effective Henry's law constant in units of M atm⁻¹. To compute the resistances representing diffusive processes, the molecular diffusivity of the species must also be provided. The user can specify two values for this: <vdmol1> and <vdmol2>. The tag <vdmol1> specifies the ratio $D_{{\text{H}}_{\mathrm{2}}\text{O}}/D_i$ of the molecular diffusivities of H₂O and the specific tracer, while <vdmol2> encodes the ratio from Eq. (9), more precisely $\mathtt{\text{vdmol2}}=\mathrm{1}/f_D$.¹

F2 Wet deposition

To activate wet deposition for a given tracer, the following tag must be set in the tracer XML file:

<iwash type="int">1</iwash>

The scavenging ratios for in-cloud and below-cloud scavenging are specified with the tags <win> and <wsub>, respectively.

Appendix G: Lifetime-based tracers

In the following equations, k depicts the destruction rates for the given tracer, P the pressure in Pascal, COL the column of the given gas, and z_p the pressure altitude in m.

• TRCH4: parameterization by ECMWF ported to ICON-ART (Schröter et al., 2018)

  $$\begin{array}{}\text{(G1)}& \begin{array}{rl}& k_{{\mathrm{CH}}_{\mathrm{4}}}\left(P\right)=\\ & \left\{\begin{array}{ll}\frac{\mathrm{1}}{\mathrm{86}\mathrm{400}\times \mathrm{100}}& \text{if }P\le \mathrm{50}\\ \\ \frac{\mathrm{1}}{\mathrm{86}\mathrm{400}\times \left(\mathrm{100}\times \left(\mathrm{1}+\frac{\mathrm{19}\ln \left(\mathrm{10}\right)}{\left(\ln \right(\mathrm{20}){)}^{\mathrm{4}}}\frac{\left(\ln \right(P/\mathrm{50}){)}^{\mathrm{4}}}{\ln (\mathrm{10}\mathrm{000}/P)}\right)\right)}& \text{if }\mathrm{50}<P<\mathrm{10}\mathrm{000}\\ \\ \mathrm{0}& \text{if }P\ge \mathrm{10}\mathrm{000}\end{array}\right.\end{array}\end{array}$$

• TRCH3CN: parameterization as described in Fischbeck (2018) and based on Harrison and Bernath (2013) and Singh et al. (2003):

  $$\begin{array}{}\text{(G2)}& \begin{array}{rl}& k_{{\mathrm{CH}}_{\mathrm{3}}\mathrm{CN}}=\\ & \left\{\begin{array}{ll}\mathrm{1}/\mathrm{2}\mathrm{207}\mathrm{520}\mathrm{000}& \text{above tropopause}\\ \mathrm{1}/\mathrm{37}\mathrm{843}\mathrm{200}& \text{below tropopause}\\ \mathrm{1}/\mathrm{1}\mathrm{814}\mathrm{400}& \text{in the lowest nine levels }\\ & \text{above the ocean}\end{array}\right.\end{array}\end{array}$$

• TRCO2: parameterization from (Schröter et al., 2018)

  $$\begin{array}{}\text{(G3)}& k_{{\mathrm{CO}}_{\mathrm{2}}}=\mathrm{2.20}\times {\mathrm{10}}^{-\mathrm{8}}\exp (-{\mathrm{10}}^{-\mathrm{20}}{\mathrm{COL}}_{{\mathrm{O}}_{\mathrm{2}}})\end{array}$$

• TRCO: parameterization from (Schröter et al., 2018)

  $$\begin{array}{}\text{(G4)}& \begin{array}{rl}{k}_{\mathrm{CO}}\left(z_p\right)& ={\displaystyle \frac{\mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{8}}}{\mathrm{1}+\exp \left(\frac{z_p-\mathrm{85}\mathrm{000}}{\mathrm{2800}}\right)}}\\ & +{\mathrm{10}}^{-\mathrm{7}}\exp \left(-{\displaystyle \frac{(z_p-\mathrm{70}\mathrm{000}{)}^{\mathrm{2}}}{{\mathrm{7000}}^{\mathrm{2}}}}\right)\\ & +\mathrm{2}\times {\mathrm{10}}^{-\mathrm{6}}\exp \left(-{\displaystyle \frac{(z_p-\mathrm{45}\mathrm{000}{)}^{\mathrm{2}}}{\mathrm{14}{\mathrm{000}}^{\mathrm{2}}}}\right)\end{array}\end{array}$$

• TRH2O: parameterization based on Brasseur and Solomon (2006)

  $$\begin{array}{}\text{(G5)}& \begin{array}{rl}& k_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}\left(P\right)=\\ & \left\{\begin{array}{ll}\frac{\mathrm{1}}{\mathrm{86}\mathrm{400}\times \mathrm{3}}& \text{if }P\le \mathrm{0.1}\\ \\ \frac{\mathrm{1}}{\begin{array}{c}\mathrm{86}\mathrm{400}\left(\frac{\mathrm{1}}{\exp }\left(\log \left(\frac{\mathrm{1}}{\mathrm{3}}+\mathrm{0.01}\right)\right.\right.\\ -\mathrm{0.5}\left(\log \left(\mathrm{100}\right)+\log \left(\frac{\mathrm{1}}{\mathrm{3}}+\mathrm{0.01}\right)\right.\\ \left.\left.\left.\left(\mathrm{1}+\cos \left(\mathit{\pi }\log \left(\frac{P}{\mathrm{20}}\right)/\log \left(\mathrm{0.005}\right)\right)\right)\right)\right)-\mathrm{0.01}\right)\end{array}}& \text{if }\mathrm{0.1}<P<\mathrm{20}\\ \\ \mathrm{0}& \text{if }P\ge \mathrm{20}\end{array}\right.\end{array}\end{array}$$