A description of the first open-source community release of MISTRA-v9.0: a 0D/1D atmospheric boundary layer chemistry model

We present MISTRA-v9.0, a one-dimensional (1D) and box (0D) atmospheric chemistry model. The model includes a detailed particle description with regards to the microphysics, gas–particle interactions, and liquid-phase chemistry within particles. Version 9.0 (v9.0) is the first release of MISTRA as an open-source community model. A major review of the code has been performed along with this public version release to improve the user friendliness and platform independence of the model. The purpose of this public release is to maximise the benefit of MISTRA to the community by making the model freely available and easier to use and develop. This paper presents a thorough description of the model characteristics and components. We show some examples of simulations reproducing previous studies with MISTRA, finding that v9.0 is consistent with previous versions.

1 Introduction

1.1 Scientific context and purpose of the model

Atmospheric aerosols are a major component of the Earth's climate system. They significantly affect the radiative balance of the atmosphere, through direct (scattering and absorption) and indirect effects (cloud properties modification) (Carslaw et al., 2010; Boucher et al., 2014; Bellouin et al., 2020). However, their concentrations, chemical, and physical properties are still insufficiently constrained, and the variability associated with their effects is dominant in the uncertainties of climate projections (Bender, 2020). Atmospheric particles also have a fundamental role in the chemistry of the atmosphere, since they offer a large surface area and volume for (photo)-chemical reactions to occur (Andreae and Crutzen, 1997; Finlayson-Pitts, 2009; George et al., 2015; Simpson et al., 2015; Seinfeld and Pandis, 2016; Kanakidou et al., 2018). Other impacts include the reduction of visibility (see Seinfeld and Pandis, 2016, chap. 15; Zhang et al., 2020, and references therein) and health effects of pollution (e.g. Pöschl, 2005; Molina et al., 2020, and references therein).

Numerical models are essential tools to help understand the relevant processes, and make projections of their evolution in a changing climate (Ervens, 2015). While global 3D models, and specifically Earth system models (ESMs), are well suited for climate simulation, the high computing cost of coupled physical–microphysical–chemical process modelling limits the space resolution of such models. Currently, a kilometre-scale resolution is already very challenging, thus preventing fully resolved approaches for subgrid-scale processes, such as turbulence. Conversely, limited-area and 1D models can reach sufficiently fine resolution for process-resolving simulations. Ultimately, box (0D) models are designed to focus only on single grid cell processes, further reducing the computing cost as compared to 1D models. In turn, the results obtained with such models can be used to develop parameterisations for use in 3D models.

Whatever the model, a crucial step is the validation based on field measurements. Balloon or flight surveys provided valuable data for this purpose, but the number of investigated parameters is necessarily limited, with many uncertainties and unknown values. Another useful tool for understanding atmospheric chemistry and physics are the atmospheric simulation chambers (see https://www.eurochamp.org, last access: 1 July 2021). These platforms enable the simultaneous measurement of a large number of chemical species and associated physical characteristics, in a constrained volume. The resulting data sets are highly valuable to validate models. In turn, numerical models are complementary tools to help understand and interpret measured results.

In this paper, we present the 1D boundary layer chemistry model MISTRA-v9.0, including size-resolved aerosol processes as well as particle–chemistry interaction. MISTRA-v9.0 also includes a box-model (0D) configuration, which can be adapted for atmospheric simulation chamber applications. In Sect. 1.2 and 1.3, we give a brief history of the MISTRA model followed by an overview on the recent developments presented in this paper. Section 2 gives a thorough description of processes implemented in the model. Section 3 presents practical and technical aspects of MISTRA-v9.0 with the main settings, while a set of example simulations reproducing settings and configurations of previous studies is presented in Sect. 4, to show the consistency of MISTRA-v9.0 with previous results.

1.2 History of MISTRA and reference publications

The MISTRA model was originally designed to study the MIcrophysics in STRAtus clouds, and was written based on a fog model (MIFOG: Bott et al., 1990; Bott and Carmichael, 1993; von Glasow and Bott, 1999). Bott et al. (1996) developed the first version of MISTRA for the simulation of cloud microphysics in the marine boundary layer (MBL). The unique feature of this model is the use of a 2D particle distribution, with the first dimension accounting for dry particle radius, and the second dimension for the total particle radius. Based on this first version of MISTRA, Bott (1997) further included typical particle distributions of urban and rural aerosols for the study of MBLs influenced by continental air masses, and assessed the radiative forcing of stratiform clouds. The radiation code used in MISTRA, called PIFM1 (Practical Improved Flux Method, developed by Zdunkowski et al., 1982), was updated by Loughlin et al. (1997) and the new radiation code, PIFM2, was evaluated. The collision–coalescence process was implemented in MISTRA by Bott (2000, 2001). Bott (1999a) adapted the chemistry module from Bott and Carmichael (1993) in MISTRA, with special emphasis on sulfur chemistry, and studied the retroaction of cloud processing over the microphysics in Bott (1999b). Meanwhile, von Glasow (2000) developed another chemistry module for MISTRA, with special emphasis on halogen chemistry, and presented the results for cloud-free (von Glasow et al., 2002a) and cloudy MBLs (von Glasow et al., 2002b). Our paper develops the branch of MISTRA based on von Glasow (2000), whose development and application until 2015 took place under the lead of Roland von Glasow.

From the early 2000s to the mid-2010s, MISTRA was regularly improved with respect to the chemistry and the gas–particle interactions. It was used in several studies, many of them with a focus on tropospheric halogen chemistry. MISTRA was used to investigate the influence of organic coating at the surface of sea salt particles over boundary layer chemistry, and especially on bromine and chlorine chemistry in the aqueous phase (Smoydzin and von Glasow, 2007). A major development was the introduction of a module for aerosol nucleation which significantly improved the iodine chemistry (Pechtl et al., 2006, 2007). The gas-phase chemical mechanism was updated by Sommariva and von Glasow (2012).

Over the years, numerous modelling studies were performed using MISTRA (von Glasow and Crutzen, 2004; Pechtl and von Glasow, 2007; Lawler et al., 2009; Jones et al., 2010; Joyce et al., 2014) including alternative model configurations where the chemistry was computed in a 0D atmospheric box-model mode (Buys et al., 2013), and a 0D chamber mode (Buxmann et al., 2015). MISTRA was adapted to model specific environments such as volcanic plumes (Aiuppa et al., 2007; Bobrowski et al., 2007, 2015) and polar conditions (Piot and von Glasow, 2008, 2009; Buys et al., 2013). MISTRA was also used to simulate the boundary layer chemistry over the Dead Sea after implementing a calculation of chemistry in this specific liquid medium and an explicit calculation of sea–air gas exchanges (Smoydzin and von Glasow, 2009). A module for firn chemistry was developed and coupled to MISTRA to specifically address the influence of chemical reactions occurring in the snowpack on the oxidative capacity of the atmosphere in snow-covered regions (Thomas et al., 2011, 2012). In this study, we present a selection of a few specific model settings reproducing previous work, to compare the original results with those obtained using MISTRA-v9.0.

1.3 Recent developments and public release

The previous (non-public) version of MISTRA (v7.4.1) included the update of the gas-phase chemical mechanism by Sommariva and von Glasow (2012). Version 8, featuring an alternative chemical bin definition, was partly developed but not completed, thus explaining the current version number. More information about past versions of MISTRA can be found in the preface of the manual (https://github.com/Mistra-UEA/Mistra/blob/master/doc/manual_v9.0.pdf, last access: 8 July 2022). Since 2015, significant efforts have been devoted to release MISTRA as an open-source model, including major technical improvements. The original code, written in Fortran77, has been updated to Fortran90 to ease future maintenance and developments. To improve robustness and portability of the code, intensive controls throughout the code have been performed to track issues, fix bugs, and conform to strict coding rules (Metcalf et al., 2004) and coding standards (see http://www.umr-cnrm.fr/gmapdoc/IMG/pdf/coding-rules.pdf, last access: 26 October 2021 and https://wiki.c2sm.ethz.ch/pub/CM/CodingAdvicesCosmo2/europ_sd.pdf​​​​​​​, last access: 9 July 2022). This was achieved with the help of the Fortran analyser Forcheck (v14.6, no longer distributed), as well as standard code check options of compilers. Computing efficiency has also been improved by factorising parts of code, and re-indexing arrays to respect column-major order in Fortran (i.e. innermost do-loops should be leftmost indexes). The chemical “Kinetic PreProcessor” (KPP; Damian et al., 2002; Sandu and Sander, 2006) has been updated to the latest version 2.2.3 (https://people.cs.vt.edu/~asandu/Software/Kpp/ last access: 23 June 2021) with minor tuning for use in MISTRA (see the “Code availability” section at the end of the paper). Overall, several technical developments have been implemented to make the model as user-friendly as possible, and easier to adopt. The model code of MISTRA-v9.0 now has improved readability, documentation, and is available under licence EUPL-v1.1 on https://github.com/Mistra-UEA (last access: 6 April 2022​​​​​​​).

2 Scientific description

2.1 Overview of the model components

MISTRA is a 1D model of the MBL. The vertical grid is separated into three regions: the lowest part is made of 100 layers with a constant thickness of 10 m, followed by 50 layers with logarithmically equidistant layers up to 2000 m height. The third region is a constant atmosphere with characteristics based on the standard atmosphere. It extends up to 50 km height and is only used for radiation calculations. These vertical grid settings (number and thickness of layers) can be configured easily if required.

Figure 1 schematically shows the most important processes that are included in the model for a cloudy MBL. The meteorological and microphysical part consists of the boundary layer model, MISTRA, described in detail by Bott et al. (1996) and Bott (1997). The most important processes are turbulent mixing, condensation, evaporation, and radiative heating. Apart from dynamics and thermodynamics, MISTRA includes a detailed microphysical module that calculates particle growth explicitly and includes feedbacks between radiation and particles. The radiative-transfer parameterisation is a standard two-stream code using 6 spectral bands for visible and 12 bands for infrared radiation (Loughlin et al., 1997). A chemistry module computes the atmospheric chemistry in the gas phase and in the particles. Gas-phase chemistry is active in all model layers; aerosol chemistry is only in layers where the relative humidity has been greater than the deliquescence humidity and not dropped below the crystallisation humidity (as discussed in Sect. 2.3.3). When a cloud forms, cloud droplet chemistry is also active. Fluxes of sea salt aerosol and gases from the ocean are included (see Sect. 2.3.6). A nucleation module is also included to account for new particles nucleated from the gas-phase species (see Sect. 2.3.7).

Figure 1 Schematic depiction of the most important processes included in the 1D boundary layer model MISTRA. The free troposphere and marine boundary layer are denoted as FT and MBL, respectively.

2.2 Meteorology, microphysics and thermodynamics

The model is 1D, thus all variables are taken to be horizontally homogeneous. The set of prognostic variables comprises the horizontal components of the wind speed u and v, the specific humidity q, and the potential temperature θ. The Boussinesq approximation is applied and the air pressure is derived from the large-scale hydrostatic equilibrium.

The set of governing equations for these prognostic variables is

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial u}{\partial t}}=-w{\displaystyle \frac{\partial u}{\partial z}}+{\displaystyle \frac{\partial }{\partial z}}\left(K_{\mathrm{m}}{\displaystyle \frac{\partial u}{\partial z}}\right)+f_{\mathrm{c}}(v-v_{\mathrm{g}})\text{(2)}& {\displaystyle \frac{\partial v}{\partial t}}=-w{\displaystyle \frac{\partial v}{\partial z}}+{\displaystyle \frac{\partial }{\partial z}}\left(K_{\mathrm{m}}{\displaystyle \frac{\partial v}{\partial z}}\right)-f_{\mathrm{c}}(u-u_{\mathrm{g}})\text{(3)}& {\displaystyle \frac{\partial q}{\partial t}}=-w{\displaystyle \frac{\partial q}{\partial z}}+{\displaystyle \frac{\partial }{\partial z}}\left(K_{\mathrm{h}}{\displaystyle \frac{\partial q}{\partial z}}\right)+{\displaystyle \frac{C}{\mathit{\rho }}}\text{(4)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial \mathit{\theta }}{\partial t}}& =-w{\displaystyle \frac{\partial \mathit{\theta }}{\partial z}}+{\displaystyle \frac{\partial }{\partial z}}\left(K_{\mathrm{h}}{\displaystyle \frac{\partial \mathit{\theta }}{\partial z}}\right)-{\left({\displaystyle \frac{p_{\mathrm{0}}}{p}}\right)}^{\frac{R_{\mathrm{a}}}{c_{\mathrm{p}}}}{\displaystyle \frac{\mathrm{1}}{c_{\mathrm{p}}\mathit{\rho }}}\\ & \times \left({\displaystyle \frac{\partial E_{\mathrm{n}}}{\partial z}}+LC\right),\end{array}\end{array}$$

where f_c is the Coriolis parameter, u_g and v_g are the geostrophic wind components, K_m and K_h are the turbulent exchange coefficients for momentum and heat, L is the latent heat of condensation, C the condensation rate, ρ the air density, p the air pressure, p₀ the air pressure at the surface, R_a the specific gas constant for dry air, c_p the specific heat of dry air at constant pressure, and E_n the net radiative flux density, respectively. The first term on the right of each equation is the large-scale subsidence. Strictly, in a 1D framework, the vertical velocity w should be zero everywhere, otherwise this implies a downward mass transport (for w<0) without lateral outflow at the bottom of the 1D model column, as would occur in the real atmosphere. Therefore the mass balance is violated if subsidence is included. However, including subsidence is essential for modelling stratiform cloud evolution (e.g. Driedonks and Duynkerke, 1989). In runs where only aerosol chemistry is studied, i.e. in runs without clouds, the vertical velocity is set to zero (w=0) in the model to avoid this problem, while for the cloud runs, subsidence is usually included.

Turbulence is treated with the level 2.5 model of Mellor and Yamada (1982) with the modifications described in Bott et al. (1996) and Bott (1997). The turbulent exchange coefficients, K_m and K_h, are calculated via stability functions $S_{\mathrm{m}/\mathrm{h}}$ and $G_{\mathrm{m}/\mathrm{h}}$, where the subscript “m” stands for shear and “h” for buoyancy production. The prognostic equation for the turbulence kinetic energy (TKE) e is

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\displaystyle \frac{\partial e}{\partial t}}& =-w{\displaystyle \frac{\partial e}{\partial z}}+{\displaystyle \frac{\partial }{\partial z}}\left(K_{\mathrm{e}}{\displaystyle \frac{\partial e}{\partial z}}\right)+{\displaystyle \frac{(\mathrm{2}e{)}^{\mathrm{3}/\mathrm{2}}}{l}}\\ & \times (S_{\mathrm{m}}{G}_{\mathrm{m}}+S_{\mathrm{h}}{G}_{\mathrm{h}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{16.6}}}),\end{array}\end{array}$$

assuming a constant dissipation ratio (last term on the right). For more details and an explanation of the calculation of the mixing length l, the exchange coefficient K_e for TKE, and the functions $S_{\mathrm{m}/\mathrm{h}}$ and $G_{\mathrm{m}/\mathrm{h}}$, see Mellor and Yamada (1982), Bott et al. (1996), and Bott (1997).

Figure 2 The 2D particle spectrum as a function of the dry particle radius a and the total particle radius r. Added are the chemical bins: sulfate aerosol bin (I), sea salt aerosol bin (II), sulfate cloud droplet bin (III), and sea salt droplet bin (IV). For simplicity, a 35 by 35 bin grid is plotted, in the model 70 × 70 bins are used.

The microphysics is treated using a joint 2D particle size distribution function f(a,r), where a is the dry particle radius the particles would have if no water were present in the particles, and r is the total particle radius. The 2D particle grid is divided into 70 logarithmically equidistant spaced dry aerosol classes. The minimum dry particle radius is generally set to 0.005 µm and the maximum radius 15 µm. Choosing these values allows one to account for all accumulation mode particles and most of the coarse particles. The minimum and maximum, as well as the number of bins for both dimensions of the particle spectrum are adjustable. Each of the 70 dry particle classes is associated with 70 total particle radius classes, ranging from the actual dry particle radius up to 60 µm (150 µm in cloud runs). See Fig. 2 for a depiction of 2D particle grid. The prognostic equation for f(a,r) is

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\displaystyle \frac{\partial f(a,r)}{\partial t}}& =-w{\displaystyle \frac{\partial f(a,r)}{\partial z}}+{\displaystyle \frac{\partial }{\partial z}}\left(K_{\mathrm{h}}\mathit{\rho }{\displaystyle \frac{\partial f(a,r)/\mathit{\rho }}{\partial z}}\right)\\ & -{\displaystyle \frac{\partial }{\partial z}}\left(w_{\mathrm{t}}f\right(a,r\left)\right)-{\displaystyle \frac{\partial }{\partial r}}\left(\dot{r}f\right(a,r\left)\right).\end{array}\end{array}$$

Again, subsidence is the first term on the right, followed by turbulent mixing, particle sedimentation (w_t is the sedimentation velocity), and changes in f due to particle growth ($\dot{r}=\mathrm{d}r/\mathrm{d}t$). The 2D particle spectrum is initialised with distribution depending on the type of aerosol chosen (see Bott, 1997, and references therein). Currently, particle distributions are provided for typical marine, rural, and urban air masses. Other distributions are available for specific studies, such as a polar distribution (see Buys et al., 2013, and the corresponding example simulation). Particles are initialised with a water coating according to the equilibrium radius of the dry nucleus at the ambient relative humidity. During the time integration, particle growth is calculated explicitly for each bin of the 2D particle spectrum using the growth equation following Davies (1985) (see also Bott et al., 1996):

$$\begin{array}{}\text{(7)}& r{\displaystyle \frac{\mathrm{d}r}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{1}}{C_{\mathrm{1}}}}\left[C_{\mathrm{2}}\left({\displaystyle \frac{S_{\mathrm{\infty }}}{S_{\mathrm{r}}}}-\mathrm{1}\right)-{\displaystyle \frac{F_{\mathrm{d}}(a,r)-m_{\mathrm{w}}(a,r)c_{\mathrm{w}}\mathrm{d}T/\mathrm{d}t}{\mathrm{4}\mathit{\pi }r}}\right],\end{array}$$

where m_w(a,r) is the liquid water mass of the particle, c_w is the specific heat of water, S_∞ is the ambient supersaturation and S_r is the supersaturation at the droplet's surface according to the Köhler equation:

$$\begin{array}{}\text{(8)}& S_{\mathrm{r}}=\exp \left[{\displaystyle \frac{A}{r}}-{\displaystyle \frac{B{a}^{\mathrm{3}}}{r^{\mathrm{3}}-a^{\mathrm{3}}}}\right],\end{array}$$

where factors A and B account for the Kelvin effect and the solute effect, respectively.

The change in particle radius is not determined by changes in water vapour saturation alone, but also by the net radiative flux at the particle's surface F_d(a,r), that leads to temperature changes and therefore to condensation or evaporation. The constants, C₁ and C₂, in Eq. (7) are

$$\begin{array}{}\text{(9)}& C_{\mathrm{1}}={\mathit{\rho }}_{\mathrm{w}}L+{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}{C}_{\mathrm{2}}}{D_{\mathrm{v}}^{\prime }{S}_{\mathrm{r}}{\mathit{\rho }}_{\mathrm{s}}}}{C}_{\mathrm{2}}=k^{\prime }T{\left[{\displaystyle \frac{L}{R_{\mathrm{v}}T}}-\mathrm{1}\right]}^{-\mathrm{1}},\end{array}$$

where ρ_w is the specific density of water, ρ_s is the saturation vapour density and R_v the specific gas constant for water vapour. The thermal conductivity k^′ of moist air and the diffusivity of water vapour $D_{\mathrm{v}}^{\prime }$ have been corrected for gas kinetic effects following Pruppacher and Klett (1997) (their Eqs. 13.20 and 13.14, respectively). For the accommodation coefficient of water (condensation coefficient), a value of α_c=0.036 is used (see Table 5.4 in Pruppacher and Klett, 1997 for a compilation of measured α_c values; in Table 13.1 they use α_c=0.036 as “best estimate”).

The condensation rate C in Eq. (4) is determined diagnostically from the particle growth Eq. (7).

Collision–coalescence processes are not included in the model because this leads to difficulties when redistributing the chemical species in the particles. A version of MISTRA including collision–coalescence without considering chemistry does exist (Bott, 2000), and this limitation of MISTRA-v9.0 is discussed in Sect. 4.1.2.

For the calculation of the radiative fluxes, a δ-two-stream approach is used (PIFM radiative code: Zdunkowski et al., 1982; Loughlin et al., 1997). The radiative fluxes are used for calculating heating rates and the effect of radiation on particle growth. The radiation field is calculated with the aerosol/cloud particle data from the microphysical part of the model, so feedbacks between radiation and particle growth are fully implemented. The calculation of photolysis frequencies is described in Sect. 2.3.5.

2.3 Chemistry

The multiphase chemistry module comprises chemical reactions in the gas phase as well as in deliquescent aerosol and cloud particles. Transfer between gas and aqueous phase and surface reactions on particles are also included. The reaction set was based on that of Sander and Crutzen (1996) plus some organic reactions from Lurmann et al. (1986). It has been updated and expanded by von Glasow and Crutzen (2004) to include a better description of the oxidation of dimethylsulfide (DMS). Iodine chemistry was significantly improved by Pechtl et al. (2006, 2007). Further updates to the chemical mechanism were done by Sommariva and von Glasow (2012). The current mechanism is provided in the model manual (tables in Appendix D). In the following, the term aqueous phase is used as generic term for sub-cloud aerosol, interstitial aerosol (i.e. non-activated aerosol particles in cloudy layers), and cloud particles. Aqueous chemistry is not computed above the top of the boundary layer (i.e. the top of clouds, if present).

2.3.1 Gas phase and uptake

The prognostic equation for the concentration of a gas-phase chemical species c_g (amount per air volume) including subsidence, turbulent exchange, deposition on the ocean surface, chemical production and destruction, emission, and exchange with the aqueous phases is

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\displaystyle \frac{\partial c_{\mathrm{g}}}{\partial t}}& =-w{\displaystyle \frac{\partial c_{\mathrm{g}}}{\partial z}}+{\displaystyle \frac{\partial }{\partial z}}\left(K_{\mathrm{h}}\mathit{\rho }{\displaystyle \frac{\partial c_{\mathrm{g}}/\mathit{\rho }}{\partial z}}\right)+P-S{c}_{\mathrm{g}}+E\\ & -D{c}_{\mathrm{g}}-\sum _{i=\mathrm{1}}^{n_{\mathrm{kc}}}\left[\stackrel{\mathrm{‾}}{k_{\mathrm{t},i}}\right(w_{\mathrm{l},i}{c}_{\mathrm{g}}-{\displaystyle \frac{c_{\mathrm{a},i}}{H_{\mathrm{s}}^{\mathrm{cc}}}}\left)\right].\end{array}\end{array}$$

Again subsidence is the first term on the right and is included only in runs with clouds, otherwise w=0. The second term on the right-hand side of Eq. (10) describes the vertical turbulent mixing. Chemical production and sink (i.e. loss) terms are denoted as P and S, respectively. The emission E as well as dry deposition D are effective only in the lowermost model layer. The calculation of the dry deposition velocity $v_{\mathrm{g}}^{\mathrm{dry}}$, that is needed for the determination of D, is explained in Sect. 2.3.6. Note that both E and D are not inserted as fluxes in Eq. (10). Instead, the actual fluxes have to be divided by the thickness of the lowermost model layer to yield D and E. The last term in Eq. (10) describes the transport between the gas phase and the aqueous phases according to the formulation by Schwartz (1986) (see also Sander, 1999). In this term, n_kc is the number of aqueous classes (see Sect. 2.3.2), $H_{\mathrm{s}}^{\mathrm{cc}}$ is the dimensionless Henry constant obtained by $H_{\mathrm{s}}^{\mathrm{cc}}=H_{\mathrm{s}}^{\mathrm{cp}}RT$, where $H_{\mathrm{s}}^{\mathrm{cp}}$ is in $\mathrm{mol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{Pa}}^{-\mathrm{1}}$, and w_l,i is the dimensionless liquid water content (V_aq $/$ V_air​​​​​​​) of bin i.

For a single particle, the mass transfer coefficient k_t is defined as

$$\begin{array}{}\text{(11)}& k_{\mathrm{t}}=({\displaystyle \frac{r^{\mathrm{2}}}{\mathrm{3}{D}_{\mathrm{g}}}}+{\displaystyle \frac{\mathrm{4}r}{\mathrm{3}\stackrel{\mathrm{‾}}{v}\mathit{\alpha }}}{)}^{-\mathrm{1}},\end{array}$$

with the particle radius r, the mean molecular speed $\stackrel{\mathrm{‾}}{v}=\sqrt{\mathrm{8}RT/\left(M\mathit{\pi }\right)}$ (M is the molar mass), the accommodation coefficient α, and the gas phase diffusion coefficient D_g. D_g is approximated using the mean free path length λ as $D_{\mathrm{g}}=\mathit{\lambda }\stackrel{\mathrm{‾}}{v}/\mathrm{3}$ (e.g. Gombosi, 1994, p. 125).

Chameides (1984) points out that the time needed to establish equilibrium between the gas and aqueous phases differs greatly for individual species and that soluble species never reach equilibrium in cloud droplets, emphasising the importance of describing phase transfer in the kinetic form that is used here. Audiffren et al. (1998) and Chaumerliac et al. (2000) point out that for reactive species like H₂O₂, the use of the Henry equilibrium assumption, instead of the detailed description of mass transfer in the kinetic form that is used here, would lead to significant errors in cloud-droplet concentrations.

Ambient particle populations are never monodisperse, i.e. one has to account for particles with different radii. The transfer coefficient $\stackrel{\mathrm{‾}}{k_{\mathrm{t}}}$ for a particle population is given by the integral

$$\begin{array}{}\text{(12)}& \stackrel{\mathrm{‾}}{k_{\mathrm{t}}}={\displaystyle \frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}{w}_l}}\underset{\lg r_{\min }}{\overset{\lg r_{\max }}{\int }}({\displaystyle \frac{r^{\mathrm{2}}}{\mathrm{3}{D}_{\mathrm{g}}}}+{\displaystyle \frac{\mathrm{4}r}{\mathrm{3}\stackrel{\mathrm{‾}}{v}\mathit{\alpha }}}{)}^{-\mathrm{1}}{r}^{\mathrm{3}}{\displaystyle \frac{\partial N}{\partial \lg r}}\mathrm{d}\lg r,\end{array}$$

where the size distribution function $\partial N/\partial \lg r$ depends on the type of aerosol chosen.

2.3.2 Aqueous phase

Aqueous chemistry is calculated in four bins (see Fig. 2): deliquescent particles with a dry radius less than 0.5 µm are included in the “sulfate aerosol” bin #1, whereas deliquescent particles with a dry radius greater than 0.5 µm are in the “sea salt aerosol” bin #2. Although the composition of the particles changes over time, the terms “sulfate” and “sea salt” aerosol are used to describe the origin of the particles. The particles are internally mixed by exchange with the gas phase but, as mentioned earlier, not by particle coagulation. Depending on the type of aerosol relevant to the study, various initial compositions of the aerosol bins may be chosen.

When the total particle radius exceeds the dry particle radius by a factor of 10, i.e. when the total particle volume is 1000 times greater than the dry aerosol volume, the particle and its associated chemical species are moved to the corresponding sea salt or sulfate-derived cloud particle class (#3 and #4, respectively). This threshold roughly coincides with the critical radius derived from the Köhler equation (see Eq. 8). Conversely, when particles shrink, they are redistributed from the droplet to the aerosol bins.

Therefore, in a cloud-free layer there are two (n_kc=2) aqueous chemistry bins (sulfate and sea salt aerosol), and in a cloudy layer two cloud droplet (sulfate and sea salt derived) and two interstitial aerosol (sulfate and sea salt) bins, giving a total of four (n_kc=4) aqueous chemistry bins. In each of these bins the following prognostic equation is solved for each chemical species c_a,i (amount per air volume), where the index i stands for the ith aqueous bin:

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{\displaystyle \frac{\partial c_{\mathrm{a},i}}{\partial t}}& =-w{\displaystyle \frac{\partial c_{\mathrm{a},i}}{\partial z}}+{\displaystyle \frac{\partial }{\partial z}}\left(K_{\mathrm{h}}\mathit{\rho }{\displaystyle \frac{\partial c_{\mathrm{a},i}/\mathit{\rho }}{\partial z}}\right)+P-S{c}_{\mathrm{a},i}\\ & +E-D+P_{\mathrm{pc}}+\stackrel{\mathrm{‾}}{k_{\mathrm{t},i}}(w_{\mathrm{l},i}{c}_{\mathrm{g}}-{\displaystyle \frac{c_{\mathrm{a},i}}{H_{\mathrm{s}}^{\mathrm{cc}}}}).\end{array}\end{array}$$

The individual terms have similar meanings as in Eq. (10). The calculation of the sedimentation velocity $v_{\mathrm{a},i}^{\mathrm{dry}}$, that is needed for the calculation of the dry deposition D, is explained in Sect. 2.3.6. The additional term P_pc accounts for the transport of chemical species from the aerosol to the cloud droplet regimes and vice versa when droplets are formed or when they evaporate, i.e. when particles move along the Köhler curve and get activated or unactivated. If only phase transfer is considered, Eq. (13) reduces in steady-state conditions ($\partial c_{\mathrm{a},i}/\partial t=\mathrm{0}$) to the Henry equilibrium $c_{\mathrm{a},i}=w_{\mathrm{l},i}{c}_{\mathrm{g}}{H}_{\mathrm{s}}^{\mathrm{cc}}$.

The concentration of H⁺ ions is calculated like any other species, i.e. no further assumptions are made. The charge balance is satisfied implicitly.

2.3.3 Hysteresis of particle activation

Cloud processing, i.e. the change of aerosol mass due to uptake of gases, is included based on the model of Bott (1999b).

It has been observed in many laboratory experiments that soluble aerosols remain in a highly concentrated metastable aqueous state when they are dried below their deliquescence humidity. Only when they reach the crystallisation humidity they can be regarded as “dry”. This effect is called the hysteresis effect. For NaCl, the crystallisation point is about 45 % relative humidity (Shaw and Rood, 1990; Tang, 1997; Pruppacher and Klett, 1997; and Lee and Hsu, 2000).

The crystallisation humidity for many mixed aerosol particles containing sulfate or nitrate is below 40 % relative humidity (Seinfeld and Pandis, 2016, and references therein), implying that aerosol particles that had already been involved in cloud cycles will also be in an aqueous metastable state. Therefore most soluble aerosol particles will be present in the atmosphere as metastable aqueous particles below their deliquescence humidity. If the humidity drops below the crystallisation humidity, these particles can only reactivate when the deliquescence humidity is reached.

2.3.4 Accounting for the chemical activity

Aerosol particles are usually highly concentrated solutions. Laboratory measurements show that NaCl molalities can be in excess of 10 mol kg⁻¹ (Tang, 1997) implying high ionic strengths. Therefore, it is necessary to account for deviations from ideal behaviour by including activity coefficients. The Pitzer formalism (Pitzer, 1991) is used to calculate the activity coefficients for the actual composition of each aqueous size bin. The implementation by Beiping Luo (Luo et al., 1995 and personal communication, 1999​​​​​​​) is used in MISTRA. It computes the activity of 7 main ions (H⁺, ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$, Na⁺, ${{\mathrm{HSO}}_{\mathrm{4}}}^{-}$, ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$, ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$, and Cl⁻). The activities of 15 other ions are scaled on the previous ones, based on the results from Liang and Jacobson (1999) and Chameides and Stelson (1992).

2.3.5 Photolysis

Here an overview of the calculation of the photolysis rates is given. For a detailed description see the model manual, chap. 5​​​​​​​. Photolysis is calculated online using the method of Landgraf and Crutzen (1998). The photolysis rate constant (or photodissociation coefficient) J_X for a gas X can be calculated from the spectral actinic flux F(λ) via the integral

$$\begin{array}{}\text{(14)}& J_{\mathrm{X}}=\underset{I}{\int }{\mathit{\sigma }}_{\mathrm{X}}\left(\mathit{\lambda }\right){\mathit{\varphi }}_{\mathrm{X}}\left(\mathit{\lambda }\right)F\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda },\end{array}$$

where λ is the wavelength, σ_X the absorption cross section, ϕ_X the quantum yield, and I the photochemically active spectral interval. If the integral in Eq. (14) was approximated with a sum, the number of wavelength intervals needed for an accurate approximation of the integral would be in the order of 100, which would lead to excessive computing times. Landgraf and Crutzen (1998) suggested a method using only 8 spectral intervals approximating Eq. (14) by

$$\begin{array}{}\text{(15)}& {\displaystyle }{J}_{\mathrm{X}}\approx \sum _{i=\mathrm{1}}^{\mathrm{8}}{J}_{i,\mathrm{X}}^{\mathrm{a}}{\mathit{\delta }}_i,\end{array}$$

where $J_{i,\mathrm{X}}^{\mathrm{a}}$ is the photolysis rate constant for a purely absorbing atmosphere. The factor δ_i

$$\begin{array}{}\text{(16)}& {\displaystyle }{\mathit{\delta }}_i={\displaystyle \frac{F\left({\mathit{\lambda }}_i\right)}{F^{\mathrm{a}}\left({\mathit{\lambda }}_i\right)}}\end{array}$$

describes the effect of scattering by air molecules, aerosol, and cloud particles. The actinic flux of a purely absorbing atmosphere is F^a(λ_i). The factor δ_i is calculated online for one wavelength for each interval, while the $J_{i,\mathrm{X}}^{\mathrm{a}}$ is precalculated with a fine spectral resolution and approximated during runtime from lookup tables or by using polynomials. The advantage of this procedure is that the fine absorption structures that are present in σ_X and ϕ_X are considered and only Rayleigh and cloud scattering, included in F(λ_i), are treated with a coarse spectral resolution, which is justified.

For the calculation of the actinic fluxes, a four-stream radiation code is used in addition to the two-stream radiation code used for the determination of the net radiative flux density E_n, because different spectral resolutions and accuracies are needed for these different purposes. Based on the findings of Ruggaber et al. (1997), photolysis rates inside aqueous particles are increased by a factor of two to account for the actinic flux enhancement inside the particles due to multiple scattering.

2.3.6 Emission and deposition

The emission of gases is accounted for in the model, either with constant emission fluxes (e.g. for DMS and NH₃ emitted from the sea surface), or with scenarios of emission variable in time (see the example run based on the study of Joyce et al., 2014).

Sea salt particles are emitted by bursting bubbles at the sea surface (e.g. Woodcock et al., 1953; Pruppacher and Klett, 1997). The parameterisations of Monahan et al. (1986) and Smith et al. (1993) are implemented to estimate the flux of particles. The former is advised for small to moderate wind speeds, while the latter has to be used for high wind speeds.

The dry deposition velocity for gases $v_{\mathrm{g}}^{\mathrm{dry}}$ at the sea surface is calculated using the resistance model described by Wesely (1989):

$$\begin{array}{}\text{(17)}& v_{\mathrm{g}}^{\mathrm{dry}}={\displaystyle \frac{\mathrm{1}}{r_{\mathrm{a}}+r_{\mathrm{b}}+r_{\mathrm{c}}}}.\end{array}$$

The aerodynamic resistance r_a is calculated using

$$\begin{array}{}\text{(18)}& r_{\mathrm{a}}={\displaystyle \frac{\mathrm{1}}{\mathit{\kappa }{u}_{\ast }}}\left[\ln \right({\displaystyle \frac{z}{z_{\mathrm{0}}}})+{\mathrm{\Phi }}_{\mathrm{s}}(z,L_{\mathrm{MO}}\left)\right],\end{array}$$

with the friction velocity u_∗, the von Kármán constant κ=0.4, and the stability function Φ_s which depends on the Monin–Obukhov length L_MO, the roughness length z₀ and a reference height z. The quasi-laminar layer resistance r_b is parameterised for gases as

$$\begin{array}{}\text{(19)}& r_{\mathrm{b}}={\displaystyle \frac{\mathrm{5}S{c}^{\mathrm{2}/\mathrm{3}}}{u_{\ast }}}.\end{array}$$

The Schmidt number Sc can be written as $Sc=\mathit{\nu }/D_{\mathrm{g}}$ with the kinematic viscosity of air ν and the gas diffusion coefficient D_g as in Eq. (11). The surface resistance r_c is calculated using the formula by Seinfeld and Pandis (2016) (their Eq. 19.30):

$$\begin{array}{}\text{(20)}& r_{\mathrm{c}}={\displaystyle \frac{\mathrm{2.54}\times {\mathrm{10}}^{\mathrm{4}}\mathrm{mol}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{Pa}}^{-\mathrm{1}}\mathrm{K}}{H^{\ast }T{u}_{\ast }}},\end{array}$$

with the effective Henry constant H^∗.

The dry deposition velocity of particles $v_{\mathrm{a},i}^{\mathrm{dry}}$ is calculated following Seinfeld and Pandis (2016):

$$\begin{array}{}\text{(21)}& v_{\mathrm{a},i}^{\mathrm{dry}}=\left\{\begin{array}{ll}\frac{\mathrm{1}}{r_{\mathrm{a}}+r_{\mathrm{b}}+r_{\mathrm{a}}{r}_{\mathrm{b}}{w}_{\mathrm{t}}}+w_{\mathrm{t}}& \text{lowest model layer}\\ w_{\mathrm{t}}& \text{rest of model domain.}\end{array}\right.,\end{array}$$

where the quasi-laminar resistance r_b is parameterised for particles as

$$\begin{array}{}\text{(22)}& r_{\mathrm{b}}={\displaystyle \frac{\mathrm{1}}{u_{\ast }(S{c}^{-\mathrm{2}/\mathrm{3}}+{\mathrm{10}}^{-\mathrm{3}/St})}}.\end{array}$$

The Stokes number St can be written as $St=w_{\mathrm{t}}{u}_{\ast }^{\mathrm{2}}/\left(g\mathit{\nu }\right)$ where g is the gravitational acceleration. The particle sedimentation velocity w_t is calculated in the microphysical module assuming Stokes flow, and taking into account the Cunningham slip-flow correction for particles with r<10 µm, and following Beard for larger particles (see Pruppacher and Klett, 1997).

Finally, the dry deposition term D is calculated as

$$\begin{array}{}\text{(23)}& D=\exp \left(-\mathrm{\Delta }t/h\times v_{\mathrm{g}}^{\mathrm{dry}}\right),\end{array}$$

where Δt is the model time step, and h is the height of the lowermost model layer.

2.3.7 Nucleation

A module computing the nucleation process was implemented in MISTRA by Pechtl et al. (2006). Only a brief overview is given here, while a comprehensive description is given in the model manual (chap. 4). The nucleation module developed by Pechtl et al. (2006) includes both ternary sulfuric acid–ammonia–water (H₂SO₄–NH₃–H₂O) nucleation, and homomolecular homogeneous OIO nucleation. The former is explicitly calculated as a function of H₂SO₄ and NH₃ concentrations, relative humidity, and temperature following the work by Napari et al. (2002). The latter is parameterised following Burkholder et al. (2004). Each process can be activated or not independently (see Table 1), and lead to the computation of “real” nucleation rates. In a second step, the “apparent” nucleation rate is computed following the work of Kerminen and Kulmala (2002) and Kerminen et al. (2004).

The nucleated particles computed in this module can then be integrated in the model, with three possible options: (i) no coupling, (ii) coupling with the microphysics without feedback on chemistry, and (iii) coupling with microphysics and chemistry (see Table 1).

3 Technical description

3.1 Namelist settings

3.1.1 General configuration switches

Table 1 presents the switches available to define the model configuration.

Table 1 General configuration switches in MISTRA-v9.0.

^* network Common Data Form; see https://www.unidata.ucar.edu/software/netcdf/ (last access: 21 November 2021).

3.1.2 Initialisation and run settings

Initial atmospheric conditions are set in the namelist with the parameters presented in Table 2. All these parameters have default values, even if most of them are expected to be redefined by the user to match the simulated atmosphere. Standard settings (timing and geography, run duration) are straightforward and are not detailed hereafter. Typical run duration covers a few hours to a few days. A longer run duration is sometimes necessary for model spin-up. The restart option of the model allows a single spin-up run to initialise the model and perform a sensitivity analysis from that stage, for instance. In addition to these, surface settings are detailed in Table 3.

Table 2 Namelist settings for model initialisation of MISTRA-v9.0.

^* 1 DU $=\mathrm{10}\mathrm{1325}/(\mathrm{273.15}\times \mathrm{1.380649}\times {\mathrm{10}}^{\mathrm{-}\mathrm{23}})$ m${}^{-\mathrm{2}}\approx \mathrm{2.687}\times {\mathrm{10}}^{\mathrm{25}}$ m⁻².

Monahan et al. (1986)Smith et al. (1993)

Table 3 Namelist settings for surface initialisation and parameterisation in MISTRA-v9.0.

^* Currently, constant albedo is used over the six solar wavelength bands. Alternative choice with varying albedo could be implemented with this namelist option.

3.1.3 Special runs setting

When a specific run requires multiple adjustments in various parts of the code that were not already including namelist options, a single general switch might be used instead of defining several new namelist entries for each parameterisation that requires special settings. An example of such a global switch for special configuration is given with the lpJoyce14bc switch, used to activate all relevant parts of code to reproduce the base case of the Joyce et al. (2014) study (particle distribution and composition, gas and particle emission scenario, special formulation of accommodation coefficient for N₂O₅, and of gas dry deposition, etc.). Similarly, switches lpBuys13_0D and lpBuxmann15alph are used to reproduce all relevant settings of the studies of Buys et al. (2013, 0D case) and Buxmann et al. (2015, alpha case), respectively.

3.2 How to run MISTRA

3.2.1 Get the code

The model is provided on GitHub on the following repository: https://github.com/Mistra-UEA/Mistra (last access: 9 June 2022​​​​​​​). It is released under the European Union Public Licence (EUPL) v1.1, which permits free commercial and private use and unrestricted distribution, but requires that future developments of MISTRA are shared under the same licence. The version of KPP adapted for MISTRA is provided along with the distribution, and is released under its own licence.

3.2.2 System requirements and installation

A Fortran compiler is required to compile the model code. During the recent development stages, MISTRA has been regularly compiled using either GNU Fortran (gfortran) or Intel Fortran (ifort). New users are advised to choose one of these compilers. The implementation of KPP output files into MISTRA is done with bash and csh scripts, thus any change in the chemical mechanism will require these shells.

Plotting scripts provided as example are written for Ferret (http://ferret.pmel.noaa.gov/Ferret/, last access: 4 November 2021) and NCL (NCAR, 2019), but neither are necessary to run the model. Only KPP needs to be installed on the user system and the instructions to do so are in the readme file of the KPP distribution package. Preprocessed files using the current chemical mechanism are provided in the distribution, so the installation of KPP can be skipped until the user needs to modify the chemical mechanism.

3.2.3 Prepare the chemical mechanism files

This section can be skipped if no change is applied to the current chemistry mechanism. All files related to the chemical mechanism are contained in the subdirectory ./src/mech. The chemical mechanism, written with the formalism of KPP, is contained in two main files: master_gas.eqn for gas phase reactions, and master_aqueous.eqn for the liquid phase mechanism. For convenience, all necessary steps to prepare the equation files for KPP, run KPP, and adapt the resulting output from KPP for MISTRA have been set up in a Makefile, so that the user simply has to run make in the ./src/mech directory to proceed. The resulting files are copied to the main source directory.

Figure 3 Contour plot of cloud water content (in g kg⁻¹) as a function of height and simulation time. Top: study from Loughlin et al. (1997, Fig. 1a). Bottom: MISTRA-v9.0. For both, the simulation settings are identical to those in Bott et al. (1996). Note the minimum contour level is set to 0.01 in both panels, but was displayed incorrectly in the original figure. Top panel reproduced with permission from Loughlin et al. (1997, Fig. 1a). © 1997 by John Wiley and Sons.

3.2.4 Compile the model

In ./src, after ensuring that the Makefile refers to the correct Fortran compiler, and links to the appropriate netCDF libraries, compile the model running make. The resulting executable file is mistra.

3.2.5 Set a namelist and initial chemical species concentration

As presented in details in Sect. 3.1, the namelist file allows the user to configure the model, by setting the main options and initialisation values. Several namelists are provided in the distribution and can be used as starting points to define new ones corresponding to the user requirements. The set of initial concentration, and emission of gas phase species can be set in a tab-separated table, whose name has to be specified in the namelist. If no file is specified, the ./src/mech/gas_species.csv file is used by default.

3.2.6 Set a param file and run

The param file allows the user to specify the namelist to use for the run, and to define the paths to the input, output and mechanism directories. For most cases, these directories will be the default ones (./input, ./output, and ./src/mech), and only the namelist name should be specified. Several param files are provided as example. This script is in charge of creating the subdirectory for output (which will be named the same as the param file name), and launch the model.

Figure 4 Contour plot of turbulent kinetic energy (TKE, in m² s⁻²) as a function of height and simulation time. The thick line shows the bottom of cloud, defined as LWC dropping below 0.01 g kg⁻¹. Top: study from Loughlin et al. (1997, Fig. 4a). Bottom: MISTRA-v9.0. For both, the simulation settings are identical to those in Bott et al. (1996). Top panel reproduced with permission from Loughlin et al. (1997, Fig. 4a). © 1997 by John Wiley and Sons. ​​​​​​​

Figure 5 Distribution of particles in the 2D particle grid and at different heights in the cloud. This graph is built identical to Bott et al. (1996, Fig. 12): for each dry radius class (on x axis), the total radius containing the maximum number of particles is marked with a filled circle; total radii containing 50 % to 99 % of the maximum are marked with open circles, and total radii containing 1 % to 50 % of the maximum are marked with plus signs. The 1:1 values (i.e. where total radius equal to dry radius) are represented with cross signs. The full line shows the activation radius as accounted for in MISTRA. From top to bottom row is top (785 m), middle (605 m), and bottom (555 m) of the cloud. Left column: version of 1996, panels are from Bott et al. (1996, Fig. 12). Right column: MISTRA-v9.0. Left column reproduced with permission from Bott et al. (1996, Fig. 12). © 1996 by John Wiley and Sons.

Figure 6 1D distribution of particle mass as a function of radius. (a) No-chem-MISTRA without collision–coalescence. (b) No-chem-MISTRA with collision–coalescence. (c) MISTRA-v9.0 (without collision–coalescence implemented). Each panel shows the distribution in the top (black line), middle (red line), and bottom of cloud (green line). (a), (b) ​​​​​​​ are from Bott (2020, Fig. 2, published under CC BY 4.0 licence). The simulation settings were taken from Bott (2020).

4 Consistency with previous versions

In this section, we present a series of example runs that have been performed to evaluate the model. All the examples provided here reproduce the settings of previous studies carried out with previous versions of MISTRA. For that purpose, several namelists have been introduced to hold all relevant parameters in order to reproduce the same simulation scenarios as in the original publications. These namelists, as well as the scripts used to produce the plots presented here, are available in the MISTRA repository.

4.1 Meteorology and microphysics

4.1.1 Comparison with 1996 version: LWC, TKE and 2D spectrum

The first example focuses on the physical and microphysical aspects of the model. For this purpose, the chemistry is switched off. The initialisation settings are identical to those of the original paper from Bott et al. (1996) and are provided in the namelist BTZ96. Some model changes have been maintained for this comparison, even if this leads to differences to the original version. For instance, the number of bins in the 2D particle spectrum is set to 70×70 in MISTRA-v9.0 while it was 40×50 in the version of Bott et al. (1996). However, we adjusted the minimum and maximum particle radius values so that the resolution is nearly identical in both versions.

This simulation reproduces conditions over the North Sea, for 3 simulated days (2 are shown, the first one is used as model spin-up) centred on 22 July. The radiative code used by Bott et al. (1996) has been updated from PIFM1 to PIFM2 by Loughlin et al. (1997), and the simulation settings were the same in both papers. For this reason, we compare figures from the study of Loughlin et al. (1997) with the current output of MISTRA.

Both Figs. 3 and 4 show very similar model output between the 1996 version of MISTRA and the current version. The runs are similar, qualitatively and quantitatively (the maximum LWC is 6 % higher, the maximum TKE is 2.5 % higher in MISTRA-v9.0 than in the 1996 version), without changing the findings and conclusions of the original study.

Bott et al. (1996) also showed the distribution of particles in the 2D grid, and we used the same graph format in Fig. 5. Qualitatively, the two simulations are similar. MISTRA-v9.0 exhibits more particle growth for the smallest dry radius bins; however, this happens only for a minority of the particles (“plus” signs in Fig. 5 denotes bins where particle concentration is less than half the maximum particle concentration, for each dry radius class). As stated previously, the 1996 version of MISTRA used by Bott et al. (1996) included the first version of the radiative code PIFM1, now updated to PIFM2. The differences between both radiation schemes are likely the reason for these slightly different particle distributions observed for the small dry particle radius. This figure also highlights the microphysical properties and dynamics of particle within the cloud, with the activation of particles occurring when the supersaturation (not shown) is high enough, which is the case in the upper part of the cloud. Conversely, in the middle and bottom parts of the cloud, most of the particles are found below the critical radius, even if some particles grow to above their respective critical radius, since the supersaturation is not high enough.

Figure 7 Contour plot of humidity, potential temperature, and temperature in the run based on Joyce et al. (2014). (a) Original study (figure published under CC BY 3.0 licence). (b) MISTRA-v9.0.

Figure 8 Contour plot of NO_x (a, b) and PM_2.5 (c, d) in the run based on the study of Joyce et al. (2014). (a, c) Original study (figure published under CC BY 3.0 licence). (b, d) MISTRA-v9.0. Scales are identical for both.

4.1.2 Impact of neglecting coalescence

As pointed out in the general presentation of the model (Sect. 2.2), the collision–coalescence process is not accounted for in MISTRA-v9.0, which is a limitation of the model. The collision–coalescence process was implemented in a version of MISTRA without chemistry (Bott, 2000), hereafter referred to as MISTRA-coal-nochem for brevity. In a recent study, Bott (2020) used MISTRA-coal-nochem and compared the results with and without activating this process (Fig. 6). He showed that accounting for the collision–coalescence process leads to significant differences in the particle distribution, with a bimodal spectrum of particles larger than 40 µm when collision–coalescence is activated (Fig. 6b) Conversely, when particles grow solely by diffusional uptake of water vapour, their size distribution remain in the 2 to 30 µm range (Fig. 6b). We defined a namelist, named Bott2020, reproducing the same settings as in Bott (2020), to perform a further evaluation of MISTRA-v9.0 against MISTRA-coal-nochem. The resulting 1D particle distribution is presented in Fig. 6c, and shows similar results as compared to MISTRA-coal-nochem (Fig. 6a).

Figure 9 Contour plot of height versus time for four gases and two aqueous phase species in the run based on Joyce et al. (2014). Left side: original study (figure published under CC BY 3.0 licence). Right side: MISTRA-v9.0. Scales are identical.

Figure 10 Gas-phase concentration of Br₂ and BrCl in a box-model run based on the settings of Buys et al. (2013). (a, b) Original study (Fig. 3 published under CC BY 3.0 licence). (c, d) MISTRA-v9.0. Scales are identical.

Figure 11 Gas-phase concentration in a chamber-model run based on the study of Buxmann et al. (2015) (alpha wl sensitivity experiment). Left side: original study. Right side: MISTRA-v9.0 with identical settings. Scales are identical. Left panel reproduced with permission from Buxmann et al. (2015, Fig. 5). © 2015 by CSIRO Publishing.

Despite the important differences in particle distribution when collision–coalescence process is included, this limitation in MISTRA-v9.0 is expected to have an insignificant effect on simulations without clouds (non-activated particles only). Conversely, cloudy runs should be restricted to conditions where collision–coalescence is less important, i.e. cases where no or little drizzle formation would be expected. According to Duynkerke (1998), drizzle formation starts to be important when the cloud depth is greater than 300 m. Future development plans with MISTRA-v9.0 include a re-evaluation of the feasibility of including the collision–coalescence process along with chemistry.

4.2 Chemistry in 1D simulations

A namelist reproducing the settings of the study by Joyce et al. (2014) is provided as namelist.Joyce14bc. In this study, MISTRA was used to simulate an urban pollution plume from Fairbanks, Alaska. The model was thus used in an alternative configuration, with the surface covered by snow (with the relevant physical properties). An emission scenario of NO_x (NO+NO₂) was defined, and the evolution of gas- and aqueous-phase species was evaluated. In such configuration, the meteorological parameters have a strong influence over the stability of the atmosphere, thus in turn over the vertical exchange of chemical species. In Fig. 7, key meteorological variables are presented for both the original study and the new runs obtained with MISTRA-v9.0. As expected, there is excellent agreement between both versions, which shows that the recent code developments did not alter the model with regards to the plotted variables.

Figures 8 and 9 show the comparison of gas- and particle-phase chemical species, between the original study of Joyce et al. (2014) and MISTRA-v9.0. Again, the plotted variables agree very well between both versions, with nearly identical results. In Fig. 9, the only exception is ammonium (${{\mathrm{NH}}_{\mathrm{4}}}^{+}$, bottom right panel on each side) whose maximum concentration is decreased by 60 % in MISTRA-v9.0. The reason for this change was investigated, and we found that in the original run, the initialisation of a variable in the routine computing the gas–particle exchange rates was missing. This is now corrected in MISTRA-v9.0, and explains the differences for ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$.

4.3 Box and chamber model configurations

We present two additional configurations of MISTRA-v9.0, as an atmospheric box (0D) model (Buys et al., 2013), and in chamber configuration following the study of Buxmann et al. (2015). In both cases, we set namelists (namelist.Buys13_0D and namelist.Buxmann15_alpha, respectively) with the same settings as in the original publications.

Figures 10 and 11 compare model output between each model version, and show that minor differences exist but are very limited, and the results agree well qualitatively. This comes as a further demonstration of the good consistency of MISTRA-v9.0 with previous results.

5 Conclusions

We have presented the current version of the 0D/1D atmospheric chemistry model MISTRA-v9.0, released for the first time as an open-source, community model. MISTRA-v9.0 is a versatile model with a range of capabilities, from the study of stratus cloud microphysics, radiative forcing and turbulence, to the multiphase atmospheric chemistry of the boundary layer. While its original purpose was only the study of cloud-free and cloudy MBL, MISTRA was successfully extended in previous studies to model other environments such as polar conditions and volcanic plumes. In this study, we updated the model code to comply with coding standards, and we compared current output of a range of test cases against previous studies with identical settings. Results obtained with MISTRA-v9.0 are consistent with the previous results, even after 20 years of development. MISTRA-v9.0 is a powerful tool for atmospheric chemistry research purposes, now easier to use, and free to use under the EUPL-v1.1 licence. Community input and development is welcome for MISTRA-v9.0.

Appendix A

Table A1 List of symbols.​​​​​​​

Appendix B: List of abbreviations

DMS	Dimethyl sulfide (CH3–S–CH3)
ESM	Earth system model
EUPL	European Union Public Licence
FT	Free troposphere
KPP	The Kinetic PreProcessor
LWC	Liquid water content
MBL	Marine boundary layer
MIFOG	MIcrophysical FOG model
netCDF	Network Common Data Form
PIFM	Practical Improved Flux Method (radiative code)
TKE	Turbulence kinetic energy

Appendix C: Subsidence profiles

Three subsidence profiles are currently implemented in MISTRA-v9.0.

• Option 1 follows the hyperbolic expression of Bott et al. (1996, Eq. 5).

• Option 2 linearly decreases from wmin at ground level to wmax at height h  = etaw1.

• Option 3 linearly decreases from wmin at ground level to wmax at height h = zinv.

Figure C1 Vertical profiles of subsidence, computed for w_min=0 cm s⁻¹, $w_{\max }=-\mathrm{0.5}$ cm s⁻¹, and z_inv=700 m.