ICON-ART 1 . 0 – a new online-coupled model system from the global to regional scale

1Institute for Meteorology and Climate Research, Karlsruhe Institute of Technology, Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany 2Deutscher Wetterdienst, Frankfurter Str. 135, 63067 Offenbach, Germany 3Steinbuch Centre for Computing, Karlsruhe Institute of Technology, Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany *now at: Deutscher Wetterdienst, Frankfurter Str. 135, 63067 Offenbach, Germany


Introduction
In recent years several global and regional-scale model systems have been developed to take into account the feedback between natural and anthropogenic gaseous compounds and aerosol particles and the state of the atmosphere.At the beginning those model systems mainly covered the global scale in general using a hydrostatic framework.Recent developments of those hydrostatic global chemistryclimate models take into account the dynamical and chemical coupling between troposphere and stratosphere and partly with the mesosphere.Examples of such model systems are ECHAM/HAMMOZ (see Glossary in Appendix A) (Pozzoli et al., 2008a, b), EMAC (see Glossary) (Roeckner et al., 2006;Jöckel et al., 2006Jöckel et al., , 2010)), and WACCM (see Glossary) (Kunz et al., 2011;Smith et al., 2011).In the meantime also regional-scale online-coupled model systems exist (Baklanov et al., 2014).Those regional-scale models are using a non-hydrostatic framework as this is required to resolve the relevant processes on this scale.Examples of such model systems are WRF-Chem (see Glossary) (Chapman et al., 2009) and COSMO-ART (Vogel et al., 2009).Both model systems Published by Copernicus Publications on behalf of the European Geosciences Union.
are based on meteorological models that are applied for operational weather forecast on timescales of a few days by national weather services.
Regional-scale models need boundary conditions for the meteorological as well as chemical and aerosol variables.Most often these boundary conditions are taken from a global-scale model.Here the problem arises that these model systems are inconsistent in model physics and the specification of air constituents and the chemistry involved.The global atmospheric model ICON (ICOsahedral Nonhydrostatic, Zängl et al., 2014) offers the possibility to overcome this inconsistency as it uses a nonhydrostatic framework already on the global scale and allows one-way and two-way nesting in regions of interest down to horizontal grid mesh sizes in the order of kilometers.
Based on ICON, the new model system ICON-ART is currently under development.ART stands for Aerosols and Reactive Trace gases.The extension ART has been previously used in COSMO-ART to study the feedback between aerosol, trace gases and the atmosphere (e.g., Bangert et al., 2012;Lundgren et al., 2013;Rieger et al., 2014).At the final stage of the model development ICON-ART will contain tropospheric and stratospheric chemistry, aerosol chemistry and aerosol dynamics.Moreover, as a fully online-coupled model system ICON-ART will account for the impact of gases and aerosols on radiation and clouds.By this, the feedback between gaseous and particulate matter and the state of the atmosphere will be realized.
This paper describes the basic equations, gives an overview of the physical parameterizations and numerical methods used in ICON-ART and shows results of the first applications.If not stated differently, physical parameterizations (e.g., radiation, microphysics) used for the simulations in this paper are the same as described in Zängl et al. (2014).Section 2 presents a short summary of ICON and the numerical methods used for the tracer transport.Section 3 gives the basic equations for the treatment of gaseous and particulate matter.Section 4 describes the handling of physical and chemical processes realized so far.Section 5 describes the temporal discretization and the methods applied for the coupling of ICON and ART.In Sect.6, results of first applications of ICON-ART are presented: a case study of the global distribution of short-lived bromocarbons, the spatial and temporal distribution of the ash cloud of the Eyjafjallajökull eruption in 2010, and finally an estimation of the global annual sea-salt emission.

ICON
ICON-ART is based on the nonhydrostatic model system ICON which was developed in a joint project between the German Weather Service (DWD) and the Max Planck Institute for Meteorology (MPI-M) as a unified next-generation global numerical weather prediction (NWP) and climate modeling system.The main goals which were reached during the development of ICON are -better conservation properties than in the existing global model systems GME (Majewski et al., 2002) and ECHAM (Stevens et al., 2013), with the obligatory requirement of exact local mass conservation and massconsistent transport; -better scalability on future massively parallel highperformance computing architectures; and -the availability of some means of static mesh refinement.This was subsequently concretized into the capability of mixing one-way nested and two-way nested grids within one model application, combined with an option for vertical nesting in order to allow the global grid to extend into the mesosphere (which greatly facilitates the assimilation of satellite data).The nested domains extend only into the lower stratosphere in order to save computing time.
The dynamical core is based on the set of prognostic variables suggested by Gassmann and Herzog (2008), using fluxform equations for the thermodynamic scalars density ρ and virtual potential temperature θ v .This allows for achieving local mass conservation in a straightforward way.Massconsistent tracer transport is obtained by passing temporally averaged mass fluxes to the transport scheme (see Sect. 2.1).Compared to the hydrostatic dynamical core developed as an intermediate step by Wan et al. (2013), several refinements of the model numerics have been implemented in order to reduce the amount of computational diffusion required for numerical stability.Most importantly, the velocity components entering into the divergence operator are averaged such as to obtain (nearly) second-order accuracy, and upwind-biased discretization are used for the advection of the thermodynamic scalars.Besides imposing some implicit damping on small-scale structures, the latter reduce the numerical dispersion errors compared to second-order centered differences.

Tracer transport
Tracer transport is accounted for in a time-split fashion, i.e., by treating vertical and horizontal transport separately.In the vertical, the finite-volume piecewise parabolic method (PPM) (Colella and Woodward, 1984) is applied, where the tracer distribution in each cell is reconstructed using 1-D parabolas.The specific formulation in ICON is able to cope with large Courant-numbers (CFL > 1), following the approach proposed by Skamarock (2006).
For horizontal transport, a simplified flux-form semi-Lagrangian (FFSL) scheme is used, similar to Miura (2007) and Lauritzen et al. (2011).The basic idea for computing the horizontal flux divergence, is to trace the area that is "swept" through an Eulerian cell edge during 1 time step.
In the current implementation, the swept area is approximated as a rhomboid and the tracer distribution in each cell is reconstructed using either 2-D linear, quadratic or cubic polynomials.The polynomial coefficients are estimated using a conservative weighted least squares reconstruction method (Ollivier-Gooch and van Altena, 2002).The performance of the scheme with first-order (linear) polynomials is documented in Lauritzen et al. (2014).
Specific care is taken to retain tracer and air mass consistency.Firstly, mass fluxes passed to the transport scheme are temporal averages of the mass fluxes computed within the dynamical core (during the solution of the mass continuity equation).Secondly, as part of the time-split approach, the mass continuity equation is diagnostically re-integrated, as proposed by Easter (1993).
The transport scheme preserves linear correlations given that a monotone flux limiter is applied.

Basic equations
In the following the basic equations for the treatment of gases and aerosols in ICON-ART are given.We will use the socalled barycentric mean (indicated by a hat) with respect to the density of air ρ where is a (mass-)specific variable and will be further described in the following sections.A variable with a bar on top is Reynolds-averaged.The fluctuation is given by = − .The total time derivative reads as where v is the barycentric mean of the velocity.The continuity equation is given by dρ dt = −ρ∇ • v. (3)

Gaseous tracers
For gaseous tracers, the scalar variable is given by the ratio of the partial density of gas l and the total density.This results in the barycentric-averaged mass mixing ratio g,l : which will be used in the following equations.Within ICON-ART, the spatiotemporal evolution of gaseous tracers is treated according to the following equation where ∇ • (ρv g,l ) is the change due to turbulent fluxes.The production rate due to chemical reactions is given by P l and the loss rate by L l .Emissions are accounted for by E l (processes are further explained within Sect.4).Applying Eqs.
(2) and (3), Eq. ( 5) can be rewritten in the so-called flux form: where the flux divergence ∇•( vρ g,l ) includes the horizontal and vertical advection of the gaseous compound l.

Monodisperse aerosol
For monodisperse aerosol, the scalar variable is given by the ratio of the mass concentration M l of aerosol l and the total density.This results in the barycentric-averaged mass mixing ratio l : which will be used in the following equations.The balance equation for a monodisperse tracer in flux form is given by where ∇ • (ρv l ) denotes the change due to turbulent fluxes, v sed,l is the sedimentation velocity, λ l the washout coefficient, and E l stands for the emissions flux of tracer l (processes are further explained within Sect.4).

Polydisperse aerosol
Based on the extended version of MADEsoot (Modal Aerosol Dynamics Model for Europe, extended by soot; Vogel et al., 2009), polydisperse aerosol particles are represented by several lognormal distributions.As first example of a polydisperse aerosol, we implemented sea-salt aerosol.
Mass mixing ratio and specific number are prognostic variables whereas the median diameter is a diagnostic variable.The standard deviation is kept constant.We use three lognormally distributed modes for sea-salt aerosol (Lundgren et al., 2013).A list of modes, the according median diameters at emission and the standard deviation is given in Table 1.During the simulation, some processes (namely the diameter dependent) can change the median diameter of a distribution.
As number and mass are prognostic, we require the barycentric averages of the (mass-)specific number 0,l as well as the mass mixing ratio 3,l (the indices 0 and 3 are chosen due to the proportionality to these moments of the www.geosci-model-dev.net/8/1659/2015/1.9 2.0 1.7 distribution).They are formed by using Eq. ( 1) with the ratio of number (mass) concentration N l (M l ) to the total density for the scalar variable .This results in (10)

Number
For lognormally distributed aerosol, the specific number ψ 0,l of mode l with diameter d p is given by where the shape parameters of the lognormal distribution of mode l are the standard deviation σ l and the median diameter d 0,l .As stated before, the median diameter is a diagnostic variable and the standard deviation is held constant during the whole simulation.
For the total specific number of mode l, 0,l , we solve the following prognostic equation where ∇ • (ρv 0,l ) is the turbulent flux of the specific number of mode l, v sed,0,l is the sedimentation velocity of the specific number of mode l, W 0,l denotes the loss of particles of mode l due to wet below-cloud scavenging, and E 0,l denotes the number emission flux of particles of mode l (processes are further explained within Sect.4).

Mass
The mass mixing ratio ψ 3,l of lognormally distributed aerosol of mode l at diameter d p is given by where d 3,l denotes the median diameter of mode l.For the emission scheme, the according median diameter of the emissions, d 3,l,E , is calculated with the following relation (Seinfeld and Pandis, 2006): The according prognostic equation that is solved for 3,l , is given by where ∇ • (ρv 3,l ) is the turbulent flux of the mass mixing ratio of mode l, v sed,3,l is the sedimentation velocity of the mass mixing ratio of mode l, W 3,l denotes the loss in the mass mixing ratio of mode l due to below-cloud scavenging, and E 3,l denotes the mass emission flux of mode l (processes are further explained within Sect.4).

Physical and chemical processes
Within this section, we want to present the physical and chemical parameterizations we use within ICON-ART.We have ordered the processes within this section in the same sequence as ICON-ART computes them.

Emission
Although emissions are rather a boundary condition than a physical process, we decided to include them at this point, as they are the source term for primary aerosol and also (besides chemical production) for gaseous compounds.

Very short-lived bromocarbons
In this study CHBr 3 and CH 2 Br 2 have been included as idealized chemical tracers.Both very short-lived bromocarbons were introduced into the model domain by prescribing a constant volume mixing ratio (vmr) globally for pressures greater than 950 hPa.The vmr values are taken from the WMO Ozone Assessment 2010 (World Meteorological Organization, 2011) and are listed in Table 2.

Volcanic ash
Usually the actual source strength during an ongoing volcanic eruption is not known.In order to overcome this problem there is a need to parameterize those emissions.In the following we will outline in which way the emissions are parameterized in ICON-ART.
The idea of the currently implemented parameterization is based on the experience we have made during the recent eruptions of the Iceland volcanoes (Eyjafjallajökull in April 2010 and Grimsvötn in May 2011).The only information that was available within a short time delay during these events was the height of the top of the plume of the volcano.We assume that this information will also be available in future events.For that reason we derived a parameterization that only depends on the top of the plume height.Stohl et al. (2011) used the method of inverse modeling with a Lagrangian model to derive vertical profiles of the emissions of volcanic ash for the Eyjafjallajökull during the 2010 eruption.Figure 1a shows the vertical profile of the emissions that was derived by Stohl et al. (2011) applying ECMWF data to drive the dispersion model.We normalized the emission values with the maximum value and normalized the height above surface with 12.3 km, which was the height of the top of the plume observed in this case.The latter gives the dimensionless height z : where z is the height in m and h the height of the top of the plume.By this procedure we obtain a dimensionless vertical emission profile that is shown in Fig. 1b.We assume that the shape of this normalized profile is universal.The profile is fitted with a Gaussian distribution.The result is shown by the blue curve in Fig. 1b.The fit-function f e (z ) is given by where a 1 = 0.0076, a 2 = 0.9724, a 3 = 0.4481, and a 4 = 0.3078.
Assuming that we know the total emission E tot in kg s −1 , we can calculate the vertical emission profile E(z ) by where an analytical solution of the integral in the denominator is given by So far we have calculated the emission terms independent of the size class.ICON-ART simulates the size distribution of volcanic ash using a sectional approach with six monodisperse size bins (1, 3, 5, 10, 15, and 30 µm).For distributing the total emission to the individual size classes, we used the observed size distributions close to the Eyjafjallajökull that are reported in Schumann et al. (2011).Based on these measurements we obtain distribution factors f l given in Table 3.
In order to specify E tot , we once more rely on the only information which is available within a short delay of time.Mastin et al. (2009) give a parameterization to calculate E tot as a function of the height of the plume top at the volcano: In this equation, h is given in [km] and E tot results in kg s −1 .With this parameterization, we have all the ingredients to calculate the emission term E l for Eq. ( 8): where V is the volume of a grid cell.The factor f lrt is the fraction of emitted volcanic ash which is available for longrange transport.

Sea-salt aerosol
Emissions of sea-salt aerosol are realized as described in Sect.2.5 of Lundgren et al. (2013).The emission diameters d 0,l,emiss , d 3,l,emiss and the standard deviation σ l of the three lognormally distributed modes are given in Table 1.
The emission fluxes E 0,l and E 3,l that are required for Eqs. ( 12) and ( 15) are calculated according to three different emission parameterizations depending on the mode (Lundgren et al., 2013).For the film mode (Mode A), the emission fluxes are calculated following Mårtensson et al. (2003).This parameterization depends on sea surface temperature and wind speed at 10 m above ground.Emission fluxes of the spume mode (Mode B) are parameterized as a function of wind speed at 10 m above ground as described by Monahan et al. (1986).For the jet mode (Mode C), the parameterization of Smith et al. (1993) is used which describes the emission fluxes also as a function of the 10 m wind speed.

Sedimentation
Within ICON-ART, the sedimentation is treated as an additional vertical advection with an always downward directed vertical velocity v sed,l .For the vertical advection, a finitevolume PPM with a quartic interpolation is used (Colella and Woodward, 1984).A monotonic flux limiter guarantees positive definiteness.The formulation of the vertical advection is Courant-number independent (Skamarock, 2006).

Monodisperse aerosol
The sedimentation velocity for the monodisperse particles v sed,l is based on Stokes law (e.g., Seinfeld and Pandis, 2006): where g is the gravitational acceleration [m s −2 ], D p,l is the particle diameter [m], ρ p is the particle density [kg m −3 ], and C d is the drag coefficient.C c,l is the dimensionless slip correction factor, which depends on the particle diameter and the mean free path of air λ air [m]:

Polydisperse aerosol
The sedimentation velocity of polydisperse aerosol is calculated as described in Binkowski and Shankar (1995) for the zeroth and third moment of the aerosol distribution.By this approach, the size dependence of the gravitational settling is considered.

Turbulence and dry deposition
The turbulent fluxes within ICON are treated by a onedimensional prognostic TKE turbulence scheme (Raschendorfer, 2001).The interface allows for the additional vertical diffusion of tracers.As a necessary boundary condition, a surface flux −w is required by the turbulence scheme: where C d h is the turbulent transfer coefficient for heat, v h the horizontal wind velocity, a the value of in the lowest model layer, and s a value of at the surface.
In the cases of gases and particles this surface flux is determined by the dry deposition process.A commonly used parameterization of this process is based on the dry deposition velocity v dep : with z R = 10 • z 0 where z 0 is the roughness length.The deposition velocity v dep of aerosols is calculated as described in Binkowski and Shankar (1995).Combining Eqs. ( 24) and ( 25) results in an artificial surface value: which can then be used to calculate the surface flux −w in Eq. ( 24).

Washout
One of the major sinks of atmospheric aerosol particles is the washout of particles by rain.We use different parameterizations depending on the description of aerosol (monodisperse, polydisperse) which are summarized in the following.So far, the loss of particles by nucleation and impaction scavenging is not considered.

Monodisperse aerosol
The washout rate of monodisperse aerosol is given by where l is the scavenging coefficient for particles in mode l.With the assumptions, that size and terminal fall velocity of aerosol particles are small compared to rain drops, l for monodisperse particles is given by where D r is the diameter of a rain drop, D p,l is the diameter of particles in mode l, w(D r ) is the terminal fall velocity of rain drops, E(D r , D p,l ) is the collision efficiency between particles of diameter D r and D p,l , and n r (D r ) is the rain drop number density size distribution.
To further simplify Eq. ( 28) we assume that the rain drops are monodisperse.Using the definition of the rain fall intensity where N r is the total number concentration of rain droplets and ρ w the density of water, we obtain Assuming large particles (D p,l > 1 µm ) and a constant representative rain drop size of D r = 5 × 10 −4 m we obtain E(D r , D p,l ) ≈ 1.This results in a scavenging coefficient of where l is given in [s −1 ] and R is given in [kg m −2 s −1 ].The scavenging coefficient calculated by Eq. ( 31) is then used in Eq. ( 27) to calculate the change of l by washout.

Polydisperse aerosol
The removal of polydisperse aerosol by wet deposition is parameterized as a function of the particle size distribution and the size distribution of rain droplets as described in Lundgren et al. (2013).In contrast to Lundgren et al. (2013), we use a precipitation rate depending on height to determine the local droplet distribution instead of using surface precipitation rate over the complete precipitation column.

Chemical reactions
The chemical degradation of the atmospheric trace gases is parameterized by a simplified chemistry scheme.Up to now only two very short-lived bromocarbos (CHBr 3 and CH 2 Br 2 ) were included in ICON-ART.Both substances are depleted due to chemical reactions with OH or by photolysis with no chemical production terms.As the total tropospheric chemical lifetimes for both species are known these chemical lifetimes are recalculated into a total loss rate being used in the balance equation Eq. ( 6).

Subgrid-scale convective transport
The contribution of transport and mixing by subgrid-scale convection on the temporal evolution of the tracer concentrations is an addition to the advection term in Eqs. ( 6), ( 8), (12), and (15) which is necessary at coarse grid mesh sizes ( 1 km).In order to account for this additional transport and mixing, we adapted the parameterization by Bechtold et al. (2008), which is used operationally in ICON and the IFS (Integrated Forecast System) model of ECMWF (European Centre for Medium-Range Weather Forecasts).The parameterization uses a bulk mass flux scheme and considers deep, shallow, and mid-level convection.For a detailed description of the convection scheme we refer to Bechtold et al. (2008).For the convective transport of tracers only shallow and deep convection are considered.

Numerical implementation
Within this section we give an overview of the numerical implementation focusing on the numerical time integration and the coupling structure of the host model ICON with the extension ART.

Temporal discretization
Numerical time integration of Eq. ( 6) and analogously Eqs. ( 8), (12), and (15) can be carried out applying different methods.Following the process splitting concept used for most processes in ICON, we carry out the numerical time integration of the individual processes step by step.That means that the tendency due to a certain process is calculated with the according prognostic state variables that are already updated by the previous processes.Within one integration over time from time level t i to time level t i+1 , several updates due to different processes may be performed sequentially.Starting with (t i ), the time integration scheme that leads to (t i+1 ) is outlined below: where t = t i+1 − t i and the process rates for emissions E l ( (t i )), advection ADV l ( 1 (t i+1 )), sedimentation SED l ( 2 (t i+1 )), turbulent diffusion DIF l ( 3 (t i+1 )), washout W l ( 4 (t i+1 )), chemistry CHEM l ( 5 (t i+1 )), and subgrid-scale convective transport CON l ( 6 (t i+1 )).The term for the subgrid-scale convective transport is an addition to the advection.It describes the temporal change caused by vertical mixing due to shallow and deep convection.This term appears in those cases where the horizontal grid spacing does not allow one to describe the process of shallow and deep convection explicitly.In contrast to the uniform t in Eqs. ( 32)-( 38), a subcycling for the sedimentation process (Eq.34) is carried out using the dynamics time step.Due to stability reasons the implicit Euler solution has been taken for the very short-lived substances (VSLS) tracers in this study.

Coupling ICON and ART
Within ICON, the additional ART modules are integrated in a way that ensures a flexible plug-in of process routines as well as an unaffected ICON simulation in case ART is not used.This is realized by interface modules containing a subroutine with calls to the ART modules.These calls are separated by preprocessor (#ifdef) structures.

MODULE mo_art_example_interface #ifdef __ICON_ART USE mo_art_example, ONLY: art_example
Interface modules are part of the ICON code (e.g., mo_art_example_interface), whereas the routines called by the interfaces (e.g., mo_art_example) are part of the ART code.
The sequence of calls to the different routines in the ICON code is illustrated in blue in Fig. 2. Within one complete integration over time in the ICON model, the dynamics is followed by the call to the tracer and hydrometeor advection calculation.Thereafter, the so-called fast and slow physics processes are called.Saturation adjustment, turbulent diffusion, and microphysics are accounted for as fast physical processes.Radiation, convection, the calculation of the cloud cover, and the gravity wave drag are referred to as slow physics.For stability reasons, a subcycling with a shorter time step is performed within the dynamics.
The ART processes are marked by orange boxes within the ICON-ART sequence in Fig. 2. Directly at the beginning, the emissions of aerosols and gaseous species are calculated.The advection of ART tracers is done within the ICON tracer framework followed by the sedimentation where a subcycling is used for stability reasons.Within the fast physics, turbulent diffusion, washout, and chemical reactions are treated.Finally, within the slow physics, vertical transport due to subgrid-scale convection is performed.Advection, turbulence, and convection are marked by orange frames to illustrate that these are processes that are extended inside the ICON code for the treatment of aerosol particles and trace gases.
The aerosol and gaseous concentrations treated by ART are updated directly within the according process routines or interfaces (see Sect. 5.1).We performed tests with different numbers of cores (powers of two between 64 and 1024) and found roughly a factor of 3 for an ICON-ART simulation compared to an ICON simulation without ART.The ART simulation for this purpose was performed with volcanic ash and sea-salt aerosol switched on.This shows that the scalability of ICON applies also to ICON-ART.

First applications
In this section we present examples of first applications of the model system ICON-ART.According to the processes implemented so far, we focus on very short-lived bromocarbons, aerosol from volcanic eruptions, and sea salt.The forcing of dynamics and transport in these simulations was done by parameterized processes of the NWP version of ICON and namelist parameters were set accordingly.A R2B06 grid (about 40 km horizontal grid spacing; for a detailed description of ICON grids see Zängl et al., 2014) with no nested domain has been chosen with 90 non-equidistant vertical levels up to 75 km together with a time step of 72 s.The vertical thickness of the lowest model layer is 20 m, the maximum thickness of about 2600 m is reached at the top of the model domain.

Simulation of very short-lived bromocarbons
Biogenic emitted VSLS have a short chemical lifetime in the atmosphere compared to tropospheric transport timescales (World Meteorological Organization, 2011).As the ocean is the main source of the most prominent VSLS, bromoform (CHBr 3 ) and dibromomethane (CH 2 Br 2 ), this leads to large concentration gradients in the troposphere.The tropospheric depletion of CHBr 3 is mainly due to photolysis, whereas for CH 2 Br 2 the loss is dominated by oxidation by the hydroxyl radical (OH) both contributing to the atmospheric inorganic bromine (Bry) budget.
Once released active bromine radicals play a significant role in tropospheric as well as stratospheric chemistry as it is in particular involved in ozone destroying catalytic cycles.Although the bromine budget in the stratosphere is dominated by release of Bry from long-lived source gases (e.g., halons) which is relatively well understood the contribution of biogenic VSLS to stratospheric bromine is still uncertain (e.g., Aschmann and Sinnhuber, 2013).
Thus, it is important to simulate the emissions, chemistry and transport of VSLS from the surface to the lower stratosphere reasonably having in mind that due to the short chemical lifetime in the troposphere the transport of CHBr 3 and CH 2 Br 2 is mainly expected in regions of deep convection taking place frequently as, e.g., the tropical western Pacific (e.g., Aschmann et al., 2009).To test the ability of ICON-ART to simulate this fast transport from the ocean surface into the lower stratosphere CHBr 3 and CH 2 Br 2 have been included as idealized chemical tracers.The boundary conditions and the chemical lifetimes are taken from the WMO Ozone assessment 2010 (World Meteorological Organization, 2011).They are recalculated into a destruction frequency for the implicit solution of the balance equation Eq. ( 6) (see Table 2).
The simulation was initialized with data from the ECMWF Integrated Forecast System (IFS) for 1 June 2012, 00:00 UTC.The sea surface temperature was initialized with the skin temperature from the IFS initialization data and kept constant during the simulation.The output was interpolated on a regular latitude-longitude grid with 0.5 • × 0.5 • resolution on pressure levels.
To estimate the ability of ICON-ART in the NWP mode to simulate longer timescales necessary for the diffusion of the very short-lived bromocarbons into the upper tropospherelower stratosphere (UTLS) the zonal mean of the simulated temperature and of the zonal wind are shown in Fig. 3 and compared to ERA-Interim data (Dee et al., 2011) at 1 October 2012.In this model setup ICON-ART is able to re-produce the main characteristics of the reanalyzed meteorology 122 days after initialization in the UTLS.For example, the simulated temperature agrees well with the ERA-Interim data in absolute temperature values in the tropical lower stratosphere as well as in the minimum temperatures within the stratospheric polar vortex in the Southern Hemisphere.A good agreement with the ERA-Interim data is also found for the wind fields revealing that ICON-ART in the NWP mode is suitable for the investigation of the tracer transport of the VSLS from the surface into the UTLS region.
In Fig. 4 the zonal mean of the simulated distribution of CHBr 3 and CH 2 Br 2 at 1 October 2012 is shown.Therein the fixed boundary condition for p ≥ 950 hPa is visible together with the fast upward transport into the UTLS in the tropics.Due to its longer lifetime CH 2 Br 2 is transported quasihorizontally in the upper troposphere into the mid-latitudes and also slightly higher up into the lower stratosphere compared to the about 5 times shorter lived CHBr 3 .
The simulated fast transport into the lower stratosphere occurs mainly in the tropical western Pacific region (see Fig. 5) which is in agreement with previous studies as the preferred region of the transport of VSLS into the lower stratosphere (e.g., Aschmann et al., 2009;Sala et al., 2014).For CHBr 3 the distribution at 150 hPa is more inhomogeneous than for CH 2 Br 2 due to its shorter life time pointing to the regions of fast vertical transport into the lower stratosphere.Consequently, the advection into the mid-latitudes is more visible in the longer lived CH 2 Br 2 compared to CHBr 3 .
Zonal mean profiles of the simulated CHBr 3 and CH 2 Br 2 averaged between 20 • S and 20 • N are compared to mean tropical observations in Fig. 6.The mean observations of CHBr 3 and CH 2 Br 2 are based on a compilation of data from different projects and campaigns on different platforms: the CARIBIC (see Glossary) project between 2009 and 2013 (Wisher et al., 2014), the campaigns TRACE-A in 1992, STRAT in 1996, PEM-Tropics in 1996and 1999, ACCENT in 1999, TRACE-P in 2001, Pre-AVE, AVE and CR-AVE in 2004and 2006, TC4 in 2007, HIPPO-1 to HIPPO-5 between 2009and 2011, SHIVA in 2011, and TACTS/ESMVal in 2012(G. Krysztofiak, personal communication, 2014).The observed mean tropical profiles have been multiplied by 1.70 for CHBr 3 and 1.15 for CH 2 Br 2 , respectively, for an easier comparison with the modeled profiles.This is justified as the boundary value of World Meteorological Organization (2011) differs significantly from the observed value and as the calculation of the concentrations of the VSLS is linearized due to the lifetime approach and therefore depends linearly on the concentration of the VSLS themselves.The ICON-ART results of both brominated substances exhibit the characteristic C-shape profile form and more pronounced for the short-lived CHBr 3 than for CH 2 Br 2 , both also being observed.The volume mixing ratios of about 1 pptv at about 200 hPa (about 11 km) for the longer-lived CH 2 Br 2 is in good agreement with the mean observations as well as other model studies which are in the range of about 0.9 pptv for CH 2 Br 2 .For the shorter-lived CHBr 3 , the observations are in the range of 0.3-1.1 with a mean of about 0.6 pptv (World Meteorological Organization, 2011), and thus slightly lower than the simulated volume mixing ratio.This discrepancy might be caused to a possible sampling bias of the highly variable CHBr 3 in that altitude region due to its short lifetime (e.g., Kerkweg et al., 2008;Liang et al., 2010;Sala et al., 2014;Wisher et al., 2014).

Volcanic ash
The forecast of volcanic ash particles is of great interest for aviation.Moreover, the impact of volcanic ash particles on radiation and cloud properties is of importance for climate change and most probably also for weather forecast.In order to test the capability of ICON-ART to simulate the spatial and temporal distribution of volcanic ash particles, we performed a simulation of the ash cloud of the Eyjafjallajökull eruption in April 2010.This eruption led to a shutdown of civil aviation over large parts of Europe.In response, high efforts have been undertaken to improve the prediction of the volcanic ash plume by deriving the volcanic eruption source strength and vertical emission profile from direct observations (e.g., Stohl et al., 2011).With ICON-ART's predecessor, COSMO-ART, time lagged ensembles have been pro- duced to assess the uncertainties of the volcanic ash forecast due to meteorology (Vogel et al., 2014).Based on both studies, we developed a new parameterization of the source term as outlined in Sect.4.1.2.We performed a 5 days continuous forecast, that means ICON-ART was initialized only at the beginning of the forecast period.The simulation starts at 14 April 2010, 00:00 UTC.We use a R2B06 grid that results in a horizontal grid mesh size of about 40 km and specified the observed emission heights based on Vogel et al. (2014).The emission fluxes for the different size bins are calculated using Eq. ( 21).It is assumed that a significant fraction of the total emitted mass is deposited close to the source due to the gravitational settling of large particles, aggregation of small particles and organized downdrafts.Common values for the ash fraction available for long-range transport lie between 1 and 10 % (e.g., Witham et al., 2012).We chose an appropriate value of f lrt = 0.04 in Eq. ( 21) which lies well in that range.
Figure 7 shows time series of the simulated and the observed number concentrations of particles with a diameter of 3 µm at the meteorological observatory Hohenpeissenberg, Germany.The arrival of the ash plume is captured quite well by the simulation keeping in mind that we have used a relatively coarse grid mesh size and that the forecast lead time is already 4 days when the plume reaches the station.
Lidar measurements of the volcanic ash plume were carried out by Gasteiger et al. (2010) at Maisach, Germany.The quantities derived from those measurements are not directly comparable to modeled variables of ICON-ART.Therefore, we calculated the cross sections of all size bins as a proxy.A comparison of the observed range-corrected signal and the simulated cross section of the ash particles is given in Fig. 8.It is very promising that the model simulations capture some of the main features of the observed plume.This includes the time when the maximum at higher levels occurs, as well as the descent of the plume.The thickness of the simulated plume differs from the observations which can be explained by the vertical grid spacing used for the simulations.In the height range where the maximum occurs the vertical grid spacing is in the order of 300 m.
In Fig. 9, the horizontal distribution of the total volcanic ash concentration (sum of all six size bins) is shown at four reference heights at 16 April 2010, 12:00 UTC.The shape of the ash plume is very characteristic.It spans a horizontally thin band across the northern coasts of France, Germany, Poland, and the Baltic.This band is partly connected to Iceland by an area of low concentrations over the North Sea.Although further investigation and validation is needed concerning the absolute values of the ash concentration, the spatial pattern of the simulated ash plume is in very good agreement with previous studies including model results and observations (e.g., Vogel et al., 2014, Fig. 9;Emeis et al., 2011, Fig. 12;Pappalardo et al., 2013, Fig. 11).

Sea-salt aerosol
Sea-salt aerosol is directly emitted into the atmosphere as a results of the wind stress at the sea surface.The earth's surface is roughly 70 % covered by oceans and sea salt is probably the key aerosol constituent over large parts of the oceanic regions.Consequently, sea salt plays a major role for atmospheric processes from weather to climate timescales.For this reason sea salt has to be taken into account in onlinecoupled weather forecast and climate models.In order to test the emission parameterization used in ICON-ART (see Sect. 4.1.3),a 1-year simulation was performed starting at 29 March 2014 using a R2B06 grid (about 40 km horizontal grid spacing).
The total global mass production rate of sea-spray aerosol varies strongly between different parameterizations as highlighted by the review paper of Grythe et al. (2014).For sea-salt particles with a diameter of less than 10 µm they found a range of 3-70 Pg yr −1 for the global annual emissions.Spada et al. (2013) estimated 3.9-8.1 Pg yr −1 for a size range of 0.1 to 15 µm.Using our parameterization as described in Sect.4.1.3and simulating 1 year, we obtain a global mass production for particles smaller than 10 µm of 7.36 Pg yr −1 .For particles with a diameter below 15 µm, we obtain 10.86 Pg yr −1 .This number is somewhat higher than the range given by Spada et al. (2013).Summing up the total emissions of all three modes contained in ICON-ART, we obtain 26.0 Pg yr −1 .The results are summarized in Table 4.
The horizontal distribution of the annual emissions in kg m −2 yr −1 is presented in Fig. 10.The strongest source regions can be found around Antarctica between 45 and 60 • S and in the northern Pacific and northern Atlantic between 30 and 60 • N.This is in line with recent studies concerning global sea-salt emissions (e.g., Grythe et al., 2014;Ovadnevaite et al., 2014).

Summary
We presented a first version of the extended modeling system ICON-ART.The goal of developing ICON-ART is to account for the interactions of atmospheric trace substances (gases and particles) and the state of the atmosphere within a numerical weather prediction model from the global to regional scale.This first version contains the numerical treatment of the balance equations for gaseous compounds, monodisperse particles and polydisperse particles.We presented these equations as well as the physical parameterizations and numerical methods we use to solve these equations.Two VSLS, CHBr 3 and CH 2 Br 2 , were simulated with ICON-ART.Their volume mixing ratio of about 1 pptv in the tropical upper troposphere as well as the regional distribution with the tropical western Pacific region as the main source region of stratospheric VSLS is in good agreement with observations and previous model studies. in dependence on the plume height which is commonly the first available information during an ongoing eruption.We tuned the emission parameterization based on observations.A preliminary comparison with observed lidar profiles shows that ICON-ART reproduces main features as plume height and temporal development.Moreover, the simulation shows the characteristic shape of the Eyjafjallajökull ash plume as seen in previous publications and observations.Sea salt is treated as a polydisperse aerosol with prognostic mass and number and diagnostic diameter.We conducted a 1-year simulation in order to compare our mass emission fluxes with those derived in other studies.We obtain a total emission flux of 26.0 Pg yr −1 and an emission flux of particles with diameter less than 10 µm of 7.36 Pg yr −1 .This is within the range given in the review paper by Grythe et al. (2014).
The first version of ICON-ART, which is presented in this publication, is the basis for a comprehensive modeling framework capable of simulating secondary aerosol formation and the impact of aerosol on clouds as well as on radiation.

Figure 1 .
Figure 1.(a) Vertical profile of the emissions given in Stohl et al. (2011).(b) Normalized emission profiles.The blue line gives the fit with a Gaussian distribution.

Figure 2 .
Figure 2. 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 (see Sect. 2) but ART tracers are treated inside this framework.The black circle indicates the sequence of the time integration.

Figure 3 .
Figure 3. Zonal mean of temperature (left column) and zonal wind (right column) at 1 October 2012 as given by ERA-Interim reanalysis (top row) and ICON-ART (bottom row).

Figure 5 .
Figure 5. Simulated distribution of CHBr 3 (top) and CH 2 Br 2 (bottom) volume mixing ratio at 150 hPa at 1 October 2012.Please note the different color scales.

Figure 6 .
Figure 6.Mean vertical profiles between 20 • S and 20 • N of CHBr 3 (red) and CH 2 Br 2 (blue) volume mixing ratio simulated by ICON-ART for 1 October 2012 (solid lines) and observed during different campaigns (dashed lines).Please note that, the observed mean tropical profiles have been multiplied by 1.70 for CHBr 3 and 1.15 for CH 2 Br 2 .See text for more details.

Figure 7 .
Figure 7. Simulated and observed number concentrations of particles with a diameter of 3 µm at Hohenpeissenberg, Germany.

Figure 8 .
Figure8.Top: logarithm of range-corrected signal of multiwavelength lidar system (MULIS) at = 1064 nm at Maisach from 16 April 2010, 17:00 UTC to 17 April 2010, 17:00 UTC and from 0 to 10 km above ground; white areas denote periods without measurements (taken fromGasteiger et al., 2010).The thick white line shows the hand drawn border of the top of the ash plume based on the observations.Bottom: simulated cross sections in µm 2 m −3 for the size bins 1, 3, 5, 10, and 15 µm.The brownish line is a copy of the white line drawn in the measurements into the graph of the model results.

Figure 10 .
Figure 10.Horizontal distribution of the mass production rate of sea-salt aerosol in kg m −2 yr −1 .Top: only modes A and B are considered to be comparable to Grythe et al. (2014).Bottom: sum of all three modes, A, B and C.

Table 1 .
Parameters for the lognormally distributed aerosol species.d 0,l,E (d 3,l,E ) is the median diameter of the specific number (mass) emission of mode l.The standard deviation of mode l, σ l , is held constant for the whole simulation.

Table 2 .
Chemical lifetime and boundary condition of CHBr 3 and CH 2 Br 2 in the idealized chemical tracer experiment.

Table 3 .
Schumann et al. (2011)epresentation of volcanic ash emissions within a sectional approach.The diameters of the size bin centers are denoted by d l .The distribution factors f l are taken fromSchumann et al. (2011).

Table 4 .
Annual sea-salt aerosol emissions obtained within this study compared to ranges given by different reviews of existing seasalt aerosol source functions.