Geoscientific Model Development a Kinetic Chemistry Tagging Technique and Its Application to Modelling the Stable Isotopic Composition of Atmospheric Trace Gases

Isotope composition, in many cases, holds unique information on the sources, chemical modification and sinks of atmospheric trace gases. Vital to the interpretation and use of an increasing number of isotope analyses is appropriate modelling. However, the exact implementation of isotopic information in chemistry-climate models is a challenge, and often studies use simplifications which limit their applicability. Here we implement a thorough isotopic extension in MECCA, a comprehensive kinetic chemistry sub-model. To this end, we devise a generic tagging technique for the kinetic chemistry mechanisms implemented as the sub-submodel MECCA-TAG. The technique is diagnostic and numerically efficient and supports the investigation of various aspects of kinetic chemistry schemes. We focus specifically on the application to the modelling of stable isotopic composition. The results of MECCA-TAG are evaluated against the reference sub-submodel MECCA-DBL, which is implicitly full-detailed, but computationally expensive and thus sub-optimal in practical applications. Furthermore, we evaluate the elaborate carbon and oxygen isotopic mechanism by simulating the multi-isotope composition of CO and other trace gases in the CAABA/MECCA box-model. The mechanism realistically simulates the oxygen isotope composition of key species, as well as the carbon isotope signature transfer. The model adequately reproduces the isotope chemistry features for CO, taking into account the limits of the modelling domain. In particular, the mass-independently fractionated (MIF) composition of CO due to reactions of ozone with unsaturated hydrocarbons (a source effect) versus its intrinsic MIF enrichment induced in the removal reaction via oxidation by OH is assessed. The simulated ozone source effect was up to +1‰ in 17 O(CO). The versatile modelling framework we employ (the Modular Earth Submodel System, MESSy) opens the way for implementation of the novel detailed isotopic chemistry treatment in the three-dimensional atmospheric-chemistry general circulation model EMAC. We therefore also present estimates of the computational gain obtained by the developed optimisations.


Introduction
The use of stable isotope information in atmospheric chemistry research has made steady progress, and is expected to continue to do so, driven not only by the upsurge in environmental research, but also due to important analytical developments in isotope measurement techniques.In contrast, modelling of isotope effects has been lagging behind to some degree and is largely ad hoc.This situation is not very satisfactory particularly because many stable isotope studies purport to the bold phrase "constraining the budget" by which is meant "to use stable isotope information for better quantifying the budget calculation of trace gases".For that purpose however, the sufficiently complete representation of stable isotope ratios in atmospheric chemistry models is required.Without proper modelling, the isotope information does not really constrain a trace gas budget, but merely adds another poorly constrained trace gas budget, namely that of the minor isotope labelled trace gas, e.g. 13 CH 4 .Besides the "practical" application of constraining trace gas budgets, stable isotope research can help to improve understanding several processes better.A case in point may be that Published by Copernicus Publications on behalf of the European Geosciences Union.
although satellites observe the tropospheric distribution of CO on a global scale, isotopic composition information allows inferences about the underlying processes and sources driving this distribution.
In this communication, we introduce a generic method for isotope treatment in a comprehensive kinetic chemistry modelling system that is, due to its detail, suitable for investigating the impacts of both, processes and sources.We will use CO as an example due to its rich isotopic information, for 13 C, 18 O and additionally 17 O whose variations in abundance do not follow those of 18 O in a straightforward mass-dependent manner.This may eventually be of considerable importance in studying not only its present budget, but also its changes over the past using ice core material.
Historically, process-and budget-level isotope-enabled atmospheric chemistry models diverged, embodying rather different aspects of the isotope chemistry in a particular modelling application.
The process-oriented studies comprise comprehensive chemistry schemes accounting for isotopes and related isotope effects in full detail.This approach resolves certain peculiarities of the isotopic chemistry, but typically induces computational strain.Firstly, one is compelled to simulate a greatly enlarged chemical mechanism (i.e.set of species and equations, whose numbers grow as the isotopologues and isotopomers are considered separately).Secondly, some of the simplifications commonly used in kinetic chemistry modelling may be not applicable for isotopes, thus limiting usage of reduced chemical schemes (see below).Consequently, the modelling domain is usually limited, e.g. to a column-or a boxmodel, frequently assisting only laboratory work.For instance, Yung et al. (1991) at an early stage applied a onedimensional photochemical model to study the heavy 18 O isotope enrichment in stratospheric CO 2 .They showed that a model calculation with an incorporated complex chemistry mechanism can account for the observed enrichment, and attributed its origin to the chemical exchange with ozone mediated by O( 1 D), supporting the hypothesis of the chemical origin of the effect.In subsequent work (Yung et al., 1997) the model was refined using an expanded chemical scheme that accounted for 17 O, the enrichment of which in O 3 does not follow that for 18 O in the expected mass-dependent manner.The refined chemical mechanism (accounting for two rare isotopes, plus selected isotopomers) grew many times larger.The authors remarked that the calculations could be improved using a 2-D model, yet greater computational accuracy was needed for isotopic species.
Nonetheless, the main complication of isotope chemistry modelling lies in the intricacies of isotope composition transfer.
Although normally not discerned in regular chemistry reaction mechanisms, different isotopes end up in different reaction products, thus forming "isotopically" different kinetic pathways.
For example, there are isotope exchange reactions that affect the isotopic ratios of reacting species, but not their concentrations.
A model study illustrating the significance of oxygen isotope exchange in the atmosphere was done by Lyons (2001), who elaborated the novel mechanism including oxygen isotope exchange reactions.His column-model simulations predicted both, tropospheric and stratospheric ozoneimpelled isotope enrichments for oxygenated short-lived radicals as well as notable mass-independent enrichment for stratospheric water through OH.The latter was shown to originate due to the isotope exchange with NO x , signifying a notable influence of the isotope exchange on the result.Zahn et al. (2006) incorporated this result in an extensive study on multiple isotopes in stratospheric H 2 O, which required a comprehensive mechanism for hydrogen and oxygen isotopes encompassing a large set of kinetic isotope effects (KIEs) and isotope exchange reactions.
They emphasised important peculiarities in isotope modelling, such as inapplicability of the "chemical families" concept.For instance, in common (non-isotopic) cases the conversion of one family member to another does not affect the abundance of the whole family; however, if in the course of the conversion chain the exchange with another isotope reservoir takes place, it alters the isotopic composition of the entire family.The HO x (= HO 2 + OH) family is an example of such, as HO 2 obtains its composition from molecular oxygen during the conversion from OH.Another aspect of isotope chemistry is that the isotopic composition is usually measured as the atomic ratio of the rare isotope to the abundant one, whereas the reaction kinetics are calculated on the basis of molecular concentrations.The removal of a particular isotope atom from the molecule during chemical reaction (for example, D or H during the oxidation of CH 4 to HCHO) alters its isotope ratio differently.Pieterse et al. (2009) refer to the importance of this counting effect in their study on hydrogen isotope chemistry.
The above mentioned complications are, nevertheless, commonly neglected in the budgeting studies.The isotopic composition here serves as an additional constraint on the budget of a species composed of contributions from a set of emission classes, each associated with a given isotopic signature.The modelling domain covers scales from zerodimensional box-models to three-dimensional atmospheric transport models with implicit (i.e.isotopic signatures and effects are approximated from a set of chemical tracer values) or explicit chemistry calculation schemes (e.g.Gros et al., 2002;Cuntz et al., 2003).A large number of atmosphere modelling studies use isotope variations, e.g. in CO 2 , CH 4 or N 2 O, to assess the interhemispheric, troposphere-stratosphere and biosphere exchange (see, for instance, Allan et al., 2001;Boering et al., 2004;Lassey et al., 1993;Liang et al., 2008;Liang and Yung, 2007;McCarthy et al., 2003;Quay et al., 1991).The key features of such studies are the use of relative timescales of the atmospheric transport/mixing processes and isotope ratio variations of a given tracer together with the facility to account for the isotope effects.To illustrate this, the isotopic composition of carbon in methane is relatively easy to parameterise since it has virtually no in situ sources (i.e. as a product of a chemical reaction).Thus, the isotopic composition is easily derived by combining direct (surface) emission, transport and removal terms.Applying this concept in a similar manner for the species which insitu sources cannot be neglected, most implicit-chemistry transport models incorporate simplified isotopic chemistry parameterisations.These account for the main chemical sources and sinks of the species in so-called "net reactions" (without consideration of the intermediates) accompanied by terminal KIEs (i.e.upon entering/exiting a "net reaction" chain).Unfortunately, such a concept is not very well-suited for species that have multiple in situ sources or for "net reactions" involving intermediates whose behaviour varies in different chemical and physical regimes.Recalling the examples above, it is difficult to accurately account for the isotopic exchange of a particular intermediate with other disjoined reservoirs, as well as for the other individual and non-chemical processes that alter the isotopic composition.These may be transport (e.g. the advection term is governed by the chemical gradient of the intermediate, but not the lump of them) or removal processes accompanied by KIEs.
In the exemplary case of CO, its photochemical production, initiated by methane and other VOCs, proceeds via shared intermediates, some of which are removed due to wet scavenging or dry deposition.The methane reaction chain, parameterised with the "net reaction", can be written as (Manning et al., 1997): where the yield λ approximates the effect of the chemistry regime and intermediates removal in the chain.
The estimations of λ are hitherto inconsistent being derived or assumed in a number of calculations (see, for instance, Lelieveld and Crutzen, 1991;Duncan et al., 2007).Bergamaschi et al. (2000) studied this parameter with a three-dimensional transport model incorporating implicit chemistry (i.e. as described above) and the inversion technique, to date the most extensive modelling study on CO carbon and oxygen isotopes.
Since methane carbon possesses a very distinct isotopic signature in the atmosphere, the isotopic composition of CO is very sensitive to its input modulated by λ.The resulting synthesised composition attested rather low yield values, although this result may not be exclusive, since the (unknown) fractionation effects in the chain may explain the CO isotopic signatures even when λ equals unity.Unfortunately, the limitations of the approximation preclude further application of the model results, whereas a detailed chemistry scheme would be of benefit here.Indeed, when the intermediates are treated individually, the yield λ becomes a diagnostic variable since its effect is simulated explicitly.Moreover, tracing the composition of the intermediates, particularly the precursors of CO acting in the methane chain, may shed light upon the role of the potential KIEs there.The presence of the fractionation effect in the chain will cause intermediates to exhibit isotopic compositions different from that of methane, and the model can indicate which species are sensitive to the fractionation, thus to be subjected to further experimental work.
At present, models explicitly calculating chemistry hold the potential of providing detailed isotope chemistry treatment.This has been shown by numerous studies, but not yet applied within atmospheric chemistry global circulation models (AC-GCM) that are capable of handling comprehensive chemistry mechanisms.The advantage of the coupled AC-GCM system is also its capability to account for physical processes affecting isotope ratios individually per species.The holdup in development of such systems can be probably attributed to the large computational demands and painstaking technical implementation of the detailed isotope chemistry.The motivation of the work presented here is to give an isotope extension to the comprehensive kinetic chemistry model that is implemented in the AC-GCM we use, to overcome the computational issues and facilitate the model configuration.The tagging technique proposed here (described in Sect.2) is generic and is well-suited for various kinetic chemistry investigations; it is intrinsically apt for stable isotope modelling as well.In Sect.3, we review the isotope-specific side of the kinetic chemistry modelling with respect to the approximations that are introduced.Some aspects of the adaptation of the tagging technique to stable isotope chemistry modelling are given in Sect. 4. The focus of Sect. 5 is the model evaluation, with illustrative examples followed by a performance analysis.Finally, model configuration aspects are sketched in Sect.6.We turn the reader's attention to the Appendix A where the summary of the used notations is given.

The tagging framework
The tagging technique described here is intended to serve as a diagnostic method applied to simulated kinetic chemistry systems only.The concept of tagging implies that a given characteristic (not necessarily measurable in reality, but an imaginary tag) is assigned to the molecules of a given subset of simulated chemical species.In other words, a "tag" is an additional property of a "tracer" describing the species (as defined, e.g., by Jöckel et al., 2008), to trace certain processes.For instance, one can add a tag of species "origin" to trace its atmospheric transport pathways and contribution of the sources to regional pollution.
Similarly, the information related to the chemical composition interchange can be obtained by tracing the distribution of the "tags" transferred in chemical reactions from one species to another.For instance, one can tag a single species to assess its We first declare the studied kinetic system as being regular.It incorporates the chemical mechanism, being a set of compartments (species) of molecules that interact via different pathways (chemical reactions) characterized by certain transfer speeds (reaction rates).The mechanism defines the resulting numerical system of ordinary differential equations (ODEs) with the species' concentrations as unknowns that must be solved (numerically integrated) to find the evolution of the species' concentrations in time from their given initial values.In this work we use MECCA (Module Efficiently Calculating the Chemistry of the Atmosphere; Sander et al., 2005) as the regular kinetic system and extend it by a flexible tagging procedure.MECCA is an atmospheric chemistry sub-model based on a kinetic pre-processor (KPP, Sandu and Sander, 2006).
KPP provides a set of robust and efficient chemistry kinetic system solvers allowing the simulation of chemical mechanisms of various complexities.MECCA is currently interfaced via the MESSy interface (Modular Earth Submodel System, Jöckel et al., 2005) to the photochemical 0-D (box) model CAABA/MECCA (Sander et al., 2010) and to the 3-D atmospheric chemistry AC-GCM EMAC (ECHAM5/MESSy atmospheric chemistry, Jöckel et al., 2006).The tagging technique for MECCA presented here is implemented as the sub-sub-model MECCA-TAG according to the MESSy standard (Jöckel et al., 2005).
Due to the nature of its application, tagging is not intended to provide additional information on the significance of the various chemical pathways in the mechanism (as, for example, in Lehmann, 2004), which forms a different diagnostic method.Nevertheless, all chemical pathways (i.e.reactions) in the studied mechanism are accounted for to trace the net species' exchange.Although the number of different "tags" as such is not limited, we limit our approach to molecules that have not more than one single "tag" attached.All species' molecules are partitioned this way to disjoint subsets we term tagging classes.Molecules with the same "tag" belong to the same and only class.This allows us to define the (tagging) modelling problem as such that it still deals with quantitative (statistical) characteristics, namely the number of molecules of a certain class, rather than considering the tagging of each individual molecule.By definition, the total abundance of a species is represented by the sum of its molecules of all classes.Consequently, the resulting quantitative information one obtains is the share of each class in the entire budget of every species.

The modelling problem
The tagging technique is designed to meet two criteria that define its applicability in complex and costly simulations.First, it should add minimum possible computational overhead to the entire simulated system.Second, it should be diagnostic (i.e.decoupled) in such a way that the regular system remains unaltered and produces (numerically) identical results while being tagged.However, as we will point out, the simultaneous realisation of these two objectives forms a challenge.For a subset of tagged species one needs to formulate an approach that parameterises the chemical interaction of the molecules of each class in addition to the regular mechanism.The natural approach to this is to form a new kinetic system whose mechanism replicates only the set of necessary reactions, and includes sets of species representing each class of tagged species.To give an example, consider a subset of N species tagged into M classes.To create the kinetic tagging system, we pick out all tagged reactions in which these N tagged species act.Then we replicate these reactions M times in a new mechanism for every set of species of each class.If the number of tagged reactions equals R, the resulting mechanism should include in total N times M species and R times M reactions.Furthermore, by imposing the initial conditions from the regular system and additionally the species' classes distributions, we integrate the tagged system simultaneously with the regular one.
Unfortunately, in most cases it is not sufficient to replicate only the tagged reactions and species.To illuminate this, let us regard two reactions of different order.A unimolecular reaction A k 1 −→ products describes the decomposition of a species A. The reaction rate v 1 in [molecules cm −3 s −1 ] is expressed as and depends only on the concentration [A] and the first-order rate coefficient describes the interaction of species A and B with the rate v 2 which depends on the concentration of both reactants and the second-order rate coefficient k 2 , in units of [molec −1 cm 3 s −1 ]: Now suppose that species A is tagged.It is clear that the rates of those reactions replicating the unimolecular reaction in the tagged system are to be determined by concentrations of the respective "class" of species A. However, a problem arises with the bimolecular reaction when species B is not tagged and hence its concentration is unknown, if only tagged species are replicated.We have to somehow add species B and all reactions it participates in to the new system.This will inevitably cause the same problem with non-tagged species paired in reactions with B, and so on.Finally, one may end up copying a large part of the regular mechanism plus enlarging the tagged part with M times more species and reactions between them.
Another perplexing case arises when both, A and B are tagged.How should one parameterise the way the educts of different classes react and contribute to the common products?Again, a natural way is to consider all possible combinations of possible reactants, class by class.This would require to add another 2 M reactions for each bimolecular reaction of two different species and some M * (M +1)/2 reactions for any self-reacting species.Finally, albeit that the resulting tagging system is diagnostic and decoupled from the regular one, its kinetic mechanism grows to a great extent.For comprehensive mechanisms (which include many reactions and species to tag) the tagged system will become unwieldy, too large, and hence computationally impracticable to be simulated in 3-D models.
Here, we propose and describe an approximate solution to the problem, which permits omitting species that are not tagged.The basic idea of the approach is to utilize the reaction rates v that are customarily calculated in course of the integration of the regular system.Indeed, in the above example, processing of the species A does not require the information about the concentration [B], if the rate v 2 is already known.However, at this stage the specific kinetic system formulation intervenes.The decoupled regular system is inaccessible while it is being integrated; one cannot retrieve the momentary value of v 2 but only its integral over a certain period of time.Nonetheless, as we will discuss, the latter can actually still be used under a number of assumptions.In addition, we will reconsider the species' composition exchange in higher order reactions on the basis of pair exchange for the purpose of dealing clearly with the mentioned transfer ambiguities.This also allows effective optimisations to the final tagging kinetic system.

Definitions and assumptions
We denote a subset of tagged species picked from the regular mechanism as and a subset of reactions between them as where N S and N R are the total number of tagged species and reactions, respectively.Any reaction of R can be unimolecular (thermal decomposition, photolysis or radioactive decay) or multi-molecular (bimolecular or termolecular reactions).For the purpose of tagging, it is important to discriminate how many tagged species are present among the educts (the reacting species).For a reaction with a single tagged educt, the composition transfer can be formulated as: The notation of Eq. ( 1) resembles the common chemical reaction notation, but we mention only the tagged educts and products and apply the reaction rate v r instead of the rate coefficient.So far, we will suppose v r to be a defined momentary value in units of [molecules m −3 s −1 ].
In essence, Eq. ( 1) signifies that in course of the reaction R r all its products P r (i.e., the set of species C p ) are being created from the educt species C e at rate v r .The coefficients s rp designate the stoichiometric coefficients of the respective products.Now consider a reaction of two tagged educts: Here, the molecules of educts C e and C n are transferred with the same rate v r and become distributed among the products.
In fact, Eq. ( 2) can be rewritten as the composition of two single-educt reactions: where each pair of branching ratios η n→p r and η e→p r regulates the transfer of the molecules to the particular product C p from C n and C e , respectively.To maintain the balance in (3) with respect to (2), the sum of the ratios must be equal to unity for each product p: The ratios η are the additional information required for the tagging.For example, one can dispense the explicit creation of a certain product C p from C n by setting η n→p r =1 and η e→p r =0, accordingly.When the exact transmission of the molecules from the educts to the products is not known, the following setting is proposed: = q e q e + q n , η n→p r = q n q e + q n (4) where q n and q e are the number of the selected isotope atoms in the educts.Thus the branching ratios in (4) define a proportional isotope content share of both educts in the products composition.

S. Gromov et al.: A tagging technique and isotope chemistry modelling
We give an illustrative example for the formulations above, using the reaction of quenching and interaction of O( 1 D) with O 2 : First, we rewrite this bimolecular reaction as a composition of two simpler reactions, i.e.
Given the proportional share of the educts composition in the products, the transfer ratios for each educt-product pair, according to Eq. ( 4), are Alternatively, to specify the explicit transfer of the composition from O( 1 D) to O( 3 P), one defines the ratios as In this elementary example, the products' stoichiometric coefficients s rp are all equal to unity.Equation (3) also shows that if C e and C n were species of different classes, each of them would create the products of their own class.As a result: if any higher order reaction of R having multiple tagged educts can be "decomposed" into a set of single educt reactions, we can consider the class to class exchange between particular educt species instead of dealing with combinations of the educts of different classes.This is important for the subsequent approximations that are introduced.
The use of formulation (3) and reaction rates v enables us to reduce the total number of reactions we need to introduce in the kinetic mechanism of the tagged system.The subset R may contain more than one reaction in which a particular pair of species undergoes exchange.Instead of replicating all of them, we can create one pair reaction: in which R e→p is the subset of reactions where species C e produces C p and v e→p is the pair reaction rate.For the comprehensive mechanisms, this step may lighten the numerical complexity of the tagged system by transforming the set of tagged reactions to the equivalent set of maximum 2•N S pair reactions (thus reducing the system of ODEs to solve).
As example for a pair reaction rate we illustrate, similar to the reaction (5), the quenching of O( 1 D) to O( 3 P) with molecular nitrogen: In ( 9) the transfer of the oxygen is unambiguous, thus η=1.The final O( 1 D) to O( 3 P) pair reaction encompassing reactions ( 5) and ( 9) plus assuming the transfer from ( 7) is: The pair reaction rate is always proportional to the concentration of the educt C e (e.g.O( 1 D) in the above example).Subsequently, if Eq. ( 8) were replicated in the tagging system, then the reaction rate of each class' replica is proportional to the share of this class.This result, together with the "class to class" concept from Eq. ( 3), founds the formulation of the tagged system kinetic mechanism that we propose below.
We assume that a subset C of species is tagged into M classes.Thus, in the tagged kinetic system mechanism we introduce sets of species termed as tagging classes m C: such that every regular species C n is represented class-byclass in the tagged system as: where lowercase c n and m c n denote the concentration of the species C n and m C n , respectively.Note that these species belong to the different (i.e.regular and tagged) simulated mechanisms.Further, for every class m and a pair of exchanging species, we define the pair reaction in the tagging mechanism: The rates m v e→p must satisfy the following condition to account properly for the total transfer: This ensures that the amounts of reacted molecules in the regular and tagged systems are equal.Recalling from (*) that the reaction rate is proportional to the concentration of Geosci.Model Dev., 3, 337-364, 2010 www.geosci-model-dev.net/3/337/2010/ the educt, for Eqs. ( 8) and ( 11) we can relate the rates v e→p and m v e→p as where k e→p plays the role of a first-order reaction rate coefficient.In fact, k e→p is a "superposed" rate coefficient of all single-and multi-tagged educt reactions of reacting C e , as resulting from the assumptions we considered.Since there is no difference between the molecules of the regular species C e and any tagging class species m C e , the coefficient k is identical for them.Rearranging Eq. ( 12) with respect to k e→p and retaining from Eq. ( 10) that the concentrations m c e are additive to the concentration c e , we compute the final m-th class rate m v e→p as: which means that we set the reaction rate of a class proportional to its share in the entire budget of the reacting species.Altogether Eqs. ( 3), ( 10), ( 11) and ( 13) define the tagging system kinetic mechanism that uses reaction rate values v r provided by the regular system and additional exchange information in the form of branching coefficients η r .
So far, we considered the reaction rates v r as the instantaneous values obtained from the regular system.However, as we mentioned, when tagging is applied to MECCA (the kinetic chemistry modelling system we employ), these momentary rates are not accessible due to the specifics of the KPP.Solving the regular kinetic chemistry mechanism, KPP can provide only the integral of the particular reaction rate over a given period, namely the integration step.During the integration, KPP advances the solution of the regular system by an internally adjusting time stepping procedure; the appropriate length of the external integration step may vary from seconds to hours, depending on the system simulated.Therefore, the final pair rates (Eq.8) have to be derived a posteriori from the integral value of v r .In a first approximation, the estimated reaction rate A r can be assumed to be constant throughout the integration time step: Here A r is an average value of the regular reaction rate v r over the time step of length τ .Such a primitive approximation ensues from the fact that the exact evolution of the reaction rates in the regular system during its internal integration procedure is not known.To obtain a better approximation or estimate, one needs to account for the specifics of the regular mechanism and kinetic solver, which is beyond of the scope of this work.We will rely on the fact that reaction rates constancy neither violates mass transfer/conservation in the system, nor affects the ratios of totals per class (i.e. total mass of the molecules of all species of a certain class related to the total mass of molecules in the tagged system).Notwithstanding, the use of ( 14) may introduce artificial species' composition (namely class to class ratios) exchanges.These may happen in the case of concurrent significant changes of reaction rates and the reacting species composition during the internal integration.Indeed, the constant reaction rates approximation holds well for cases in which either the species composition or the rates of reactions they are removed by, are approximately constant.To clarify the reasoning underlying this statement, we consider the following example sketched in Fig. 1.
A set of species is tagged into two classes with "blue" and "red" tags.Initially, species B and R possess exclusively "blue" and "red" molecules respectively, whilst species P and E are empty.Now assume that during the integration step the chemical regime changes radically in the middle of it.It happens that B reacts to P during the first half of the integration step, filling it up it with "blue" molecules.Then B stops, but R starts to move to P, colouring it to a certain "purple" mixture of "red" and "blue" molecules.Concomitantly, in the first half of the integration step, P spills into E. What will be the resulting colour of E by the time the integration is finished?It is completely blue.Now suppose that reaction rates do not change, they are assumed to be constant averages.During the integration, B and R simultaneously react to P colouring it straight with purple.What would be the resulting colour of E? In this case, it is purple.Because P changes its colour, the constancy of the rate between P and E makes the difference.Consequently, the constant rate approximation is correct for the transfer of either equal portions of a variable composition or equal composition at variable rates.

S. Gromov et al.: A tagging technique and isotope chemistry modelling
The magnitude of deviation of the approximated system is proportional to the length of the regular system integration step τ .In the ideal case, differences disappear as the integration step becomes so small that the rates in the regular system become virtually constant.However, for the range of applied KPP integration steps (from seconds within the box-model to tens of minutes within the AC-GCM), the corresponding rate changes appear to be still admissible to achieve a necessary satisfactory precision.We derive the actual precision by comparing the tagged system to the reference "doubled" system described below.The crucial point is the adequate choice of the regular system's integration step that should be short enough to adequately circumvent the mentioned difficulties for each particular mechanism dealt with.

Integration of the tagged system
The following sketches the numerical implementation of the tagged system originating from the assumptions formulated above.The systems are integrated in sequence, i.e. the tagged system is being integrated over the same time step in immediate succession to the regular system's integration.During each integration step, the regular kinetic solver advances the species' concentrations in the regular system to the new values and provides the integrals of the reaction rates v r which are used to derive the approximated rates A r .Further, by using A r , the concentrations of the class species m C are advanced within the tagged system according to: We use bold-faced symbols for matrices and vectors.Here m c is the class concentration vector, m f n is the class fraction (share) in the total budget of C n .The indices e, p refer to the educt and product species, r is the reaction index, and m is the class index.For the sake of notation convenience, we use the vector product of the "rate matrix" J of the size N S by N S elements by a vector of species' m-th class fractions m f .Each diagonal element J e,e contains the net reaction rate for a certain species C e , whereas each off-diagonal element J p,e contains the pair reaction rate (Eq.8) for species C e to C p .The matrix J resembles the Jacobian matrix with respect to its structure, since the position and value of each element determines the direction and rate of the molecules' transfer between the species, or its removal.
The resulting ODE system defined by Eq. ( 15) requires considerably less computational operations for its integration compared to a straightforwardly replicated mechanism.The rate matrix J within the current formulation consists of elements that remain constant while the solution is being advanced; hence, during the integration iteration the amount of transferred molecules of each class m is determined only by its fraction m f n .The rate matrix is evaluated only once before the integration starts and used to advance every set of class species.This is particularly useful for tagging into a large number of classes.To compare, the straightforward replication of the regular mechanism (see Sect. 2.2) would require the calculation of at least N R • N S • M reaction rates to advance all classes of tagged species, whereas in the formulation (15) this number reduces to N S • N S , and usually N S <N R .In consequence, the greater the complexity of the simulated mechanism, the greater the expected gain in speed.
Despite its simplicity, approach (15) as a rule inherits numerical complications from the regular system, particularly its stiffness.The implementation of the integration scheme for Eq. ( 15) in the MECCA-TAG includes a selection of different solvers (including the opportunity to plug-in a user-defined integrator) incorporating specifically those schemes that can cope with stiff systems.The latter are based on several approaches (Bloch, 1991;Eriksson et al., 2003;Press et al., 1992) adapted for the needs of tagging, but will not be detailed here.Generally, every solver employed for the tagging integration ensures that the solution of the ODE (15) is asymptotically stable.Ordinarily, different solving methods lead to slightly different equilibrium solutions for tagging, i.e. classes fractions, depending on the precision chosen.That is similar to what one discovers simulating the species concentrations in regular chemistry with different solvers.The regular chemistry (hence, the approximated reaction rates as well) is calculated by KPP for instance with the Rosenbrock solvers that have excellent stability properties (Sandu and Sander, 2006).The tagging integration intrinsically produces (again, within the given precision) the same total concentrations as in the regular mechanism.We find also, that errors in simulated class fractions due to the tagging integration itself are negligibly small compared to those introduced by the rates approximation in Eq. ( 14).Hence, only the regular mechanism integration time step length τ is decisive for the tagging precision.
To recapitulate, our method comprises an optimised kinetic chemistry mechanism tagging technique that is: 1) diagnostic and decoupled from the regular system, 2) simplifies calculations using the reaction rates calculated by the regular system solver, and 3) uses a computationally lighter integration method (constant coefficients ODE system).This method can be applied to a comprehensive chemical mechanism for various tagging purposes.The tagging is diagnostic in the sense that it helps to diagnose the regular kinetic chemistry and does not affect the solution of the regular kinetic system.However, this method itself is based on prognostic equations.

Modelling of isotope-enabled kinetic systems
Modelling of isotope chemistry-enabled systems requires a special treatment of the composition exchange between species as compared to the ordinary molecular case.The specific characteristic of any isotopically differentiated compartment is its isotopic ratio Z R, which is not always equal to the molecular ratio of the isotopologues representing it.Simulating isotopic chemistry requires that the isotopic ratios are correctly altered by the different kinetic processes and transfer between species.
Concerning tagging for simulating isotope ratios one has, apart from the standard chemical interaction, to account for isotope effects (kinetic and equilibrium isotope fractionation).Moreover, it is often the case that for ensuring the correct isotopic signature to be passed on from educts to products in the regular chemical reactions, the transfer of the rare and abundant isotopes has to be parameterised explicitly.To our knowledge, in many approaches, essential peculiarities are inadequately treated in the isotope-modelling literature.Below we discuss some shortcomings along with our approach towards the isotopeenabled chemistry mechanism modelling and its realisation in MECCA-TAG.For the concepts and formulation related to isotopes and isotopic effects, we refer to Appendix B and reviews referenced therein.

Isotope tagging assumptions
We next elucidate the case of isotopic tagging of the chemistry mechanism that involves the separation of all species of interest (those containing the isotope element of the atomic number Z) into a number of classes equal to the number of possible isotopologues considered.Although a more detailed, complete parameterisation is possible (e.g.consideration of isotopomers and multiple substituted isotopologues), the tagging presented here is restricted to the existence of only one isotopologue of a certain isotopic substitution of the element Z; the extension to a larger number of isotopologues (plus isotopomers) can be achieved analogously.Rejection of isotopologues with more than one minor isotope substitution, as well as isotopomeric differences is a simplification applied partly based on the following arguments: -The abundance of the rare isotopes for the abundant elements (H, C, N, O) is sufficiently low to neglect the presence of doubly (or more) substituted molecules.
-Given low abundances, the differences introduced by isotope effects have no significant effect on the chemical system as a whole.In contrast, for chlorine, 37 Cl and -Non-stochastic internal isotopic distribution is neglected (thus excluding special cases such as ozone where the central O atom has a different abundance of 17 O and of 18 O compared to the terminal position of oxygen atoms) since it has a minor effect on the species isotopic composition in terms of kinetic exchange.Nevertheless, the separation of certain important species into different isotopomers (for instance, symmetric and asymmetric O 3 ) can be introduced.On the other hand, the general reduction of isotopomeric differences greatly reduces the number of species (and reactions they act in) to be simulated, and hence computational costs.
Resulting from these assumptions, we set the number of tagging classes M to be equal to the number of considered isotopic substitutions of the element Z, or the number of considered isotopologues.The molecules of different classes contain either only abundant isotopes (giving the major isotopologue) or containing a singly substituted rare isotope (giving the minor isotopologues): where C is a subset of species containing element Z.
The superscripts of C and c indicate the major or any minor isotopologues.Incidentally, any minor isotopologue possessing more than one atom of the selected element Z is a molecule containing always one rare isotope of Z, with the remainder of the isotopes being the abundant ones.Thus, for species of the same isotopic composition (i.e.isotopic ratio Z R) the molecular fractions of the minor isotopologues are proportional to the number of the atoms q n of the element Z these species contain.For example, for isotopic carbon in CO (q CO =1) and C 5 H 8 (q C 5 H 8 =5) of the same ratio 13C R, the fraction of molecules containing rare 13 C differs by a factor of five.

Isotope composition transfer
In the kinetic system, the transfer of the composition between species in course of the reaction must occur without violation of the mass conservation.Normally, reaction stoichiometry for Eq.(1) satisfies the following condition: where q e and q p are the number of the atoms of a given isotope Z in the educt and products used to define the For the correct mass transfer, Eq. ( 2) must hold, if the weightings q are substituted by the number of atoms of any element Z in the respective species.Moreover, for the correct transfer of isotopic composition, the respective quantity of the rare and abundant atoms must be preserved as well, in order to preserve the isotopic ratio of the transferred portion of atoms.Recalling the definition of the major and minor isotopologues, we state the correct system representing Eq. ( 1) for each isotopologue as: Here, superscripts of C indicate major and i-th minor isotopologues, respectively; i is the minor isotopologue index (1≤i≤M−1).Similarly to Eqs. ( 2) and ( 3), any high-order reaction can be represented as a composition of multiple reactions of the kind of Eq. ( 17).So far, we have not yet accounted for any (kinetic) isotope effect in Eq. ( 17).
The first term in squared brackets describes the probability of the rare isotope to be transferred to the current product, while the second term describes the creation of the major isotopologues from the abundant isotopes of the educt, eliminating discrepancies caused by neglecting of the doubly (or more) substituted isotopologues.Indeed, in case q p is smaller than q e the term (1−q p /q e ) is positive and accounts for the major isotopologue of C p constructed from excessive abundant atoms left after the decomposition of C e and redistribution of its atoms to the minor isotopologues.When q p equals q e the term cancels, thus describing the straight transfer of the composition.In the case of q p being larger than q e the mentioned term has a negative sign, physically describing the amount of molecules that has to be taken from the major isotopologue pool to complete at least one valid minor isotopologue of the C p species.The "probability" q p /q e in all cases is greater than zero, meaning that the minor isotopologue is created unconditionally.
The system (17) preserves both the mass and the isotopic ratio in course of the reaction: the number of transferred rare isotopes is conserved, whereas the deficiency or excess of the abundant atoms is eliminated at the expense of the major isotopologue pool (see Appendix C for the calculation).Virtually, we perform the repartitioning of the transferred reacted atoms to the valid isotopologues; such scheme satisfies the requirements introduced in Sect.3.1.It is essential (for any isotope kinetic chemistry modelling system) to verify that such repartitioning of the composition is parameterised correctly, otherwise the system either starts to violate mass conservation or introduces inadmissible conversion of the rare isotope atoms into abundant or vice versa.Obviously, formalism (17) relies on the assumption that the probability of each reactant atom to constitute every product is equal; diverging cases should be considered individually (e.g. by applying a concept analogous to that introduced for Eqs. 2 and 3).
We provide some examples to clarify the above.Consider an oxygen isotopic tagging case (using 16 O and 18 O) in a reaction of a peroxy molecule (HO 2 radical) to give two single oxygen-bearing molecules (OH radicals): The transfer of the oxygen composition from the HO 2 to OH reservoir is: Thus, each minor HO 2 isotopologue creates both, major 16 OH and minor 18 OH.Similarly, consider a more complex case of the isotope tagging of the reaction of O( 1 D) with O 2 (Eq.5, decomposition (6) and transfer conditions (7) from Sect.2.3): At first glance, the stoichiometric coefficients in Eq. ( 18) look odd; nevertheless, they are a result of the combination of the branching ratios η and the mass/isotopic ratio balance condition (Eq.17).One can ascertain that the number of rare and abundant isotope atoms on the left and right side are equal.The negative sign of the term in the last equation means that we take two 16 O isotopes from the major isotopologue pool in form of one 16 O 2 molecule in order to build every two 18 O 16 O isotopologues.

Kinetic isotope effect
Due to changes in physical properties caused by isotopic substitution, the reaction rates of the minor isotopologues differ (generally only slightly) from those of the major isotopologues.These kinetic isotope effects (KIEs) in irreversible reactions can be attributed to various causes, like differences in zero point energy etc., and are often quantitatively given as fractionation factors α (see Appendix B).Note that the definition we use deviates from the IUPAC-recommendation as it is intuitively preferable and Geosci.Model Dev., 3, 337-364, 2010 www.geosci-model-dev.net/3/337/2010/uses the rate of the minor isotopologue over that of the major one.Accounting for KIEs, Eq. ( 17) is modified for the minor isotopologue: where the e,i α r designate the fractionation factors for the i-th minor isotopologue of C e , respectively.Each minor isotopologue in Eq. ( 19) may create both, minor and major isotopologues which belong to the different tagging classes.Such would require essential changes to the parameterisation (Eq.15) since previously disjoint classes start to exchange their composition.Moreover, the KIE alters the reaction rates such that they are no longer strictly proportional to the respective class' share, as we derived in Eqs. ( 12) and ( 13).To resolve these complications, we further reformulate Eq. ( 15) in terms of exclusive atomic tagging (i.e. the consideration of atomic fractions of the rare and abundant isotope atom in the species composition), which makes calculations easier.Nonetheless, the rate of the reaction of each isotopologue is determined by its molecular fraction.

Isotope exchange reactions
A minor isotopologue of one species can exchange its rare isotope with the abundant isotope of another species major isotopologue (e.g.H 18 2 O + 16 OH − − H 16 2 O + 18 OH).When this exchange reaction is sufficiently fast, isotopic equilibrium is obtained.In terms of isotopic tagging, the onward (i.e.left to right in the example reaction above) isotope exchange is the conversion of the minor (H 18 2 O) to major (H 16  2 O) isotopologues of one reagent and vice versa (i.e.major 16 OH to minor 18 OH) for another at a given rate.Noteworthy, the backward (conformably, right to left) isotope exchange may have a different rate (compared to the onward reaction), owing to the KIEs.Therefore, the general description of an isotope exchange reaction is: The first equation in (20) introduces the reference reaction rate v r to which the pair of isotope exchange reaction rates is related to; a, b are the indices of exchanging species and a,i β r , b,i β r denote the kinetic fractionation factors for the i-th minor substitution of C a and C b in the pair.The ratios λ incorporated in the reaction rates are analogous to those introduced in Eqs. ( 3) and ( 4), but describe the probability of one species' minor isotopologue to transfer its rare atom thus creating the minor isotopologue of the partner.The particular derivation of λ a→b r shown above relies on the assumption that the isotopomeric differences do not play a role in the exchange; every rare isotope of C a has equal probability to be exchanged with each abundant isotope of C b .In some studies (for example, by Johnston et al., 2000) the onward reaction rate is explicitly reported; then for the backward reaction, λ and β are expressed relatively to the onward reaction rate with λ r = λ b→a r /λ a→b r and β r = b,i β r − a,i β r , respectively, while for the onward reaction they both equal to unity.
The rates of each reaction in Eq. ( 20) are determined by the abundance of the different (major and minor) isotopologues of both species and thus cannot be correctly approximated in the tagged system as being linearly proportional to the reference reaction rate v r .The use of Eq. ( 13) would slightly (since maj c a min,i c a ) change the reaction rates, thus introducing an additional kinetic effect when product λ • β differs for onward and backward directions.In fact, the linearisation (Eq.13) can be used in case the mentioned artificial KIE is insignificant compared to the effect of the species composition exchange in the ordinary reactions.Nevertheless, further we apply proper adjustments to the tagging scheme to treat isotopic exchange in a non-linear way.
To illustrate the formulation (20), consider an example of the oxygen isotope exchange between OH and NO 2 : The probability of 2/3 in the onward exchange reaction rate of Eq. ( 21) signifies the statistical expectation of the 18 O isotope exchange between OH and NO 2 .Conformably, the backward reaction probability is by a factor of two smaller (1/3), because when minor N 18 O 16 O swaps one of its oxygen atoms with 16 OH, we expect that only in half of the cases the rare 18 O isotope becomes exchanged.

Conservative characteristics of an isotopic system
To assess whether the parameterisation of the described effects in the chemistry isotope modelling system is adequate, we consider some of its invariable characteristics.
The following features must be observed in a closed system (i.e. when the mechanism is isolated so that no removal or introduction of the molecules from outside of it is allowed, and reactions have no misbalanced stoichiometry): 1. the total atomic mass and also the ratio of the rare to abundant isotope(s) masses is preserved; 2. the previous characteristic is insensitive to the presence of KIEs and isotope exchange reactions, changes of macro parameters (e.g.pressure or temperature), and changes of the species abundances; 3. if there are no kinetic isotope effects introduced in the system, and all compartments have the same initial isotopic composition, then no isotope ratio change is expected for any species despite alterations in its abundance.
If these characteristics are not reproduced, the system obviously introduces improper conversion of the rare isotope atoms into abundant ones or vice versa.The third condition points directly to incorrect isotopic transfer in the regular reactions.
Note that usually the isotope ratio in the system is a very sensitive characteristic.For example, our reference simulations with isotope-enabled MECCA mechanisms show that for a typical ratio of the stable carbon isotopes corresponding to 0‰ (with respect to the carbon V-PDB 1 standard ratio), a deviation of the δ 13 C value of the total carbon in the system (characteristic 1) larger than merely 10 −10 permil (equivalent to the concentrations calculation precision of absolute and relative tolerance values of 10 molecules/cm 3 and 10 −3 , correspondingly) is a sign of potential error.

MECCA-TAG isotopic modelling specifics
The following section is devoted to extending the molecular tagging technique we introduce to account for the isotope chemistry specifics outlined above.We describe also some additional processing of the species budgets required for diagnostic tagging when it is applied to modelling systems like MECCA.Unless stated otherwise, the notation follows the one used in the preceding sections.

Alterations to the tagged system formulation due to isotopic effects
The changes to the tagged system deal with the additional processing of the classes representing minor isotopologues and their interaction with the major isotopologues class.
To formulate the isotopically tagged system, we recur the considerations described in Sects.2.3 and 2.4 and the specifics of the isotopic transfer from Eq. ( 17).For the major isotopologue containing only abundant isotope atoms, the integration proceeds in the same way as in Eq. ( 15), using its molecular fraction maj f n : 1 The Pee-Dee Belemnite (abbreviated as PDB or V-PDB) standard ratio of 1123.72×10−5 is used for δ 13 C reporting (see Appendix B).
The rate matrix terms include the branching ratios (Eq.4) derived with respect to the atomic content of the selected isotope and composition transfer.
Using the already calculated "major" rate matrix J, we form the "minor" matrices whose elements include the correction terms that account for the kinetic isotope effects for the minor isotopologues: min,i J p,e = J p,e , C p ,C e / ∈ C KIE , else The C KIE and R KIE refer to the subsets of species and reactions experiencing kinetic isotope fractionation, respectively.The correction terms for the production and loss elements are of opposite signs to properly account for the absolute value of the current reaction rate contributing to the corresponding production and loss "one-way" rates.The correction terms for each reaction are calculated just once before the integration.
The elements of the rate matrices are referring to the molecular exchange rates, i.e. the number of molecules of the particular class transferred per unit of time.Recalling the isotopic composition transfer rules (Sect.3.2) we need to consider the transfer of the rare and abundant isotopes from the minor isotopologues separately.For that, we introduce rare and abundant atom fractions χ of the minor isotopologue that are proportional to its molecular fraction and atomic content: where q n is the number of atoms of the selected isotope in the molecule C n .Thus, for instance, for the oxygen composition of NO 3 (q n =3) we account for one rare ( rare χ=1/3) and two abundant ( abu χ =2/3) oxygen isotopes in the minor isotopologue.According to the definitions (Sect.3.1), one transferred rare isotope implies that one minor isotopologue has ended up in the product, a substitution of the molecular fraction of the minor isotopologues with the atomic one properly accounts for the transfer of minor isotopologue molecules: The bracketed term in Eq. ( 25) accounts for the redistribution of the abundant atoms in the products in order to create valid isotopologues: the positive part equals the total amount of abundant isotopes reacted via minor isotopologues, while the negative denotes the fraction that is taken to complete the minor isotopologue molecules.If the set of tagged species consists only of those having one isotope, the mentioned term cancels out and Eq. ( 25) becomes identical to the one for the molecular system (Eq.15).Finally, we need to account for the isotopic exchange as it is the pertinent process inducing interchange between major and minor classes.From Eq. ( 20) one derives the net rate of the major-to-minor isotopologues conversion of the exchanging species: For the pair of reagents C a , C b exchanging in the subset of reactions R IEX , the conversion rate I has the same magnitude; indeed, the net effect of the isotope exchange reaction is equal to the effect of the process when one reagent swaps its rare atoms with abundant ones of the partner until species are in isotopic disequilibrium (i.e.I is non-zero).
Summing the contribution of all isotope exchange reactions in which a certain species participates, we form isotope exchange component vectors i γ : where R n IEX is the subset of isotope exchange reactions of the species C n , p is the index of a particular reaction r exchange partner.Every element of i γ n equals the net rate of onward (or backward, when negative) "conversion" of the major isotopologue maj C n into the minor min,i C n due to all isotope exchange reactions in which species C n takes part.
Further, i γ is added to Eq. ( 25): The resulting system (28) is integrated in the same way as (15) inheriting its optimisation features.Some extra work is required to process the minor isotopologues' abundant isotopes' transfer term and the isotope exchange (i.e.calculation of the vectors i γ in every iteration) although the latter is much less expensive compared to the processing of the ordinary kinetic exchange.

Budget corrections
In case there are kinetic isotope effects introduced in the mechanism, different reaction rates for certain isotopologues will result in a slight but systematic deviation of the tagged species budgets in the tagged system compared to those in the regular mechanism.To avoid integration over this propagating divergence the correction of the species budgets is necessary.It is performed by scaling the class species budgets to the total budget of the regular mechanism, while the newly obtained isotopic composition remains unaffected: The arrows indicate that the left-hand-side is re-defined.
The correction is performed after the integration of both, regular and tagged, systems is finished at the time (t 0 +τ ).The fractions maj f n and min,i f n are calculated as given in ( 22) and ( 24 ) from not yet corrected values of maj c n and min,i c n .One could reformulate the KIE fractionation factors to achieve the same products' composition as from Eq. ( 29) without the correction.This may be done by assigning the following fractionation factors to each (including the major) isotopologue reaction rate: Here, the notation follows that of Eq. ( 19) and e,maj α r and e,min,i α r denote new fractionation factors for the major and minor isotopologue(s) reactions.The amount of produced molecules equals to those in the regular mechanism, whilst the effect on the reagents' ratios change is of the same KIE magnitude of e,i α r .Unfortunately, a reformulation like Eq. ( 30) involves recalculation of the corresponding www.geosci-model-dev.net/3/337/2010/Geosci.Model Dev., 3, 337-364, 2010 pair rates (elements of J) for every isotopologue class in every iteration of the integration of ( 28), which may increase computational effort considerably.Since for typical integration step lengths (see Sect. 2.3) the discrepancies to be eliminated by Eq. ( 29) do not significantly change the chemistry, we choose the a posteriori option (Eq.29) which is obviously faster and does not interfere with the main function of the isotopic tagging, namely to transfer and preserve the isotopic ratios.

Accounting for a specific composition source
It is common practice in kinetic chemistry modelling, that complex reactions involving several intermediate steps are approximated (or "lumped") using one net reaction of a certain rate by omitting (usually short-lived) intermediates.
Nevertheless, the latter can exchange their composition with other species in course of the reaction chain, which is of importance for the isotope chemistry and thus needs to be accounted for.For example, the reaction of methane oxidation may be parameterised in the regular mechanism as (O 2 is omitted among the reactants): which is a superposed sequence of two specific reactions, namely: Virtually, k net equals k 1 since the rate-determining step of this sequence is the reaction of methane with the hydroxyl radical (k 1 k 2 at typical atmospheric concentrations of OH and O 2 ).The methyl radical reacts with molecular oxygen, which gives its signature to the methylperoxy radical, whereas the water molecule incorporates its oxygen atom from OH.As a result, atmospheric oxygen resets the CH 3 O 2 signature, thus determining the composition of all its subsequent products, whilst OH contributes to the composition of the atmospheric water (Dubey et al., 1997).
To account for such a particular sourcing in the tagged system (Eq.28), a simple change of the corresponding production term in the rate matrix J in Eq. ( 22) is required.It is noteworthy that the source specification violates the isotopic composition transfer balance; it is introduced to account for the correct isotopic source signatures in those comprehensive mechanisms that include reactions of high complexity.At atmospheric conditions in the presence of highly abundant species, like atmospheric water or molecular oxygen, this approach is a valid approximation.For any precise setup (e.g. a laboratory experiment reproduction within a box-model) it is adamant that all intermediates are accounted for.

Doubling technique
As an alternative to the diagnostic tagging, we introduce in addition a reference method (hereafter referred to as "doubling") implemented in the MECCA-DBL subsubmodel.The doubling is -alike tagging -a diagnostic, but not a decoupled method.It embodies the expansion of the regular mechanism by additional doubled species for each tagged regular species together with additional reactions they act in.In the explicit doubling setup, the original species and reactions are replaced by the doubled ones.In the implicit setup, additional reactions become virtual.In other words they include only the doubled species thus leaving the regular mechanism reactions unchanged.Other, non-tagged, educts in virtual reactions are replicated as reaction products thus becoming catalysts: only the doubles are produced or consumed.The high-order reactions are decomposed using Eq. ( 3) to a set of single-order reactions.To account optionally for a particular composition source (see Sect. 4.3) in the doubled system, a set of additional reactions with the rate proportional to the source species' class fraction is introduced.In the case of isotopic doubling, the isotopespecific transfer and exchange reactions are formulated exactly as in Eqs. ( 17) to ( 20).Similar to the tagging case, the implicitly doubled mechanism requires the doubled species budget correction (see Sect. 4.2) when kinetic isotope effects are introduced.The explicitly doubled mechanism replaces the original one and thus does not need to be corrected for KIE influence.
Both doubling setups yield the same result (within the numerical precision), whereas the implicit doubling introduces no changes to the regular reactions and allows simultaneous doubling of several different configurations (i.e.given sets of tagged species, reactions and number of classes).
In addition, the implicit setup enables a comparison of the results obtained by doubling and tagging methods implemented within the same simulated chemical mechanism (which is more preferable due to the independence of the result on the numerical precision of the simulated system).Tables 1a and 1b present examples of implicit and explicit doubling, respectively, of the regular mechanism for the isotopic tagging described above.
The doubling usually enlarges a comprehensive mechanism extensively.
We use this method (which is computationally very expensive) as the reference, which is used to assess the adequacy and accuracy of the tagging approach.Nevertheless, for less resource-demanding applications (e.g.box-or column-model), doubling extends MECCA to a fully fledged isotope-enabled kinetic chemistry model.

Simulations performed with isotope tagging setups
To verify whether both, the tagging and the doubling approaches, can correctly reproduce the isotopic composition, several simulations have been performed with the CAABA/MECCA box-model for the carbon-12/13 and oxygen-16/17/18 tagging configurations.As criterion of proper reproduction by tagging, the absolute differences of the simulated isotopic ratios of selected species for tagging and doubling were required to be within Z R=10 −9 .For carbon, this is equivalent to a change of less than 10 −4 permil in δ 13 C, for oxygen it is equivalent to deviations of 5×10 −4 and 2.5×10 −3 permil for δ 18 O and δ 17 O, respectively2 .This limit is a trade-off between the precision and speed of the tagged system, and is at least orders of magnitude below experimental precision.
In the next sections, we describe the results of the simulations.They all show a good run-time performance and further reveal that tagging offers a remarkable speedup compared to the doubling, due to the optimisations used.

Test case with a carbon isotope prototype mechanism
A simple, but representative test was performed on the prototype mechanism derived from the complete MECCA chemistry mechanism (similar to that applied by Jöckel et al., 2006 on a global scale) but excluding NMHC chemistry.The resulting mechanism is schematically depicted in Fig. 2 and has the following properties: -The mechanism contains intermediate and nonintermediate species (CH 4 as a source and CO 2 as reservoir) and two artificial source species (C X , C Y ).
-The species exchange their isotopic composition via different reactions (the rates are taken from the original MECCA reactions list), some species are recycled.
-A kinetic isotope effect for the CO+OH reaction is used equivalent to an enrichment of +6.5‰ in δ 13 C(CO) at 1 bar (Röckmann et al., 1998b).
The prototype mechanism was configured for the CAABA box-model.A set of six simulations, in which the species' isotopic composition was traced, was carried out to check the isotopic composition exchange between the species as follows: -All species were initialised with the same isotopic composition (setups 1, 2).
-The dominating source of carbon in the system was either methane (3, 4) or artificially added C X and C Y species (5, 6) with corresponding values of δ 13 C(CH 4 ) = -52‰ and δ 13 C(C X ,C Y ) = -32‰.
-The system was closed, except for setups 5 and 6, where methane was absent and C X and C Y were added.
-Simulations were performed with a parameterised seasonal cycle determining the photolysis rates and the OH concentration.
www.geosci-model-dev.net/3/337/2010/Geosci.Model Dev., 3, 337-364, 2010 • s rp q p q e min,i C p + 1 − q p q e maj C p (reaction with the specified particular source) b Notation as before.Bracketed species are those from the regular mechanism used as virtual catalysts.Here we feature the reaction rate constant k r instead of reaction rate v r .a For this bimolecular reaction only one pair of equations in which C n is used as a catalyst is shown, another pair for C e is to be written in analogous way; thus one reaction of two tagged species of M isotopologues requires 2•M reactions to be added.b The isotopic composition of the specific source C s is accounted via its isotopologues molecular fractions f s incorporated in the reaction rate coefficient.Each particular source requires an additional set of M equations.The educts are destroyed separately in an additional set of virtual equations.
The simulated isotopic composition of CO (Fig. 3) in all cases adequately reflects the dominating source signature as well as the kinetic isotope effect.Interestingly, in setup 6 the system exhibits a noticeable KIE feedback via the OH channel.More OH availability (summer) results in CO enriched in 13 C due to the oxidation reaction KIE, while less OH (winter) leads to the takeover of the transfer of lighter composition of C X and C Y through HCHO (i.e. the KIE "competes" with the precursor signature).
A noticeable magnitude of the isotopic composition variation can be explained by the low CO abundance (and hence the more pronounced oxidation KIE signal) caused by the much smaller (compared to methane, see setup 4) source of carbon from C X and C Y , whose contribution to the intermediates (methanol and formaldehyde) results also from oxidation by OH.Excluding the non-closed system setups (5 and 6), the total carbon content and isotopic ratio of the system were found to be conserved.

The carbon and oxygen isotopic composition of CO
Another simulation was performed with the carbon-12/13 and oxygen-16/17/18 isotope tagging of the comprehensive mechanism (see Jöckel et al., 2006, Appendix B) with the CAABA/MECCA box-model.This mechanism includes ozone-and non-methane hydrocarbon (NMHC) related tropospheric chemistry.The main purpose of this study is the reproduction of the oxygen isotopic composition in view of the different factors (see Fig. 4) that determine the massindependent fractionation (MIF) in CO (see also Röckmann et al., 1998aRöckmann et al., ,b, 2002)).The simulated mechanism enables us to study the oxygen isotopic transfer through the various Notation as before.The isotope exchange and particular composition source reactions are parameterised identical to the implicit doubling case.Similarly, the doubled species are used as the catalysts in virtual reactions.
a Shown example of bimolecular reaction of two doubled species presents only reactions of major and ith minor isotopologues and a pair of i-th and j -th minor isotopologues.
Considering all possible combinations of reactants' isotopologues, the total number of replicated reactions for M isotopologues amounts to 2 M .Fig. 3. Carbon monoxide isotopic composition evolution (solid lines) for different simulation setups (1-6) with the prototype mechanism (Fig. 2).Dashed and dash-dotted lines denote the formaldehyde composition for setups 5,6 and 3,4, respectively.
intermediates together with the direct contributing factors, such as alkene ozonolysis and MIF during the CO oxidation by OH.
The tagging setup was simulated within a tropospheric box with prescribed monthly averaged trace gas surface emissions of the corresponding isotopic composition (see Table 2).
A particular carbon and oxygen transfer mechanism was specified (referring to Sander et al., 2005;Weston, 2001, listed in the supplement, and oxygen isotope exchange reactions were added, following Lyons, 2001;Zahn et al., 2006 and references therein).In addition, the anomalous (MIF) kinetic isotope effect was added for CO+OH corresponding to depletions of about 9.4‰ in δ 18 O(CO) and 0.2‰ in δ 17 O(CO), which results in an enrichment of +4.7‰ in 17 O(CO).The anomalous KIE was also introduced for the ozone formation reaction (O 2 +O) yielding enrichments of δ 18 O(O 3 ) of 90‰ and δ 18 O(O 3 ) of 78‰ correspondingly, resulting in the value of 17 O(O 3 ) of about +30‰.The non-random distribution of oxygen isotopes within the ozone molecule was not considered.According to estimates (Brenninkmeijer and Röckmann, 1997;Weston, 2001), the kinetic isotope effect leading to the  Röckmann et al. (1998a), Röckmann et al. (1998b, 2002) for references.c Listed are the enrichments (depletions) in CO due to the KIE.
production of mass-dependently fractionated (MDF) carbon monoxide of δ 18 O(CO) ∼0‰ from the methane source chain was added.The ozone KIE was set to be either anomalous (reference setup) or mass-dependent (MDF-O 3 setup) in the different simulations.
The isotopic carbon tagging setup included the KIE for the CO+OH reaction of a magnitude equivalent to an enrichment of +6.5‰ in δ 13 C(CO).The isotopic signature of methane was constrained to the value of δ 13 C(CH 4 ) = -47.2‰.In combination with the introduced KIEs for methane oxidation by OH and O( 1 D) (Saueressig et al., 2001) this source conclusively contributes to CO with a very negative signature of -52‰ compared to the other primary sources of about -27‰.In the spring and summertime available OH triggers the KIE-escorted sink and light methane-based refilling of the CO burden, thus the competition of these two processes mainly determines the variability of the 13 C/ 12 C ratio in CO.To assess the sensitivity of the mechanism to the OH level, the box was simulated in two distinctive regimes of OH seasonality emerging in the box positioned at the NH latitudes of 70 • N (high-latitude setup, HL) and 30 • N (low-latitude setup, LL).We further introduced a subsidiary tagging configuration of the mechanism to estimate the fraction of the carbon in CO that has originated from CH 4 .For this, all carbonaceous species were tagged into two classes, namely "M" (methane) and "N" (nonmethane) and all molecules of CH 4 were assigned to the class "M", whilst the other species' molecules were initialised Generally the box-model reproduces CO (Fig. 5, upper left) and its isotopic composition well.
The average δ 13 C(CO) signature (Fig. 5, lower left) for the low-latitude (LL) box is lower due to greater input from methane (from 19% in the late winter to 37% in summer) compared to the high-latitude (HL) box (10% to 31%, correspondingly).Since the latter has much less OH, the maxima and minima in δ 13 C(CO) are less pronounced and shifted in time towards the summer.During the fall-winter season the HL box reflects the signature of the primary sources in a characteristic "shoulder" at -27‰.The observed composition reflects the effect of mixing of CO from lower latitudes which is longer exposed to the OH during the transport.The low-latitude CO is less influenced by the transport effect and the LL box captures this.
The oxygen composition of CO is more influenced by the source signatures change during the year.The simulated δ 18 O(CO) signature in both boxes (Fig. 5, upper right) is overestimated by 5-7‰.We attribute this to the underestimated signatures of the main sources (Table 2).One reason may be that the δ 18 O value of CO from CH 4 or NMHC oxidation is not 0‰ as we assumed, but lower.The variability of δ 18 O(CO) in the LL box is a factor of two lower than that of the HL box due to the increased input from these in-situ sources in the spring-summer season.In the fall-winter season CO is refilled with heavier surface sources.Analogously to the 13 C case, due to the absence of transport, δ 18 O(CO) in the HL box starts to decrease with a delay, as OH appears later.(Gros et al., 2001;Röckmann et al., 2002) for CO for high (80 • N, filled symbols) and low latitudes ( 30• N, open symbols).The lower right panel presents the mass-independent fractionation signal ( 17 O) in CO and additional results of the mass-dependently fractionated O 3 setup for high latitude (dash-dotted line) and low latitude (dotted line) simulations.
The essence of the oxygen isotope simulation result is shown in the triple-isotope plot (Fig. 6, reference setup for HL).The enrichments in 17 O and 18 O are found to be massindependent (the magnitude is proportional to the vertical deviation from the MDF-line) reflecting a contribution of the large MIF-signal of ozone.The chemical "mixing" via the intermediates and the isotope exchange reactions cause the HO x and NO x groups to "move" along mixing lines of "dilution" with the end members being ozone and the two main oxygen reservoirs of air O 2 and H 2 O.
The CO anomaly is reproduced and can be attributed dominantly to the reaction with OH.In CO+OH the isotope effect for 18 O is inverse and the concomitant fractionation for 17 O is uncommonly not of about a half of that of 18 O, but is slightly positive.Thus, 17 O accumulates larger than expected for MDF composition during the OH based oxidation of CO.The 17 O(CO) maximum falls on the spring (LL) and mid-summer (HL) period and is weakened during the fall-winter seasons (Fig. 5, lower right) because of surface and in-situ sources dominating.The seasonal cycle in the HL box is delayed and the maxima are weak due to the low OH regime.In the fall-autumn seasons, 17 O(CO) in the HL box is also overestimated which can be explained by the imbalance between in situ MIF sources (ozonolysis of unsaturated hydrocarbons) and MDF in-situ and surface sources.Comparing the reference case with the MDF-O 3 setup, it becomes clear that the large ozone MIF signal is found to contribute to the CO anomaly with additional 1‰ for the presented regime.This in contrast to the nitrogenated compounds whose 17 O values are directly related to that of O 3 .For the low-latitude box there are no 17 O(CO) observations available.
Obviously, important processes such as transport and mixing cannot be represented by a box-model, it is especially noted with the high latitude observations, which are difficult to reproduce because of the "arctic haze" conditions ruling in the winter.Our further studies will focus on a better estimation of emission strengths and signatures, as well as the realisation of consistent AC-GCM simulations.Nevertheless, this comprehensive isotopic chemistry mechanism is able to capture the main features of seasonal changes in species compositions as well as the "chemical mixing" feature.Brenninkmeijer et al. (2003) and Thiemens (2006).e Includes sulphur chemistry.

Performance analysis: tagging versus doubling
To estimate the benefit of the optimised tagging technique, we compare the performance of both methods (tagging and doubling) in a series of box-model simulations with kinetic chemistry mechanisms of various complexity.Table 3 lists some MECCA mechanisms, which we "tagged" and "doubled" to simulate a stable carbon isotope chemistry configuration (thus accounting for two additional isotopologue classes).Each selected mechanism was repeatedly simulated in box-model setups for regular (MECCA), doubled (MECCA-DBL) and tagged (MECCA-TAG) cases.
We compare the wall-clock time spent for a one month simulation period.For the integration of the regular and doubled mechanisms we used the 3rd order Rosenbrock solver (ROS3) with adaptive time stepping; it is provided by the KPP package in MECCA and is very powerful in handling stiff problems efficiently.The integration of the tagging mechanism (Eq.28) was performed with a simple solver on a modified Runge-Kutta algorithm with adaptive time-step (RKA).The latter is not very well suited for systems of high stiffness, yet in other cases usually shows a decent performance due to its robustness (Press et al., 1992).Although MECCA-TAG incorporates other more complex solvers, this analysis is to show that even the RKA implemented in the tagged system offers a remarkable speed gain in the competition with the highly efficient ROS3.
To make the comparison more representative, two distinct simulation regimes, namely "free" and "forced", were chosen.In the "free" regime the initial conditions (species abundances and classes ratios) are set once before the simulation starts while the system remains closed (i.e.no introduction or removal of the molecules is allowed) during the simulation.As a consequence, the system proceeds towards an equilibrium state in both, abundance and isotopic signature.In the "forced" regime some species are permanently emitted into the box thus keeping the system in a perturbed state.These two regimes are the extremes of conditions of the chemistry grid-boxes processed in an AC-GCM, e.g., ranging from grid-boxes positioned at the surface layer (with emissions as part of the boundary conditions) and at the levels where the perturbation is smaller due to smoother species gradients and weaker mixing terms.For the presented simulations, the "forced" regime integration required more computational time because of the larger stiffness of the ODE system to be solved in contrast to the "free" regime.Comparing the overall results, we get a rough estimate of the chemistry mechanism performance to be expected in 3-D applications.
In Fig. 7 (left) the simulation times for the different mechanisms and regimes are shown; the average complexity of the mechanisms increases along the abscissa.Figure 7 (right) also depicts the estimated average additional time (expressed as a fraction of the regular mechanism integration time) required to process the doubled and tagged systems.In general, the tagged system outperforms doubling, as expected.Speed gains of a factor 3 and larger are achieved when simulating the M6, M7 and AR mechanisms; for the former two the time saved is to a certain extent proportional to the fractions of the tagged species and reactions.The least speed gain was achieved with M9, which includes sulphur chemistry accounting additionally for carbonaceous species involved in the reactive chemistry, which is inefficiently simulated with RKA in the "forced" regime.The same difficulty appears if higher hydrocarbons were included.For the EVAL mechanism the required CPU-time increased by 30% for tagging and by 70% (80% for M7) for doubling.The mechanism including all gas-phase reactions (AR) requires comparable CPU-times as EVAL, with slightly better results for the tagging, potentially due to a more efficient handling of the halogen chemistry by RKA.
Evidently, each mechanism has its own individual characteristics that may be considered individually for a better optimisation strategy.We note that for the molecular tagging case the integration time is expected to be ∼30% less compared to the isotopic tagging; the latter needs one kinetic exchange derivative calculation more, operating with atomic fractions of the species.The overall performance of the tagging technique, however, is very promising, especially in perspective of the tagging of the complex chemistry mechanisms in 3-D model environments.

MECCA-TAG/DBL implementation in MESSy
Our realisation of the MECCA-TAG and MECCA-DBL subsubmodels allows the simultaneous simulation of an arbitrary number of tagging configurations (i.e.sets of tagged species and reactions).For the tagging, this concept is of a great use when more than one tagging configuration is needed to be simulated with the same regular chemistry mechanism.A particular example is a single AC-GCM simulation with www.geosci-model-dev.net/3/337/2010/Geosci.Model Dev., 3, 337-364, 2010 various separate stable isotope configurations that tag the same incorporated regular chemistry mechanism.We stress again that none of each configuration interferes with the others or with the regular mechanism, thus remaining purely diagnostic.
The setup of a particular configuration is controlled by the configuration file (the details are presented in the electronic supplement), which contains all necessary information to perform the full (isotopic) tagging of a given mechanism of desired complexity.The configuration files are used for both, the tagging and the doubling sub-submodels that utilise the same pre-processing routines.The latter parse the selected MECCA chemistry mechanism (i.e. the species and reactions files) and tagging configuration files and automatically generate the sub-submodels code plugged into the destination simulation code.Parsing routines significantly facilitate the model setup, additionally checking for possible overlooked inconsistencies, like missing species or mass imbalance in any of the tagged reactions.The latter are also diagnosed for unaccounted sources of the isotopologues in the mechanism (with respect to the approximations described in Sect.4.3).The parsing routine codes are written in Pascal (it offers greater flexibility in handling string and text information compared to Fortran) and built with Free Pascal Compiler (http://www.freepascal.org/).The execution of the code is controlled (switched on/off) via the MECCA conventional Fortran95 namelist.
The implementation of the MECCA-TAG/DBL is in conformity with the MESSy standard (Jöckel et al., 2005).This implies the relocation of the code parts to the dedicated layers of the MESSy interface structure according to their functionality.The interface part of MECCA-TAG is responsible for passing the kinetics information from MECCA to the MECCA-TAG kinetics driver.The latter prepares the pair rates for all tagged species and invokes the integration and processing of each particular configuration.Each of the tagging configurations resides (see Fig. 8) in the submodel core layer (SMCL) thus being independent from the basemodel.The driver is called from within the submodel interface layer (SMIL) after the call to the regular MECCA integration.The interface layer routines (which are basemodel-dependent) at the BMIL layers contain subroutines processing I/O and other maintenance and provide the interface to other MESSy submodels and base models.Being in fact an extension of MECCA, MECCA-DBL utilises its infrastructure routines.MECCA-TAG and MECCA-DBL are being included in the latest releases of the MESSybased box-model CAABA/MECCA (Sander et al., 2010) thus providing diagnostic and isotope modelling extensions to this model.

Summary
We introduce the kinetic chemistry mechanism tagging technique intended for various applications in our modelling system.Assumptions and approximations are derived and discussed; optimisations of the technique as required for its application in the comprehensive kinetic chemistry module MECCA are presented.The optimised variant is implemented in the sub-submodel MECCA-TAG and intended for the application in complex simulations (specifically 3-D) to reduce the computational demands but with reasonable precision.
As a particular application, our approach for the modelling of the isotope chemistry-enabled mechanisms is described.We review the kinetic chemistry isotopic specifics and peculiarities that are required to be accounted for.In particular, a specific "isotope" transfer approach has to be applied, which differs from the original "molecular" exchange perspective.The evaluation simulations performed with the box-model, focusing on stable isotope tagging of the comprehensive mechanism, result in a reasonable agreement with observational data.
Additionally, as a reference, we introduce the subsubmodel MECCA-DBL which comprises the implementation of a non-optimised straightforward doubling of the existing MECCA chemical mechanism.
Taking over from MECCA its thorough capabilities and precision, MECCA-DBL can be used as a full isotope-chemistry enabled module intended for less resource-demanding applications (e.g.simulations with the box-model).Meanwhile, at the rather high precision/speed trade-off, MECCA-TAG shows a remarkably better performance.
MECCA-TAG/DBL allow flexible multi-configuration tagging of the comprehensive chemistry mechanisms additionally supporting isotope kinetic specifics.The submodels are conform to the Modular Earth Submodel System (MESSy), thus being automatically ready to be coupled to Earth System Models by means of the MESSy framework.abu,i χ n , respective shares of abundant and rare isotope atoms abu,i χ n = (q n −1) q n • min,i f n rare,i χ n (of the same atomic number Z) in min,i C n isotopologue rare,i χ n = 1 q n • min,i f n min,i J Copy of the rate matrix J with some changed elements accounting for KIE for i-th minor isotopologue in the isotopically tagged system [molec cm −3 s −1 ] i R The isotopic ratio of the i-th rare isotope (see also Appendix -Isotopomers: from "isotopic isomers", are the set of different isotope position configurations of a selected isotopologue.Hence, 15 N 14 N 16 O, 14 N 15 N 16 O are the isotopomers, so are CH 2 DCH=O and CH 3 CD=O.

(Conventional definitions)
-Isotopic ratio R is the atomic abundance ratio of any minor to a major abundant isotope: 13C R(CO) =

C
12 C in CO -The delta-notation is used to express the difference of a given sample isotopic ratio relative to the reference material sample, usually reported in permil: -The isotopic fractionation constant ε expresses the difference of one compartment isotopic ratio relative to another, often used to describe isotopic composition changes due to various processes (e.g.soil uptake fractionation, KIE, etc.)For example, the fractionation of 18 O isotopes in water due to the vapour pressure isotope effect is defined as: -The isotopic fractionation factor α determines how much faster/slower the reaction with the minor isotopologue goes relative to the reaction with the major isotopologue: where i k CH 4 +OH is the CH 4 +OH→ reaction rate corresponding to the i-th isotopologue.Note that this definition differs from the UIPAC-recommended definition of α as the ratio of reaction rate of the isotopically light (lower atomic mass) species to the rate of the heavy (higher atomic mass) ones (Coplen et al., 2002).
-The capital delta notation is used to express the composition deviation from general, largely observed correlation between minor isotopologues enrichments, called mass-independent isotopic fractionation (MIF): 17 O ≡ δ 17 O + 1 where β determines the slope of the line on which the data points of the mass-dependently fractionated composition fall in the three-isotope plot.The value β is close to 0.5 but varies individually for a particular species and fractionation process.Here, for universality, we use the value of β=0.528 derived for the meteoric waters (Barkan and Luz, 2005) for all species.With relation to the atmospheric reservoir (where the magnitude of isotope ratios deviation is not large) the relation above can be linearly approximated as: We refer the reader to the literature (Brenninkmeijer et al., 2003;Johnson et al., 2002;Kaye, 1987;Kohen and Limbach, 2006;Schauble, 2004) for an extensive review of the isotopic effects and their application in atmospheric studies.

Isotope transfer calculation
Below we elucidate the isotope transfer in the system (17).Consider the reaction of the minor isotopologue (line 2) in the Eq. ( 17): min,i C e k r −→ p∈P r s rp q p q e min,i C p + 1 − q p q e maj C p At the left hand side (LHS), an isotopologue possessing q e isotopes of the element of interest is reacting.We consider the case when q e is greater than unity; hence, the reacting isotopologue carries one rare and (q e −1) abundant isotopes.For example, for the rare 13 C-isotopologue of propene (C 3 H 6 ), q e =3 and the isotopologue composition is 13 C 12 C 2 H 6 .Now, consider a case when the right hand side (RHS) has a single product C p , with the number of the isotopes q p .The stoichiometry condition (16) requires that s rp =q e /q p .Three cases are possible then: 1. q e = q p , i.e. number of isotopes is equal in the educt and the product, and s rp =1; 2. q e > q p , i.e. the product species pool receives one rare and (q p −1) abundant isotopologues; 3. q e < q p , i.e. the educt is incorporated into the product.
For the case (1), the transfer of the rare and abundant isotopes is unambiguous.In the cases 2) and 3), the stoichiometric coefficient is not equal to unity to create the same number of isotopes redistributed in the product molecules.Noteworthy, the LHS has only one rare isotope, thus only one rare isotope can be created at the RHS.The additional coefficient q p /q e next to the product minor isotopologue in the RHS ensures that; thus not more than one rare isotope is transferred as a part of one minor isotopologue.As an artificial example for the case 2), if a 13 C-propene isotopologue produces formaldehyde molecules (HCHO), then s HCHO =3 and only one of these HCHO carries 13 C isotope, the rest two acquire the abundant 12 C from propene: The case 3) is somewhat intricate.Imagine, for the example only, that now HCHO produces propene.The stoichiometric coefficient in this case is s C 3 H 6 =1/3, i.e. one needs to assemble three carbons from the formaldehyde molecules to create C 3 H 6 .However, we allow the minor isotopologue of propane to carry not more than one 13 C isotope, thus the rest isotopes have to be 12 C.These can be taken only in form of propene major isotopologues, so that the overall isotopic ratio of the propene pool is preserved.Conclusively, for this example of case 3), it follows: Again, the coefficients require that only one minor isotopologue is being created, and it is completed with the abundant isotopes from the major isotopologues.A more realistic example of the case 3) is the bimolecular reaction of the production of the ozone isotopologues from minor isotopologues of O( 3 P) and O 2 , considering 16 O and 18 O isotopes: Here, the bimolecular reaction decomposition (Eq. 3) is used with the branching ratios (Eq.4), so that effective stoichiometric coefficients become s O( 3 P) =1/3 and s O 2 =2/3, respectively.This argumentation can be straightforwardly extended to the cases when the RHS has more than one product.Finally, we calculate the number of rare and abundant isotopes at the LHS and RHS in the general formulation (17).First, we count the number of RHS rare isotopes only, i.e. n rare = p∈P r s rp q p q e = 1 q e p∈P r s rp q p = 1 assuming one rare isotope in each minor isotopologue and recalling the stoichiometry condition (16).Similarly, the number of the abundant atoms at the RHS is n abu = p∈P r s rp q p q e q p − 1 + 1 − q p q e q p = (q e − 1) q e p∈P r s rp q p = (q e − 1) taking into account that the major and minor isotopologue of the products have q p and (q p − 1) abundant isotopes, respectively.Thus, the total number of transferred isotopes is q e , of which one is rare.

Fig. 2 .
Fig. 2. Schematic representation of carbon interchanges in the prototype mechanism.Intermediates are shown in blue; methane, carbon dioxide and two artificial species (C X , C Y ) are non-transient in-situ C source species.Pathway captions refer to the internal MECCA reaction number and second reactant.

Fig. 4 .
Fig. 4. Main sources and sinks determining the 17 O signal in carbon monoxide.

Fig. 5 .
Fig.5.Simulated CO mixing ratio and isotopic composition in the box.Solid and dashed lines refer to the HL and LL box setups, respectively (see text).Symbols denote the observations(Gros et al., 2001;Röckmann et al., 2002) for CO for high (80 • N, filled symbols) and low latitudes (30• N, open symbols).The lower right panel presents the mass-independent fractionation signal ( 17 O) in CO and additional results of the mass-dependently fractionated O 3 setup for high latitude (dash-dotted line) and low latitude (dotted line) simulations.

Fig. 6 .
Fig.6.Triple-isotope plot of the main oxygen-bearing gas-phase species composition simulated in the reference setup.The shaded areas indicate typical compositions observed for certain species reported byBrenninkmeijer et al. (2003) andThiemens (2006).

Fig. 7 .
Fig. 7.Left: total time spent to integrate the regular, doubled and tagged systems.Solid and dashed lines denote "free" and "forced" regimes, respectively (see text for details).Right: the average additional time required for the doubled and tagged systems as a fraction of the regular mechanism integration time.The vertical lines on the top of the bars represent the variability due to the different regimes.

Fig. 8 .
Fig.8.Schematic of the MECCA-TAG and MECCA-DBL principal structure (SMCL level is detailed).The regular simulated mechanism (R) provides one-way information on the tagged reaction rates and species budgets to the MECCA-TAG driver.Further each configuration (T1, T2) uses approximated rates to advance the tagged species.Additional species and reactions of the doubling configurations (D1, D2) are added to the regular mechanism, thus MECCA-DBL is just an extension of MECCA.

S
. Gromov et al.: A tagging technique and isotope chemistry modelling m C n tagged mechanism n-th species m-th class species m C = m C 1 ,..., m C N S , C n = m C n m c n tagged mechanism n-th species m-th class species concentration c share of m-th class species in total n-th species molecular abundance i minor isotopologue index i=1...M−1 maj C n , isotopically tagged mechanism n-th species major maj C n is equivalent to 1 C n min,i C n and i-th minor isotopologues min,i C n is equivalent to 1+i C n , i=1...M−1 maj c n , isotopically tagged mechanism n-th species major if C n refers to CO, for oxygen isotopes min,i c n and i-th minor isotopologues concentration major isotopologue: maj c n = [C 16 O], minor isotopologues: min,1 c n = [C 17 O], min,2 c n = [C 18 O] maj f n , respective shares of the major and i-th minor maj f n = f n isotopologue molecules in n-th species total molecular abundance; maj f n is equivalent to 1 f n min,i f n is equivalent to 1+i f n , i=1...M−1 min,i f n = min,i c n maj c n + KIE fractionation factor for min,i C e isotopologue educt in r-th reaction 13 CH 4 α CH 4 +OH = k13 CH 4 +OH k12 CH 4 +OH a,i β r isotope exchange reaction fractionation factor for min,i C a isotopologue λ a→b r probability of rare isotope transfer from C a to C b species in the isotope exchange i-th rare isotope atom transfer from C a to C b species in the isotope exchange reaction r [molec cm −3 s −1 ]i γ a net rate of gain (loss) of i-th rare isotope atoms i γ a = − species due to all isotope exchange reactions it participates in [molec cm −3 s −1 ] a example reaction: CH 3 OOH + OH −→ 0.7• CH 3 O 2 + 0.3 • H 2 CO + 0.3 • OH + H 2 O b example reaction: O 2 + O( 3 P) −→ O 3 Appendix BUsed terms and definitions for the isotopesNot to be confused with different terms and definitions used to express the isotopic composition and effects, the following definitions are used throughout this document.Those are: (IUPAC definitions, McNaught and Wilkinson, 1997, http://goldbook.iupac.org)-Isotopologue: a molecular entity that differs only in isotopic composition (number of isotopic substitutions), e.g.CH 4 , CH 3 D, CH 2 D 2 are isotopologues.
Gromov et al.:A tagging technique and isotope chemistry modelling individual contribution to the budget of the others.Another prominent application is the isotopic tagging, when the presence of the rare isotope in the molecule constitutes indeed a natural intrinsic "tag".This isotope tagging is discussed below.

Table 1a .
Example for implicit isotopic doubling of the regular chemical mechanism.

Table 1b .
Example for explicit isotopic doubling of the regular chemical mechanism.

Table 2 .
Emission and in-situ sources and their signatures used in the box-model simulations a,b .
a c Includes stratospheric reactions set.d Includes halogen chemistry (excluded in other mechanisms).