The Chemical Mechanism Integrator Cminor v1.0: a stand-alone Fortran environment for the particle-based simulation of chemical multiphase mechanisms

We present version 1.0 of the Chemical Mechanism Integrator (Cminor), a fully modularized modern Fortran software package for the computational simulation of skeletal and detailed chemical kinetic systems derived from atmospheric and combustion chemistry. Cminor aims for the efficient simulation of complex chemical mechanisms by using various mathematical techniques. These are tailored to systems of ordinary differential equations (ODEs), having the specific structure arising from chemical reaction systems. Additionally, a high-speed mechanism parser allows the user to interchange reactions or their parameters in an ASCII format text file and immediately start a new simulation without recompiling, enabling fast and numerous simulations. Cminor's solver technique is based on Rosenbrock methods. Different measures of local errors and an analytical Jacobian matrix approach are implemented, where efficiency is obtained by exploiting the sparsity structure of the Jacobian.

Cminor can be run in one of three configurations:

• A box-model framework for either pure gas-phase mechanisms or a multi-modal aerosol distribution dissolved in mono-dispersed cloud droplets.

• A rising adiabatic parcel, in which the activation of multi-modal aerosols is represented by solving the droplet condensation equation.

• A constant volume environment, where thermodynamic properties are evaluated by polynomial functions of temperature according to the standards of the Chemkin thermodynamic data base.

The software package is evaluated by applying seven different chemical mechanisms. Three of them are from the field of air-quality modeling and three are from the area of combustion kinetics, ranging from 7 species and 10 reactions to 10 196 species and 23 098 reactions. The last mechanism describes sulfur accumulation in clouds, which is tested along with a rising parcel and condensating cloud droplets.

1 Introduction

The understanding of physical processes that determine atmospheric states of climatological, ecological, social, political or agricultural relevance has undergone significant improvement in theoretical and modeling aspects. While chemical processes are omnipresent in any state of the atmosphere, they have drawn relatively minor attention when physical aspects of the atmosphere are considered. Nevertheless, the problem of considering the atmosphere as a physicochemical system, rather than an entirely physical system, cannot be avoided. Therefore, models of atmospheric processes usually work with assumptions about the underlying chemical state. Major aspects, e.g., sensitivity to the chemical state, chemical processing, or inter-dependencies, are mostly ignored. In contrast to this, the community of atmospheric chemistry has shown a number of important possible influences of chemical processes. These concern some of the most relevant aspects of atmospheric physics, like radiative properties, aerosol size spectra, cloud field dynamics and precipitation (see, e.g., Murphy et al., 1998; Kulmala et al., 2000; Kreidenweis et al., 2003; Christner et al., 2008; Zaveri et al., 2010; Kazil et al., 2011; Goel et al., 2020; Kupc et al., 2020; Zhang et al., 2024). Atmospheric chemical processes are described by a collection of reactions, constituting a specific chemical mechanism, possibly including multiple phases (gaseous, aqueous, solid) and phase changes, i.e., nucleation or dissolution. Timescales of chemical reactions might depend on temperature, pressure, insolation, or chemical composition, but generally vary largely from reaction to reaction.

By stating an ordinary differential equation (ODE) for the time evolution of every chemical species' concentration in a mechanism, large and stiff ODE systems arise. Therefore, the simulation of comprehensive chemical reaction mechanisms requires the development of efficient algorithms, particularly suited for the characteristics of ODE systems of chemical kinetics. The software package Cminor (Chemical Mechanism Integrator), written in modern Fortran, was developed as an environment for the simulation of atmospheric chemistry mechanisms of any size and type.

Several software packages have been developed for the analysis of kinetic systems, e.g., ChemKin (Kee et al., 1996), SPACCIM (Wolke et al., 2005), KPP (Sandu and Sander, 2006), CAABA/MECCA (Sander et al., 2011), and SpeedChem (Perini et al., 2012). As a matter of fact, the ongoing further development and creation of new chemical mechanisms shows that a tool providing an efficient solver and high flexibility is absolutely necessary. This goal is achieved by Cminor.

This paper is organized as follows. The next section introduces how a Cminor simulation is set up, including various input options. The third section assesses the equations solved by Cminor, followed by a section on the numerical methods used by Cminor. Afterwards, simulation results are shown. The paper is concluded by a last section summarizing results and perspectives of this work. All variables not explicitly defined in the text are listed in the Supplement, Sect. S11.

2 Model Input

The Cminor software package is a two-in-one simulation tool for chemical kinetics problems in atmospheric and combustion chemistry. Cminor reads a reaction system provided as an ASCII file, generates a set of stiff ODEs, and solves it directly. One of the advantages of Cminor is the flexible high-speed ASCII mechanism parser where the generation time of internally used structures grows linearly with the mechanism size. This allows the user to quickly add/remove reactions or manipulate reaction rate constants within the reaction mechanism file. The ability to quickly redefine these parameters is crucial for investigating newly developed chemical mechanisms, in contrast to KPP (Sandu and Sander, 2006), where a fixed Fortran code for every new chemical mechanism has to be generated, compiled, and executed.

2.1 General setup

For a given model scenario with initial values and prescribed emission and deposition rates (i.e., constant sources and sinks for gas phase species), the user is able to simulate the evolution of species concentrations and other variables in a box-model framework over a defined time interval. The Cminor box-model is depicted in Fig. 1.

The syntax rules for combustion mechanisms can be taken from the ChemKin documentation (Kee et al., 1996), which is widely used for combustion kinetics. The necessary thermodynamic data base, which includes the coefficients for fitting polynomials, follows the same ChemKin syntax rules (Kee et al., 1990). For atmospheric chemical mechanisms, a detailed description of the syntax rules is provided in the Supplement. In any case, the input data has to be provided in readable ASCII format. Assistance in the syntax rules is also given by the examples provided in the repository, as well as on request by the authors. A short summary shall be given in the following.

Running a simulation requires four files beside the actual model code:

• the *.sys-file,

• the *.dat-file,

• the *.ini-file, and

• the *.run-file.

The *.sys-file is solely to list all reactions of the mechanism to simulate. The description of a reaction consists of at least three lines. The first line specifies the reaction class (gas phase, aqueous phase, forward/backward reactions, phase transfer), the second line declares the reaction itself, and the third line contains an identifier for the reaction rate constant type, with parameter values. The third line determines the specific formula of the rate constant, i.e., the speed, of the reaction. There are various reaction rate constant types, like

• Arrhenius laws (temperature dependencies),

• photolytic processes (insolation dependencies),

• Troe reactions (pressure dependencies),

• pH-dependent formulas,

• humidity-dependent formulas,

and others. Each type comes in a multitude of slightly varying formulas, depending on parameters like activation energy of the reaction. All available rate constant formulas are listed in the Supplement, Sect. S4. The description of reactions is exemplified in Table 1. The optional fourth line, starting with “FACTOR:”, sets a passive species to be included as an educt for the reaction. Passive species are mostly primary components of air, like O₂ or N₂, which are not included as a separate equation. Instead, their concentration is treated as constant since they are abundant. The available FACTOR species are as well listed in the Supplement, Sect. S3.5.

Cminor supports use of the SMILES notation for species throughout all files (Weininger, 1988). A species can have any name following these rules:

• the name of the species must be in capital letters,

• square brackets enclosing the whole name indicate passive species (e.g., [O2]),

• an “a” in front of a species indicates presence in droplets, i.e., an aqueous species (e.g., aO2),

• one or multiple “p” or “m” after the species indicate ions and their charge (e.g., CO3mm means CO${}_{\mathrm{3}}^{\mathrm{2}-}$), and

• allowed special characters in species names are [, ], (, ), =, +, and *.

Table 1 Overview and examples of the structure describing a specific reaction in the *.sys-file. The repertoire of rate constants and the meaning of the rate parameters (A, B, E/R, …) can be found in the Supplement, Sect. S3.5.

The *.dat-file lists molar masses, accommodation coefficients, and diffusion coefficients for gaseous species, as well as molar masses and charges for aqueous species. Moreover, the names of organic peroxy radicals (RO₂ species), gaseous or aqueous, have to be listed here if they are needed as (non-constant) factors in the mechanism. If specified as a factor for a reaction, the respective rate constant will be multiplied with the sum of the current concentrations of all RO₂ species specified in this file.

The *.ini-file contains the species' initial values, default is zero, and constant emission and deposition rates for gas phase species. Also, the aerosol composition and distribution are specified here.

The *.run-file is a collection of Fortran NAMELISTs to specify the remaining simulation parameters, e.g., temperature, a trigger for adiabatic parcel modeling, output-file name, time-integration method, and more. All NAMELISTs and their variables are listed in the Supplement, Sect. S9, including their default values. Multiple examples are given in the repository Rug et al. (2025).

Figure 1 Graphical depiction of Cminor. Initially, the box contains a fixed number of cloud condensation nuclei that can activate to cloud droplets. Next to gas phase reactions, aqueous phase reactions take place within the cloud droplets. The letters A_(g) and B_(g) denote the concentrations in the gas phase and A_(a) and B_(a) the concentrations in the aqueous phase, respectively. The magnification shows the processes near and inside of the cloud droplets. The interchange between the gas and aqueous phase is specified according to the Schwartz approach (Schwartz, 1986).

2.2 Liquid water content

The liquid water content (LWC) describes the mass of condensed water per volume of air. The LWC for a mono-dispersed droplet population is given by LWC =ρ_wVN, where ρ_w=1000 kg m⁻³ is the mass density of liquid water, V the droplet volume, and N the number concentration of droplets. Within this simplified model, the time dependent LWC(t) is given by a piecewise linear function mimicking the condensation process, where t is the independent time variable, $b_{\mathrm{1}},\mathrm{\dots },b_{\mathrm{6}}$ are the bounds of the linear function, and LWC_min and LWC_max the respective minimum and maximum LWC. If the simulation time exceeds the bound b₆, the function will be continued periodically, such that $b_i^{\prime }=b_i+(b_{\mathrm{6}}-b_{\mathrm{1}})$ for $i=\mathrm{1},\mathrm{\dots },\mathrm{6}$, as indicated in Fig. 2. From LWC(t), the droplet wet radii are calculated and used for aqueous chemistry and phase transfer. In the case of an adiabatic parcel, the LWC is determined by the droplet condensation equation, which will be described in Sect. 3.3, below. For prescribed LWC(t), all droplets contain the same amount of liquid water while different amounts of dissolved aerosol are possible.

Figure 2 Pseudo liquid water content function LWC(t). All boundaries can be defined by the user.

2.3 Aerosol modes

The composition of aerosol is initialized at the beginning of the simulation. The syntax for the setup of the aerosol modes is described in the Supplement, Sect. S6.3. We assume that the aerosol is completely dissolved in the droplets. The aerosol is assumed to be log-normally distributed with an arbitrary number of modes. Each mode can have a specific composition, resembling, e.g., coexisting sea-salt and ammonium-salt particles, or any internally and/or externally mixed aerosol distribution. The aerosol is assigned to a user-specified number of droplet classes, each representing a number of equal hydrometeors with a specific amount of liquid water and dissolved aerosol mass. The aerosol mass is determined by discretizing the aerosol distribution using logarithmically equidistant bins, one for each droplet class. To constrain the aerosol masses to a reasonable range, the aerosol distribution of each mode is bounded by determining the 0.5 %-radius-quantile and the 99.5 %-mass-quantile. If all modes contain the same aerosol type, this procedure is applied for the sum of these modes, using the respective smallest and largest quantiles of the individual modes. If the modes differ in aerosol type, the procedure is applied to each mode individually. In addition to the initial aerosol mass, general initial values for any aqueous species can be prescribed.

If the adiabatic parcel option is used, aerosol activation is determined by Köhler theory (Sect. 3.3). Alternatively, when the change in LWC is prescribed (Sect. 2.2), only aerosols with dry radii greater than a user-specified threshold are assumed to be activated, and the LWC is distributed equally among the activated particles.

3 Model Equations

All equations solved by Cminor are discussed here. Any variable, meteorological parameter or value not defined here is listed in the Supplement, Sect. S11.

3.1 Description of chemical kinetic systems

Consider a system of n_R reactions involving n_S chemical species that can be represented in the general form

$$\begin{array}{}\text{(1)}& {\displaystyle }\sum _{j=\mathrm{1}}^{n_{\mathrm{S}}}{\mathit{\nu }}_{ij}^{\mathrm{e}}{S}_j\stackrel{k_i}{\to }\sum _{j=\mathrm{1}}^{n_{\mathrm{S}}}{\mathit{\nu }}_{ij}^{\mathrm{p}}{S}_j,i=\mathrm{1},\mathrm{\dots },n_{\mathrm{R}},\end{array}$$

where the stoichiometric coefficients of educts ${\mathit{\nu }}_{ij}^{\mathrm{e}}$ and products ${\mathit{\nu }}_{ij}^{\mathrm{p}}$ are stored as sparse matrices. S_j denotes the name of the jth chemical species in the mechanism and k_i is the rate constant. Table 2 summarizes the submodules that provide the rate constant k_i for the case of atmospheric chemical mechanisms, a detailed description of all implemented types of rate constants can be seen in the Supplement, Sect. S4.

Table 2 Tabular overview of possible combinations of reaction types used in simulations of atmospheric chemistry.

All functions are implemented in a way that minimizes numerical operations, such as $(+,-,*,/,**)$, and calls to intrinsic functions, e.g., exp (⋅), $\log (\cdot ),$ sin (⋅) and cos (⋅), while performing these efficiently. Next to hard coded functions, the user is also able to declare own functions that do not match any of the functions provided, see Supplement, Sect. S4.1, to introduce custom rate constants.

The general reaction kinetics formulation Eq. (1) leads to a set of n_S ODEs expressing mass conservation in the system as

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\text{d}{c}_j}{\text{d}t}}=\sum _{i=\mathrm{1}}^{n_{\mathrm{R}}}{\mathit{\nu }}_{ij}{r}_i(\mathit{c},T,\mathit{\varphi })+c_j^{\text{emis}},j=\mathrm{1},\mathrm{\dots },n_{\mathrm{S}}\end{array}$$

where d$c_j/$dt is the change in concentration of species j, ${\mathit{\nu }}_{ij}=({\mathit{\nu }}_{ij}^{\mathrm{p}}-{\mathit{\nu }}_{ij}^{\mathrm{e}})$ is the difference between stoichiometric coefficients of products and educts of species j in reaction i, and $r_i(\mathit{c},T,\mathit{\varphi })$ is the reaction rate, depending on the concentrations of all reactants in reaction i, expressed by c, temperature T, and possibly other parameters ϕ. In addition, one can specify emission and deposition rates for gas phase species, which are summed up in the vector c^emis. While the emission rates are specified as actual rate in the *.ini-file, deposition is determined by a given rate constant from the *.ini-file. The sink/source term $c_j^{\text{emis}}$ is calculated as

$${\displaystyle }{c}_j^{\text{emis}}=c_j^{\text{in}}-c_j^{\text{out}}{c}_j,j=\mathrm{1},\mathrm{\dots },n_S,$$

where $c_j^{\text{in}}$ is the given emission rate and $c_j^{\text{out}}$ the given deposition rate constant. The reaction rate is calculated by

$$\begin{array}{}\text{(3)}& {\displaystyle }{r}_i(\mathit{c},T,\mathit{\varphi })=k_i(T,\mathit{\varphi })\cdot \prod _{j=\mathrm{1}}^{n_{\mathrm{S}}}{c}_j^{{\mathit{\nu }}_{ij}^{\mathrm{e}}},i=\mathrm{1},\mathrm{\dots },n_{\mathrm{R}}.\end{array}$$

The rate constants k_i may depend on temperature and other parameters like insolation (introducing a time dependency), or third body collisions (concentration dependency). The temperature plays a major role in combustion, and a respective equation is included in the calculations, as described later. While rate constants with non-temperature dependencies are implemented and can be used, they are not treated as non-autonomous parts of the system in the numerical procedure, beside being evaluated whenever needed. This is a common approximation and assumes not too rapid changes (e.g., a constant solar zenith angle, which is a good approximation for typical time steps below a few minutes, or third body collisions, that depend on species that are basically passive).

The right-hand side of Eq. (3) represents the law of mass action. In general, several other concentrations influence the production rate of each species, therefore the ODE system is highly coupled.

Because Cminor is capable of simulating gas phase combustion mechanisms, the initial value problem needs closure through energy conservation. A perfectly adiabatic constant-volume reactor is considered, as this condition is the one usually occurring when incorporating the solution of detailed chemistry into multidimensional computational fluid dynamics (CFD) codes, where the species and internal energy source terms for each cell are computed as part of a sub-cycling strategy for the whole code (Perini et al., 2012). In particular, the internal energy source term gives a change in the average reactor temperature

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\text{d}T}{\text{d}t}}=-{\displaystyle \frac{\mathrm{1}}{{\overline{c}}_{\mathrm{v}}(\mathit{c},T)\mathit{\rho }}}\sum _{j=\mathrm{1}}^{n_{\mathrm{S}}}\sum _{i=\mathrm{1}}^{n_{\mathrm{R}}}{U}_j\left(T\right)\left[{\mathit{\nu }}_{ij}{r}_i\right(\mathit{c},T,\mathit{\varphi })+c_j^{\text{emis}}],\end{array}$$

containing the temperature dependent polynomial fitting function U_j(T) describing the molar internal energy values of species j (Kee et al., 1996; Perini et al., 2012), the mass-averaged specific heat at constant volume ${\overline{c}}_{\mathrm{v}}(\mathit{c},T)={\sum }_{j=\mathrm{1}}^{n_{\mathrm{S}}}(c_j\partial U_j(T)/\partial T)$, and the average system density ρ.

3.2 Mass transfer between phases

The gaseous-aqueous mass transfer in multi-phase systems has to be defined by a pseudo first-order equilibrium reaction A_(g)⇌A_(a), where subscript g denotes the concentration of species A in the gas phase and subscript a the concentration of A in the aqueous phase, respectively. The forward and backward rate coefficients are calculated using the accommodation coefficient α, the gas phase diffusion constant D_g, and Henry's law constant. The correct syntax for species names is denoted in the Supplement, Sect. S3.3.

With Henry's law coefficient $A:=K_{\mathrm{H}}^{\mathrm{⊝}}$ in [$\mathrm{mol}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{atm}}^{-\mathrm{1}}$] and $B:=-{△}_{\text{soln}}H/R$ in [K], where ${△}_{\text{soln}}H=$ molar enthalpy of dissolution in [J mol⁻¹] and R is the gas constant. Then, Henry's law constant k_H is computed using the temperature-dependent equation

$$\begin{array}{}\text{(5)}& {\displaystyle }{k}_{\mathrm{H}}\left(T\right)=A\cdot \exp \left[B\cdot \left(T^{-\mathrm{1}}-T_{\text{ref}}^{-\mathrm{1}}\right)\right],\end{array}$$

where A and B are the parameter given in the third line of the reaction structure (Table 1), T is the actual temperature, and T_ref=298.15 K is a reference temperature. The droplet wet radius r is then used to calculate the mass transfer coefficient according to Schwartz (1986),

$$\begin{array}{}\text{(6)}& {\displaystyle }{k}_{\text{mt}}={\left(k_{\text{diff}}\cdot r^{\mathrm{2}}+k_{\text{acc}}\cdot r\right)}^{-\mathrm{1}},\end{array}$$

with $k_{\text{diff}}=\mathrm{1}/\left(\mathrm{3}{D}_{\mathrm{g}}\right)$, where D_g is the gas phase diffusion coefficient, and $k_{\text{acc}}=\mathrm{4}/\left(\mathrm{3}\widehat{\mathit{\alpha }}v\right)$, where $\widehat{\mathit{\alpha }}$ is the mass accommodation coefficient, $v=\sqrt{\mathrm{8}RT/\left(\mathit{\pi }{m}_{\text{mol}}\right)}$ is the mean molecular speed, and m_mol the molar mass of the gas. The reaction rates for gas uptake by droplets is

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\text{d}{A}_{\left(\mathrm{a}\right)}}{\text{d}t}}{|}_{\text{gas uptake}}=-{\displaystyle \frac{\text{d}{A}_{\left(\mathrm{g}\right)}}{\text{d}t}}{|}_{\text{gas uptake}}=k_{\mathrm{mt}}{\displaystyle \frac{{\text{LWC}}_i}{{\mathit{\rho }}_{\mathrm{w}}}}\left[A_{\left(\mathrm{g}\right)}\right],\end{array}$$

and

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{\displaystyle \frac{\text{d}{A}_{\left(\mathrm{g}\right)}}{\text{d}t}}{|}_{\text{gas escape}}& =-{\displaystyle \frac{\text{d}{A}_{\left(\mathrm{a}\right)}}{\text{d}t}}{|}_{\text{gas escape}}\\ & =k_{\mathrm{mt}}\left[A_{\left(\mathrm{a}\right)}\right]/\left(k_{\mathrm{H}}\right(T\left)R_{\text{atm}}T\right)\end{array}\end{array}$$

for the escape of gas from droplets, with R_atm=0.0820574 $\mathrm{L}\mathrm{atm}{\mathrm{mol}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$. The division by ρ_w converts the LWC from [kg m⁻³] to [m³ m⁻³], ensuring unit independent rates. The underscore i indicates that if there are multiple droplet classes, only the water volume of a single droplet class is determining the uptake of gases by this class.

3.3 Adiabatic parcel setup and condensational growth

In terms of the liquid water mixing ratio of droplet class i, q_l,i, the equation for condensational growth of water mass reads

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{\text{d}{q}_{\mathrm{l},i}}{\text{d}t}}=n_i\mathrm{4}\mathit{\pi }{r}_{\mathrm{d},i}^{\mathrm{2}}{\displaystyle \frac{S-S_{\text{eq},i}}{F_k{S}_{\text{eq},i}(r_{\mathrm{d},i}+r_{\mathit{\beta }})+F_d(r_{\mathrm{d},i}+r_{\mathit{\alpha }})}},\end{array}$$

where r_d,i is the radius of the droplets of this class and n_i the corresponding number mixing ratio of droplets.

$${\displaystyle }S={\displaystyle \frac{q_{\mathrm{v}}}{q_{\mathrm{S}}\left(T\right)}}$$

is the ambient relative humidity with q_v and q_S(p,T) being current and saturation vapor mixing ratio, respectively.

$$\begin{array}{}\text{(10)}& {\displaystyle }{S}_{\text{eq},i}={\displaystyle \frac{n_{\mathrm{w},i}}{n_{\mathrm{w},i}+n_{\mathrm{s},i}}}\exp \left({\displaystyle \frac{\mathrm{2}\mathit{\sigma }}{r_{\mathrm{d},i}{\mathit{\rho }}_{\mathrm{w}}{R}_{\mathrm{v}}T}}\right)\end{array}$$

is the value of the Köhler curve for the particles' current radius, where n_w,i and n_s,i are the concentrations of liquid water and dissolved species in moles, respectively, σ is the surface tension between liquid water and air, and R_v is the gas constant of water vapor.

$$\begin{array}{}\text{(11)}& {\displaystyle }{F}_k={\displaystyle \frac{L_{\mathrm{v}}}{k_{T}T}}\left({\displaystyle \frac{L_{\mathrm{v}}}{R_{\mathrm{v}}T}}-\mathrm{1}\right)\text{and}{F}_d={\displaystyle \frac{R_{\mathrm{v}}T}{D_{\mathrm{v}}{e}_s\left(T\right)}}\end{array}$$

are variables accounting for latent heat release and diffusion of vapor (Lamb and Verlinde, 2011). Lastly,

$$\begin{array}{}\text{(12)}& {\displaystyle }{r}_{\mathit{\beta }}={\displaystyle \frac{k_T}{{\mathit{\beta }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}p}}{\displaystyle \frac{(\mathrm{2}\mathit{\pi }{R}_{\mathrm{a}}T{)}^{\frac{\mathrm{1}}{\mathrm{2}}}}{c_{\mathrm{va}}+R_{\mathrm{a}}/\mathrm{2}}}\text{and}{r}_{\mathit{\alpha }}={\left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{R_{\mathrm{v}}T}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{D_{\mathrm{v}}}{{\mathit{\alpha }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}}}\end{array}$$

are physical relaxation parameters (Fukuta and Walter, 1970). To follow current common nomenclature, we changed the names from the original reference, because they refer to what is now called “accommodation coefficient“ as “condensation coefficient“. We kept the letter α referring to what is now usually called the accommodation coefficient, as for the gases in Sect. 3.2. By recommendations in Pruppacher and Klett (1978) and Fukuta and Walter (1970), we set ${\mathit{\beta }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}=\mathrm{1.0}$ and ${\mathit{\alpha }}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}=\mathrm{0.0415}$. These values can be changed to the user's preferences without recompilation, as described in the Supplement, Sect. S9.4.

We have to include a few other equations to account for an adiabatically rising parcel. The water vapor mixing ratio q_v is only affected by condensation and evaporation, i.e., the change in liquid water mixing ratio q_l, such that

$${\displaystyle \frac{\text{d}{q}_{\mathrm{v}}}{\text{d}t}}=-{\displaystyle \frac{\text{d}{q}_{\mathrm{l}}}{\text{d}t}}$$

with

$${\displaystyle }{q}_{\mathrm{l}}=\sum _{i=\mathrm{1}}^{n_{\mathrm{D}}}{q}_{\mathrm{l},i},$$

where n_D is the number of droplet classes. A user-specified updraft velocity

$${\displaystyle \frac{\text{d}z}{\text{d}t}}=w$$

determines the rate of pressure change, where we assume instantaneous adjustment of the parcel's pressure to the ambient conditions. From the ideal gas law,

$$\begin{array}{}\text{(13)}& {\displaystyle }{\mathit{\rho }}_{\mathrm{a}}={\displaystyle \frac{p}{R_{\mathrm{a}}T}},\end{array}$$

where ρ_a is the mass density of the air in the parcel and R_a the specific gas constant of dry air. Differentiation in time gives

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\text{d}{\mathit{\rho }}_{\mathrm{a}}}{\text{d}t}}{\displaystyle }=-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{a}}g}{R_{\mathrm{a}}T}}{\displaystyle \frac{\text{d}z}{\text{d}t}}-{\displaystyle \frac{p}{R_{\mathrm{a}}{T}^{\mathrm{2}}}}{\displaystyle \frac{\text{d}T}{\text{d}t}}.\end{array}$$

Here, $\text{d}T/\text{d}t$ considers the expanding parcel (dry adiabatic lapse rate) and latent heat release following Eq. (9), such that

$$\begin{array}{}\text{(15)}& {\displaystyle \frac{\text{d}T}{\text{d}t}}=-{\displaystyle \frac{g}{c_{\mathrm{pa}}}}{\displaystyle \frac{\text{d}z}{\text{d}t}}+{\displaystyle \frac{L_{\mathrm{v}}}{c_{\mathrm{pa}}}}{\displaystyle \frac{\text{d}{q}_{\mathrm{l}}}{\text{d}t}}.\end{array}$$

4 Numerical Integration

The variety of time constants inherent to most sets of kinetic equations has long been recognized as a barrier to their efficient integration by traditional methods, as stability requirements dictate an integration step size limited by the smallest time constant. From the numerical point of view, atmospheric chemistry is challenging due to the coexistence of very stable (e.g., CH₄) and very reactive (e.g., O(¹D)) species. Mathematically speaking, this phenomenon is better known as stiffness. Therefore, a major task is the integration of stiff systems of ODEs in reasonable time (Sandu and Sander, 2006; Sandu et al., 1997). Since we are dealing with extremely stiff ODE systems, the solution usually relies on robust solvers.

For this purpose, Rosenbrock methods are employed (Hairer et al., 1991). Cminor's implementation allows the user to change the method coefficients via Fortran NAMELIST in order to find the most effective integration scheme for each specific reaction mechanism. A variety of Rosenbrock methods (i.e., sets of coefficients) to choose of comes with the Cminor code already, but the implementation allows the user to easily incorporate different methods or additional solvers.

The novelty in this section lies in the compilation of all of the numerical efforts in one consistent formulation of the Rosenbrock method and describing the system in terms of linear algebra operations. Cminor provides a framework for efficient simulation of pure gas-phase mechanisms, multi-phase mechanisms, polydisperse droplets and aerosol, the droplet condensation equation including chemical effects, the temperature equation for combustion (Wolke and Knoth, 2002; Wolke et al., 2005; Perini et al., 2012), along with a fast mechanism parser, common as well as more complicated rate constants, and the possibility of quick simulations of different scenarios.

4.1 Rosenbrock methods

This section will give a brief introduction to Rosenbrock methods in order to clarify the notation. The starting point is the autonomous initial value problem

$$\begin{array}{}\text{(16)}& {\displaystyle }\dot{\mathit{y}}=\mathit{f}\left(\mathit{y}\right),t>t_{\mathrm{0}},\mathit{y}\left(t_{\mathrm{0}}\right)=(y_{\mathrm{1},\mathrm{0}},\mathrm{\dots },y_{n,\mathrm{0}}{)}^T.\end{array}$$

This assumption does not restrict the problem for non-autonomous systems $\dot{\mathit{y}}=f(t,y)$, because it can easily be transformed by treating the time variable t as dependent, such that $\dot{t}=\mathrm{1}$. The ODE system contains only first order derivatives d$y/$dt, which are usually non-linear functions f of concentration. Most stiff solvers benefit from using implicit formulations such as Rosenbrock methods which arise from diagonal-implicit Runge-Kutta (DIRK) methods. The simplest scheme of this kind is the backward Euler method

$$\begin{array}{}\text{(17)}& {\displaystyle }{\mathit{y}}_{n+\mathrm{1}}{\displaystyle }={\mathit{y}}_n+hf\left[{\mathit{y}}_{n+\mathrm{1}}\right],\end{array}$$

where h is the step size and y_n the current state vector. To avoid implicit systems and iterative solvers, Eq. (17) can be linearized, leading to

$$\begin{array}{}\text{(18)}& {\displaystyle }{\mathit{y}}_{n+\mathrm{1}}{\displaystyle }={\mathit{y}}_n+hf\left[{\mathit{y}}_n\right]+h\left[{\mathit{y}}_{n+\mathrm{1}}-{\mathit{y}}_n\right]f^{\prime }\left[{\mathit{y}}_n\right].\end{array}$$

This is commonly referred to as a linear-implicit method. It can be seen as restricting the iterative solver to one step of Newton's method, see Hairer et al. (1991) for elaborate derivations. To achieve a higher order of consistency, Rosenbrock (1963) proposed to generalize this linearly implicit approach to methods using more stages. The general form of a non-autonomous s-stage Rosenbrock method produces a next step solution as

$$\begin{array}{}\text{(19)}& \begin{array}{rl}{\displaystyle }& {\displaystyle }{\mathit{y}}_{n+\mathrm{1}}={\mathit{y}}_n+{\sum }_{i=\mathrm{1}}^s{b}_i{\mathit{k}}^{\left(i\right)},\\ {\displaystyle }\text{with}& {\displaystyle }{\mathit{k}}^{\left(i\right)}=hf\left[t_n+{\mathit{\alpha }}_{i}h,{\mathit{y}}_n+{\sum }_{j=\mathrm{1}}^{i-\mathrm{1}}{\mathit{\alpha }}_{ij}{\mathit{k}}^{\left(j\right)}\right]\\ {\displaystyle }& {\displaystyle }+{\mathit{\gamma }}_i{h}^{\mathrm{2}}{f}_t\left[t_n,{\mathit{y}}_n\right]+hJ\left[t_n,{\mathit{y}}_n\right]\cdot {\sum }_{j=\mathrm{1}}^i{\mathit{\gamma }}_{ij}{\mathit{k}}^{\left(j\right)},\end{array}\end{array}$$

where the coefficients b_i, α_ij, and γ_ij define the particular method, and are chosen such that a certain accuracy and stability is granted (Zhang et al., 2011). J refers to the Jacobian of the system. For linear-implicit methods the method coefficients have the additional feature that α_ij=0 for all j≥i, γ_ij=0 for all $j\ge i+\mathrm{1}$, and

$$\begin{array}{}\text{(20)}& {\displaystyle }{\mathit{\alpha }}_i=\sum _{j=\mathrm{1}}^{i-\mathrm{1}}{\mathit{\alpha }}_{ij},{\mathit{\gamma }}_i=\sum _{j=\mathrm{1}}^i{\mathit{\gamma }}_{ij}.\end{array}$$

The partial derivative f_t in Eq. (19) is equal to zero for autonomous systems. Reducing the computation cost is crucial for large systems because each time step requires not only s−1 matrix-vector products, but also s matrix decompositions. It is therefore advantageous to choose the same diagonal elements for γ_ii=γ, $i=\mathrm{1},\mathrm{\dots },s$. Thus, only one decomposition is needed per time step. In order to avoid the matrix-vector products on the right-hand side of Eq. (19), the stage vectors k⁽ⁱ⁾ are substituted by ${\mathit{u}}^{\left(i\right)}={\sum }_j^i{\mathit{\gamma }}_{ij}{\mathit{k}}^{\left(j\right)}$ for $i=\mathrm{1},\mathrm{\dots },s$ (see Hairer et al., 1991). The modified Rosenbrock method with s internal stages is then given by

$$\begin{array}{}\text{(21)}& \begin{array}{rl}{\displaystyle }{\mathit{y}}_{n+\mathrm{1}}& {\displaystyle }={\mathit{y}}_n+{\sum }_{j=\mathrm{1}}^s{m}_j{\mathit{u}}^{\left(j\right)},\\ {\displaystyle }\text{with}& {\displaystyle }(\mathbf{I}-h\mathit{\gamma }J){\mathit{u}}^{\left(i\right)}=hf(t_n+{\mathit{\alpha }}_{i}h,{\mathit{y}}_n+{\sum }_{j=\mathrm{1}}^{i-\mathrm{1}}{a}_{ij}{\mathit{u}}^{\left(j\right)})\\ {\displaystyle }& {\displaystyle }+{\mathit{\gamma }}_i{h}^{\mathrm{2}}{f}_t(t_n,{\mathit{y}}_n)+{\sum }_{j=\mathrm{1}}^{i-\mathrm{1}}{d}_{ij}{\mathit{u}}^{\left(j\right)},i=\mathrm{1},\mathrm{\dots },s,\end{array}\end{array}$$

where I is the identity matrix, and Γ=(γ_ii), $\left(a_{ij}\right)=\left({\mathit{\alpha }}_{ij}\right){\mathrm{\Gamma }}^{-\mathrm{1}}$, $\left(d_{ij}\right)=\text{diag}({\mathit{\gamma }}_{\mathrm{11}}^{-\mathrm{1}},\mathrm{\dots },{\mathit{\gamma }}_{ss}^{-\mathrm{1}})-{\mathrm{\Gamma }}^{-\mathrm{1}}$, $\left(m_j\right)=\left(b_j\right){\mathrm{\Gamma }}^{-\mathrm{1}}$, and γ the method parameters.

4.2 Adaptive step size control

Usually, the purpose of adaptive step size control is to achieve some predetermined accuracy in the solution with minimum computational effort. The use of Rosenbrock methods allows us to calculate the local error of the current state vector y. This task can be carried out by computing a second state vector $\widehat{\mathit{y}}$, which is given by an embedded formula where the order of the embedded method $\widehat{p}$ is usually lower than order p of the actual formula ($\widehat{p}=p-\mathrm{1}$). The step size control is implemented by estimating the local error ${\widehat{e}}_j=y_j-{\widehat{y}}_j$, and is then scaled by the denominator ${\text{scal}}_j={\text{tol}}_j^{\mathrm{A}}+max\mathit{\left\{}\right|y_j|,|{\widehat{y}}_j\left|\mathit{\right\}}\cdot {\text{tol}}^{\mathrm{R}}$, where ${\text{tol}}_j^{\mathrm{A}}$ is the absolute and tol^R the relative tolerance given by the user. Note that ${\text{tol}}_j^{\mathrm{A}}$ can be set for each species separately, where tol^R is equal for all species. Cminor's implementation allows to specify two different absolute tolerances, for gaseous and aqueous species. The actual measurement of the error induced by the Rosenbrock method is computed either via the maximum norm or Euclidean norm, denoted by $\parallel \cdot \parallel$. The new step size is then calculated by the following expression:

$$\begin{array}{}\text{(22)}& {\displaystyle }{h}^{\text{new}}=\left\{\begin{array}{l}max\mathit{\{}{h}_{min},h^{\text{old}}max\mathit{\{}\mathrm{0.1},\mathrm{0.8}\parallel \widehat{\mathit{e}}{\parallel }^{-\mathit{\delta }}\mathit{\left\}}\mathit{\right\}},\\ \text{ if }\parallel \widehat{\mathit{e}}\parallel >\mathrm{1},\\ min\mathit{\{}{h}_{max},\mathrm{0.8}{h}^{\text{old}}\parallel \widehat{\mathit{e}}{\parallel }^{-\mathit{\delta }},\mathrm{2}{h}_{\text{old}}\mathit{\}},\\ \text{ otherwise. }\end{array}\right.\end{array}$$

The exponent δ is defined as $\mathit{\delta }=\mathrm{1}/(p+\mathrm{1})$. A new step size h^new is estimated that will yield an error of tol^R on the next step or the next try at taking this step, as the case may be. The step size is multiplied by 0.8 to avoid failures. In addition, the step size h^new is bounded by user-given values for a minimal and a maximal step size, h_min and h_max, respectively.

4.3 Analytical Jacobian matrix formulation

The initial value problem is described by Eqs. (2) and (4), where the dimension of the ODE system is $n=n_{\mathrm{S}}+\mathrm{1}$, with $\dot{\mathit{c}}\in {\mathbb{R}}^{n_{\mathrm{S}}}$ and $\dot{T}\in \mathbb{R}$. The ODE system can then be represented in a very short manner as matrix-vector products,

$$\begin{array}{}\text{(23)}& \begin{array}{rl}{y}^{\prime }& =f(\mathit{c},T)=\left(\begin{array}{c}\dot{\mathit{c}}\\ \dot{T}\end{array}\right)=\left(\begin{array}{c}{\mathit{\nu }}^{T}r+{\mathit{c}}^{\text{emis}}\\ -\frac{\mathrm{1}}{{\overline{c}}_{\mathrm{v}}\mathit{\rho }}{U}^T[{\mathit{\nu }}^{T}r+{\mathit{c}}^{\text{emis}}]\end{array}\right),\\ & y\left(t_{\mathrm{0}}\right)=(c_{\mathrm{1}}^{\mathrm{0}},\mathrm{\dots },c_{n_{\mathrm{S}}}^{\mathrm{0}},T^{\mathrm{0}}{)}^T.\end{array}\end{array}$$

The use of Rosenbrock methods Eq. (21) requires the calculation of the Jacobian matrix of the ODE system Eq. (23),

$$\begin{array}{}\text{(24)}& \begin{array}{rl}J& ={\displaystyle \frac{\partial f(\mathit{c},T)}{\partial (\mathit{c},T)}}=\left[\begin{array}{cccc}\mathrm{\ddots }& & \mathrm{⋰}& \mathrm{⋮}\\ & \frac{\partial {\dot{c}}_{\mathrm{l}}}{\partial c_j}& & \frac{\partial {\dot{c}}_{\mathrm{l}}}{\partial T}\\ \mathrm{⋰}& & \mathrm{\ddots }& \mathrm{⋮}\\ \mathrm{\cdots }& \frac{\partial \dot{T}}{\partial c_j}& \mathrm{\cdots }& \frac{\partial \dot{T}}{\partial T}\end{array}\right],\\ & j=\mathrm{1},\mathrm{\dots },n_{\mathrm{S}}l=\mathrm{1},\mathrm{\dots },n_{\mathrm{S}},\end{array}\end{array}$$

where the first block $\partial {\dot{c}}_{\mathrm{l}}/\partial c_j$ is a square sparse matrix of dimension n_S containing the derivatives of the change of c_l with respect to c_j. The simplification of this part of the Jacobian matrix is related to the assumption that k=k(T) in Eq. (3) depends only on temperature and not on other variables, e.g., time, pressure, or concentrations of catalytic species. This simplification is justified by the assumed slow change of these parameters, such that

$$\begin{array}{}\text{(25)}& \begin{array}{rl}{J}_{\mathrm{cc}}& ={\displaystyle \frac{\partial }{\partial c_j}}\left({\displaystyle \frac{\text{d}{c}_{\mathrm{l}}}{\text{d}t}}\right)=\sum _{i=\mathrm{1}}^{n_{\mathrm{R}}}\left[{\mathit{\nu }}_{il}{\displaystyle \frac{\partial r_i}{\partial c_j}}\right]={\mathit{\nu }}^T{D}_{\mathrm{r}}{\mathit{\nu }}^{\mathrm{e}}{D}_{\mathrm{c}}^{-\mathrm{1}},\\ & j=\mathrm{1},\mathrm{\dots },n_{\mathrm{S}},l=\mathrm{1},\mathrm{\dots },n_{\mathrm{S}},\end{array}\end{array}$$

where $D_{\mathrm{r}}=\text{diag}(r_{\mathrm{1}},\mathrm{\dots },r_{n_{\mathrm{R}}})$ and $D_{\mathrm{c}}^{-\mathrm{1}}=\text{diag}(\mathrm{1}/c_{\mathrm{1}},\mathrm{\dots },\mathrm{1}/c_{n_{\mathrm{S}}})$ are diagonal matrices containing the reactions rates and the inverse concentrations, respectively. The second block $\partial {\dot{c}}_{\mathrm{l}}/\partial T$ consists of a full column vector and represents the derivatives of change in c_l with respect to T. These are expressed as

$$\begin{array}{}\text{(26)}& \begin{array}{rl}{J}_{\mathrm{c}T}& ={\displaystyle \frac{\partial }{\partial T}}\left({\displaystyle \frac{\text{d}{c}_{\mathrm{l}}}{\text{d}t}}\right)=\sum _{i=\mathrm{1}}^{n_{\mathrm{R}}}\left[{\mathit{\nu }}_{il}{\displaystyle \frac{\partial r_i}{\partial T}}\right]={\mathit{\nu }}^T{D}_{\mathrm{r}}\mathbf{\text{K}},\\ & l=\mathrm{1},\mathrm{\dots },n_{\mathrm{S}},\end{array}\end{array}$$

where ${\mathbf{\text{K}}}_i=\partial k_i/\partial T\cdot k_i^{-\mathrm{1}}$ contains the derivatives of the ith rate constants with respect to T, multiplied by the inverse rate constant of reaction i. As mentioned before, k_i is assumed to be only temperature-dependent. The third part of the Jacobian matrix $\partial \dot{T}/\partial c_j$ is a full row vector and contains the derivatives of change in T with respect to c_j. It is expressed as

$$\begin{array}{}\text{(27)}& \begin{array}{rl}{J}_{T\mathrm{c}}& ={\displaystyle \frac{\partial }{\partial c_j}}\left({\displaystyle \frac{\text{d}T}{\text{d}t}}\right)=-{\displaystyle \frac{\mathrm{1}}{{\overline{c}}_{\mathrm{v}}\mathit{\rho }}}[C_{\mathrm{v},j}\dot{T}+\sum _{l=\mathrm{1}}^{n_{\mathrm{S}}}{U}_j{\displaystyle \frac{\partial {\dot{c}}_l}{\partial c_j}}],\\ & j=\mathrm{1},\mathrm{\dots },n_{\mathrm{S}}\\ & =-{\displaystyle \frac{\mathrm{1}}{{\overline{c}}_{\mathrm{v}}\mathit{\rho }}}[{\mathbf{\text{C}}}_{\mathrm{v}}\dot{T}+U^T{\mathit{\nu }}^T{D}_{\mathrm{r}}{\mathit{\nu }}^{\mathrm{e}}{D}_{\mathrm{c}}^{-\mathrm{1}}],\end{array}\end{array}$$

where $C_{\mathrm{v},j}=\partial U_j/\partial T$ is the constant volume specific heat of the jth species. The first item in the Jacobian matrix, i.e., the entry in position $(n_{\mathrm{S}}+\mathrm{1},n_{\mathrm{S}}+\mathrm{1})$, is the derivative of the average mixture temperature change rate with respect to the system temperature itself,

$$\begin{array}{}\text{(28)}& \begin{array}{rl}{J}_{TT}& ={\displaystyle \frac{\partial }{\partial T}}\left({\displaystyle \frac{\text{d}T}{\text{d}t}}\right)=\\ & -{\displaystyle \frac{\mathrm{1}}{{\overline{c}}_{\mathrm{v}}\mathit{\rho }}}\mathit{\{}\mathit{\rho }\dot{T}{\displaystyle \frac{\partial {\overline{c}}_{\mathrm{v}}}{\partial T}}+\sum _{l=\mathrm{1}}^{n_{\mathrm{S}}}[C_{\mathrm{v},l}{\dot{c}}_l+U_l{\displaystyle \frac{\partial {\dot{c}}_l}{\partial T}}\left]\mathit{\right\}}=\\ & -{\displaystyle \frac{\mathrm{1}}{{\overline{c}}_{\mathrm{v}}\mathit{\rho }}}\left[\mathit{\rho }\dot{T}{\displaystyle \frac{\partial {\overline{c}}_{\mathrm{v}}}{\partial T}}+{\mathbf{\text{C}}}_{\mathrm{v}}\dot{c}+U^T{\mathit{\nu }}^T{D}_{\mathrm{r}}\mathbf{\text{K}}\right].\end{array}\end{array}$$

4.4 Integration scheme for box-model chemistry

In this section, the resulting numerical method using the classical linear algebra formulation is shown. The ODE system Eq. (23) is assumed to be autonomous. Substituting the values of J, Eq. (24), into the Rosenbrock method Eq. (21) yields the classic linear algebra formulation for computing a new state vector. Note that the partial derivative f_t, see Sect. 4.1, will be dropped. Also, still, for box-model atmospheric chemistry systems, the equation for T will be dropped. Consequently, for atmospheric systems, no dependencies of the rate constants k_i are considered in the numerical procedure, i.e., the system is assumed to be fully autonomous from the numerical point of view. Still, rate constants of photolytic reactions are updated during integration, according to current time and the specified location, i.e., zenith angle of the sun. Accordingly, this gives

$$\begin{array}{}\text{(29)}& \begin{array}{rl}& {\mathit{y}}_{n+\mathrm{1}}={\mathit{y}}_n+\sum _{i=\mathrm{1}}^s{m}_i\left[\begin{array}{c}{\mathit{u}}^{\left(i\right)}\\ {\stackrel{\mathrm{̃}}{u}}^{\left(i\right)}\end{array}\right],\text{with}\\ & \left[\begin{array}{cc}I-h\mathit{\gamma }{\mathit{\nu }}^T{D}_{\mathrm{r}}{\mathit{\nu }}^{\mathrm{e}}{D}_{\mathrm{c}}^{-\mathrm{1}}& -h\mathit{\gamma }{\mathit{\nu }}^T{D}_{\mathrm{r}}\mathbf{\text{K}}\\ -\frac{h\mathit{\gamma }}{{\overline{c}}_{\mathrm{v}}\mathit{\rho }}({\mathbf{\text{C}}}_{\mathrm{v}}\dot{T}+U^T{J}_{\mathrm{cc}})& \mathrm{1}+\frac{h\mathit{\gamma }}{{\overline{c}}_{\mathrm{v}}\mathit{\rho }}(\mathit{\rho }\frac{\partial {\overline{c}}_{\mathrm{v}}}{\partial T}\dot{T}+{\mathbf{\text{C}}}_{\mathrm{v}}\dot{c}+U^T{J}_{\mathrm{c}T})\end{array}\right]\\ & \left[\begin{array}{c}{\mathit{u}}^{\left(i\right)}\\ {\stackrel{\mathrm{̃}}{u}}^{\left(i\right)}\end{array}\right]=h\left[\begin{array}{c}{\mathit{\nu }}^T{r}^{\left(i\right)}+{\mathit{c}}^{\text{emis}}\\ -\frac{\mathrm{1}}{{\overline{c}}_{\mathrm{v}}^{\left(i\right)}\mathit{\rho }}\left(U^{\left(i\right)}{)}^T\right({\mathit{\nu }}^T{r}^{\left(i\right)}+{\mathit{c}}^{\text{emis}})\end{array}\right]\\ & +\sum _{j=\mathrm{1}}^{i-\mathrm{1}}{d}_{ij}\left[\begin{array}{c}{\mathit{u}}^{\left(j\right)}\\ {\stackrel{\mathrm{̃}}{u}}^{\left(j\right)}\end{array}\right],i=\mathrm{1},\mathrm{\dots },s.\end{array}\end{array}$$

The values of both diagonal matrices D_r and D_c are calculated at the beginning of each time step (i.e., at y⁽ⁿ⁾), meaning that they do not form intermediates. In contrast, the reaction rates r_i, the mass-averaged specific heat at constant volume ${\overline{c}}_{\mathrm{v},i}$, and the internal energy values $U_i^T$ on the right-hand-side form the intermediate solution for each stage

$$\begin{array}{ll}{\displaystyle }{r}^{\left(i\right)}& {\displaystyle }=r^{\left(i\right)}\left(t_n+{\mathit{\alpha }}_{i}h,{\mathit{y}}_n+\sum _{j=\mathrm{1}}^{i-\mathrm{1}}{a}_{ij}\left[\begin{array}{c}{\mathit{u}}^{\left(i\right)}\\ {\stackrel{\mathrm{̃}}{u}}^{\left(i\right)}\end{array}\right]\right),\\ {\displaystyle }& {\displaystyle }{\overline{c}}_{\mathrm{v}}^{\left(i\right)}={\overline{c}}_{\mathrm{v}}^{\left(i\right)}\left({\mathit{y}}_n+\sum _{j=\mathrm{1}}^{i-\mathrm{1}}{a}_{ij}\left[\begin{array}{c}{\mathit{u}}^{\left(i\right)}\\ {\stackrel{\mathrm{̃}}{u}}^{\left(i\right)}\end{array}\right]\right),\\ {\displaystyle }& {\displaystyle }{U}^{\left(i\right)}=U^{\left(i\right)}\left(T+\sum _{j=\mathrm{1}}^{i-\mathrm{1}}{a}_{ij}{\stackrel{\mathrm{̃}}{u}}^{\left(i\right)}\right)\end{array}$$

Note that the value for T is located in the last entry in the state vector y⁽ⁿ⁾, see Eq. (23).

4.5 Direct sparse linear solver

In this work, direct linear solvers are considered to calculate the solution of the linear systems in the integration procedure. For reasons of simplicity, the linear systems are denoted as Ax=b in this section. A direct solver strategy requires the factorization of the coefficient matrix A into a lower and an upper triangular matrix L and U, where $\mathbf{A}=\mathbf{L}\cdot \mathbf{U}$. The implemented sparse solver consists of three parts, (i) the symbolic phase, where the actual factor matrix LU(A) is built and stored as separate sparse matrix as it has the same structure (i.e., non-zero pattern) in every time step, (ii) the numerical factorization within the Rosenbrock method, and (iii) a solution phase for the triangular systems. While factorizing sparse matrices, generally fill-in occurs and the factors of LU(A) become more dense than the original matrix A. However, to take advantage of the sparseness of the coefficient matrix A, the equations must be arranged in special order. To find a good ordering for the sparse matrix A, we have to determine a permutation P that minimizes the fill-in in the factors L and U. Since the problem of finding a permutation P that minimizes the amount of fill-in can be reduced to an NP-complete task (Yannakakis, 1981), finding the optimal solution becomes computationally infeasible for larger systems due to the exponential growth of required computational resources. Yet, Cminor uses the minimum-degree ordering heuristic proposed by Markowitz (1957) in a symmetric manner, meaning the pivot always lies on the diagonal to preserve existing stability. Thus, if rows i and j are swapped, then columns i and j are swapped as well. The (symbolic) elimination process is denoted in Algorithm 1. A is assumed to be unsymmetric and regular, which does generally not ensure that a solution exists. Nevertheless, numerical pivoting is not implemented in this approach. The numerical phase (ii) takes place once at each time step within the Rosenbrock method. By using the Markowitz (1957) strategy, the growth of the number of non-zero elements is roughly bounded by the factor 2.

Algorithm 1 Gaussian Elimination with a Minimum-Degree heuristic: Minimizing fill-in by selecting pivots with the smallest row-column non-zero product, guided by a restriction vector Restr to prioritize certain species for permutation to the final rows and columns.

1:

procedure Eliminate_MD(A₀)

2:

for $(i=\mathrm{1}:$ dim(A₀)) do

3:

# count non-zero elements of every row and column in sub-matrix A_i

4:

# choose pivot element: $j_{\mathrm{p}}=min\mathit{\left\{}\right[$row$\left(k\right)-\mathrm{1}]\cdot [$col(k)−1]} for $\left\{\begin{array}{ll}k\ge i,k\notin \text{Restr}& \text{if }\exists k\ge i:k\notin \text{Restr}\\ k\ge i& \text{else}\end{array}\right.$

5:

# swap row and column i↔j_p

6:

# elimination of column j_p in sub-matrix A_i

7:

end for

8:

end procedure

4.6 Vectorized LU decomposition

For multi-phase systems in a box-model, droplets of the same size will behave equally as long as the initial concentrations are equal in each droplet. If this is not the case, as is for specifying log-normally distributed aerosol, each aqueous species has to be considered once for each droplet class. A droplet class refers to a number of droplets of equal size and chemical composition. This corresponds to the concept of superdroplets commonly used in cloud microphysics (Shima et al., 2009; Hoffmann et al., 2015).

When considering multiple droplet classes, the number of equations in atmospheric systems grows from

$${\displaystyle }{n}_{\text{Gas}}+n_{\text{Henry}}+n_{\text{Aq}}\text{ to }{n}_{\text{Gas}}+n_{\text{Henry}}+n_D\cdot n_{\text{Aq}},$$

where n_Gas is the number of purely gaseous species, i.e., not dissolving, n_Henry the number of gaseous species that may dissolve in droplets, n_Aq the number of aqueous species in the mechanism and n_D the number of droplet classes. This is a significant increase for mechanisms with a prominent aqueous phase. For instance, a number of 50 droplet classes is common in super-droplet models. Since the equations for each droplet class are equal in structure, as the same reactions occur in every droplet class, the matrix A has a specific, repeating pattern. Examples are shown in Fig. 3, displaying the sparsity pattern of the iteration matrix of a sulfur oxidation mechanism in clouds for one (a) and five droplet classes (b, c). The repeating pattern, in turn, can be represented by considering the matrix for one droplet class, with vectors representing the entries of aqueous species and mass transfer terms.

Figure 3 Sparsity patterns: iteration matrix of sulphate oxidation mechanism in clouds; one droplet class (a), five droplet classes ordered by droplet classes (b), five droplet classes ordered by species (c).

By doing so, the (symbolic) matrix and its (symbolic) LU decomposition stay the same for any number of droplet classes. Therefore, only the solution phase (iii) changes from scalar operations to vector operations. These are fast and minimize indirect accesses to the respective arrays, which make up a major portion of the total solution process, and therefore also of the total simulation. Of course, since there are five droplet classes in this example, the matrix to be decomposed is (b) or (c), in Fig. 3, which is just a matter of ordering.

In the following, we explore whether the decomposition can be determined using only (a). As shown in (c), the vector structure can be represented by arranging all entries corresponding to aqueous species into diagonal matrix blocks. If every zero entry modification (fill-in), i.e., change of a zero entry to a non-zero entry while decomposing the matrix, is part of a new diagonal matrix block, then it is possible to work with (a) to find the decomposition of (c). Consequently, when a diagonal block matrix transforms into a fully populated matrix during the decomposition process, the vector representation becomes impractical. In such iteration matrices, the vectorized decomposition approach is no longer viable. Also, any single entry appearing off-diagonal in a diagonal matrix block prevents the vector representation of the matrix as its position would have to be stored.

Every row in the iteration matrix corresponds to an equation of a specific species, while a column resembles the influence of a specific species on the other species. This means, every row and every column relate to one specific species. By permuting the rows and columns representing Henry species, i.e., species that are gaseous but may dissolve in droplets, to the last rows and columns, the vector structure will be preserved during decomposition. Ordered like this, the iteration matrix for multi-phase mechanisms exhibits the four regions

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathbf{A}=\\ {\displaystyle }& {\displaystyle }\left(\begin{array}{cc}\text{in-phase reactions}& \begin{array}{c}\mathrm{\text{Henry species}}\\ \to \text{aqueous or purely gaseous species}\end{array}\\ \begin{array}{c}\mathrm{\text{aqueous or purely gaseous species}}\\ \to \text{Henry species}\end{array}& \begin{array}{c}\mathrm{\text{Henry species}}\\ \to \text{Henry species}\end{array}\end{array}\right),\end{array}$$

where purely gaseous species means not dissolving in droplets.

In the following, the preservation of the vector structure by this ordering is reasoned by an induction approach. It is shown that, having a matrix which can be represented with vector entries and which is ordered in the proposed way, a new fill-in will always occur vector-wise, preserving the possibility to represent the matrix with vector entries. The following Algorithm 2 generates the LU decomposition of a matrix A.

Algorithm 2 LU Decomposition: Generating a lower triangular and an upper triangular matrix, whose product yields A. The matrices L and U are stored in the upper and lower part of the over-written matrix A.

1:

for  i = 1 : nSpc  do

2:

for  j = 1 : i-1  do

3:

A[i,j] = A[i,j]/A[j,j]

4:

for  k = j+1 : nSpc  do

5:

A[i,k] = A[i,k] - A[i,j] ⋅ A[j,k]

6:

end for

7:

end for

8:

end for

It can be seen that a necessary condition for fill-in at entry A[i,k] is

$$\begin{array}{}\text{(30)}& {\displaystyle }\exists j>\mathrm{0},j<i,j<k:\mathbf{A}[i,j]\cdot \mathbf{A}[j,k]\ne \mathrm{0}.\end{array}$$

With this necessary condition, the logical path of the following reasoning is

fill-in → Eq. (30) → appropriate fill-in (vector-wise where needed).

The bottom left and upper right regions of A can in any case be represented by vector entries as the vector representation only prevents storing zeros in the upper left region. The lower right region has scalar entries anyhow, representing gaseous reactions of dissolving species (not shown in Fig. 3). Therefore, it is sufficient to show that if fill-in occurs for the case of i and k both representing aqueous species, fill-in in the form of a diagonal matrix block follows.

Let $i,j,k>\mathrm{0}$ be a triple of indices, for which $\mathbf{A}[i,j]\cdot \mathbf{A}[j,k]\ne \mathrm{0}$ holds for a given matrix A, where the gaseous, dissolving species are permuted to the last rows and columns. If $i,k\le n_{\text{Aq}}\cdot n_{\text{D}}+n_{\text{Gas}}$, the entry A[i,k] is in the region of either diagonal matrix entries, i.e., aqueous species, or single-value entries, i.e., purely gaseous species. As mentioned above, this is the only important case for which preservation of the vector structure of A needs to be shown. Any entry $\mathbf{A}[i,j]\ne \mathrm{0}$ comes from a reaction of species j to species i (or an equivalent previous fill-in). Since there are no reactions between purely gaseous species and aqueous species, it follows that $i,j,$ and k all belong to the same type of species, either purely gaseous or aqueous. It also follows that fill-in can never occur between aqueous and purely gaseous species. If they are gaseous, no vector structure has to be preserved, in this case the fill-in is irrelevant. If they are aqueous, $i,j,$ and k can be written as

$$\begin{array}{}\text{(31)}& \begin{array}{rl}{\displaystyle }i& {\displaystyle }=n_{\text{Gas}<i}+(S_i-\mathrm{1})\cdot n_{\text{D}}+d,\\ {\displaystyle }j& {\displaystyle }=n_{\text{Gas}<j}+(S_j-\mathrm{1})\cdot n_{\text{D}}+d,\\ {\displaystyle }k& {\displaystyle }=n_{\text{Gas}<k}+(S_k-\mathrm{1})\cdot n_{\text{D}}+d,\end{array}\end{array}$$

where n_Gas<i is the number of purely gaseous species that were permuted before species i, $S_i\in \mathit{\{}\mathrm{1},\mathrm{\dots },n_{\text{Aq}}\mathit{\}}$ is the aqueous species number disregarding droplet classes, and $d\in \mathit{\{}\mathrm{1},\mathrm{\dots },n_{\text{D}}\mathit{\}}$ is the droplet class in which species i is dissolved. The analagous nomenclature is used for the variables containing j and k. Note that d is the same for $i,j,$ and k. Before fill-in, the matrix is assumed to have proper vector entries, and we re-iterate that there is no reaction from one droplet class to another. Thus, d is the same for $i,j,$ and k. Now, w.l.o.g., it is assumed that d=1. For other values of d, analogous calculations yield the same result. The vector structure is preserved if

$$\begin{array}{}\text{(32)}& \begin{array}{rl}& \forall l_{\mathrm{1}},l_{\mathrm{2}}\in \mathit{\{}\mathrm{0},\mathrm{\dots },n_{\text{D}}-\mathrm{1}\mathit{\}},l_{\mathrm{1}}\ne l_{\mathrm{2}},\forall j^{*}>\mathrm{0},j^{*}<i^{*}\\ & :=i+l_{\mathrm{1}},j^{*}<k^{*}:=k+l_{\mathrm{2}}:\mathbf{A}[i^{*},j^{*}]\cdot \mathbf{A}[j^{*},k^{*}]=\mathrm{0}.\end{array}\end{array}$$

Under this condition, the entry $\mathbf{A}[i^{*},k^{*}]$ will not undergo fill-in, i.e., every off-diagonal entry in the submatrix

$${\displaystyle }\mathbf{A}[i:i+n_{\text{D}}-\mathrm{1},k:k+n_{\text{D}}-\mathrm{1}]$$

stays equal to zero. Then, the possible fill-in is a diagonal matrix, which is what needs to be shown. The updated representation of Eq. (31), which is

$$\begin{array}{rl}{i}^{*}=i+l_{\mathrm{1}}& =n_{\text{Gas}<i}+(S_i-\mathrm{1})\cdot n_{\text{D}}+d+l_{\mathrm{1}},\\ k^{*}=k+l_{\mathrm{2}}& =n_{\text{Gas}<k}+(S_k-\mathrm{1})\cdot n_{\text{D}}+d+l_{\mathrm{2}},\end{array}$$

helps recognizing that, since l₁≠l₂ by Eq. (32), species i^* and k^* now belong to separate droplet classes, namely d+l₁ and d+l₂, respectively. Since the matrix is assumed to have proper vector entries before the examined fill-in, it can only have non-zero values at entries which belong to the same droplet class. Any $j^{*}>\mathrm{0}$ with $j^{*}<i^{*}$ and $j^{*}<k^{*}$ may correspond to an aqueous species of a specific droplet class, where either $\mathbf{A}[i^{*},j^{*}]$ or $\mathbf{A}[j^{*},k^{*}]$ or both are equal to zero because of the aforementioned reason, or it may correspond to a purely gaseous species, where the respective matrix entries are zero, too, because purely gaseous species do not interact with the aqueous species i^* and k^*. Thus, Eq. (32) holds.

Hence, the vector structure of A, e.g., the rightmost matrix in Fig. 3, is preserved during decomposition for the proposed ordering, which is to permute the gaseous dissolving species to the last rows and columns. In opposition to this, if i represents a Henry species, which is permuted between aqueous and purely gaseous ones, Eq. (32) is violated for l₁=1, l₂=0 and $j^{*}=j$. In this case, the vector structure is not preserved.

4.7 Numerical treatment of parcel equations

For the Jacobian J, we need derivatives of Eqs. (9), (14), and (15) with respect to the rest of the variables. To minimize floating point operations, we neglect these derivatives, resulting in a sort of explicit scheme for these equations, which we assume is reasonable as they behave not as stiff as the chemical species' equations. The only exception is for the condensational equation, where we calculate the derivative with respect to the water masses, where $i,j\in \mathit{\{}\mathrm{1},\mathrm{\dots },n_{\mathrm{D}}\mathit{\}}$ and n_D is again the number of droplet classes. Here,

$$\begin{array}{ll}{\displaystyle \frac{\partial }{\partial q_{\mathrm{l},j}}}\left({\displaystyle \frac{\text{d}{q}_{\mathrm{l},i}}{\text{d}t}}\right)& {\displaystyle }={\displaystyle \frac{\partial }{\partial q_{\mathrm{l},j}}}\left[D_i\cdot (S-S_{\mathrm{eq},i})\right]\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\partial D_i}{\partial q_{\mathrm{l},j}}}\cdot (S-S_{\mathrm{eq},i})+D_i\cdot {\displaystyle \frac{\partial (S-S_{\mathrm{eq}})}{\partial q_{\mathrm{l},j}}}\\ \text{(33)}& {\displaystyle }& {\displaystyle }\approx D_i\cdot {\displaystyle \frac{\partial (S-S_{\mathrm{eq},i})}{\partial q_{\mathrm{l},j}}}.\end{array}$$

Since supersaturations in liquid clouds are small, we neglect $(\partial D_i/\partial q_{\mathrm{l},j})\cdot (S-S_{\mathrm{eq},i})$, and D_i summarizes all remaining factors in Eq. (9). The remaining term splits up in the derivatives

$$\begin{array}{}\text{(34)}& {\displaystyle \frac{\partial S}{\partial q_{\mathrm{l},j}}}\text{ and }-{\displaystyle \frac{\partial S_{\mathrm{eq},i}}{\partial q_{\mathrm{l},j}}}.\end{array}$$

The first term is zero since S has no direct relation to the liquid water mixing ratio. For the second term, we omit the index indicating the droplet class because for i≠j the term is zero. Also, we again summarize factors to R and K, the Raoult and the Kelvin term, respectively. Thus,

$$\begin{array}{}\text{(35)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial q_{\mathrm{l}}}}{S}_{\mathrm{eq}}& {\displaystyle }={\displaystyle \frac{\partial }{\partial q_{\mathrm{l}}}}(R\cdot K)\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\partial }{\partial q_{\mathrm{l}}}}\left[{\displaystyle \frac{q_{\mathrm{l}}/M_{\mathrm{w}}}{q_{\mathrm{l}}/M_{\mathrm{w}}+n_{\mathrm{s}}}}\cdot \exp \left({\displaystyle \frac{a}{r_{\mathrm{d}}}}\right)\right]\\ {\displaystyle }& {\displaystyle }=R^{\prime }\cdot K+R\cdot K^{\prime }\end{array}\end{array}$$

The derivative of the Raoult term with respect to water mass is easily calculated as

$$\begin{array}{}\text{(36)}& {\displaystyle }{R}^{\prime }={\displaystyle \frac{n_{\mathrm{s}}}{M_{\mathrm{w}}(q_{\mathrm{l}}/M_{\mathrm{w}}+n_{\mathrm{s}}{)}^{\mathrm{2}}}}.\end{array}$$

For the derivative of the Kelvin term, we note that the radius r_d might indicate the radius of a droplet not only consisting of water, but also soluble and insoluble material. However, assuming ideal solution behavior for droplets of a size larger than a few nanometers, the volume of a single droplet might be calculated as $V=V_{\mathrm{w}}+V_{\text{rest}}$ with $V_{\mathrm{w}}=q_{\mathrm{l}}/\left({\mathit{\rho }}_{\mathrm{w}}n\right)$, n is the number of droplets (in the droplet class) per kg of air, and V_rest the volume of dissolved species. Doing so, we get

$$\begin{array}{}\text{(37)}& \begin{array}{rl}{\displaystyle }{K}^{\prime }& {\displaystyle }={\displaystyle \frac{\partial \exp \left(\frac{a}{r_{\mathrm{d}}}\right)}{\partial r_{\mathrm{d}}}}\cdot {\displaystyle \frac{\partial r_{\mathrm{d}}}{\partial V_{\mathrm{w}}}}\cdot {\displaystyle \frac{\partial V_{\mathrm{w}}}{\partial q_{\mathrm{l}}}}\\ {\displaystyle }& {\displaystyle }=\left[-{\displaystyle \frac{a}{r_{\mathrm{d}}^{\mathrm{2}}}}\exp \left({\displaystyle \frac{a}{r_{\mathrm{d}}}}\right)\right]\cdot {\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}{\left[{\displaystyle \frac{\mathrm{3}}{\mathrm{4}\mathit{\pi }}}(V_{\mathrm{w}}+V_{\text{rest}})\right]}^{-\frac{\mathrm{2}}{\mathrm{3}}}\cdot {\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{w}}n}}\\ {\displaystyle }& {\displaystyle }=-{\displaystyle \frac{a}{\mathrm{4}\mathit{\pi }{\mathit{\rho }}_{\mathrm{w}}n{r}_{\mathrm{d}}^{\mathrm{4}}}}\exp \left({\displaystyle \frac{a}{r_{\mathrm{d}}}}\right).\end{array}\end{array}$$

The resulting classical linear algebra scheme is

$$\begin{array}{}\text{(38)}& \begin{array}{rl}& {\mathit{y}}_{n+\mathrm{1}}={\mathit{y}}_n+\sum _{i=\mathrm{1}}^s{m}_i\left[\begin{array}{c}{\mathit{u}}^{\left(i\right)}\\ {\stackrel{\mathrm{̃}}{u}}^{\left(i\right)}\end{array}\right],\text{with}\\ & \left[\begin{array}{cc}I-h\mathit{\gamma }{\mathit{\nu }}^T{D}_{\mathrm{r}}{\mathit{\nu }}^{\mathrm{e}}{D}_{\mathrm{c}}^{-\mathrm{1}}& \mathbf{0}\\ \mathbf{0}& I-h\mathit{\gamma }{J}_{\text{parcel}}\end{array}\right]\left[\begin{array}{c}{\mathit{u}}^{\left(i\right)}\\ {\stackrel{\mathrm{̃}}{\mathit{u}}}^{\left(i\right)}\end{array}\right]=\\ & h\left[\begin{array}{c}{\mathit{\nu }}^T{r}^{\left(i\right)}+{\mathit{c}}^{\text{emis}}+\mathit{c}\frac{\dot{\mathit{\rho }}}{\mathit{\rho }}\\ {\dot{\mathit{q}}}_{\mathrm{l}}^{\left(i\right)}\\ -\frac{g}{c_{\mathrm{pa}}}w+\frac{L_{\mathrm{v}}}{c_{\mathrm{pa}}}{\dot{q}}_{\mathrm{l}}^{\left(i\right)}\\ -{\dot{q}}_{\mathrm{l}}^{\left(i\right)}\\ -\frac{{\mathit{\rho }}_{\mathrm{a}}g}{R_{\mathrm{a}}T}w-\frac{p}{R_{\mathrm{a}}{T}^{\mathrm{2}}}\left(-\frac{g}{c_{\mathrm{pa}}}w+\frac{L_{\mathrm{v}}}{c_{\mathrm{pa}}}{\dot{q}}_{\mathrm{l}}^{\left(i\right)}\right)\\ w\end{array}\right]\\ & +\sum _{j=\mathrm{1}}^{i-\mathrm{1}}{d}_{ij}\left[\begin{array}{c}{\mathit{u}}^{\left(j\right)}\\ {\stackrel{\mathrm{̃}}{\mathit{u}}}^{\left(j\right)}\end{array}\right],i=\mathrm{1},\mathrm{\dots },s,\end{array}\end{array}$$

where ${\stackrel{\mathrm{̃}}{\mathit{u}}}^{\left(j\right)}$ denotes the jth intermediate of the transformed parcel variables,

$$\begin{array}{}\text{(39)}& {\displaystyle }{\dot{\mathit{q}}}_{\mathrm{l}}^{\left(i\right)}=\text{diag}\left({\mathbf{D}}^{\left(i\right)}\right)(S^{\left(i\right)}\cdot \mathbf{1}-{\mathbf{S}}_{\mathrm{eq}}^{\left(i\right)})\end{array}$$

and

$$\begin{array}{}\text{(40)}& {\displaystyle }{J}_{\text{parcel}}=\left(\begin{array}{ccccc}{\mathbf{\text{J}}}_{mm}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -{\mathbf{\text{J}}}_{mm}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right),\end{array}$$

with

$$\begin{array}{}\text{(41)}& {\mathbf{J}}_{mm}=-\text{diag}\left(\mathbf{D}\right)\left[\text{diag}\left({\mathbf{R}}^{\prime }\right)\mathbf{K}+\text{diag}\left(\mathbf{R}\right){\mathbf{K}}^{\prime }\right].\end{array}$$

The bold letters indicate vectors containing values of all droplet classes. This represents an approximate Jacobian matrix. Note that the concept of vector entries of Sect. 4.6 is applied here (if n_D>1). Wherever no index to indicate the stage value is written, the value of the first stage is taken, i.e., i=1, as the Jacobian is only calculated at the beginning of a time step. Again, the reaction rates r⁽ⁱ⁾ form the intermediate solution for each stage as

$${\displaystyle }{r}^{\left(i\right)}=r\left(t_n+{\mathit{\alpha }}_{i}h,{\mathit{y}}_n+\sum _{j=\mathrm{1}}^{i-\mathrm{1}}{a}_{ij}\left[\begin{array}{c}{\mathit{u}}^{\left(i\right)}\\ {\stackrel{\mathrm{̃}}{\mathit{u}}}^{\left(i\right)}\end{array}\right]\right).$$

5 Evaluation and First Results

Cminor is evaluated by the numerical simulation of six different chemical mechanisms in a box-model setup, followed by a comparison to two other chemical solvers with a box-model gas-phase simulation, and a simulation of a rising adiabatic parcel. Five reaction mechanisms are derived from atmospheric chemistry modeling, and three mechanisms from the field of combustion chemistry, which are used in internal combustion engine CFD, autoignition and flame simulations. These simulations are meant to present the functionality of the solver. All results can be reproduced by a python script, provided in the repository Rug et al. (2025). Tables 3 and 4 provide an overview of the conditions simulated and the references to the chemical mechanisms. All simulations are carried out on a MacBook Pro 2015 with a 2.5 GHz Quad-Core Intel Core i7 processor.

The first atmospheric mechanism is the Chapman mechanism (Sandu and Sander, 2006), using 7 species and 10 reactions to explain the presence of the ozone layer in earth's stratosphere. The results are shown in Fig. 4a.

Sandu and Sander (2006)Stockwell et al. (1997)Jenkin et al. (2012)Kreidenweis et al. (2003)Ervens et al. (2003)Bräuer (2015)

Table 3 Mechanism features for the four atmospheric (multi-phase) systems. For multiple droplet classes, aqueous species and reactions need to be considered once for every droplet class in the ODE system, due to different concentrations. The resulting number of species and reactions in the ODE system is shown in parenthesis, if multiple droplet classes are present.

Patel et al. (2004)Seiser et al. (2000)Herbinet et al. (2008)

Table 4 Mechanism features for the three constant-volume reactor combustion systems.

Figure 4 Numerical simulation results of three atmospheric (multi-phase) chemistry mechanisms for a one day time interval, with panel (a) presenting the Chapman, (b) the RACM+CAPRAMv2.4, and (c) the MCMv3.2+CAPRAMv4.0α mechanism. To visualize all species' concentrations in a single plot, the concentrations have been normalized by individual scaling factors, as indicated in the legend. The blue columns indicate the formation of cloud droplets, which is required for aqueous-phase reactions. CPU time is the total time spent for integration with a relatively strict relative error tolerance of $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{4}}$.

Figure 5 Numerical simulation results of three combustion mechanisms: (a) ERC n-heptane, (b) LLNL n-heptane, and (c) LLNL methyl-decanoate. CPU time is the total time spent for integration with a relatively strict relative error tolerance of $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{4}}$.

The second atmospheric mechanism, RACM+CAPRAMv2.4, consists of the gas phase Regional Atmospheric Chemistry Model (RACM) by Stockwell et al. (1997), coupled with the Chemical Aqueous Phase Radical Mechanism version 2.4 (CAPRAM2.4) by Ervens et al. (2003). RACM is an upgrade of previously developed mechanisms for acid deposition in the troposphere, and intends to be valid for remote to polluted conditions in the whole troposphere. It is still a condensed mechanism, i.e., not a near-explicit description of all resolved processes, but it resembles thoroughly revised inorganic as well as organic chemistry. The aqueous phase mechanism CAPRAM2.4 contains inorganic, extended organic, and transition metal chemistry. Simulation results have been analyzed earlier by Tilgner et al. (2008), including time resolved source and sink studies focusing particularly on multiphase phase processing of radical oxidants and of C2–C4 organic compounds. Results are presented in Fig. 4b.

The third atmospheric mechanism, MCMv3.2+CAPRAMv4.0α, consists of the detailed Master Chemical Mechanism version 3.2 (MCMv3.2) developed by Jenkin et al. (2012), which is a near-explicit, pure gas phase mechanism of the troposphere. Coupled with version 4.0α of the Chemical Aqueous Phase Radical Mechanism (CAPRAM4.0α) (Bräuer, 2015), it represents a detailed documentation of chemical processes in the gaseous and aqueous phase of the troposphere. The Master Chemical Mechanism was developed to investigate the degradation of emitted volatile organic compounds (VOCs), which have a major influence on the chemistry of the troposphere, contributing to the formation of ozone, secondary organic aerosol (SOA), and other secondary pollutants. CAPRAM4.0α is the most comprehensive version of the CAPRAM mechanism and is also used to study oxidant budgets and SOA formation in the aqueous phase. Figure 4c shows the corresponding results.

In homogeneous charge compression ignition engines, the mixture of fuel and oxidizer is compressed to the point of auto-ignition, creating heat to be transformed into work by the engine. Investigating and optimizing such engines requires deep understanding of the ongoing chemistry and intermediate chemical states of the fuel as well as of physical properties, especially temperature. The following three combustion mechanisms were chosen because they were used as representatives for combustion chemistry in Perini et al. (2012).

Figure 6 Comparison of eight species of a simulation of the MCMv3.2 mechanism with KPP and Cminor. Black lines show simulation results of Cminor, yellow dots indicate the results of the simulation using KPP. Red dots indicate a 2 % relative difference between the results of KPP and Cminor, blue dots indicate 5 % relative difference, and cyan dots 10 % relative difference between the two models.

Being used as a reference fuel, n-heptane studies are crucial and a number of mechanisms exist to simulate ignition delay times, heat release rates, and more, some being more reduced versions and some containing more detailed descriptions of the chemical processes. Two are shown here, a strongly reduced version from the Engine Research Center (ERC) at the University of Wisconsin-Madison (Patel et al., 2004), and a slightly larger mechanism developed at the Lawrence Livermore National Laboratory (LLNL) in California (Seiser et al., 2000). The results are shown in Fig. 5a and b. The third combustion mechanism, also developed at LLNL (Herbinet et al., 2008), describes the consumption of an alternative fuel, driven by methyl decanoate. Such large methyl esters are characteristic for rapeseed or soybean derived biodiesels and such fuels show different behavior, e.g., early formation of carbon dioxide, which this larger mechanism is able to reproduce. It serves as an example for state-of-the-art comprehensive combustion mechanisms and results are shown in Fig. 5c.

The last two mechanisms are used to validate the correctness of the solver. The MCMv3.2 (Jenkin et al., 2012) is used to analyze Cminor’s ability to solve pure chemistry. The mechanism is simulated using Cminor and KPP with the same initial and environmental conditions, numerical methods, and tolerances. The results match, while the times to perform a simulation differ, demonstrating different usages of the two solvers. Also, KPP has been maintained and optimized over many years and creates a source code perfectly fitting one specific mechanism. Cminor has unused optimization potential in matrix-vector operations and possibly elsewhere. The simulation times are 4.55 and 3.08 s for Cminor and KPP, respectively (a factor of 1.48). This is primarily because KPP needs less computation time per numerical time step. While KPP does 1358 time steps in total, Cminor does 1394. Reading and preparing the simulation, on the other hand, takes 5 h and 52 min for KPP, and 2.2 s for Cminor. While this showcases that Cminor enables fast and flexible modification and analysis of mechanisms, sensitivites, and different scenarios, an optimization effort should be made to reach closer to KPP’s simulation times.

For the parcel model, we chose to reproduce the case used in the intercomparison study by Kreidenweis et al. (2003) and in Jaruga and Pawlowska (2018). The sulfur oxidation mechanism deployed there describes the irreversible accumulation of sulfur in cloud droplets via oxidation pathways with hydrogen peroxide and ozone. This models the aging of aerosol, in this case the growth, by existence of the aqueous phase, i.e., a cloud. Results are shown in Fig. 7, where yellow dots indicate the results presented in Jaruga and Pawlowska (2018). Most distinctly, Cminor predicts a cloud droplet concentration of 397 cm⁻³, while Jaruga and Pawlowska (2018) estimated a concentration of 269 cm⁻³. This difference is due to the lower water accommodation coefficient of 0.0415 used by Cminor, while Jaruga and Pawlowska (2018) used a value of 1.0. Using a water accommodation coefficient of 1.0, Cminor predicts a cloud droplet number concentration of 338 cm⁻³. All values are well within the range presented in the intercomparison study Kreidenweis et al. (2003), the results of Jaruga and Pawlowska (2018) being at the lower end, while Cminor is located more in the upper region.

Figure 7 Simulated conditions in the adiabatic parcel setup (Kreidenweis et al., 2003). Subplot (a) shows the evolution of the aerosol size distribution with the number concentration of activated aerosol, denoted by N, (b) the liquid water mixing ratio, (c) the SO₂ concentration (both gaseous and dissolved), and (d) the negative logarithm of the volumetric mean H⁺ concentration in the droplets. Yellow dots indicate the results presented in Jaruga and Pawlowska (2018).

Finally, Fig. 8 shows a collection of statistics of the mechanisms presented here. It gives a glimpse on how the computational effort and simulation times scale with the number of species and reactions, how the sparsity grows in the respective matrices, and of the robustness of the symmetric Markowitz ordering strategy (Markowitz, 1957), used to prevent sparse matrices from becoming dense matrices upon decomposition. This can only be seen as an empirical demonstration since the computational effort depends heavily on aspects of the chemical mechanism and cannot be predicted in terms of species number or other simple measures. The computation time for one integration time step shows linear behavior on the log-log scale, but this cannot be used as a general estimate as there are more factors determining computation times. For total simulation time, especially the number of steps needed is important. The total simulation times can be seen in Figs. 4 and 5.

Figure 8 Overview of some statistics of the presented mechanisms. The abscissa is the number of species, determining where the mechanism statistics are placed, to preserve any scaling behavior with the number of species (which is the number of equations in the ODE system). The ordinates do not have units, unless indicated in the legends. The bars indicate n_R, the number of reactions, the number of non-zeros in the iteration matrix I−hγJ, and the number of non-zero entries in the LU decomposition of the iteration matrix. The lines with circles indicate the time to prepare the mechanism, i.e., everything what has to be done before integration, namely, reading reactions and symbolic decomposition. Lines with squares show the time spent for one Rosenbrock step. Times for atmospheric and combustion mechanisms are shown separately, because the temperature equation for combustion increases the computational effort, as clearly visible in the figure. The grey line with triangles shows the decreasing density of the decomposed iteration matrix, i.e. the number of non-zero entries divided by the matrix size $n_{\mathrm{S}}^{\mathrm{2}}$.

6 Summary and Outlook

In this work, we presented Cminor – a highly efficient stand-alone solver environment designed to simulate chemical kinetic systems, from skeletal mechanisms to detailed descriptions of the atmosphere's chemical processes in the gaseous and aqueous phase. The optimized implementation allows for fast and accurate simulations of these chemical systems. Additionally, the lifting of an adiabatic parcel and cloud droplet formation can be simulated, which are influenced by chemical processes and vice versa. Furthermore, Cminor can solve the energy conservation equation needed for a perfectly adiabatic constant-volume reactor to investigate combustion chemistry instead of atmospheric systems.

The solver of Cminor uses linear-implicit Rosenbrock numerical schemes, which are well-suited for the integration of stiff systems of chemical kinetics. Moreover, Cminor’s solver is implemented in a way that is tailored to specifically consider and exploit the characteristics of ODE systems of chemical kinetics. The following list of key features illustrates the improvements of the solver. Using sparse linear algebra is essential, but also exploiting matrix block structures, the sparsity of the analytical Jacobian matrix, only considering relevant derivatives which corresponds to an approximated Jacobian matrix, efficient evaluation of the rate constants, a symbolic and droplet-wise vectorized algebra and LU decomposition, and automated adaptive step size with error control.

The input syntax of atmospheric mechanisms for Cminor is outlined in the text, and elaborated on in the Supplement. A collection of rate constants can be chosen of to simulate Arrhenius-type laws, photolytic reactions, third-body interactions, pressure-dependencies, humidity-dependencies, gas-aqua phase transfer, dissociation, pH-independence, and special reaction types, including custom formulas. The setup of a simulation is designed to be intuitive and only requires use of readable ASCII files to list the chemical reactions, initial values, and simulation parameters such as temperature, output frequency, and error tolerances. For combustion systems, the system description follows the widely used ChemKin syntax rules and JANAF standards.

Flexibility is achieved by the stand-alone nature of Cminor, which reads the mechanism and other inputs to construct the required ODE system on the fly with a high-speed mechanism parser. Chemical mechanisms, ambient conditions, and all simulation parameters can be changed and studied without recompiling, by typing a new parameter in the respective text file and re-executing the Cminor run command. Output is written in NetCDF format, or binary files for more abundant data.

The code is fully modularized and does not require the use of additional 3rd party software, except the output data format (NetCDF) and the basic Linear Algebra Pack (LAPack). Cminor v1.0 is open source and available under the terms of the GNU General Public License version 3.0.

We showed that Cminor is applicable to smaller as well as modern detailed descriptions of chemical processes in atmospheric and combustion chemistry. State-of-the-art multiphase mechanisms, e.g., the Master Chemical Mechanism version 3.2 (Jenkin et al., 2012) coupled with the Chemical Aqueous Phase RAdical Mechanism version 4.0α (Herrmann et al., 2005), with more than 10 000 species and 23 000 reactions, or a methyl decanoate mechanism (Herbinet et al., 2008), derived from bio-fuel chemistry, with nearly 3000 species and 17 000 reactions can be simulated. Also, Cminor is able to accurately reproduce previously studied effects of chemical processes on aerosol size distributions (Kreidenweis et al., 2003). We found that it takes 3.08 s to simulate the Master Chemical Mechanism version 3.2 with a state-of-the-art KPP code, while Cminor required 4.55 s, which is slower by a factor of 1.48. Reading the mechanism and preparing the simulation, however, took 5 h and 52 min for KPP, and 2.2 s for Cminor. This overhead is required any time a reaction, a parameter, or an environmental variable is changed, demonstrating the flexibility and handiness of Cminor. Nonetheless, further optimization is desirable to reach closer to KPP’s simulation time, for example in Cminor's linear algebra routines.

Cminor's current box-model setup can be used to study sensitivites, new mechanisms, specific interactions or other detailed research interests. The workflow of Cminor also enables ensemble simulations and analysis. Moreover, Cminor can easily be coupled with three-dimensional models, enabling the integration with comprehensive modeling frameworks, such as chemical transport or physical processes of the atmosphere. The droplet condensation equation resembles an important interaction of chemical and physical processes, which is accurately modeled by Cminor.

In future work, we aim to couple Cminor to a large-eddy simulation (LES) model with Lagrangian cloud microphysics. The droplet classes of Cminor match the concept of superdroplets in Lagrangian cloud models, while the LES simulates the transport of chemical constituents of the atmosphere. This will enable the investigation of the interplay of atmospheric chemistry with dynamics and cloud physics on larger scales. Due to its modular structure, Cminor is easily extendable. Resolving processes including a particulate phase would be a valuable addition, as well as ionic strength effects and other peculiarities for small water contents. Adding the possibility to include particle nucleation will be subject of future work. In this way, Cminor will contribute to a more systematic understanding of the influences of atmospheric chemistry on physical, meteorological, and climatological processes.