The atmospheric chemistry box model CAABA/MECCA-3.0

. We present version 3.0 of the atmospheric chemistry box model CAABA/MECCA. In addition to a complete update of the rate coefﬁcients to the most recent rec-ommendations, a number of new features have been added: chemistry in multiple aerosol size bins; automatic multiple simulations reaching steady-state conditions; Monte-Carlo simulations with randomly varied rate coefﬁcients within their experimental uncertainties; calculations along Lagrangian trajectories; mercury chemistry; more detailed isoprene chemistry; tagging of isotopically labeled species. Further changes have been implemented to make the code more user-friendly and to facilitate the analysis of the model results. Like earlier versions, CAABA/MECCA-3.0 is a community model published under the GNU General Public License.


Introduction
In a previous publication by Sander et al. (2005), we presented the multi-purpose atmospheric chemistry model MECCA (Module Efficiently Calculating the Chemistry of the Atmosphere). Since its release, the code development has branched with different modelers adding and using different new features (Stickler et al., 2006;Xie et al., 2008;Taraborrelli et al., 2009;Riede et al., 2009;Keene et al., 2009;Gromov et al., 2010;Kubistin et al., 2010). We have now merged these additions into the most recent model version Illustration showing that the MECCA chemistry submodel can also be plugged into a global, 3-dimensional model, e.g., a general circulation model (GCM). Likewise, the photolysis submodel JVAL can also be plugged into a global model.
The CAABA/MECCA software has been written for UNIX/Linux operating systems, it does not run natively under MS-Windows. However, it is now possible to run the model inside a virtual Ubuntu Linux machine with the VMware Player software. After the model has run on the virtual Linux machine, the user can access the output under Windows and analyse it there. Information how to obtain the code on a virtual machine is given on the MECCA web page.

CAABA model description
CAABA represents an air parcel (i.e., a box) in the atmosphere. Several scenarios are available for the environmental conditions of the simulation. As an example, the marine boundary layer scenario is shown schematically in Fig. 2. The model setup can be changed via the xcaaba shell script and the caaba.nml namelist file.

Submodels
In addition to the MECCA chemistry, the model calculates exchange with air masses outside of the box and the photolysis of chemical species. These processes are coded as individual submodels, which are connected to CAABA via the Modular Earth Submodel System (MESSy) interface by Jöckel et al. (2005), as shown in Fig. 1. The processes of entrainment, emission, detrainment, and deposition are calculated by the submodel SEMIDEP (Simplified EMIssion and DEPosition). For box model calculations, there is no difference between entrainment from above and emissions from below the air parcel. Both processes can therefore be treated as fluxes into the box. Likewise, detrainment and deposition are both treated as fluxes out of the box. Assuming a well-mixed marine boundary layer of height z mbl (cm), the change in concentration c (cm −3 s −1 ) can be calculated (e.g., Kerkweg et al., 2006b) from the emission or entrainment flux J e (cm −2 s −1 ) as: The concentration change c due to deposition (e.g., Kerkweg et al., 2006a) depends on the concentration c (cm −3 ) and the deposition velocity v d (cm s −1 ): CAABA contains 3 submodels which can provide rate coefficients J (s −1 ), also called "J -values", for photolytical reactions: -The submodel SAPPHO (Simplified And Parameterized PHOtolysis rates) provides a simple function: Here, ϑ is the solar zenith angle. The three parameters a (s −1 ), b (dimensionless), and c (dimensionless) are defined for each photolysis reaction. This method has also been used in previous model versions .
-The submodel JVAL calculates J -values using the method of Landgraf and Crutzen (1998). Its implementation into the MESSy system will be described in a separate publication (Sander et al., 2011).
-Finally, the submodel READJ reads J -values from a file. This option is useful for either applying measured J -values, or offline-calculated J -values obtained from a different model.

Multiple simulations and the steady-state mode
In the base configuration, CAABA/MECCA calculates the temporal evolution of the chemistry inside an air parcel. This is ideal for sensitivity studies analyzing the effect of individual reactions inside a large chemical mechanism. To make measurements that can be compared directly with these model results, it would be necessary to follow an air parcel, e.g. by attaching instruments to a balloon that is drifting with the wind. Usually, however, measurements consist of data originating from different air masses. To permit a comparison of these measurements with model results, we have added two features to CAABA. One is the "steady-state mode", which automatically terminates a simulation when a selected species reaches a steady-state, i.e. when its concentration change is less than a user-defined threshold. A useful application of this mode is to fix the physical boundary conditions and the concentrations of long-lived species (e.g., alkanes, alkenes) to their measured values and then calculate steady-state values of short-lived species (e.g., OH, HO 2 ). A comparison of the simulated short-lived species to their measured values yields insight into how well the reaction mechanism is understood. The second feature that has been added to CAABA provides a user-friendly way to perform multiple simulations. The code loops over several sets of input data, and performs one model run for each of them. For example, the data sets could contain measurements made at different times. For each data set, first the input data for the simulation are read (e.g., measured concentrations and offline-calculated photolysis rate coefficients), and then a model simulation is performed, terminating when a steady state is reached. This is illustrated in Fig. 3a where a box simulation is made for each data set of air-borne measurements. This approach had already been implemented in previous model versions (de Reus et al., 2005;Stickler et al., 2006;Kubistin et al., 2010) and is now available in the current version of CAABA. Although it usually makes sense to combine multiple simulations with the steady-state mode, each of them could also be used independently: For example, multiple runs can be performed for a fixed simulation time (not until reaching steady state) to compare the ozone destruction rates d [O 3 ]/dt at the end of the individual model simulations.

Lagrangian trajectories
With the submodel TRAJECT by Riede et al. (2009), CAABA is turned into a trajectory-box model, simulating atmospheric chemistry along pre-calculated Lagrangian trajectories. In contrast to the steady-state mode (Sect. 2.2) with fixed physical boundary conditions for each simulation, TRAJECT does not aim for steady state and varies longitude, latitude, pressure, temperature, relative humidity, and (optionally) radiation during simulation time along each trajectory ( Fig. 3b). As the trajectory information is provided off-line via an input file to CAABA, its source is variable. It could be a kinetic trajectory model using ECMWF analyses, or even a simple, manually composed data set mimicking the temperature and pressure gradient of a laboratory experiment.
The preference for steady-state mode or trajectory mode depends wholly on the user's intent for the model study.
The steady-state mode may for instance help to generate information about short-lived substances in chemical equilibrium when only longer-lived substances were measured and supplied to the model. In trajectory mode, chemical processes can be studied under varying physical conditions as they occur for instance along atmospheric trajectories.
In addition to the general functionality of a trajectory-box model, CAABA/TRAJECT simulations also offer a method to analyse results from global models. For that purpose, the trajectory-box model is configured for maximum consistency with the global model: -Trajectories are based on the global-model wind fields.
-The global-model chemical mechanism is used.
-Initial tracer concentrations for CAABA/TRAJECT are sampled from the global model.
-Photolysis rate coefficients in CAABA/TRAJECT are scaled to the global-model values.
CAABA in trajectory mode then re-simulates the atmospheric chemistry of the global model during transport. Note that mixing with surrounding air masses, such as diffusion, scavenging, deposition, or convective transport, is currently not considered in CAABA/TRAJECT and is thus the only systematic cause for a deviation in chemistry between trajectory-box model results and global-model results. When comparing the output of the two models along the trajectories, any deviation can then be attributed to mixing (and its indirect effect on chemistry) during trajectory travel time.
Consequently, after defining the concentration at the beginning of the global-model trajectory as transport contribution, the subsequent contributions of mixing and chemistry onto the concentrations of chemical species in the global model can be separated and quantified by the trajectory-box model. For further details and an application example, see Riede et al. (2009).

Stratospheric model simulations
The chemistry mechanism contains reactions for the whole atmosphere, including the stratosphere and mesosphere. These are routinely used when MECCA is coupled to a global 3-D model (e.g., Jöckel et al., 2006). The CAABA/MECCA box model, however, has so far mainly been used for simulations focusing on the lower troposphere. As a new feature, CAABA now offers an easy way to switch to a stratospheric scenario. Typical values for initializing the chemistry at 10 or 20 hPa are available. Further values for other conditions can be added if needed. Photolysis rate coefficients are calculated in the JVAL submodel, based on the vertical profiles of ozone, water vapor, pressure and temperature. These profiles have been updated and extended up to a level of 1 Pa (≈ 80 km altitude) using results from the ECHAM5/MESSy1 model by Jöckel et al. (2006). An application of CAABA/MECCA in the stratosphere and mesosphere has recently been presented by Baumgaertner et al. (2010).

MECCA model description
MECCA is a chemistry submodel that contains a comprehensive atmospheric reaction mechanism. The MECCA version presented by Sander et al. (2005) included: (1) the basic O 3 , CH 4 , HO x , and NO x chemistry, (2) non-methane hydrocarbon (NMHC) chemistry, (3) halogen (Cl, Br, I) chemistry, and (4) sulfur chemistry. A suitable subset of the chemistry mechanism can be selected via the xmecca shell script and the mecca.nml namelist file. Recent extensions of MECCA are presented in the following sections.

Multiple aerosol size bins
Previous versions of MECCA were able to calculate aerosol chemistry in one or two bins of aerosol radius, hereafter referred to as size bins. Modeling the marine boundary layer, these were normally assigned to a submicrometer sulfate mode and a supermicrometer sea-salt mode, respectively. Considering the high demand of CPU time for aqueous-phase chemistry, this is already at the limit of complexity for global models (e.g., Kerkweg et al., 2008). However, to study the interaction of gas-phase chemistry with aqueous-phase chemistry in aerosols of different sizes in detail, more than just two aerosol bins are required. To achieve this, we now allow up to 99 aerosol size bins in MECCA. Each size bin is described by its radius (m) and its liquid water content (m 3 m −3 ). These two parameters determine the number concentration of each bin (m −3 ) in the box model. A lifetime is assigned to each size bin, representing the size-dependent loss rate of particles via sedimentation and dry deposition. Assuming that aerosol production and loss are in steady state, an appropriate production term keeps the aerosol concentration constant during the simulation. Each size bin may be initialized with an individual composition of cations (NH + 4 ) and anions (SO 2− 4 , HSO − 4 , Cl − , Br − , I − , IO − 3 , HCO − 3 ). In case there is a difference between the concentrations of cations and anions, it is assumed to consist of chemically inert ions, e.g., Na + for sea salt. Transfer between the gas phase and each aerosol size bin is calculated using the mass transfer coefficient k mt (Schwartz, 1986;Sander, 1999).
Although there are no direct reactions between aqueousphase species in different aerosol bins, they are coupled indirectly due to the interaction with the gas phase. MECCA with multiple aerosol size bins has recently been used for a comparison to ship-borne measurements of size-resolved marine aerosol (Keene et al., 2009).

Monte-Carlo simulations
The value of each rate coefficient in the chemical mechanism is associated with an uncertainty. Monte-Carlo simulations are an efficient method to analyse how sensitively the model reacts to these uncertainties. A MECCA Monte-Carlo simulation is an ensemble of simulations, in each of which all rate coefficients are varied randomly. The range of model results can be related to the range of the varied input parameters. Such an analysis has been performed with a previous MECCA version by Kubistin et al. (2010) to study the effect of rate-coefficient uncertainties on the production of OH radicals. To give a further example, Fig. 4 shows how O 3 mixing ratios change when rate coefficients are varied.
In each individual Monte-Carlo simulation j , all rate coefficients k i are varied by a Monte-Carlo factor: Here, k MC i,j is the rate coefficient of reaction i used in the Monte-Carlo simulation j . It is defined as the product of the recommended value k i and the Monte-Carlo factor f Here, x i,j is a normally-distributed pseudo-random number centered around zero, and f i is the measurement uncertainty which can be derived from data compilations, e.g. Atkinson et al. (2004Atkinson et al. ( , 2006Atkinson et al. ( , 2007. If not known, a value of f i = 1.25 is assumed. The Marsaglia polar method (see Marsaglia and Bray (1964) and http://en.wikipedia.org/wiki/Marsaglia polar method) is used to transform uniformly-distributed into the normally-distributed pseudo-random numbers x i,j . MECCA uses two options from the submodel MAIN RND to generate uniformly distributed pseudo-random numbers: either the Fortran95 standard subroutine RANDOM NUMBER or the Mersenne-Twister algorithm (Matsumoto and Nishimura, 1998). A disadvantage of the Fortran95 standard subroutine is that different compilers will produce different pseudo-random number sequences of different quality. The Mersenne-Twister algorithm should be used in case identical sequences are required with different compilers.

Mercury chemistry
Atmospheric gas and aqueous-phase chemistry of mercury has been included by Xie et al. (2008). The mechanism development focused on the oxidation of Hg(0) by O 3 , OH, and halogen radicals. These reactions are not only important to explain mercury-depletion events in polar regions, but may also be of significance with respect to the global distribution of mercury (e.g., Steffen et al., 2008;Jöckel et al., 2010). However, it should be noted that the uncertainties of current kinetic laboratory measurements are quite large. Thus, regular updates of the rate coefficients in future versions of MECCA will be necessary.

Chemistry of organic trace gases
The oxidation mechanism for organic trace gases is significantly different from that in previous versions of MECCA. The main developments are briefly described here.
The self-and cross-reactions of organic peroxy radicals R i O 2 are treated following the permutation reaction formalism (Jenkin et al., 1997;Saunders et al., 2003), which is a simplification of the formalism used by Madronich and Calvert (1990). Each individual peroxy radical reacts with a generic RO 2 , which represents the sum of all peroxy radicals: RO 2 = i R i O 2 . A pseudo-unimolecular reaction is then assigned to each individual peroxy radical: with an averaged rate coefficient k i that includes the concentration of the generic RO 2 .
A new condensed degradation scheme for isoprene (2methyl-1,3-butadiene, C 5 H 8 ) has been implemented by Taraborrelli et al. (2009). Replacing the Mainz Isoprene Mechanism (MIM) by Pöschl et al. (2000), the new scheme is called MIM2. It closely reproduces the results of the explicit isoprene mechanism of the Master Chemical Mechanism (MCM, http://mcm.leeds.ac.uk/MCM) by Saunders et al. (2003). The buffering capacity of MIM2 against changes in NO emissions is good and avoids the large biases for nitrogen oxides which are characteristic for MIM. Compared to MCM, MIM2 computes small relative biases (<6 %) for most intermediates and has already been used for global atmospheric simulations Butler et al., 2008). MIM2 is characterized by a low degree of species lumping. Thus, modifications are easy to implement and unlikely to introduce unintentional side effects on other parts of the mechanism, which can occur when changing heavily parameterized schemes. This feature is particularly important since a poor knowledge of the low-NO x isoprene oxidation has been linked to the large model bias for the HO x measurements over a pristine tropical forest . Recent theoretical and experimental advancements in isoprene oxidation chemistry may explain those discrepancies (Peeters et al., 2009, and others), and MIM2 is suitable for an implementation of these results (e.g., Stavrakou et al., 2010).
The chemistry of propane (C 3 H 8 ) and n-butane (C 4 H 10 ) has been taken from the MCM and simplified. The OHinitiated oxidation of ethene is taken from Orlando et al. (1998), fully representing the temperature and NO x dependence of the product distribution. Ozonolysis of ethene is taken from the MCM and simplified according to the principles described by Taraborrelli et al. (2009).
Finally, the whole mechanism is now mass-conserving with respect to carbon, including CO and CO 2 . It is thus suitable for assessments of carbon-isotopic trace-gas composition and CO production from biogenic sources.

NH 3 chemistry
A mechanism for the chemical degradation of gaseous NH 3 has been added. Even though this sink is of minor importance for NH 3 (the chemical loss is only a few percent compared to deposition and heterogeneous chemistry), the degradation products contribute to the budgets of N 2 O and NO x . The reaction pathways towards either N 2 O or NO depend on the concentrations of NO 2 and O 3 : under low NO x conditions with mixing ratios below 1 nmol mol −1 , NH 3 acts as a source for NO x , otherwise as a sink. The implemented reactions are taken from the condensed mechanism of Kohlmann and Poppe (1999).

Kinetic diagnostics and isotope chemistry
A tagging technique for the MECCA reaction mechanism has been implemented by Gromov et al. (2010). The new sub-submodel MECCA-TAG allows the explicit on-line calculation of the contributions of selected species to the budgets of other species of interest within a complex chemical mechanism. For example, one can estimate the contributions of methane and isoprene to the production of CO, as well as to the budgets of any associated intermediate like HCHO or CH 3 O 2 . In addition, tagging allows the direct calculation of the species' yields, e.g., number of molecules of CO produced per initially reacted methane or isoprene. Considering the complex chemistry, MECCA-TAG substantially facilitates the analysis of the model results, considering the complex chemistry of isoprene oxidation in MECCA (MIM2, see Sect. 3.4).
An extended application of tagging is isotope chemistry modelling, which is now supported in MECCA. Modeling of isotope chemistry requires careful consideration of isotopic composition transfer, specific isotope exchange reactions and kinetic isotope effects. MECCA-TAG offers a comfortable isotope chemistry parameterisation with all necessary input contained in a single user-friendly configuration file (suffix cfg). Only the list of species and isotopes of interest must be specified, without adding or changing any reaction: the parsing routines automatically process the equation (suffix eqn) and species (suffix spc) files of the selected chemical mechanism and generate additional Fortran95 code for CAABA/MECCA. The execution of the sub-submodel is controlled via a Fortran95 namelist. To date, we have implemented stable carbon ( 12 C/ 13 C) and oxygen ( 16 O/ 17 O/ 18 O) isotope configurations in MECCA that have been evaluated for CO and relevant species . Isotope chemistry of other elements, e.g. hydrogen, sulfur, nitrogen, is to be added in the future.

Further changes
To improve the up-to-dateness and the user-friedliness, and to facilitate the analysis of the model results, several minor changes have been applied to MECCA: -Rate coefficients have been updated to recent recommendations and laboratory studies Atkinson et al., 2004Atkinson et al., , 2006Atkinson et al., , 2007. A complete list of chemical reactions, rate coefficients, and references is available in the file meccanism.pdf in the Supplement.
-During the generation of a new chemical mechanism, the user needs to answer several questions interactively about the selected reactions, the numerical solver, and other features. It is now preferred to collect all answers in a batch file (suffix bat) and then produce a new chemical mechanism non-interactively. Similarly, the new replacement files (suffix rpl) allow easy switching between the standard mechanism and a mechanism modified for testing purposes. Both features are explained in more detail in the CAABA/MECCA user manual (see Supplement).
-The kinetic preprocessor KPP, which performs the numerical integration of the chemical reaction mechanism, has been updated to version 2.1 by Sandu and Sander (2006), which provides a greater number of different numerical integrators. To facilitate the installation and to ensure compatibility, the KPP software has been included into the MECCA distribution.
-The new "groundhog day" feature endlessly repeats the solar cycle of the first model day. This can be used for a model spin-up on a fixed day of the year, i.e. neglecting the seasonal cycle.
-The implementation of diagnostic tracers allows the output of reaction rates associated with the selected reactions. Alternatively, the output of all reaction rates can be switched on.

Summary
We have presented the current version of the versatile atmospheric chemistry box model CAABA/MECCA-3.0. Its chemical and numerical flexibility, its modularity, and its portability make it an ideal tool for many applications. The code is publicly available for atmospheric chemistry and climate research. Our future plans for model development include to: -update of the photolysis rate coefficients according to recently measured spectra, -develop statistical methods to analyse the results from Monte-Carlo simulations automatically, and -extend the MIM2 reaction scheme for isoprene and monoterpenes further.