Awesome

Welcome to the PyCHAM (CHemistry with Aerosol Microphysics in Python Box Model) software for modelling of indoor environments, including aerosol chambers. Funding has been provided by the EUROCHAMP-2020 research project and the National Centre for Atmospheric Science (NCAS). Please open an issue on the GitHub repository or contact Simon O'Meara (simon.omeara@manchester.ac.uk) with any issues, comments or suggestions.

PyCHAM is an open-source computer code (written in Python) for simulating aerosol chambers. It is supplied under the GNU General Public License v3.0. The license document is provided with the software (LICENSE) and contains information around modification and conveyancing.

PyCHAM was written by Simon P. O'Meara, Shuxuan Xu and Ademipo Onanuga.

Table of Content

1. Documentation
2. Installation
3. Running
4. Testing
5. Inputs
6. Outputs
7. Photochemistry
8. Gas-particle Partitioning
9. Numerical Considerations
10. Quick Plotting Tab
11. Flow Mode
12. Indoor Air Quality Modelling
13. Major Version Updates
14. Ambient Modelling
15. Frequently Asked Questions
16. Acknowledgements

Documentation

The README file you are now viewing serves as the PyCHAM manual, explaining how to setup the software and use it. As an additional resource, we also provide an introductory video.

The article published in the Journal for Open Source Software explains the underlying mechanisms of PyCHAM and its purpose. This article was reviewed using v0.2.4 of PyCHAM. Additionally, the article published in Geophysical Model Development provides a detailed introduction of PyCHAM and its use. This article was reviewed using v2.1.1 of PyCHAM The DOI for all PyCHAM releases is: 10.5281/zenodo.3752676.

Version numbers of PyCHAM try to adhere to the semantic described by semver, which is: MAJOR.MINOR.PATCH.

Installation

There are two options for installing, via conda and via pip. The pip method takes longer as the openbabel package has to be installed separately. The instructions below for the pip method currently apply only to linux and macOS, whilst the conda instructions apply to windows, linux and macOS.

Install via conda

1. Download the PyCHAM repository from github.com/simonom/PyCHAM

2. Download and install the package manager Miniconda (Anaconda is also suitable but takes more memory and takes longer to install) using the following address and selecting the appropriate operating system version: https://docs.conda.io/en/latest/miniconda.html.

3. Ensure conda is operating correctly, the method varies between operating systems and is explained here

The following steps are at the command line:

4. cd into the directory where the PyCHAM package is stored.

5. Use the following command to install: conda env create -f PyCHAM_OSenv.yml -n PyCHAM, where OS is replaced by your operating system name (win (for Windows), lin (for Linux), mac (for macOS)).

6. Activate the environment with the command: conda activate PyCHAM

6a. For Windows only, now also install openbabel using the following at the command line: conda install openbabel. In addition, it has been noted that when following the above steps in Windows, openbabel does not always link to its data directory. You can test this by looking for open babel warning messages when running the PyCHAM examples (to run PyCHAM please see Running). A message such as 'the aromatic.txt file could not be found' indicates the data library is not set. To set, search for Control Panel in the Windows search bar then User Accounts and User Accounts again. Then select 'Change my environment variables' for User variables select New then create the new variable with Variable name 'BABEL_DATADIR' and with Variable value set as the path to the relevant open babel data directory. To identify your open babel data directory you should look for the open babel folder containing aromatic.txt at a location such as: C:\Users\your_account_name\.conda\pkgs\openbabel-2.4.1-py37_5\share\openbabel. Once the path is suppled to the Variable value box select OK to create the variable then open a new anaconda prompt to restart PyCHAM.

6b. On Mac (please skip this step if on Linux or Windows) it has been noted that when following the above steps and then running PyCHAM (to run PyCHAM please see Running), an error is given stating that Developer Tools have not been installed. If this error displays please follow the on-screen instructions to install Developer Tools.

7. Install is complete for Linux, Mac and Windows; to run PyCHAM please see Running.

Install via pip

1. Ensure that swig is installed on your system. For example, for macOS, a command at the command line like: brew install swig , and for linux a command at the command line like: sudo apt install swig

2. Ensure that eigen is available on your system. For example, for macOS, a command at the command line like: brew install eigen, and for linux, download the latest stable eigen release here, then unzip and move the unzipped folder into the /usr/local/include folder using a command at the command line like: sudo mv eigen-3.3.9 /usr/local/include

3. Create a virtual environment in a suitable location running with at least python3.6. For example, if your command line recognises python3.6 as python3.6, the command to make a virtual environment called 3env on macOS is: python3.6 -m venv 3env. Note that python3.6 in this example should be replaced with the appropriate command for recognising python on your machine (often just python).

4. Activate this environment. For example, for a virtual environment called 3env on macOS and linux the command at the command line is: source 3env/bin/activate. For example, for a virtual environment called 3env on Windows the command at the command line is: .\3env\Scripts\activate

5. openbabel must be installed separately to PyCHAM, begin by downloading and unzipping the tar file (file name containing .tar.) for the latest version of openbabel on github. Note that the unzipped version can be stored in the Downloads folder for the installation process.

The following steps are at the command line:

6. At the command line cd into the unzipped openbabel folder, for example: cd openbabel-3.1.1

7. Create a build directory: mkdir build

8. Change into the build directory: cd build

9. Using cmake, prepare the openbabel build files. This requires several specifications, the -DCMAKE_INSTALL_PREFIX specification should be the path to the site-packages folder of the virtual environment created above: cmake -DRUN_SWIG=ON -DPYTHON_BINDINGS=ON -DCMAKE_INSTALL_PREFIX=~/3env/lib/python3.6/site-packages ..

9a) If the above command causes a message at the command prompt that eigen cannot be found this can be fixed by stating its location, for example: cmake .. -DRUN_SWIG=ON -DPYTHON_BINDINGS=ON -DCMAKE_INSTALL_PREFIX=~/3env/lib/python3.6/site-packages -DEIGEN3_INCLUDE_DIR=/usr/local/Cellar/eigen/3.3.9/include/eigen3

10. Complete installation with: make install

11. Test that openbabel is functioning with: python

11a) Then inside the python interpreter use: import openbabel

11b) If this works fine (no error message) continue to step 12 If this returns: ModuleNotFoundError: No module named 'openbabel', then quit the python interpreter: quit()

11c) If error seen in step above, add the openbabel path to the python path, for example: export PYTHONPATH=~/3env/lib/python3.6/site-packages/lib/python3.6/site-packages/openbabel:$PYTHONPATH

12. Test that pybel is functioning with: python

12a) Then inside the python interpreter use: import pybel

12b) If no error message seen continue to next step. During testing this threw a relative import error which was corrected by changing the relevant line (given in the error message) in the pybel.py file (file location given in the error message) from "from . import openbabel as ob" to "import openbabel as ob"

13. Ensure pip and wheel up to date: pip install --upgrade pip wheel

14. Install PyCHAM and its dependencies in the virtual environment: python -m pip install --upgrade PyCHAM

Install is complete, to run PyCHAM please see Running.

Running

1. For model inputs, ensure you have: a .txt file chemical reaction scheme, a .xml file for converting the component names used in the chemical reaction scheme file to SMILE strings and a .txt file stating values of model variables (e.g. temperature) - see details about these three files below and note that example files are available in PyCHAM/input

2. Once Installation is complete and the appropriate environment has been activated (see Installation), use the command line to change into the top level directory PyCHAM (the directory above the PyCHAM main file).

3. There are two choices for starting up the programme. If you are new to PyCHAM, and/or have simulations that can be practically selected manually, then begin the programme from the command line to use PyCHAM via the graphical user interface: python PyCHAM. Alternatively, if you wish to automate simulation setup and run, then edit the module 'automated_setup_and_call.py' to your needs and call this from the command line: python automated_setup_and_call

4. The PyCHAM graphical user interface (GUI) should now display on your screen. Using the 'Simulate' tab, one can select the folder containing all input files using the 'Select Folder Containing Input Files' button. This will search the selected folder for the input files (chemical reaction scheme, xml and model variables). For the chemical scheme, files with filenames including 'chem' will be identified. For the xml, files with filenames including 'xml' will be identified. For the model variables, files with filenames including 'var' will be identified. Any identified files will then be displayed in the GUI (see below for details on the contents of the chemical scheme, xml and model variables input files).

5. To select any of the input files individually, one can use the corresponding GUI button.

6. The GUI will display the found inputs provided in the selected model variables file. For inputs not stated in this file, the displayed variables are default.

7. Problems with the input files will be displayed in the GUI - this functionality is under development, meaning that not all problems are currently captured.

8. Once the first simulation is ready (through selection of the desired combination of correct input files described above), the user chooses between a single simulation or adding to batch, with the latter allowing multiple simulations to be queued.

9a. If the user chooses a single simulation to run, a progress bar will show, which represents the time through the experiment as a fraction of the total experiment time.

9b. If the user chooses to add to batch, then further simulations can be chosen by repeating steps 4-7 above. When ready, the batch can be run with the start series of simulations button. The progress bar then represents individual experiments and the current simulation is shown in the GUI. Note that when adding to batch input files should be located in different folders (rather than changing the inputs inside a folder already selected for batch between adding to batch).

10. The 'Plot' tab allows multiple plotting options. The Standard Results Plot produces two sub-plots in one figure: one with the particle number distribution, secondary aerosol mass, and particle number concentration against time, and another plot that shows the gas-phase concentrations of specified components with time (the specified components are those with initial concentrations given in the model variables file).

11. The 'Quit' button will stop the programme. If it does not work, the ctrl+z key combination in the console window can cease operations safely. In both cases Python will release all memory associated with the simulation.

Testing

Unit tests for PyCHAM modules can be found in the PyCHAM/unit_tests folder. Call these tests from the home folder for PyCHAM, with: python test_module.py with module replaced by the name of the PyCHAM module to be tested. For some unit tests example inputs are required, the chemical scheme files for these are stored in unit_tests/input with file names beginning with test_ ..., therefore we recommend users do not use chemical schemes with the same naming convention to prevent confusion. Where required, model variables for unit tests either use the default values or those given in the unit test script and use the xml file provided in PyCHAM/input.

Continuous integration testing can be completed using the '.travis.yml' (home folder) and 'test_TravisCI.py' (unit_tests folder) files at the Travis CI website.

Various example run outputs are saved in the PyCHAM/output/ folder. To reproduce these, select from PyCHAM/input for the relevant example. The ind_AQ_ex and ambient_constrained_ex are further described in the Indoor Air Quality Modelling and Ambient Modelling sections, respectively (found below). Note that the example output may vary between releases so please check correspondence.

Inputs

Chemical Scheme file

The chemical scheme file states the reactions and their rate coefficients in the gas- and aqueous-phases.

An example chemical scheme file is given in the PyCHAM/input folder, called 'example_scheme.txt', which has been obtained from the Master Chemical Mechanism (MCM) website (KPP version) and modified.

Results are automatically saved in PyCHAM/output/name_of_chemical_scheme_file/name_given_in_model_variables_input_file_for_saving.

The unit tests described above save results with the prefix 'test_', therefore we recommend using a different convention for chemical scheme names to prevent confusion.

Markers are required to recognise different sections of the chemical scheme. The default markers are for the MCM KPP format, however, others can be specified using the chem_scheme_markers input in the model variables input file. A guide to chem_scheme_markers is given in the Model Variables .txt file section below. This includes how to distinguish between gas- and aqueous-phase reactions.

Reaction rate coefficients for chemical reactions and generic rate coefficients must adhere to the following rules: The expression for the rate coefficient can use Fortran type scientific notation or python type; acceptable math functions: EXP, exp, dsqrt, dlog, LOG, dabs, LOG10, numpy.exp, numpy.sqrt, numpy.log, numpy.abs, numpy.log10; rate coefficients may be functions of TEMP, RH, M, N2, O2 where TEMP is temperature (K), RH is relative humidity (0-1), M, N2 and O2 are the concentrations of third body, nitrogen and oxygen, respectively (# molecules/cm3 (air)).

Inside the chemical scheme file, the expression for the reaction rate coefficient of a chemical reaction and the reaction itself must be contained on the same line of the file, with some delimiter (described above with chem_scheme_markers) separating them.

Chemical Scheme .xml file

An example is given in the inputs folder (of the Github repository), called 'examples_xml.xml'. It has a two line header, the first states that the mechanism is beginning (<mechanism>) and the second states that the species definition is beginning (<species_defs>). The end of the species list must be marked (</species_defs>) and finally, the end of the mechanism must be marked (</mechanism>).

Beneath this, every component included in the reactions of the chemical scheme must have its SMILES string given. To add new components, use this three line example per new component:

<species species_number="s6058" species_name="O2">

<smiles>O=O</smiles>

</species>

Here the first line states that the species definition is beginning and gives a unique code, s6058 in this case, and its chemical scheme name, in this case O2. The second line provides the SMILES string, in this case O=O. The third line states that the definition is finished. For information on SMILES please see: SMILES website.

For Master Chemical Mechanism chemical schemes the associated xml file can be acquired from the MCM website.

Model Variables .txt file

An example is provided in the inputs folder (of the Github repository), called 'example_model_var.txt' , this can include the following variables separated by a return (so one line per variable), note that if a variable is irrelevant for your simulation, it can be omitted and will be replaced by the default. In addition, you can find more information around photochemistry and flow mode instruments in the corresponding section of this README file.

Outputs

Model results are saved to the folder specified by the user in the model variables file, using the model variable res_file_name (described above). PyCHAM automatically places this folder inside the PyCHAM/output/name of chemical scheme/ folder. Several files are stored in this output folder. The concentrations of components is stored in the file with name beginning 'concentrations_all_components_all_times'. This is a comma separated value (csv) file. If you would like suggestions for code to open and view results, please see below, along with the plotter_gp.py file and the retr_out.py files.

A minimum working example (for plotting the time profile of the gas-phase concentration of a given component) is (note that you may need to activate the PyCHAM environment in order to have the necessary packages available (numpy and matplotlib)):

# state path to output folder (for Windows Operating System use \\ to separate folders rather than /))
output_by_sim = 'path to your output folder'

# combine folder path with specific file name for component names
fname = str(output_by_sim + '/model_and_component_constants')

# open file containing component name
const_in = open(fname)
# create empty dictionary to hold component names
const = {}

# loop through lines of file containing component names
for line in const_in.readlines():
		
	dlist = [] # empty list to hold values
	for i in line.split(',')[1::]:
			
		if (str(line.split(',')[0]) == 'chem_scheme_names') or (str(line.split(',')[0]) == 'SMILES') or (str(line.split(',')[0]) == 'space_mode'):
			i = i.strip('\n')
			i = i.strip('[')
			i = i.strip(']')
			i = i.strip(' ')
			i = i.strip('\'')
			dlist.append(str(i))
			
	const[str(line.split(',')[0])] = dlist

# close file with component names
const_in.close()

# isolate component names from dictionary
comp_names = const['chem_scheme_names']

# withdraw times (s) -----------------
fname = str(output_by_sim+'/time')
# import numpy package
import numpy as np
# load times
t_array = np.loadtxt(fname, delimiter=',', skiprows=1)
timehr = t_array/3600.0 # convert from s to hr
# ------------------------------------------

# combine folder path with specific file name for component concentrations
fname = str(output_by_sim + '/concentrations_all_components_all_times_gas_particle_wall')
# load file, omitting headers
y = np.loadtxt(fname, delimiter=',', skiprows=1)

# state name of component you want to plot
comp_names_to_plot = 'APINENE'

# get index of this component
indx_plot = [comp_names.index(comp_names_to_plot.strip())]

# import packages for plotting
import matplotlib.pyplot as plt

# make plot
plt.plot(timehr, y[:, indx_plot])

# show plot
plt.show()

Photochemistry

Chemical schemes may include photochemical reactions where the rate of reaction is dependent on light intensity. Several of the model variables described here in the Model Variables .txt file section are relevant to correct modelling of photochemistry and these will be further detailed here.

The input variables light_status and light_time determine when the chamber is illuminated or dark.

The input variable act_flux_path states the actinic flux (photon/cm2/nm/s) as a function of wavelength (nm). For chambers with artificial light (lamps) it is necessary to supply this file so that PyCHAM knows the light intensity spectrum. PyCHAM will automatically interpolate the wavelengths and corresponding actinic fluxes given in act_flux_path to unit wavelength resolution (every 1 nm) to ensure correct integration of photolysis rate across the spectrum. Inside the photofiles folder are examples of act_flux_path (e.g. Example_act_flux.csv), including the file for Manchester Aerosol Chamber (MAC) (MAC_Actinic_Flux_Spectrum.csv). The required format is a comma separated value file with wavelength (nm) in the first column and the corresponding actinic flux (photon/cm2/nm/s) in the second column. No headers are allowed.

For chambers with natural light (open roof), users may also supply an act_flux_path representing the relevant solar light intensity spectrum. However, if natural light is present and the chemical scheme is derived from the Master Chemical Mechanism then PyCHAM will use the parameterisation of Hayman (1997), described in Saunders et al. (2003), to estimate the photolysis rates of the Master Chemical Mechanism. In this model setup users may also supply the day number of the year (# days) (DayOfYear model variabled) time of day (Greenwich Mean Time (GMT)/Coordinated Universal Time (UTC) in seconds (not hours:minutes:seconds)) that the experiment starts (daytime_start model variable), the latitude (lat model variable) (degrees) and longitude (lon model variable) (degrees). These inputs allow the solar zenith angle to be calculated according to the first chapter of Environmental UV Photobiology (1993): "The Atmosphere and UV-B Radiation at Ground Level" by S. Madronich (Environmental UV Photobiology, 1993). This setting of solar radiation and deriving the photolysis rates for MCM is the default setting when light_status is set to illuminated.

Photolysis rate is the product of actinic flux, component absorption cross-sections (wavelength dependent) and quantum yield (wavelength dependent) integrated over the relevant range of the light spectrum. By default PyCHAM assumes the Master Chemical Mechanism photolysis reaction rate coefficients require estimation. For this reason, the PyCHAM software comes with the component absorption cross-sections and quantum yields as recommended by the Master Chemical Mechanism v3.3.1 website: http://mcm.leeds.ac.uk/MCMv3.3.1/parameters/photolysis.htt.

For photolysis reactions in chemical schemes other than Master Chemical Mechanism, the user can supply their own file for component absorption cross-section and quantum yield (the values are to be contained in the same file). The file name should be stated in the photo_par_file model variable and the file should be stored in PyCHAM/photofiles. A short example is given in the photofiles folder, called example_inputs.txt. File must be of .txt format with the formatting: <br> J_n_axs <br> wv_m, axs_m <br> J_n_qy <br> wv_M, qy_m <br> J_end <br> where n is the photochemical reaction number, axs represents the absorption cross-section (cm2/molecule), wv is wavelength (nm), _m is the wavelength number, and qy represents quantum yield (fraction). J_end marks the end of the photolysis file.

If a lamp power is given in terms of watts (J (kg m^2/s^2)/s), and a spectrum per unit wavelength is provided then it is possible to convert to actinic flux in units of photons/s. To do this use the photon energy formula: E (J (kg m^2/s^2)) = h (J (kg m^2/s^2) s)*c (m/s)/lamda (m), where h is Planck's constant (6.626e-34 Js), c is the speed of light (3e8 m/s) and lambda is the wavelength (m), to get the energy of one photon. Then, divide the provided Watts (J/s) per unit wavelength value by the result to get the actinic flux (photons/s) (note that by definition the actinic flux corresponds to a unit wavelength and a unit area).

Gas-particle Partitioning

Whilst the numerical treatment of gas-particle partitioning is described in the Geophysical Model Development (GMD) paper, here the link between model variables and gas-particle partitioning is further described.

Setting the seed particle composition through the model variables seed_name and seedx has no automatic effect on the activity coefficient of components (and therefore on gas-particle partitioning via activity), since PyCHAM does not automatically estimate activity coefficients.

To estimate gas-particle partitiong, PyCHAM uses the difference in concentration of a component between the gas and particle phase, which depends on the mole fraction of that component in the particle phase. Particle-phase mole fractions are a function of the concentrations of all components present in a particle and, for seed particle components, depends on the product of their particle-phase concentration and their dissociation constants (seed_diss model variable). The dissociation constant is the number of ions an inorganic component dissociates into at infinite dilution per molecule of that component, e.g. for ammonium sulphate the constant is 3 and for sodium chloride it is 2. Gas-particle partitioning is also a function of the activity coefficient of a component with respect to particles. The activity coefficient is not estimated automatically in PyCHAM, but maybe set by the user through the act_user and act_comp model variables. This is currently limited to one value per component for all size bins throughout a whole simulation. More information on activity coefficients can be found on page 407 of Seinfeld and Pandis 2016.

The default model treatment for gas-particle partitioning of water is the same as for other components, including the ability of the user to set the activity coefficient for water. The user may additionally (for water only) use the drh_ft, erh_ft and H2O_hist model variables to set whether water is able to partition to particles or not at a given relative humidity and water partitioning history. If partitioning is not allowed, particles are modelled as completely dry (zero mole fraction of water), whilst if partitioning is allowed, water is allowed to partition as described above in this section and in the [Geophysical Model Development (GMD) paper] (https://doi.org/10.5194/gmd-14-675-2021). PyCHAM cannot automatically determine deliquescence or efflorescence relative humidities, nor activity coefficients of water.

Numerical Considerations

In this section aspects of PyCHAM affecting numerical stability and computation time speed-up are discussed.

Speed-up of computation time can be achieved through solving gas-particle partitioning of water separately to other processes (the other processes are gas-wall partitioning of water, partitioning of non-water components between gas-particle and gas-wall, and chemical reactions). For a system with ~1000 chemical reactions and 32 particle size bins a speed-up of a factor ~500 was seen when this separation was introduced. Separation is done by default but can be turned off by setting the ser_H2O model variable to 0. To the best of our knowledge, solving water gas-particle partitioning separately has negligible effect on integration estimates.

The Ordinary Differential Equation (ODE) solver package is solve_ivp from Scipy. For the integration of the vapour-particle partitioning of water problem the Radau integration method is used as testing indicates this gives least computation time. For integration of other processes (vapour-particle partitioning of non-water components and chemistry) problems, the backward differentiation formula (BDF) method is used as it is well suited to stiff problems.

The user can supply their own integration tolerances (int_tol) in the model variables file. By default PyCHAM uses tolerances that were found to suit the problems presented in the GMD software decription paper (cited above). However, non-stiff problems can be solved with less computation using higher integration tolerances, whilst stiffer problems may become unstable unless lower tolerances are used.

The ODE solver can fail to find a stable solution. Sometimes this is due to a problem being too stiff. If this occurs a message is printed to the PyCHAM GUI. The initial message describes that instability has been noted and a smaller integration time step is being tried to improve stability. Below a certain small time step the message changes and the simulation stops. It has been decided that a stable solution cannot be found. In this situation, to help the user identify the cause of instability, a file called ODE_solver_break_relevant_fluxes is produced which contains the fluxes of components with negative concentrations output by the ODE solver (where negative concentrations are physically unrealistic and therefore indicative of causing/resulting from integration instability). The user could identify small or large fluxes (relative to other fluxes in this file) to gain indication of the processes causing instability and therefore determine whether a change of inputs is possible or whether PyCHAM is unsuitable for the problem in question.

Stability of the ODE solver may be compromised by a change in gas-particle partitioning of water that is too great for the solver. We recommend that instantaneous changes to relative humidity therefore do not surpass 0.01. In the example input titled rh_change_example is an example of relativly large changes to relative humidity that PyCHAM can cope with, however much greater changes would cause instability.

If the user-supplied model variables includes a seed component that is volatile then it may evaporate. This can cause instability in the ODE solver. Because this scenario is (in our experimental experience) unlikely in reality, PyCHAM automatically makes seed components non-volatile. However, this can be over-ridden by specifying a pure component saturation vapour pressure for seed components in the user-supplied model variables.

Computation time speed-up and ODE solver stability may be gained by increasing the value of the user-supplied model variable z_prt_coeff. This variable determines for which size bins gas-particle partitioning of components (including water) is allowed. If the gas-particle partitioning coefficient of a particle size bin is sufficiently low that it comprises less than this fraction of the total gas-particle partitioning coefficient (sum of coefficients across all size bins), then for this size bin partitioning is treated as zero. If a size bin represents a relatively very small fraction of total partitioning the partitioning problem can become too stiff and prevent the ODE solver from concluding. To the best of our knowledge the default value of 1.e-9 for z_prt_coeff has negligible effect on model estimates. However, if users state a higher number, we recommend a sensitivity test for their simulation between this user-supplied number and one an order of magnitude higher to check whether a significant effect is seen - if not then the user-supplied number is shown to be conservative as it does not affect results - if it is then the user must decide whether the loss of accuracy is acceptable.

To reduce stiffness when gas-particle partitioning components reduce to very small concentrations (less than the fraction given by the z_prt_coeff model variable when excluding water and seed, but summing over gas and particle abundances), gas-particle partitioning of these species is now turned off by default in partit_var.py. Users can remove this assumption by setting z_prt_coeff model variable to 0. The fractions are recalculated after every integration step, so that any increases in abundance above the z_prt_coeff fraction mean that gas-particle partitioning restarts.

A change in boundary conditions, such as lighting, can cause concentrations to tend toward zero. In turn, this can cause negative concentrations of components. Therefore, limited negativity is allowed when checking for instability in the solution. This can be seen inside the ode_updater module below the call to the ode_solv module. An example input folder called neg_conc_example can invoke the allowable negative concentrations. Gas-phase (and possibly particle-phase) output for the components with MCM names: O, O1D and PROPALO show a tendency toward zero after 50 minutes of experiment, with fluctuations above and below zero.

Quick Plotting Tab

Particle mass concentration (whether total of all components or excluding certain components) is calculated based on the concentration of components and their molecular weight. At a given time, for a given component in a given particle size bin (note that the concentration of a component represents the sum over all particles present in a size bin) the formula is: ((# molecules/cm3)/(Avogadro's Number (# molecules/mol)))*(g/mol)*1.e12 = ug/m3. This equation can be extended over certain size bins and certain groupings of components to attain the desired particle-phase mass concentration.

Flow Mode

When preparing model variable inputs for an instrument in flow mode (e.g. a flow tube or a chamber in flow reactor mode), the dilution factor and continuous influx of component model variables can be used. To simulate removal of a constant fraction of the chamber's volume per second, set the dil_fac model variable accordingly. For example, if the residence time in the instrument is 10 seconds, then 0.1 of the volume is removed per second, so use dil_fac = 0.1. If components of interest (including potentially water) are injected to the instrument to replace the components lost through chamber air being extracted, then use the model variables: const_infl, const_infl_t and Cinfl to describe their continuous influx.

Indoor Air Quality Modelling

For simulating atmospheres inside a building, an example set of input files is provided at PyCHAM/input/ind_AQ_ex. When simulating buildings for first time, it is recommended to first read the notes above regarding Flow Mode. This is because ventilation of the building will extract a certain fraction of air per second, which will be replaced by outdoor air.

Other features that are relevant to indoor air simulation include: the wavelength-dependent attenuation of solar radiation (see the model variable trans_fac above); the difference rates of deposition of gases and vapours to surfaces (possibly due to differing reaction rates on surfaces) (see the vol_Comp model variable above); the emission of gases/vapours into the building either from outdoors, or from indoor surfaces can be dealt with through the const_infl model variable (described above), which has the option of specifying a path to a spreadsheet containing emissions into the gas phase (e.g. if the number of components or/and the number of time points is sufficiently great that it is more neatly contained in a spreadsheet).

PyCHAM has been tested against observations from homes with the purpose of verifying that it can include the key aspects of atmospheric science that are particularly important for the indoor environment: surface-gas-particle parititioning of organics (verified against [Lunderberg et al. (2020)] (https://dx.doi.org/10.1021/acs.est.0c00966)), ozone reaction on surfaces to generate oxidised organics (verified against [Morrison et al. (2022)] (https://pubs.rsc.org/en/content/articlelanding/2022/EM/D2EM00307D)), the effect of ventilation; the effect of indoor emissions on indoor aerosol (gas and particle), infiltration of outdoor pollution indoors, infiltration of indoor pollution outdoors.

The [Lunderberg et al. (2020)] (https://dx.doi.org/10.1021/acs.est.0c00966) study observes that components with volatilities comparable to C13-C22 alkanes show increased mixing ratios in the gas-phase+particle-phase with increased temperature, with decreased dependency as volatility decreases. They conclude that increased temperature decreases the fraction of these components on surfaces. Furthermore, components with volatilities comparable to C23-C31 alkanes show increased mixing ratios in the gas-phase+particle-phase with increased mass concentrations of particulate matter with diameter below 2.5 um. They conclude that a reservoir for these components on surfaces acts to modulate the particle-phase concentration. The PyCHAM inputs to generate the Lunderberg et al. (2020) results and the code to generate the PyCHAM results below are stored in PyCHAM/input/ind_AQ_ex/Lunderberg2020. Note that to run this code for your system, you will need to set the dir_path and ret_path variables to the path for your system. The plots below show results from PyCHAM simulations and reproduce the Lunderberg et al. (2020) observations.

Major Version Updates

Major version updates are updates that are not backward compatible, meaning that using the same model inputs in the updated version may give results that are different to those created using the same inputs in previous versions. Major version updates are represented by a change to the first number in the verison number. The version number represents MAJOR.MINOR.PATCH changes. In 2024, work with Yarê Baker (then of Forschungszentrum Jülich GmbH) recognised that PyCHAM version 4 and below had been using the incorrect calculation for the mean free path (mfp) of trace components in air inside the partit_var_prep module (which is then applied in the partit_var module to estimate Knudsen number for the transiton regime correction factor for the diffusion coefficient of components (required to estimate rate of gas-particle partitioning)). In PyCHAM version 4 and below equation 15.24 of Jacobson 2005 (Fundamentals of Atmospheric Modeling) was used, but this equation is intended for the mean free path of the main (relatively small) components of air, such as N2 and O2, through air, whereas in PyCHAM version 5 and above, Eq. 9.15 of Seinfeld and Pandis 2016 (Atmospheric Chemistry and Physics, 3rd edition, print ISBN 9781118947401) is used, which is better suited to the mean free path of relatively large trace components of air through air. The effect on resulting particulate matter is demonstrated below, where results from v4.8.0 (using the old mean free path calculation) are compared against results from v5.0.0 (using the new mean free path calculation) for the secondary particulate matter formed during the simulation described in Figure 1 of the PyCHAM GMD paper (2021). It should be noted that although the relative effect can reach a factor of around three, this relates to a relatively small absolute difference, due to a relatively small absolute mass concentration.

Ambient Modelling

Several considerations may need to be made when accurately simulating the ambient (outdoor) environment. The ambient_constrained_ex demonstrates these aspects, including use of the model variables: dil_fac, cont_infl, pconc, pconct and pcont model variables to represent possible mixing of box air with air from adjoining spaces, or direct emission into box air. Furthermore, obs_file demonstrates constraint against observations, and spin_up demonstrates a spin up of the gas-phase chemistry and gas-particle partitioning.

Frequently Asked Questions

Why does PyCHAM crash without an error message? This has been observed using the conda install on a Windows 10 operating system with a Intel(R) Core(TM) i7-8500y CPU @ 1.50Ghz processor with processor speed 1400 MHz. Checking the Event Viewer application on Windows, under the Windows Logs/Application tab, showed that libblas was crashing. To correct, the PyCHAM virtual environment was activated in the command line, then from the command line conda was used to uninstall libblas and its dependents, then conda was used to install libblas again followed by its dependents. This solved the issue.

How are seed component concentrations initialised? The concentration of seed particles is based on the seed particle properties supplied by the user in the model variable file. The molar volume (cm3/mol) of the seed component is calculated using: (g/mol)/(g/cm3) or (molar mass)/(component density), where density is estimated from the Girolami method of UManSysProp. The volume of seed particles per size bin is calculated from the number size distribution stated in the model variable file. Seed particle volumes per size bin are then divided by the molar volume of seed components to estimate the concentration of the seed components per size bin: # molecules/cm3 (seed components per size bin) = (cm3/cm3 (seed particle per size bin)/(cm3/mol) (seed components))*Avogadro's constant.

What if I need an automated method for generating the xml file? A suggested code is contained in SMILES_generator.py - which allows xml file generation for the PRAM chemical mechanism, as well as accretion product formation from two organic peroxy radicals.

Acknowledgements

This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 730997. Simon O'Meara received funding support from the Natural Environment Research Council through the National Centre for Atmospheric Science.