Difference between revisions of "Photochemistry"
| Line 9: | Line 9: | ||
In order to use photochemistry, you simply need to have the following line in '''callphys.def''': | In order to use photochemistry, you simply need to have the following line in '''callphys.def''': | ||
photochem = .true. | photochem = .true. | ||
| + | Other flags related to photochemistry are: | ||
| + | * '''photoheat''' to turn on or off heating by UV absorption. | ||
| + | * '''jonline''' to turn on or off online photolysis. | ||
| + | * '''stellarflux''' to indicate name of the stellar flux file. This file is located in datadir/stellar_spectra/photochem_stellar_spectra/, has two column: first the wavelength in nm, and second the flux in photons/cm²/s/nm, and has a 2-lines header. | ||
| + | * '''depos''' to turn on or off deposition. | ||
=== traceur.def === | === traceur.def === | ||
| Line 22: | Line 27: | ||
=== Photolysis === | === Photolysis === | ||
| + | ==== Online photolysis ==== | ||
| + | In the baseline use of the photochemistry model, photolysis are computed at the runtime (online), by solving the UV radiative transfer using the TUV solver, considering species-specific UV cross section. o do so, first make sure that the apropriate flag is set to true in the '''callphys.def''': | ||
| + | jonline = .true. | ||
Photolysis reactions are identified by a 0 at column 101 of the '''reactfile'''. They need two types of files (which need to be placed in '''datadir/cross_sections/'''): | Photolysis reactions are identified by a 0 at column 101 of the '''reactfile'''. They need two types of files (which need to be placed in '''datadir/cross_sections/'''): | ||
# cross sections (syntax: first line: wavelength in nm, second line cross section in cm<sup>2</sup>, first line is ignored) | # cross sections (syntax: first line: wavelength in nm, second line cross section in cm<sup>2</sup>, first line is ignored) | ||
| Line 29: | Line 37: | ||
o2 hv o o 0 jo2_o 4 150 200 250 300 O2_150.abs O2_200.abs O2_250.abs O2_300.abs o2_o.yld | o2 hv o o 0 jo2_o 4 150 200 250 300 O2_150.abs O2_200.abs O2_250.abs O2_300.abs o2_o.yld | ||
In this case, 4 different cross section files exist, corresponding to the temperature range 150 K-200 K, 200 K-250 K, 250 K-300 K, >300 K. | In this case, 4 different cross section files exist, corresponding to the temperature range 150 K-200 K, 200 K-250 K, 250 K-300 K, >300 K. | ||
| + | ==== Offline photolysis ==== | ||
| + | ==== Constant photodissociation rates ==== | ||
| + | Sometimes photodissociation rates are simply given as a constant throughout the atmosphere. One way to reproduce that is to treat these reactions as type 1 reaction (see below). | ||
=== Reaction rate formulae === | === Reaction rate formulae === | ||
| Line 38: | Line 49: | ||
For example, the following line describes the 3-body reaction CO + O + M = CO<sub>2</sub> + M, whose reaction rate is <math>k=2.2\times10^{-33}e^{-\frac{1780}{T} }[M]</math> (be careful, here <math>[M]</math> corresponds to the molecular density of the atmosphere, not that of the third body): | For example, the following line describes the 3-body reaction CO + O + M = CO<sub>2</sub> + M, whose reaction rate is <math>k=2.2\times10^{-33}e^{-\frac{1780}{T} }[M]</math> (be careful, here <math>[M]</math> corresponds to the molecular density of the atmosphere, not that of the third body): | ||
co o M co2 1 2.2e-33 -1780.0 0.0 1 1.0 | co o M co2 1 2.2e-33 -1780.0 0.0 1 1.0 | ||
| + | |||
| + | == Boundary conditions == | ||
| + | === Deposition velocity === | ||
| + | === Production rate === | ||
| + | === Escape === | ||
== Using hard-coded rate constants == | == Using hard-coded rate constants == | ||
| − | You might sometimes want to track a reaction whose rate constant is not capture by either of the three above parametrizations. In this case you'll need to get your hands dirty and hard-code the rate constant parametrization in the code. Notice that the reaction should not be included in the | + | You might sometimes want to track a reaction whose rate constant is not capture by either of the three above parametrizations. In this case you'll need to get your hands dirty and hard-code the rate constant parametrization in the code. Notice that the reaction should not be included in the reactfile in this case (it is however a good practice to indicate it as a commented line and write "hard-coded', both to keep track of your network and because the post-processing tools will notify it by reading the reactfile). Here is how to proceed: |
First, identify the type of your reaction: currently, the model distinguishes between 1) bimolecular reactions (A + B -> C (+ D)), 2) quadratic reactions (A + A -> B (+C)), 3) quenching reactions (A + C -> B + C), and 4) heterogeneous reactions (e.g. A + ice -> B + C). Bimolecular reactions are given type '''v4''', quadratic reactions are given type '''v3''', and both quenching an heterogeneous reactions are given type '''vphot''' (the same as photodissociation reactions). | First, identify the type of your reaction: currently, the model distinguishes between 1) bimolecular reactions (A + B -> C (+ D)), 2) quadratic reactions (A + A -> B (+C)), 3) quenching reactions (A + C -> B + C), and 4) heterogeneous reactions (e.g. A + ice -> B + C). Bimolecular reactions are given type '''v4''', quadratic reactions are given type '''v3''', and both quenching an heterogeneous reactions are given type '''vphot''' (the same as photodissociation reactions). | ||
| Line 77: | Line 93: | ||
!end do | !end do | ||
(except the first one) to use the parametrization from the JPL 2015 collection (linked above). Notice that several parametrizations are implemented; to use another one simply uncomment another block, or write your own. | (except the first one) to use the parametrization from the JPL 2015 collection (linked above). Notice that several parametrizations are implemented; to use another one simply uncomment another block, or write your own. | ||
| + | |||
| + | === Hard-coding reactions with 3 products === | ||
| + | The model does not natively handle reactions with 3 products. The trick is to instead split it into 2 "semi-reactions", for instance instead of HNO3 + OH -> H2O + O + NO2, add 0.5 HNO3 + 0.5 OH -> H2O + 0.5 NO2 and 0.5 HNO3 + 0.5 OH -> O + 0.5 NO2 (or any stoichiometrically-equivalent couple). Each semi-reaction then has the same constant rate as the total reaction. Don't forget to count them as two distinct reactions. | ||
Revision as of 15:02, 23 October 2024
Generic PCM
General Use
Compilation
Photochemistry in the generic model is largely independent of the and does not require to recompile it (just make sure that you have linked the LAPACK and BLAS libraries).
callphys.def
In order to use photochemistry, you simply need to have the following line in callphys.def:
photochem = .true.
Other flags related to photochemistry are:
- photoheat to turn on or off heating by UV absorption.
- jonline to turn on or off online photolysis.
- stellarflux to indicate name of the stellar flux file. This file is located in datadir/stellar_spectra/photochem_stellar_spectra/, has two column: first the wavelength in nm, and second the flux in photons/cm²/s/nm, and has a 2-lines header.
- depos to turn on or off deposition.
traceur.def
All chemical species must be listed in traceurs.def (which needs to follow the modern layout).
Reaction network
Besides the species- and reaction-specific data files (see below), you need one file: reactfile (located in the chemnetwork folder) for the reaction network. The first line of the file is ignored, then for each line, the syntax of the reactfile is as follows:
- Columns 1-50 are reserved for reactants. Their names must match those given traceurs.def, in particular, if you have both water vapour and water ice as tracers, named h2o_vap and h2o_ice respectively, the name to use for water vapour in chemical reaction is h2o_vap. Reactant species must be separated by any number of blank spaces.
- Columns 51-100 are reserved for products, with the same rules applying. Notice that in the case of a third body reaction, the third body appears only with the reactants and is omitted on the products side.
- Column 101 is an integer designating the type of reaction: 0 for photolysis, 1, 2 or 3 for other reactions, depending on their reaction rate formula. The four next spaces are blank.
- Column 106 and after are dedicated to parameters specific to each reaction type.
Note that lines can be commented out by starting them with the character "!".
Photolysis
Online photolysis
In the baseline use of the photochemistry model, photolysis are computed at the runtime (online), by solving the UV radiative transfer using the TUV solver, considering species-specific UV cross section. o do so, first make sure that the apropriate flag is set to true in the callphys.def:
jonline = .true.
Photolysis reactions are identified by a 0 at column 101 of the reactfile. They need two types of files (which need to be placed in datadir/cross_sections/):
- cross sections (syntax: first line: wavelength in nm, second line cross section in cm2, first line is ignored)
- yield rates (syntax: first line: wavelength in nm, second line: photodissociation yield, first line is ignored)
These data can be found in benchmark cases for common species (e.g. water), or on the following pages for more exotic ones: Mainz spectral atlas, phidrates, Leiden University. The parameters for photolysis reaction in the reactfile are given as follows, from column 106: an identifier character string starting with "j", an integer corresponding to the number of cross section files (one for each temperature range), the lower bound of each temperature range, the name of each cross section file, the name of the yield rate file. For example, the following line describes the photolysis of O2 into two O atoms in the ground state:
o2 hv o o 0 jo2_o 4 150 200 250 300 O2_150.abs O2_200.abs O2_250.abs O2_300.abs o2_o.yld
In this case, 4 different cross section files exist, corresponding to the temperature range 150 K-200 K, 200 K-250 K, 250 K-300 K, >300 K.
Offline photolysis
Constant photodissociation rates
Sometimes photodissociation rates are simply given as a constant throughout the atmosphere. One way to reproduce that is to treat these reactions as type 1 reaction (see below).
Reaction rate formulae
Many reaction rate parametrizations can be found in here. The model currently accounts for three different formula for the reaction rate (more exotic ones can be hard-coded in the chimiedata_h.F90 file, as indicated in the TRAPPIST-1e example). The chosen formula is by 1, 2, or 3 at column 51 or the reactfile.They read respectively:
- The simple form of the reaction rate is
\[k=a\left(\frac{T}{T_0}\right)^ce^{\frac{-b}{T} }[M]^d\]Its parameters (to give in the reactfile in the same order, starting from column 106) are \(a, b, c, T_0, d\).
- Termolecular reactions have low- and high-pressure-limiting rate constant (\(k_0\) and \(k_\infty\), respectively). Their effective rate constant is then given by
\[k=ge^{-\frac{h}{T} }+\frac{k_0\left(\frac{T}{T_0}\right)^ne^{\frac{-a_0}{T} }[M]^{d_{\text up} } }{1+\frac{k_0\left(\frac{T}{T_0}\right)^ne^{\frac{-a_0}{T} } }{k_\infty\left(\frac{T}{T_0}\right)^me^{\frac{-b_0}{T} } }[M]^{d_{\text down} } }fc^{\frac{1}{1+\left[\log_{10}\left(\frac{k_0\left(\frac{T}{T_0}\right)^ne^{\frac{-a_0}{T} } }{k_\infty\left(\frac{T}{T_0}\right)^me^{\frac{-b_0}{T} } }[M] \right) \right]^2 } }\]Its parameters (to give in the reactfile in the same order, starting from column 106) are \(k_0, n, a_0, k_\infty, m, b_0, T_0, f_c, g, h, d_{\text up}, d_{\text down}\).
- Finally, a last form of the reaction rate is
\[\frac{k_0\left(\frac{T}{T_0}\right)^ne^{\frac{-a_0}{T} }[M]^{d_{\text up} } }{1+\frac{k_0\left(\frac{T}{T_0}\right)^ne^{\frac{-a_0}{T} } }{k_\infty\left(\frac{T}{T_0}\right)^me^{\frac{-b_0}{T} } }[M]^{d_{\text down} } }[F_{\text cent}]^X,\]with \(F_{\text cent}=1-ae^{-\frac{T}{T^{***} } }+ae^{-\frac{T}{T^* } }+ae^{-\frac{T^{**}}{T } }\), and \[X=\left[1+\left[\frac{\log_{10}\left(\frac{k_0\left(\frac{T}{T_0}\right)^ne^{\frac{-a_0}{T} } }{k_\infty\left(\frac{T}{T_0}\right)^me^{\frac{-b_0}{T} } }[M] \right)+c }{N-d\left(\log_{10}\left(\frac{k_0\left(\frac{T}{T_0}\right)^ne^{\frac{-a_0}{T} } }{k_\infty\left(\frac{T}{T_0}\right)^me^{\frac{-b_0}{T} } }[M] \right) \right)+c } \right]^2 \right]^{-1},\]with \(c=-0.4-0.67\log_{10}(F_{\text cent})\), \(N=0.75-1.27\log_{10}(F_{\text cent})\), and \(d=0.14\).
- Its parameters (to give in the reactfile in the same order, starting from column 106) are \(k_0, n, a_0, k_\infty, m, b_0, T_0, a, T^{***}, T^*, T^{**}, d_{\text up}, d_{\text down}\).
For example, the following line describes the 3-body reaction CO + O + M = CO2 + M, whose reaction rate is \(k=2.2\times10^{-33}e^{-\frac{1780}{T} }[M]\) (be careful, here \([M]\) corresponds to the molecular density of the atmosphere, not that of the third body):
co o M co2 1 2.2e-33 -1780.0 0.0 1 1.0
Boundary conditions
Deposition velocity
Production rate
Escape
Using hard-coded rate constants
You might sometimes want to track a reaction whose rate constant is not capture by either of the three above parametrizations. In this case you'll need to get your hands dirty and hard-code the rate constant parametrization in the code. Notice that the reaction should not be included in the reactfile in this case (it is however a good practice to indicate it as a commented line and write "hard-coded', both to keep track of your network and because the post-processing tools will notify it by reading the reactfile). Here is how to proceed:
First, identify the type of your reaction: currently, the model distinguishes between 1) bimolecular reactions (A + B -> C (+ D)), 2) quadratic reactions (A + A -> B (+C)), 3) quenching reactions (A + C -> B + C), and 4) heterogeneous reactions (e.g. A + ice -> B + C). Bimolecular reactions are given type v4, quadratic reactions are given type v3, and both quenching an heterogeneous reactions are given type vphot (the same as photodissociation reactions).
Then, open the file LMDZ.GENERIC/libf/aeronostd/chimiedata_h.F90. This file has a dedicated section for hard-coded reactions. Some of them are already there, commented out. We'll use one of them (the reaction CO + OH -> CO2 + H) as a example on how to proceed. This is a bimolecular reaction, so we first need to increment the counter of hard-coded reactions of type v_4:
integer, parameter :: n4_hard_coding = 1
In the master branch version of the code, the variable n4_hard_coding is set to 0. If you're version of the code already has hard-coded reactions, then simply add one to the current value of n4_hard_coding.
In the same file, you will find a subroutine called indice_HC for registering hard-coded reactions. You can create a space for your reaction, taking example of what's already there. For instance our example reaction is already there:
!===========================================================
! e001 : CO + OH -> CO2 + H
!===========================================================
! nb_reaction_4 = nb_reaction_4 + 1
! indice_4(nb_reaction_4) = z4spec(1.0, indexchim('co'), 1.0, indexchim('oh'), 1.0, indexchim('co2'), 1.0, indexchim('h'))
To use it, you must uncomment the two last lines. The top one increments the number of reactions of type v4, and the second one registers the reacting species for your new reaction.
Then go to the subroutine reactionrates_HC. This is where you will write down the parametrization formula into the reaction rates table (the variable v_4 for bimolecular reactions, v_3 for quadratic reactions, or v_phot for all others). Again, you can take example of what's already there. If you want to use our example reaction, find the dedicated space, opening with the line:
!--- e001: oh + co -> co2 + h
and simply uncomment the following lines:
!nb_reaction_4 = nb_reaction_4 + 1
to update the index at which the rate is written in table v_4, and then (for instance):
! jpl 2015 !do ilev = 1,nlayer !!oh + co -> h + co2 ! rate1 = 1.5e-13*(t(ilev)/300.)**(0.0) !!oh + co + m -> hoco ! ak0 = 5.9e-33*(t(ilev)/300.)**(-1.0) ! ak1 = 1.1e-12*(t(ilev)/300.)**(1.3) ! rate2 = (ak0*dens(ilev))/(1. + ak0*dens(ilev)/ak1) ! xpo2 = 1./(1. + alog10((ak0*dens(ilev))/ak1)**2) ! ! v_4(ilev,nb_reaction_4) = rate1 + rate2*0.6**xpo2 !end do
(except the first one) to use the parametrization from the JPL 2015 collection (linked above). Notice that several parametrizations are implemented; to use another one simply uncomment another block, or write your own.
Hard-coding reactions with 3 products
The model does not natively handle reactions with 3 products. The trick is to instead split it into 2 "semi-reactions", for instance instead of HNO3 + OH -> H2O + O + NO2, add 0.5 HNO3 + 0.5 OH -> H2O + 0.5 NO2 and 0.5 HNO3 + 0.5 OH -> O + 0.5 NO2 (or any stoichiometrically-equivalent couple). Each semi-reaction then has the same constant rate as the total reaction. Don't forget to count them as two distinct reactions.
