PyCHAM (v1.3.4): a Python box model for simulating aerosol chambers

In this paper the CHemistry with Aerosol Microphysics in Python (PyCHAM) box model software for aerosol chambers is described and assessed against benchmark simulations for accuracy. The model solves the coupled system of ordinary differential equations for gas-phase chemistry, gas-particle partitioning and gas-wall partitioning. Additionally, it can solve for coagulation, nucleation and particle loss to walls. PyCHAM is open source, whilst the graphical user interface, modular structure, manual and suite of tests for troubleshooting and tracking the effect of modifications to individual modules have been 5 designed for optimal usability. In this paper, the modelled processes are individually assessed against benchmark simulations, and key parameters described. Examples of output when processes are coupled are also provided. Sensitivity of individual processes to relevant parameters is illustrated along with convergence of model output with increasing temporal and spatial resolution. The latter sensitivity analysis informs our recommendations for model setup. Where appropriate, parameterisations for specific processes have been chosen for their general applicability with their rationale detailed here. It is intended that 10 PyCHAM aids the design and analysis of aerosol chamber experiments, with comparison of simulations against observations allowing improvement of process understanding that can be transferred to ambient atmosphere simulations.

processes represented in PyCHAM through comparison with benchmark simulations.
Demonstration of PyCHAM application in which all major processes are influential is provided by simulating an experiment based on the role of nitrate radical (NO 3 ) oxidation of limonene in SOA evolution (Fig. 1). Such an experiment has implications for indoor air quality at night time when the photolysis of NO 3 ceases (Waring and Wells, 2015), therefore lights were turned off for this simulation. Following a similar approach to the experiment of Fry et al. (2011), the effect of NO 3 in the presence 70 of ozone (O 3 ) can be replicated through introduction of O 3 and nitrogen dioxide (NO 2 ) to the chamber whilst removing the effect of the hydroxyl radical (OH) through addition of excess carbon monoxide (CO). At sufficient concentrations, this mixture initiates particle nucleation (Fry et al., 2011), which we simulate here through the nucleation parameterisation described below.
O 3 , NO 2 and limonene are injected again later but with the addition of seed aerosol for reproduction of indoor environments with substantial existing particulate matter. 75 In Fig. 1 the particle organic nitrate curve demonstrates that around 15 % of SOA comprises organonitrates that result from NO 3 reaction. The particle inorganic nitrate curve represents the contribution of dinitrogen pentoxide (N 2 O 5 ) and nitric acid (HNO 3 ) to particles. Here we simulate the hydrolysis of N 2 O 5 into the aqueous phase of particles and wall by setting its activity coefficient to zero and its accommodation coefficient according to Lowe et al. (2015). Studies have recently revealed the role of highly oxidised molecules (HOM) (Ehn et al., 2014), with the Peroxy Radical Autoxidation Mechanism (PRAM) 80 simulating their chemistry (Roldin et al., 2019). For the results in Fig. 1, the PRAM scheme has been coupled with that of the Master Chemical Mechanism (MCM) (Jenkin et al., 1997;Saunders et al., 2003). It is intended that simulations such as The primary difference between multiphase processes in simulation chambers and the real atmosphere is the presence of 85 walls. Accurate representation of these processes in chamber models requires reasonable and realistic representation of the chamber walls (e.g Matsunaga and Ziemann, 2010;Zhang et al., 2015). Clearly PyCHAM must reasonably capture the partitioning of components and deposition of particles to chamber walls as incorrect reproduction makes comparison against Figure 1. Limonene oxidation in the dark with and without seed particles. In (a), the gas-phase concentrations of key components and in (b), the particle properties. At the start 10 ppb limonene, 22 ppb NO2 and 500 ppm CO are introduced. At 1.5 hours 38 ppb O3 and 8 ppb NO2 are injected (injection 1). At 4 hours a further injection of O3 (45 ppb), limonene (10 ppb) and NO2 (19 ppb) is coincident with an injection of seed aerosol (10 µg m −3 with a mean diameter of 0.2 µm) (injection 2). In (b), the total particle number concentration ([N]) corresponds to the the first of the right axes, mass concentrations of: secondary organic aerosol mass concentration (SOA), components deposited to walls (wall deposit), sum of particulate organic components with a nitrate functional group (particle organic nitrate) and sum of particulate inorganic components with a nitrate functional group (particle inorganic nitrate) correspond to the second of the right axes.
Number size concentrations correspond to the filled contours, colour bar and left axis. Whilst [N] includes seed particles, SOA excludes the seed material and all mass concentrations exclude water. measurements misleading or impossible. However, with correct wall loss constraint, simulations such as Fig. 1 can be compared against observations for verifying process understanding. 90 3 PyCHAM structure For ease of navigation, PyCHAM has a modular structure with each key physicochemical process assigned an individual module. Unit tests are provided for key modules, which allow the user to check a particular process is functioning correctly. It is intended that these tests will be useful for troubleshooting and for analysing the effects of modifications to modules.
At the core of PyCHAM lies simultaneous numerical integration of three coupled processes: gas-phase photochemistry, 95 vapour-particle partitioning and vapour-wall partitioning. The ordinary differential equations (ODEs) for these processes are solved by the backward differentiation formula, for which utility has been demonstrated in similar atmospheric applications (Jacobson, 2005), from the CVODE Sundials software (Hindmarsh et al., 2005). We use Assimulo (Andersson et al., 2015), a python wrapper for sundials allowing communication between the solver and Python code. The model structure is outlined in the schematic of Fig. 2, where we introduce the user-defined 'update time interval', which is the time between updates to 100 4 https://doi.org/10.5194/gmd-2020-234 Preprint. Discussion started: 20 August 2020 c Author(s) 2020. CC BY 4.0 License. On pressing the 'Run Model' button PyCHAM loops over the 'update time interval' until the experiment end time is reached. On each loop PyCHAM first checks whether any discontinuous changes to the chamber condition occur during the proposed time interval, and automatically reduces the interval if confirmed, such that the change can be implemented at the correct time at the start of the subsequent interval. Depending on problem stiffness the integrator sets sub-time steps and generates results at each. The final stage of a loop where the update time interval has been reached is solution of coagulation, particle deposition to walls and nucleation, which update the particle number concentration. The right column notes these temporal resolution aspects, with the mentioned sensitivity of key system properties investigated in Section 11. natural light flux and particle number concentration due to coagulation, particle deposition to wall and nucleation (the former applicable during open roof experiments and the latter applicable for experiments with particles). Particle number concentration and photolysis rates are constants in the ODEs for gas-particle partitioning and gas-phase chemistry, respectively. The update time interval is passed to the integrator, which adaptively sets sub-time steps depending on problem stiffness. For discontinuous changes to chamber conditions, such as injection of gas-phase components or seed particles, automatic adaption ensures they 105 occur at the start of a time interval.
Whilst coagulation and particle loss to wall have timescales of minutes to hours, nucleation can cause substantial changes within seconds (Section 11). Users may increase the update time interval, which has the advantage of decreased simulation time, and the disadvantage of divergence from high resolution estimates. Simulation sensitivity to temporal resolution is investigated in Section 11, including recommendations for maximum time intervals. 110 5 https://doi.org/10.5194/gmd-2020-234 Preprint. Discussion started: 20 August 2020 c Author(s) 2020. CC BY 4.0 License. In the example model output of Fig. 1 several features of PyCHAM are demonstrated. First, the coupling of gas-phase chemistry and resulting partitioning of vapours with sufficiently low volatility to particles and walls. Nucleation has been simulated prior to the introduction of seed particles with approximate values given for the tuned nucleation parameters (Section 10) that ensure nucleation begins at the introduction of ozone, has a duration of thirty minutes and produces a peak number concentration similar to that of the seed particle. Finally, coagulation and particle wall loss (the latter using the model of McMurry and 115 Rader (1985)), contribute to the decay in particle number concentration.
The PyCHAM software is initiated via the command line to generate a graphical user interface (GUI). Via the GUI, users select three files (Fig. 2) representing: i) the chemical scheme, ii) a file associating the chemical identifiers inside the chemical scheme to their Simplified Molecular Input Line Entry System (SMILE) strings (Weininger, 1988), and iii) a model variables file. A fourth button on the GUI starts the simulation.

120
A parsing module interprets the chemical scheme and uses the chemical identifier conversion file to match component identifiers to their SMILE strings. Additionally, modules are automatically created that will calculate chemical tendencies and track the process tendencies (change due to chemistry and partitioning) of specified components. A gas-phase initiation module interprets the user-defined starting concentrations of components, whilst a particle-phase initiation module establishes any seed particles at experiment start. The integration module is then called where the ODEs for gas-phase photochemistry and 125 partitioning are solved (Fig. 2). A saving module stores results by default for gas and wall concentrations, corresponding time and constants such as component molecular weight. If the user has setup the simulation to include particles, then component particulate concentration and particle number size distributions (with and without water) are also saved, furthermore, if the user has defined components to track, then the process tendencies of these are saved.
The fifth and final button on the GUI will display and save graphs of the temporal profiles of number size distribution, 130 secondary aerosol mass concentration, total particle number concentration, and the gas-phase concentrations of the components whose initial concentrations are user-defined. Besides these default plots, users can find additional examples of plotting scripts (those used for figures here) in the PyCHAM repository. The programme can be stopped via the terminal when in integration mode, or outside this mode it can be terminated by closing the GUI.
Below we describe and verify the processes described above as coded in PyCHAM. Necessarily each process is examined 135 in isolation, however, Fig. 1 and its associated text illustrate the coupling of mechanisms for a real world application.

Model variables and component properties
As described above, to initiate PyCHAM the user selects a completed model variables input file (an example is provided with the software). The available variables are extensive to allow adaptability to a range of experiments, consequently for a given experiment, many of the variables in this file may be left empty Here we introduce the available variables, whilst details such 140 as default values and units are provided in the appendix Table A. Several model variables (PyCHAM names of variables given in brackets) are purely functional, these include the name of the output file (res_file_name), whether to update the component property estimation files (umansysprop_update), which are described below in this section, the markers used to separate sec-tions of the chemical scheme (chem_scheme_markers), names of files containing actinic flux (act_flux_file) and absorption cross-sections and quantum yields for photochemistry (photo_par_file). Total simulation time (total_model_time), the time 145 interval for updating integration constants (update_step) and the time interval for recording results (recording_time_step) are also available.
Chamber temperature (temperature) can change during a simulation by stating the corresponding time (tempt), whilst pressure (p_init) and relative humidity (rh) are also input. For simulations involving natural light, latitude (lat), longitude (lon), day of year (DayOfYear) and start time (daytime_start) are required inputs. Whether artificial or natural, users specify when 150 (light_time) light is on or off (light_status). Any dilution rate (dil_fac) should be stated, or else the default is zero.
Initial concentrations (C0) of specified trace gases (Comp0) are stated separately to the concentrations (Ct) of specified trace gases (Compt) injected effectively instantaneously at set times (injectt) during the experiment. Specified components (const_comp) will have a constant concentration for the entire experiment. The final option for introducing named components (const_infl) is to state the rate of their influx (Cinfl) during a set period of the experiment (const_infl_t). The process tendencies 155 of certain components (tracked_comp) can be recorded, which is helpful for analysis and troubleshooting.
To simulate gas-wall partitioning, the mass transfer coefficient (kgwt) and effective absorbing wall mass concentrations (eff_abs_wall_massC) can be set.

160
Whether the experiment is seeded or involves nucleation, users state the number of size bins (number_size_bins), size at lowermost size bin boundary (lower_part_size) and at uppermost boundary (upper_part_size), and whether to have linear or logarithmic spacing of size bins (space_mode). Setting size bin number to zero turns off particle considerations.
For seeded experiments, the component comprising the seed (seed_name), its molecular weight (seed_mw), density (seed_dens) and dissociation constant (core_diss) can be input. Either the particle concentration per size bin or the total particle concentra-165 tion can be input (pconc), along with the time of particle injection (pconct). If the total particle concentration is given, this can be distributed across size bins by stating the mean radius (mean_rad) and standard deviation (std).
For nucleation experiments, the nucleating component (nuc_comp) can be changed from the default, as can the radius of newly nucleated particles (new_partr). To specify the temporal profile of nucleation, three parameters (detailed in Section 10) are input (nucv1, nucv2, nucv3).

170
Coagulation (coag_on) and particle loss to wall (McMurry_flag) can be turned off and on. If the latter is turned on, users can specify the size-dependent loss to walls (inflectDp, Grad_pre_inflect, Grad_post_inflect and Rate_at_inflect), or can invoke the McMurry and Rader (1985) model by also inputting the chamber wall surface area (Cham_SA), the charge per particle (part_charge_num) and the chamber electric field (elec_field), which are detailed in Section 9.
The components included in the user-defined chemical scheme are automatically allocated three properties by the PyCHAM  (Topping et al., 2016) which is updated on the first run of PyCHAM and 180 at the request of the user (via the model variables file) thereafter (requires internet connection). By default the UMansSysProp module applies the liquid density estimation method of Girolami (1994) (recommended by Barley et al. (2013)) and the liquid saturation vapour pressure estimation method of Nannoolal et al. (2008) (recommended by O'Meara et al. (2014). Component vapour pressures have a first order effect on absorptive partitioning between phases, however estimates for certain components, particularly those with relatively low vapour pressures as these are most difficult to measure experimentally and therefore 185 inform estimation methods, are associated with considerable uncertainty (O'Meara et al., 2014). Consequently, users can also specify the vapour pressures of certain components. Similarly, although the default activity and accommodation coefficient for all components partitioning to particles and wall is unity, users may set an alternative value for specific components. At present, activity coefficient calculations are not incorporated into PyCHAM.

190
For a chamber experiment including injection of reactive components, chemical reactions in the gas-phase drive the disequilibria that can affect the composition of gas, particle and wall. As mentioned above, schemes such as the MCM provide near-explicit gas-phase chemistry mechanisms for numerous organic precursors, and developments such as PRAM (Roldin et al., 2019) can be used to provide supplementary detailed updates to our understanding of atmospheric chemistry. PyCHAM is designed to accommodate any such detailed chemical schemes whilst also accepting very simplified or even empty (e.g.

195
for a control simulation comprising only seed particles) chemical equation files. Whilst the software manual details the requirements for input chemical schemes and chemical identifier conversion files, here we describe how PyCHAM deals with chemistry. Equations of the general form: where s represents stoichiometric number, r reactants and p products, are expressed as the ODEs: where n is the total number of reactants and r j is a given reactant for a given reaction. k r is the reaction rate coefficient.
For simulations involving gas-phase chemistry, users must therefore provide a reaction(s) of the form in Eq. R1 and an associated reaction rate coefficient inside a chemical scheme file. Naming of chemical components inside the chemical scheme 205 is unrestricted, however, the software must be able to convert names to SMILES (Weininger, 1988). Therefore, users must provide a separate file stating a unique SMILES string for every component (Fig. 2). Examples of both the chemical scheme and SMILES string conversion file are included in the software.
Inside the parsing module, reaction rate coefficients, reactant and product identities and their stoichiometric numbers are established from the chemical scheme file. To separate these properties either default formatting may be used, or a variant, so 210 long as the appropriate changes are made inside the model variables file. By default, MCM Kinetic PreProcessor (Sander and Sandu, 2006) formatting is used (Jenkin et al., 1997;Saunders et al., 2003) and PyCHAM has been rigorously tested using schemes and SMILE conversion files from the MCM website (Rickard and Young, 2020).
Reaction rate coefficients can be functions of temperature, relative humidity, pressure and concentrations of: third body, nitrogen, oxygen and peroxy radicals. Third body, nitrogen and oxygen concentrations are calculated by the ideal gas law with 215 the user-set temperature and pressure. As in the MCM, the chemical scheme file can include generic reaction rate coefficients (those that have an identifier which is used as the reaction rate coefficient for one or more reactions).
Photochemistry is controlled through stating light on/off times inside the model variables file. The treatment of photochemistry is determined by the user and depends on the chemical scheme employed. In the case of the MCM scheme and natural sunlight, the scattering model based on Hayman (1997) and described in Saunders et al. (2003) is invoked by stating the rel-220 evant spatial and temporal coordinates in the model variables file. For artificial lights, users must provide a file stating the wavelength-dependent actinic flux (as described in the manual). The model then calls on either the absorption cross-section and quantum yield estimates of MCM v3.3.1 or of a user-defined file.

Assessment of gas-phase chemistry accuracy
To assess the accuracy of the photochemistry section of PyCHAM, gas-particle partitioning and gas-wall partitioning were 225 turned off, leaving only gas-phase chemistry to be solved. Here we compare against AtChem2 (Sommariva et al., 2018) as a model benchmark, with both using MCM chemical schemes. Figure  whilst for NO x the initial concentration was 9.8 ppb in Fig. 3a and 0 ppb in Fig. 3b. Latitude was set to 51.51, longitude to 0.13 (London, UK) and the date to 1 July. The deviation between PyCHAM and AtChem2 was calculated using: is the AtChem2 maximum for a given component during the simulation which is the chosen scaling factor for deviations as it 235 means any difference between model estimates is referenced against a reasonable value for that component (in contrast scaling by b i,t when b i,t ∨(b i ) may introduce a very large percentage deviation for a relatively very small difference between model estimates). . Gas-phase photochemistry verified; simulations of photochemistry in an aerosol chamber exposed to natural light, where deviation is defined in Eq. 3. α−pinene ozonolysis is simulated in both plots, with α−pinene and O3 given the same initial concentrations of 21.1 ppb, and initial NOx concentration in (a) 9.8 ppb and in (b) 0 ppb. For both simulations, the environmental chamber is transparent and exposed to daylight without cloud interference, with dawn at approximately 4:00 hours. The particle-phase and vapour losses to walls are turned off in PyCHAM to be consistent with the AtChem2 model. Whilst Fig. 3 indicates that PyCHAM performs well for components with both relatively short (e.g. OH) and long (e.g. α−pinene) lifetimes, it is necessary to ascertain that agreement is gained through the correct mechanism. The chemical ten-240 dencies of formaldehyde were tracked in both PyCHAM and AtChem2 for the α−pinene ozonolysis simulations used for   in the legend, and the deviation is that from results for a time interval of 6x10 1 s. 52 s using a 2.5 GHz Intel Core i5 processor, with factor increases of 7 and 53 for 6x10 2 and 6x10 1 s resolutions, respectively.

255
Users should conduct a similar test if their chemical scheme or environmental conditions vary significantly from those here.
6 Gas-wall partitioning The partitioning of gases to the chamber wall is often termed wall loss as the net movement is from the gas phase to the wall (for an initially clean chamber wall). Traditionally this process has been viewed as an inconvenience since chamber results often depend on the concentration of gas-and particle-phases of certain components, whilst the fraction of these components 260 lost to walls is poorly constrained. Several studies have focussed on partitioning to Teflon walls, which are frequently employed (Matsunaga and Ziemann, 2010;Zhang et al., 2015;Zhao et al., 2018), however, the process remains poorly modelled across the wide range of chamber materials, relative humidities, gas-phase loading, component volatilities and activity coefficients present in chamber experiments (e.g. Krechmer et al., 2017;Stefenelli et al., 2018). It is therefore preferable to allow the user to fit vapour losses to walls through the tuning of two wall loss parameters, one primarily determining equilibrium, called the 265 effective wall mass concentration (C w ), and one determining rate of partitioning, the mass transfer coefficient (k w ). These influence gas-wall partitioning through an equation of the same framework as gas-particle partitioning (which is described in Section 7 and in Zaveri et al. (2008)):

270
where p 0 i is the liquid (sub-cooled if necessary) saturation vapour pressure of component i and γ i is its activity coefficient on the wall. Following the conclusions of Matsunaga and Ziemann (2010) and Zhang et al. (2015), k w represents factors such as gas-and wall-phase diffusion, turbulence, accommodation coefficient and the chamber surface area to volume ratio, whilst C w reflects the adsorbing and/or absorbing properties of the wall, including effects of relative humidity, surface area, diffusivity and porosity. We recommend the iterative fitting of k w and C w to observations through minimising observation-model residuals.

275
C w in PyCHAM does not vary with the mass transferred to the wall, which is consistent with the findings of Matsunaga and Ziemann (2010) and Zhang et al. (2015) that indicate the effective mass concentration of the wall is much larger than the mass concentration of transferred material.

Tuning gas-wall partitioning parameters
Next we illustrate the sensitivity to k w and C w in Eq. 4. The same simulation setup described above for Fig. 3 was used 280 though with α−pinene replaced by isoprene with a concentration at experiment start of 63.4 ppb. Seed particles comprised of ammonium sulphate with mean diameter 0.5 µm and number concentration 6x10 2 cm −3 were introduced at experiment start. Pure component liquid saturation vapour pressures were estimated by the Nannoolal et al. (2008) method and activity coefficients for all components were assumed to be unity. To begin, both k w and C w were set sufficiently low to effectively eliminate gas-wall partitioning. Second, C w was set equal to the mass concentration of seed particle (70 µg m −3 ) and k w raised 285 to 1x10 −1 s −1 at which a notable decrease in [SOA] was observed. Third, k w was held whilst C w was raised three orders of magnitude greater than the seed mass concentration. Fourth, C w was held at 70 µg m −3 and k w was raised by three orders of magnitude. The effect on SOA mass concentration is given in Fig. 6 and demonstrates that at sufficiently large values of C w , SOA production can be effectively suppressed through competitive uptake of vapours to chamber walls. However, for a given C w , there is a limit on suppression of SOA formation due to k w increase as it affects only the rate of partitioning with walls 290 rather than the condensable fraction.
To guide constraint for wall loss parameters, we follow the example of Matsunaga and Ziemann (2010), with a control experiment comprising a single semi-volatile component introduced to the chamber at the start of the simulation at 50 ppb. 2methylglyceric acid is selected as it has an estimated particle mass concentration saturation vapour pressure (C * ) of 1.15x10 2 µg m −3 at 298.15 K (the simulation temperature) and is an observed oxidation product of isoprene (Surratt et al., 2006). No 295 other components or particles are introduced. With regards to designing a control experiment for tuning C w and k w , the results shown in Fig. 6b demonstrate that a component with a C * close to the C w value has large sensitivity to C w , thereby allowing greatest ease of tuning. Note, that this sensitivity can be utilised through varying chamber temperature (and therefore the C * of a component), or through using a component with different volatility. Furthermore, to discern the effect of k w a component with substantial partitioning to walls is required. When quantifying k w it is worthwhile considering the required precision, 300 because as Fig. 6a demonstrates, above a certain value, no further effect on SOA concentration results. Particle number size concentration can change as a direct consequence of three processes modelled by PyCHAM: gasparticle partitioning, coagulation and nucleation. Whilst coagulation and nucleation are discussed below, readers are referred to Zaveri et al. (2008) for a thorough explanation of gas-particle partitioning. The transition regime correction factor required for the partitioning estimation in PyCHAM is from Fuchs and Sutugin (1971). Furthermore, the unit test test_kimt_calc is available  (2005)).
In this section we focus on redistribution of particles between size bins following a change in size due predominately to gas-particle partitioning. PyCHAM does this by adopting the moving-centre structure (Jacobson, 2005). The moving-centre approach has the advantage of minimal numerical diffusion and the ability to accommodate populations of particles of varying 315 modes (e.g. a nucleation event in the presence of pre-existing particles). However, it suffers from loss of accuracy due to averaging of particles originally from different size bins that have grown, shrunk or coagulated to a given size bin (Zhang et al., 1999). In contrast the full-moving structure does not average particles of different size bin together, and can therefore exactly model certain chamber scenarios, one of which we use below to verify the moving-centre approach against. Full-moving is not used in PyCHAM however, as it lacks the general applicability of moving-centre, since it cannot account for additions to the 320 particle population after experiment start (such as through injection of seed or nucleation).
To maintain stability in numerical solutions and improve model accuracy, a condition inside the moving-centre module of PyCHAM iteratively reduces the boundary condition time interval if particles in a size bin change volume sufficiently to be allocated to a size bin beyond the adjacent one, or if particles have an unrealistic negative volume (possibly due to evaporation).
Here we assess the moving-centre method through analysis of output during two periods of relatively substantial (and there-325 fore testing) condensational growth and compare to benchmark simulations. The simulations also illustrate two further means of component influx to chambers using PyCHAM in addition to the simulations above where components were introduced with an initial pulse. In the first case a constant flux of sulphuric acid is added to a chamber with seed aerosol typical of hazy conditions following the benchmark simulation of Zhang et al. (1999). For consistency with the benchmark, gas and particle partitioning to walls was turned off and sulphuric acid was assumed to be non-volatile. The analysis section of Zhang et al.

330
(1999) notes that to resolve the growth of smallest particles in this scenario, spatial resolution must be at least 100 size bins, therefore we use this value and set the update time interval to 90 s for a total 12 hour simulation.
The exact solution to this condensational growth problem is given by the full-moving structure in Fig. 7a (taken from Fig.   3 of Zhang et al. (1999)). When PyCHAM is compared against the exact solution, the tri-modal distribution is present with mean values at the correct particle size though with some disagreement in the peak height and spread. The degree of agreement 335 is significantly better than for the 13 size bin moving-centre simulation presented in Zhang et al. (1999) and indicates that PyCHAM is operating as intended. Our results in Fig. 7a are a two-point moving average which is often necessary for the moving-centre structure because its requirement that all particles in a size bin be transferred to the adjacent bin means that some bins will intermittently have zero particles. Another case of relatively intense vapour-particle partitioning is provided by the example of cloud condensation nuclei 340 experiencing varying degrees of water vapour supersaturation. Chamber experiments may involve injections of a component at specific times and the model variables input file can accommodate such a scenario. Making use of this function we reproduce the benchmark simulation of Jacobson (2005) (Fig. 13.8) where relative humidity is increased to 100.002 % every minute (including at simulation start) for nine minutes, with results analysed after ten minutes. Seed particles are assumed non-volatile and wall interactions are turned off. The parameters: temperature, seed component dissociation constant, molecular weight and 345 density are not disclosed by the reference simulation, therefore we set these as: 318.15 K, 1.0, 200 g mol −1 and 1 g cm −3 , respectively. The comparison between the Jacobson (2005) result in Fig. 7b and PyCHAM certainly shows agreement in the main feature of this simulation, which is the initially larger particles out competing smaller particles for water condensation to grow to water droplet size (D p > 10 µm). It should be noted that this is a very much more stringent test of the representation of partitioning than is ever intended for PyCHAM, which will not generally be used for the huge mass flux of condensing 350 material experienced under water supersaturated conditions. Nevertheless, the PyCHAM result gives reasonable agreement considering that key parameters (such as seed component dissociation) may vary between simulations and taken together with Fig. 7a verifies the operation of gas-particle partitioning and the moving-centre structure. turbulent inertial motion, turbulent shear and Van der Waals collision were taken from Jacobson (2005). The unit test test_coag produces a plot of coagulation kernels that can be compared to Fig. 15.7 of Jacobson (2005) to verify accuracy. Once the coagulation kernel for each pair of particle size bins (β) has been found, the combinations of size bins (denoted j and z) whose coagulation produces a particle of size bin k at time t are identified: where j max is the largest size bin that can undergo coagulation to produce a particle in size bin k. This equation is not 365 mass-conserving. However, mass transfer between size bins is estimated using the number fraction of particles coagulating.
For example, for a size bin k, the gain in molecular concentration of component i due to coagulation between size bins j and z (according to Eq. 6 and Eq. 7) is: and the loss of molecular concentration from size bin k due to coagulation is: Equations 7-9 imply that coagulation directly influences the particle number-and mass-size distributions. We now asses the sensitivity of the number size distribution and mass conservation to temporal (represented by the time interval for updating coagulation) and spatial (represented by number of size bins) resolution. A relatively complex initial distribution with four number modes is taken from ambient observations at Claremont, California on August 27, 1987 (Jacobson, 2005) and assumed 375 to comprise non-volatile material. Results are presented for a six hour simulation in Fig. 8 where particle wall loss was turned off to allow clearer assessment of the coagulation sensitivity. In the top row of Fig. 8 no gas-phase chemistry was allowed, whilst in the bottom row, a single chemical reaction with reaction rate 5.6x10 −17 molec −1 s −1 between α−pinene and O 3 (both with initial concentrations 100 ppb) was modelled to produce a single low volatility product with saturation vapour Figure 8. Sensitivity of the coagulation process to changes in temporal resolution (given in the legend) and spatial resolution (given in column titles). In the top row no chemistry occurred whilst in the bottom row a semi-volatile species was produced, as detailed in the main text. Results are for the end of a simulated six hour experiment. The ∆nv temporal resolution value given in the inset text is the percentage change in total non-volatile particle-phase material from the start to finish of the experiment, demonstrating mass conservation in the model. pressure of 1x10 −10 Pa, whilst gas-wall partitioning was turned off. For the chemistry case, approximately 500 µg m −3 of 380 secondary material was formed, compared to 90 µg m −3 of seed material. Columns in Fig. 8 are distinguished by the size bin resolution as presented in the column titles, and within each plot temporal resolution is varied.
The inset text of Fig. 8 (∆nv) gives the fractional change in non-volatile material from the start to end (six hours) of the simulation for the three temporal resolutions. It is clear that the coagulation equations introduce negligible error to mass conservation. Two features are present in the top row (no chemistry) of Fig. 8: first, that in terms of number concentration, 385 coagulation overwhelmingly affects the number concentration of smaller particles -note that such particles are sufficiently small in volume that they may coagulate with a larger particle without causing it to grow a size bin; second, that only for the smallest particles (below a diameter of 3x10 −2 µm in this case) is a sensitivity to temporal resolution clear across all size bin resolutions. For the no chemistry case there is demonstrable coupling of spatial and temporal resolution, with an increase in the former indicating greater sensitivity to the latter. However, the resolution considerations change substantially when we 390 consider the case with gas-particle partitioning (bottom row of Fig. 8). In this instance, the effect of partitioning dominates the change in number-size distribution and no sensitivity of coagulation to spatial or temporal resolution is discernible. We recommend users consider these examples in addition to the nature of their simulation and objective when deciding whether temporal or spatial resolution will significantly impact results.
9 Particle deposition to walls 395 As with gas-wall partitioning, the loss of particles to chamber walls can significantly invalidate chamber results if unaccounted for and has been detailed in previous publications (Crump and Seinfeld, 1981;McMurry and Rader, 1985;Nah et al., 2017;Wang et al., 2018). During control experiments the deposition rate of particles to walls can be inferred through observations of the rate of decay of particles of varying size (with coagulation accounted for) (Charan et al., 2019). Several studies have published results from such experiments (McMurry and Rader, 1985;Wang et al., 2018), including a relatively large dataset 400 from the EUROCHAMP2020 project (Oliveri, 2018). Comparison of inferred wall loss rates indicate that diffusion and settling enhance the loss rates of relatively small and relatively large particles, respectively (Crump and Seinfeld, 1981), however the absolute values and size-dependent gradient of the loss rates vary significantly between control experiments. Even for a given chamber, significant variations appear with changes to relative humidity, disturbance to walls due to air conditioning, and, for teflon chambers, with time since the chamber walls experienced frictional force to create electrostatic charge (Wang 405 et al., 2018). Currently no method is available to measure the required inputs that a particle deposition model would need to satisfactorily reproduce observations across all chambers and conditions, therefore in PyCHAM users have three options to estimate particle wall deposition. Here we describe the options and provide examples of their use.

Users select wall loss treatment with the McMurry_flag option in the model variables input file. The default (if left empty)
is no loss of particles to wall, which can be used for estimating wall loss corrected values such as aerosol yield. If set to one, 410 the model of McMurry and Rader (1985) is used, which is based on the particle deposition model of Crump and Seinfeld (1981) but with electrostatic effects. Studies have found the Crump and Seinfeld (1981) and McMurry and Rader (1985) approach to reproduce measured particle wall losses well (Chen et al., 1992;Kim et al., 2001). Selecting McMurry and Rader (1985) requires the user to also input the chamber surface area, the average charge per particle and the average electric field inside the chamber, where the latter two may be set to zero for nullifying electrostatic effects. With the test_wallloss module 415 users can confirm that PyCHAM accurately reproduces Fig. 2 of McMurry and Rader (1985), as shown here in Fig. 9, which demonstrates the effect of changing the charge number per particle.
If user sets the McMurry_flag option to zero then a customised particle deposition rate dependence on particle size is available. This option allows application of known or best estimate deposition rates (β) to the model, as recommended by Wang et al. (2018). Four further inputs are required for this option: the particle diameter at which the inflection in deposition 420 rates occurs (D p,f lec ) (where the inflection point marks a change in dependance of deposition rate with particle size), the rate of particle deposition to wall at the inflection point (β f lec ), and the gradients of the deposition rate with respect to particle diameter before (∇ pre ) and after the inflection (∇ pro ), where a linear dependence in log-log space is assumed, consistent with observations (Charan et al., 2019). The equations for deposition rate in this instance are given in Eq. 10, and example dependencies of rate with particle size provided by Fig. 9.

425
18 https://doi.org/10.5194/gmd-2020-234 Preprint. Discussion started: 20 August 2020 c Author(s) 2020. CC BY 4.0 License. Figure 9. Example dependencies of the particle deposition to wall rate using the model of McMurry and Rader (1985) in the solid lines, where the charge per particle is given by n and other inputs given by inset text (R is spherical-equivalent chamber radius, E is the average electric field in the chamber and ke is the coefficient of eddy diffusion). The dashed lines demonstrate the observation-based deposition rate utility of PyCHAM given in Eq. 10, with inputs at the top of the plot.

430
The simulation of nucleation to produce newly-formed suspended particles is one of the most active areas of ongoing atmospheric research and many important advances in observing the nucleation process have been, and will continue to be, made through appropriate measurements in chamber experiments and their interpretation (Dada et al., 2020). PyCHAM is not intended to interpret and examine chamber experiments designed to resolve the mechanisms involved in molecular clustering, Figure 10. Effect of varying the nucleation parameters on simulated particle number concentration when considering only nucleation.
nucleation and early growth in particle formation and there are tools much better suited to these processes. However, the use of 435 PyCHAM in simulating chamber processes in the presence of new particle formation necessitates a phenomenological accommodation of the process. Users are therefore able to provide parameters to a Gompertz function for cumulative new particle number, allowing them to fit to observed number size distributions without inferring mechanistic insight: where P 1 is the number concentration of new particles after time t that enter the smallest size bin, and nuc vn are the user-440 defined parameters. The resulting function forms an asymmetrical sigmoidal curve with time, whilst the parameters allow the amplitude (nuc v1 ), onset (nuc v2 ), and duration (nuc v3 ) of the curve to be adjusted, as shown in Fig. 10. The Gompertz function provides a sigmoidal form with faster increase in new particle number prior to peak rate of production than after (Fig. 10).

445
As with gas-wall partitioning parameters, nucleation parameters should be fitted to measurements by minimising modelobservation residuals. For this process the total particle number concentration may be used, however, the greater amount of data in number size distributions introduces stronger constraint, making it the preferred observation for fitting.
With the shape, size, composition and growth mechanism of the clusters that act as the nucleus of particles subject to ongoing research, in PyCHAM default properties are currently assigned, with a view to advance representation as understanding to enable simulations of coupled photochemistry and aerosol microphysics in seeded and unseeded experiments. However, a more rigorous mechanistic representation of nucleation and early growth should be readily accommodated and will be required 455 before PyCHAM is suitable for investigating new particle formation.

Sensitivity to temporal and spatial resolution
In PyCHAM, temporal resolution is represented by the time interval for updating ODE constants, as described in Section. 3 and spatial resolution is determined by the number of size bins. Whilst both resolutions can be decreased to decrease the time required for simulation, it can also introduce inaccuracies because PyCHAM processes are sensitive to changes to number 460 size distributions (updated after each time interval) and particle size. Although it is beyond the scope of this paper to assess resolution sensitivity across all possible PyCHAM parameter space, in this section we compare the divergence of outputs from simulations with decreasing temporal and spatial resolution against a high resolution reference for extremes of the relevant parameter space: seeded experiments with no gas-particle partitioning and both seeded and nucleation experiments with relatively large condensational growth of particles. As in Section 6, two-methylglyceric acid is used in the simulations with partitioning 465 as its vapour pressure at simulation temperature (298 K) makes it semi-volatile. Results here determine the recommended temporal and spatial resolution, provide a useful illustration of sensitivity and may help users perform sensitivity tests for their individual model inputs.
For the simulations without partitioning, the effect of resolution on particle number size distribution and total number concentration is considered, whilst for the partitioning simulations, concentration of secondary material is also relevant. To 470 allow comparison of low resolution simulations with the high simulation reference, output from the former is interpolated to the resolution of the latter. Divergence between a low resolution simulation (g) and the high resolution reference is represented by a single absolute percentage deviation (σ g ). For number size distribution, divergence is averaged over size bins containing particles and time steps, with Y of the former and Z of the latter:

475
where n g,ti,k is the particle number concentration at time step t i in size bin k from a lower resolution simulation andn ti,k is the output from the reference maximum resolution simulation. Where the two agree exactly the contribution to σ is zero, and where one output is zero and the other is greater, σ is at a maximum of 100. The term ∨(n g,ti,k ,n ti,k ) means the greater of n g,ti,k andn ti,k is used as denominator. Where the outputs from the two resolutions are similar this choice of denominator makes negligible difference, however, where one is much greater than the other it limits the divergence to a helpful (for For total number concentration and total secondary material concentration, divergence is calculated as the percentage deviation averaged over time steps: where N represents either total number concentration or total secondary material concentration.

485
For simulations assessing sensitivity to temporal resolution, 128 logarithmically spaced size bins are used, for which Fig. 8 indicates no limitation to accuracy due to spatial resolution. For seeded simulations, we use the same initial number size distribution as in Fig. 8, as this gives a relatively broad range of particle sizes, which is necessary to fully appreciate the sizedependent effects of coagulation, particle loss to wall and nucleation. All simulations were run for 24 hours and the reference simulation had an update time interval of 6 s.

490
Results for temporal resolution sensitivity are shown across three plots. The first, given in Fig. 11a represents the no partitioning case, with sensitivity assessed for two setups: only coagulation, and both coagulation and wall loss turned on. Coagulation proceeds as described in Section 8, whilst wall loss is described in Section 9, with the following inputs to recreate a size-dependent wall loss profile similar to n=3 in Fig. 9: D p,f lec = 1.0 µm, β f lec = 1.0x10 −4 s −1 , ∇ pre = ∇ pro = 1.5. Consequently, the particle loss to wall is relatively large and the sensitivity results are conservative. Fig. 11a indicates that under 495 this scenario, particle loss to walls considerably increases sensitivity to temporal resolution compared to the only coagulation case, with average deviations of 10 % for both total particle number concentration and number size distribution occurring at resolutions around two orders of magnitude finer than for coagulation alone. Fig. 11 show orders of magnitude variation in simulation times with changing resolutions. All simulations in this section were using a 2.5 GHz Intel Core i5 processor.

500
For Fig. 11b, two-methylglyceric acid is introduced at a rate of 1.0x10 −2 ppb s −1 and increases sensitivity to temporal resolution compared to Fig. 11a. A resolution of around 60 s is required to attain divergences less than 10 % for total number concentrations and secondary material, whilst for number size distribution, not even the lowest temporal resolution of 12 s can produce average divergence below 10 % when compared against the reference case of 6 s. This reflects the steep gradients in particle size-number space that are generated during intense condensational growth periods (e.g. Fig. 7), since small changes 505 to the update time interval can vary the size bins that particles concentrate in. This effect is even further pronounced when an extremely low volatility organic component is injected at the simulation start at 1 ppb to act as a nucleating agent in an unseeded simulation, with results given in Fig. 11c. We use the nucleation parameters for Eq. 11 of: nuc v1 = 1x10 4 , nuc v2 = −1x10 1 and nuc v3 = 1x10 2 , for a relatively rapid nucleation period (lasting only ten minutes) and therefore conservative assessment of sensitivity. Fig. 11c indicates that although particles do not grow to the same size bin(s) as the reference case, 510 the total concentration of number and secondary material diverges from the reference case by around 10 % for an update time interval of 60 s. Given the conservative nature of these simulations we therefore recommend a maximum update time interval of 60 s. However, with the coagulation case effectively representing zero wall loss and showing considerably less divergence, if users can demonstrate relatively low particle wall loss, a coarser resolution could be applied. Figure 11. Sensitivity of number size distribution (NSD), total particle number concentration (N ) and total concentration of secondary material ([SOA]) to temporal resolution (a-c) and to spatial resolution (d-f). The effects of coagulation alone (coag.) and combined with particle loss to wall (coag. & wall) are probed. In (a) and (d) no gas-particle partitioning is allowed, whilst in (b),(c), (e) and (f) it is, as twomethylglyceric acid is continuously injected at a rate of 1.0x10 −2 ppb s −1 . Also in (c) and (f) an extremely low volatility organic component is present at simulation start at 1 ppb, and set as the nucleating component for an unseeded simulation. The legend for lines in all plots is given in the upper right. Contours represent the time taken for simulation.
In Fig. 11d-f, the same simulation scenarios are applied as for temporal resolution sensitivity, but now we investigate spatial 515 resolution sensitivity. Using a fixed temporal resolution of 60 s, results show the divergence of 8, 32 and 64 size bins compared against results for 128 size bins (all logarithmically spaced). In Fig. 11d, when coagulation alone is effective, reasonable agreement is seen in total number concentrations across size bin resolutions. However, when wall loss is also considered, the relatively high loss rate of small particles to the wall leads to a strong dependence of divergence on resolution. The number size distribution divergence is poor across all scenarios, indicating that if this is an important output for users (e.g. when fitting 520 nucleation parameters, comparison against number size distribution is very useful), users should employ 128 size bins. Results for the partitioning cases in Fig. 11e and f show that without wall loss, total number concentration and secondary material concentration gives reasonable agreement of 10 % or less when 8 size bins are used. However, when partitioning is active, 32 size bins is the minimum resolution for divergence of approximately 10 % or less. Consequently, whilst recognising substantial differences between scenarios and user requirements, we recommend a size bin number of 32.

525
Whilst the contours in Fig. 11 provide the simulation time taken using a 1 reaction chemical scheme (ensuring that twomethylglyceric acid is recognised), Table 1 demonstrates simulation times using the α−pinene ozonolysis scheme of the MCM, which comprises approximately 1x10 3 reactions. Relevant combinations of spatial and temporal resolution are provided for a 6 hour unseeded experiment with α−pinene and ozone introduced at the start to generate a nucleation episode. The PyCHAM (CHemistry with Aerosol Microphysics in Python) software for aerosol chambers has been described. Its open source repository is given in Section 2. PyCHAM has been designed for optimal ease of use (from online access to output) whilst being broadly able to address scientific problems of current relevance across a range of aerosol chamber and experimental configurations (Section 2). We have provided a model output for the dark oxidation of limonene to illustrate the coupling of modelled processes: gas-phase chemistry, gas-particle partitioning, gas-wall partitioning, redistribution of particles between 535 size bins, particle loss to wall, coagulation and nucleation (Sections 2 and 3).
The steps to run a simulation using the software's GUI were described in Section 3 and the methods for estimating or setting component properties explained in Section 4. The setting up and solution of gas-phase photochemical reactions is detailed in Section 5, including comparison against the AtChem2 model (Sommariva et al., 2018) for verification and illustration of the effect of varying temporal resolution on model output for a system subject to varying natural light intensity.

540
For gas-wall partitioning this paper details (Section 6) a parameterisation that aims to satisfy the breadth of chamber characteristics and recommends a method for tuning to observations. In Section 7, gas-particle partitioning and the moving-centre structure for redistributing particles between size bins was introduced and assessed against benchmark simulations.
Coagulation was detailed in Section 8 and shown to introduce negligible loss of mass conservation. With Section 9, the three options for treating particle losses to walls were detailed and the resulting deposition rates as a function of particle diameter 545 were exemplified, including assessment against the benchmark of McMurry and Rader (1985). Similar to gas-wall partitioning, nucleation in PyCHAM is treated with a parameterisation that aims to optimise model versatility, with examples of parameter effects provided (Section 10).
In Section 11 the sensitivity of key outputs to temporal and spatial resolution were illustrated and informed our recommendations of a minimum update time interval of 60 s and a minimum of 32 size bins. We also show a high sensitivity of model 550 accuracy to the rate of particle wall loss and note that users could use lower resolutions if wall loss is lower than used in our tests.
Papers in preparation demonstrate further the utility of PyCHAM and its evaluation when assessed against observations.
These papers include both phenomenological and mechanistic approaches to coupled photochemistry and aerosol microphysics, both of which the model readily accommodates. new_partr Radius of newly nucleated particles (cm), if empty defaults to 2.0e-7 cm. inflectDp The particle diameter (m) at the inflection point of the size-dependent wall deposition rate.

Grad_pre_inflect
Negative log10 of the gradient of particle wall deposition rate against the log10 of particle diameter before inflection (/s). For example, for the rate to decrease by an order of magnitude every order of magnitude increase in particle diameter, set to 1.

Grad_post_inflect
Log10 of the gradient of particle wall deposition rate against the log10 of particle diameter after inflection (/s). For example, for the rate to increase by an order of magnitude for every order of magnitude increase in particle diameter, set to 1.

Rate_at_inflect
Particle deposition rate to wall at the inflection point for size-dependent particle loss to walls For example, for a two size bin simulation with 10 and 5 particles/cc in the first and second size bin respectively introduced at time 0 s, and later at time 120 s seed particles of concentration 6 and 0 particles/cc in the first and second size bin respectively are introduced, input is: pconct = 0; 120 (and the number_size_bins variable above = 2).
pconc Either total particle concentration, in which case should be a scalar, or particle concentration per size bin, in which case length should equal number of particle size bins (# particles/cc (air)). If an array of numbers, then separate numbers by a comma. If a scalar, the particles will be spread across size bins based on the values in the std and mean_rad variables below. To turn off particle considerations leave empty. If seed aerosol introduced at multiple times during the simulation, separate times using a semicolon. For example, for a two size bin simulation with 10 and 5 particles/cc in the first and second size bin respectively introduced at time 0 s, and later at time 120 s seed particles of concentration 6 and 0 particles/cc in the first and second size bin respectively are introduced, the input is: pconc = 10, 5; 6, 0 (and the number_size_bins variable above = 2).

Input Name
Description seed_name Name of component comprising the seed particles, can either be core for a component not present in the chemical scheme file, a name from this file, or H2O for water, note no quotation marks needed seed_mw Molecular weight of seed component (g/mol), if empty defaults to that of ammonium sulphate -132.14 g/mol seed_dens Density of seed material (g/cc), defaults to 1.0 g/cc if left empty mean_rad Mean radius of particles (um), defaults to a flag that tells software to estimate mean radius from the particle size bin radius bounds given by lower_part_size and upper_part_size variables above. If more than one size bin the default is the mid-point of each. If the lognormal size distribution is being found (using the std input below), mean_rad should be a scalar representing the mean radius of the lognormal size distribution. If seed particles are introduced at more than one time, then mean_rad for the different times should be separated by a semicolon. For example, if seed particle with a mean_rad of 1.0e-2 um introduced at start and with mean_rad of 1.0e-1 um introduced after 120 s, the input is: mean_rad = 1.0e-2; 1.0e-1 and the pconct input is pconct = 0; 120.