A description of the first open source community release of the 1D atmospheric chemistry model MISTRA-v9.0

We present MISTRA-v9.0, a one dimensional (1D) 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 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. In the past 20 years, MISTRA has been used in over 25 studies to address a wide range of scientific questions. The purpose of this public 5 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 version 9.0 is consistent with previous versions.

introduction of a module for aerosol nucleation which improved the iodine chemistry (Pechtl et al., 2006. 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;Lawler et al., 2009;Jones et al., 2010;Joyce et al., 2014) including alternative model configurations where the chemistry was computed in a zero-dimension (0D) atmospheric box mode (Buys et al., 2013), and a 0D chamber mode 60 (Buxmann et al., 2015). MISTRA was also adapted to model specific environments such as volcanic plumes Bobrowski et al., 2007Bobrowski et al., , 2015 polar conditions von Glasow, 2008, 2009;Buys et al., 2013), and was coupled with an ocean model for studies over the Dead Sea von Glasow, 2007, 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(Thomas et al., , 2012. In this study, we present a 65 selection of a few specific model settings reproducing previous work, to compare the original results with those obtained with MISTRA-v9.0.

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 . A version 8 featuring an alternative chemical bin definition was partly developed but not 70 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). 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, 75 and conform to strict coding rules (Metcalf et al., 2004). 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 doloops should be leftmost indexes). The chemical "Kinetic PreProcessor" (KPP) has been updated to the latest version 2.2.3 (https://people.cs.vt.edu/~asandu/Software/Kpp/ last accessed 23 June 2021). Overall, several technical developments have 80 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.  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 90 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 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. A nucleation module is also included to account for new particles nucleated from the gas phase species.

Meteorology, microphysics and thermodynamics
The model is one-dimensional, 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: 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 105 coefficients for momentum and heat, L is the latent heat of condensation, C the condensation rate, ρ the air density, p the air pressure, p 0 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 one-dimensional 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 110 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. 115 (1996) and Bott (1997). The turbulent exchange coefficients K m and K h are calculated via stability functions S m/h and G m/h , where the subscript m stands for shear and h for buoyancy production. The prognostic equation for the turbulence kinetic energy (TKE) e is: assuming a constant dissipation ratio (last term on the right). For more details and an explanation of the calculation of the 120 mixing length l, the exchange coefficient K e for the TKE, and the functions S m/h and G m/h see Mellor and Yamada (1982), Bott et al. (1996), andBott (1997).
The microphysics is treated using a joint two-dimensional particle size distribution function f (a, r) where a is the dry aerosol radius the particles would have if no water is present in the particles, and r is the total particle radius. The two-dimensional particle grid is divided into 70 logarithmically equidistant spaced dry aerosol classes. The minimum aerosol radius is generally 125 set to 0.005 µm and the maximum radius 15 µm. Choosing these values allows 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 aerosol classes is associated with 70 total particle radius classes, ranging from the actual dry aerosol radius up to 60 µm (150 µm in cloud runs). See Figure 2 for a depiction of 2D particle grid. The prognostic equation for f (a, r) is: 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 (ṙ = dr/dt). The 2D particle spectrum is initialised with distribution depending on the type of aerosol chosen (see Bott, 1997, and ref. 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 135 (see for instance 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 integration, particle growth is calculated explicitly for each bin of the 2D particle spectrum using the growth equation after Davies (1985) Bott et al. (see also 1996): 140 with the ambient supersaturation S ∞ and the supersaturation at the droplet's surface S r according to the Köhler equation: 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 145 C 1 and C 2 in equation (7) are: m w (a, r) is the liquid water mass of the particle, c w and ρ w are the specific heat and 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 v have been corrected for gas kinetic effects following Pruppacher and Klett (1997) (their equations 13.20 150 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 equation (4) is determined diagnostically from the particle growth equation (7).
Collision-coalescence processes are not included in the model because this leads to difficulties when redistributing the 155 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 (PIFM2 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 160 and particle growth are fully implemented. The calculation of photolysis frequencies is described in Sect. 2.3.5.

Chemistry
The multiphase chemistry module comprises chemical reactions in the gas phase as well as in 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 165 by von Glasow and Crutzen (2004) to include a better description of the oxidation of DMS. Iodine chemistry was significantly improved by Pechtl et al. (2006. 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.

Gas phase and uptake
The prognostic equation for the concentration of a gas phase chemical species c g (in mol m −3 ) including subsidence, turbulent exchange, deposition on the ocean surface, chemical production and destruction, emission and exchange with the aqueous phases is: Again subsidence is the first term on the right and is included only in runs with clouds, otherwise w = 0. P and S are chemical production and loss terms, respectively, H cc s is the dimensionless Henry constant obtained by H cc  (10) describes the transport from the gas phase into the aqueous phases according to the formulation by Schwartz (1986) (see also Sander, 1999), n kc is the number of aqueous classes as explained below. For a single particle, the mass transfer coefficient k t is defined as: with the particle radius r, the mean molecular speedv = 8RT /(M π) (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 g = λv/3 (e.g. Gombosi (1994), p. 125). Chameides (1984) points out that the time needed to establish equilibrium between the gas and aqueous phase differs greatly for individual species and that soluble species never reach equilibrium in cloud droplets, emphasizing the importance 190 of describing phase transfer in the kinetic form that is used here. Audriffen et al. (1998) and Chaumerliac et al. (2000) point out that for reactive species like H 2 O 2 , 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 particle with different radii. The mean transfer coefficient k t for a particle population is given by the integral: where the size distribution function ∂N/∂ lg r depends on the type of aerosol chosen.

Aqueous phase
Aqueous chemistry is calculated in four bins (see Figure 2): deliquescent aerosol particles with a dry radius less than 0.5 µm are included in the "sulfate aerosol" bin #1, whereas deliquescent particles with a dry aerosol radius greater than 0.5 µm are in 200 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 get internally mixed by exchange with the gas phase but, as mentioned earlier, not by particle collisions. 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 205 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 210 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 (in mol m −3 ), where the index i stands for the i-th aqueous bin: The individual terms have similar meanings as in equation (10). The calculation of the sedimentation velocity v dry a,i , that is 215 needed for the calculation of the dry deposition D, is explained at the end of this section. 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, equation (13) reduces in steady state conditions (∂c a,i /∂t = 0) to the Henry equilibrium c a,i = w l,i c g K cc h . The concentration of H + ions is calculated like any other species, i.e. no further assumptions are made. The charge balance 220 is satisfied implicitly.

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 aerosol remains in a highly concentrated metastable aqueous state when they are dried below their deliquescence humidity. Only when they reach the crystallisation humidity they 225 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 (1998) and references therein), implying that aerosol particles that already had been involved in cloud cycles will also be in an aqueous metastable state. Therefore most soluble aerosol particles will be present in the atmosphere 230 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.

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 235 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 Luo (Luo et al. (1995)

Photolysis
Here an overview of the calculation of the photolysis rates is given, for a detailed description see the model manual, Chapter 5.
Photolysis is calculated online using the method of Landgraf and Crutzen (1998). The photolysis rate (or photo dissociation coefficient) J X for a gas X can be calculated from the spectral actinic flux F (λ) via the integral: where λ is the wavelength, σ X the absorption cross section, φ X the quantum yield and I the photochemically active spectral interval. If the integral in equation (14) were 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 (14) by: where J a i,X is the photolysis rate for a purely absorbing atmosphere. The factor δ i : describes the effect of scattering by air molecules, aerosol and cloud particles. F a (λ i ) is the actinic flux of a purely absorbing atmosphere. The factor δ i is calculated online for one wavelength for each interval, while the J a i,X are pre-calculated with a fine spectral resolution and are 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 260 by a factor of two to account for the actinic flux enhancement inside the particles due to multiple scattering.

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

265
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. (1993a) are implemented to estimates 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 dry g at the sea surface is calculated using the resistance model described by Wesely (1989): The aerodynamic resistance r a is calculated using: 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 0 and a reference height z. The quasi-laminar layer resistance r b is parameterised 275 for gases as: The Schmidt number Sc can be written as Sc = ν/D 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 (1998) 280 with the effective Henry constant H * .
The dry deposition velocity of particles v dry a,i is calculated after Seinfeld and Pandis (1998): where the quasi-laminar resistance r b is parameterised for particles as:

285
The Stokes number St can be written as St = w t u 2 * /(gν). 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 after Beard for larger particles (see Pruppacher and Klett (1997)).  Table 1 presents the switches available to define the model configuration.

Initialisation and run settings
Initial atmospheric conditions are set in the namelist with the parameters presented in Table 2. Note that default values are 295 set for all of them, however they should be systematically redefined by the user to match the simulated atmosphere. Standard settings (timing and geography, run duration) are straightforward and are not detailed hereafter. In addition to these, surface settings are detailed in Table 3.

Special runs setting 300
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 require special settings. An example of such 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 Joyce et al. (2014) study (particle distribution and composition, gas and particle emission scenario, special formulation of accommodation coefficient for N 2 O 5 , and of gas dry deposition,   netcdf request netCDF output files (true, default setting) or not (false).

Get the code
The model is provided on GitHub on the following repository: https://github.com/Mistra-UEA/Mistra. It is released under the 310 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.

System requirements and installation
A Fortran compiler is mandatory to compile the model code. During the recent development stages, MISTRA has been regularly 315 compiled using either GNU Fortran (gfortran) or Intel Fortran (ifort). New users are advised to choose one of those compilers.   µm minimum dry particle radius in the 2D particle spectrum. rnw1 µm maximum dry particle radius in the 2D particle spectrum. rw0 µm minimum total particle radius in the 2D particle spectrum. rw1 µm maximum total particle radius in the 2D particle spectrum.
zinv m inversion height. It must be lower than the highest prognostic layer. rhMaxBL 1 maximum relative humidity (default=1) in the MBL (additional constrain to xm1w parameter).
rhMaxFT 1 maximum relative humidity in the FT (additional constrain to xm1i parameter).
cGasListFile names of user gas files for non radical species.
cRadListFile names of user gas files for radical species.
scaleO3_m D.U. ozone column, to scale the photolysis rates computed. It has no effect on the radiation calculation.
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 and NCL (NCL), 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 320 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.

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 325 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.

Compile the model 330
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.

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 335 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.

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 340 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.

Consistency with previous versions
In this section, we present a series of example runs that have been performed to evaluate the model. 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 of 40 × 50 in the 355 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) centered on the 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 360 study of Loughlin et al. (1997) with the current output of MISTRA.
Both Fig. 3 and Fig. 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 shown the distribution of particles in the 2D grid, and we used the same graph format in Fig. 5.

365
Qualitatively, the two simulations are similar. MISTRA-v9.0 exhibits more particle growth for the smallest dry radius bins;  the simulation settings are identical to those in Bott et al. (1996). Top panel reproduced with permission from Loughlin et al. (1997, Fig. 4a).
Copyright (1997)    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 this slightly different particle distributions observed for 370 small dry particle radius. This figure also highlights the microphysics 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.

Impact of neglecting coalescence 375
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-coalescencce process leads to significant differences in the particle distribution, with a bimodal spectrum 380 with 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 one-dimensional particle distribution is presented in Fig. 6c, and shows similar results as compared to MISTRA-coal-nochem. Despite the important differences in particle distribution when collision-coalescence process is included, this limitation in MISTRA-v9.0 is expected to have insignificant effect for 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 390 collision-coalescence process along with chemistry.

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 surface covered by snow (with the relevant physical properties). An emission scenario of NO x was defined, 395 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. 400 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 (NH 4 + , bottom right panel on each side) whose maximum concentration is decreased by 60 % in MISTRA-v9.0. The reason of 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, 405 and explains the differences for NH 4 + .

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.    -Option 2 linearly decreases from wmin at ground level to wmax at height = etaw1.