Introducing atmospheric photochemical isotopic processes to the PATMO atmospheric code

Abstract

The study of production and depletion of chemical species is vital for the understanding of composition and evolution of planetary atmospheres. We present an implementation of photoinduced isotopic effects into the PATMO code (Ávila et al., 2021) code designed for the study of stable isotopes and photo-induced isotopic effects. With respect to the original code PATMO, where the photochemistry was not included, this report extends capability of the model to set photochemical processes for stable isotopes and thus enhancing its applicability. The PATMO code is flexible and allows the edition of new chemical reactions without need for hard code them. We also test how changes in spectral resolution affects the calculation of isotopic effects during the photodissociation of oxygen. We found that for a highly structured spectrum such as the Schumann-Runge band a spectral resolution larger than 0.005 nm is necessary for accurately modeling these isotopic effects. We also show that SO₂ and SO photodissociation couples in a complex shielding fashion and significantly affects the photo-induced isotopic effects. The model was also benchmarked against today’s Earth atmosphere, where the solar UV flux and the ozone profile of the US Standard atmosphere of 1976 was reproduced with the simple Chapman mechanism and improved with the implementation of NOx and HOx radicals.

Introduction

The chemical composition of planetary atmospheres is the outcome of a complex combination of physico-chemical processes. The atmospheric composition of today’s planet Earth is the most thoroughly studied system (e.g., Kawamiya et al., 2020; Morgenstern et al., 2017). The level of accuracy and complexity of these models is desirable for other atmospheric systems such the early Earth and other planets. For the study of isotopic effects associated to photochemical processes one-dimensional (1D) models are currently the most effective and at the same time accurate tools.

The study of the composition of atmospheres for planets such as the early Earth or exoplanets present a challenging case since there are limited observational data that can be used to validate the model. For the case of Earth’s past atmospheres, the geological record is the only source of information about atmospheric conditions. The lack of observational data forces modelers to rely on first principles and empirical parametrization determined by observations on present-day Earth and its neighboring planets.

The early stages of the study of atmospheric chemistry for systems outside present-day Earth was first summarized by Yung and DeMore (1999) and lately by Catling and Kasting (2017). The study of the atmosphere of rocky planets requires an understanding of photochemistry since the composition and distribution of trace constituents are governed by photochemical reactions. This situation is also significantly affected by surface emission and sink processes such as chemistry and deposition. The key processes required to estimate atmospheric composition are photolysis, kinetics, molecular diffusion, atmospheric escape, dry and wet deposition, and condensation and sedimentation (Catling and Kasting, 2017).

1D models have proven to be a reliable tool for studying several aspects of planetary atmospheres. The literature is extensive, and a review is beyond the scope of this paper. Seminal work includes photochemical models for Venus (Yung and Demore, 1982), Mars (Nair et al., 1994), Titan (Yung et al., 1994), several stages of the early Earth atmosphere (Kasting, 1993; Pavlov et al., 2000), and exoplanets (Kasting et al., 1993). A similar methodology has been implemented to study atmospheric chemical disequilibrium and how this parameter can be used for detecting exoplanet life (Simoncini et al., 2013). The development and application of these models face the problem of coping with large numbers of parameters that are difficult to calculate from first principles and therefore need to be included as “free parameters”. Physical parameters such as condensation, vertical diffusion, and rates of dry and wet deposition vary from system to system, but they can be estimated within some degree of accuracy (Giorgi and Chameides, 1985). Fluxes of surface emissions are specific to each planet and epoch, and they need to be treated as assumed arbitrary conditions for many systems. In order to avoid time consuming and unimportant calculations like photodissociation rate constants or pressure dependent third-body reactions, chemical reaction networks are reduced to a minimal number of reactions needed to properly describe a given system. Therefore, reducing (H₂-rich), oxidized (N₂-rich, or CO₂-rich), and oxidizing (O₂-rich) atmospheres have been studied with different chemical networks (Catling and Kasting, 2017). To this aim in recent years, flexible methods of coding have made possible modifications of chemical networks without the need to hard code all variables and hence easing the benchmarking process, as for example KROME (Grassi et al., 2014).

Numerous physical and chemical phenomena are accounted for in modern one-dimensional models and have also been applied to many systems. The increment in computational power in the past decades has made possible the use of accurate numerical schemes. Such improvement includes the calculation of photodissociation rate constants at different altitudes (pressures and temperatures), implementation of high-resolution spectra and how change in concentrations from initial conditions affects the UV light opacity properties and their associated rate constants.

These novel techniques are fundamental to explore complicated chemical models, as for example isotopologues, a powerful tool to reveal geochemical mechanisms. Implementing them to chemical networks can double, triple or quadruple (such is the case of ^32/33/34/36S) the number of reactions when single substitution is considered, and this number grows exponentially when double, triple, and more substitutions are considered (clumped isotopes such as ^32/33/34/36S^16/17/18O). There are few studies in the literature that add stable isotopes to their numerical schemes (Claire et al., 2014; Izon et al., 2017; Pavlov and Kasting, 2002). However, there are no attempts to account for isotopic effects caused by photodissociation when high-resolution spectra are considered, and there are no successful attempts to account for mass balance models that include isotopic effects.

In this work we present a photochemical model constructed from basic physical and chemical equations that have a long standing in the literature. The code employed allows to model of a wide range of atmospheres with a significant variability in composition and in redox states. The model also treats as free parameters surface gravity, solar (electromagnetic radiation) flux and planetary specific parameters such as vertical temperature profile, mean molecular mass of the atmosphere, coefficients of vertical diffusion, fluxes for surface emissions, and dry and wet deposition. Stable isotopes are treated as independent chemical species, so that isotopic ratios and non-mass dependent coefficients can be calculated. The model is also capable to calculate chemical disequilibrium; this property has been proposed has a signature of living systems and is therefore applicable to the search of life in exoplanets (Barge et al., 2017).

This paper is structured as follows; in Section “Model Description” we describe in detail the equations of the photochemical model. Section “Study Cases” presents benchmark calculations of present Earth. Conclusions and future directions are discussed in Section “Summary”.

Model Description

This section presents the basic equations included in the model. A detailed description of the chemical networks, chemical species, optically active species, photoreactions and boundary conditions are presented in the study case reactions. This report is a follow up work by Ávila et al. (2021) and focused on chemical networks, photochemistry and isotopic effects.

A code for studying planetary atmospheres (PATMO)

The evolution of each chemical species starting from a given initial condition is the key feature of the code presented in this paper. The Planetary ATMOspheres (PATMO) code is a flexible tool to compute the chemical evolution of planet atmospheres. PATMO uses a multi-layered vertical 1D scheme including (i) time-dependent chemistry (including reverse chemistry and disequilibrium), (ii) molecular/eddy diffusion, (iii) multi-frequency ultraviolet absorption, (iv) photochemistry, (v) deposition associated mechanisms (wet, dry, and gravitational), (vi) surface emissions, and (vii) stable isotopes.

Similarly to KROME (Grassi et al., 2014) PATMO is a PYTHON preprocessor that produces optimized FORTRAN code for specific chemical networks and options for specific planetary conditions. PATMO is a self-consistent and stand-alone code, but it can be also employed as an external library and can be called at every time-step from any framework code, as for example models of a planet embedded in a protoplanetary disk. PATMO follows the approach (Hu et al., 2012; Rimmer and Helling, 2016) of assuming operator splitting between (i) chemistry and diffusion, and (ii) ultraviolet light opacities. Sophisticated radiative transfer models account for absorption, emission, and scattering processes in the UV, Visible, NIR and IR. Since PATMO focusses on chemical processes, absorption in the UV is the only radiative transfer associated process relevant to this model. For this reason, we use the term UV absorption instead of radiative transfer. However, since PATMO uses a line-by-line approach to radiative transfer Visible, NIR and IR spectra can be easily added to the model.

PATMO calculates both (i) chemistry and (ii) molecular and Eddy diffusion equations within the same numerical scheme, i.e., Eqn. 9 in (Hu et al., 2012), but solved time-dependently with the ordinary differential equations (ODEs) solver DLSODES (Hindmarsh, 1983). We take advantage from the fact that in each layer the chemical network is the same, but with different abundances and rate coefficients: this allows a very efficient, vectorization of the ODE system achieving a fast code even if the solver itself is not vectorized. The code includes a module aimed at computing reverse reactions following NASA polynomials (McBride et al., 1993) employing the method described in the next section and the coefficients from Burcat’s thermochemistry database (http://garfield.chem.elte.hu/Burcat/burcat.html). The system of ODEs that includes forward and backward kinetics is particularly complex and numerically inefficient to be solved, therefore the use of fast solvers plays a key role in the numerical stability of the problem. In the next sections we describe in detail the physical components of the model.

PATMO is an open-source code fully available at the Bitbucket repository. All chemical networks, absorption spectra and boundary conditions required to reproduce the results presented in this paper are available at the same repository (https://bitbucket.org/tgrassi/planatmo/src/master/, https://github.com/SophiaAtmo/PATMO).

Fundamental equations

The model solves chemical kinetics and diffusion by using the methodology proposed in the works of Hu and Rimmer (Hu et al., 2012; Rimmer and Helling, 2016), we report hereafter their equations for completeness. The operator that evolves the system over time is described by the ODEs system

  

$$\frac{\partial n_{i,j}}{\partial t}=P_{i,j}-n_{i,j}{L}_{i,j}-\frac{\partial {\mathrm{\Phi }}_{i,j}}{\partial z}$$(1)

where n is the number density [molecules cm^–3], P is the chemical production rate [molecules cm^–3 s^–1], nL the chemical depletion rate [molecules cm^–3 s^–1], and Φ is the vertical diffusion flux [molecules cm^–2 s^–1] for i species and j layers. Further details are presented in Hu et al. (2012) and a description of the discretization procedure is presented in Appendix A.

Photochemistry and radiative transfer

PATMO is capable of handling multi-frequency bin-based photochemistry (Hu et al., 2012; Rimmer and Helling, 2016), where the photochemical cross-sections is introduced and the impinging radiation modeled accordingly. The code used by the Python builder allows an arbitrary number of bins (spectral resolution) as required by specific systems. The photochemical rate coefficients are computed as

  

$$J_{i,z}={\int }_{{\lambda }_1}^{{\lambda }_2}{\varphi }_{\lambda }\left(\lambda \right){\sigma }_{\lambda }\left(\lambda \right)I_{\lambda }\left(z\right)d\lambda ,$$(2)

Where [λ₁, λ₂] is the spectral range, ϕ_i(λ) the quantum yield for the species i [molecules photon^–1], σ_i(λ) its UV absorption cross-sections [cm² molecules^–1], I_λ(z) the solar flux at altitude z [photon cm^–2 s^–1 nm^–1]. The solar flux represents the number of photons that pass through a horizontal plane at altitude z per units of time and surface. The calculation of the solar flux at a given altitude for multiple molecules with UV absorbing properties is as follows,

  

$$I_{\lambda }\left(z\right)=I_{\lambda }\left(\infty \right)e^{-\frac{{\tau }_z}{cos\theta }},$$(3)

where I_λ(∞) is the solar flux at the top of the atmosphere, θ the zenith angle respect to the Equator, and τ_z is the optical depth,

  

$${\tau }_{z,\lambda }=(\infty -z)\sum _i{C}_i{\sigma }_{i,\lambda },$$(4)

The term (∞ – z) represents the total concentration from the top of the atmosphere to altitude z [cm] of molecule i, C_i is the concentration of molecule i [molecules cm^–3]. Photodissociation rate constants are dependent on factors that affect UV absorption cross-sections such as pressure and temperature. These rate constants are also sensitive to atmospheric composition of optically active molecules and their concentrations and therefore functional to altitude.

The properties of equation Eqn. 2 are studied in Subsection “Present Earth”, where a modern Earth is simulated with a 21% oxygen composition (as initial conditions). Under these initial conditions it cannot reproduce the current solar flux until ozone is formed by the Chapman cycle (Figs. A. 8a–c).

Kinetics and reverse chemistry

The reverse kinetic constants k_r are calculated from forward rate coefficients k_f considering local chemical equilibrium at each time step,

  

$$k_r=\frac{k_f}{k_{eq}}{\left(k_{b}T\right)}^{\delta n}$$(5)

where –Δ_RG is the free energy of the reaction, T the temperature, R the universal gas constant, k_b is the Boltzmann constant and δn the difference between the number of products and the number of reactants, where the correction factor is in cgs units. For more detailed description of how reverse reaction rate constants are calculated see 3.6 of Grassi et al. (2014).

Code structure

The PATMO code follows a similar structure as KROME (Grassi et al., 2014) where a PYTHON pre-processor creates the necessary FORTRAN files to compute optical opacities, photodissociation rates, chemical evolution and transport for a user defined list of chemical reactions. Several user function that can be called to facilitate the analysis of the constructed model are automatically added to the FORTRAN code. The source code and a detailed explanation of the installation procedure and available user functions is available as an open-source project in a Bitbucket repository (https://bitbucket.org/tgrassi/planatmo/src/master/).

Study Cases

This section presents a set of specific models constructed with PATMO. Each system uses different chemical networks, photo-reactions, upper and lower boundary conditions, emission, and deposition settings.

Benchmark calculation

Basic benchmark calculations have been discussed at Ávila et al. (2021). In this section the core features of the code a 0-dimentional system with a simple chemical network was tested against KROME. The purpose of this single box with no transport and no photo-chemistry is to verify the ability of the production-depletion equations for a simple system. Table 1 summarizes the chemical network used in this case while Table 2 presents the initial conditions for the chemical species. In this simulation O(³P), O(¹D), O₂, O₃, and N₂ concentrations were kept fixed, while H, OH and OH₂ were monitored through time until they reached equilibrium. The calculation endpoint is defined as d[x]/dt < threshold (x = H, OH, and OH₂) where the threshold value is set to satisfy this specific case. These results (Figs. 1a, b) show that the basic equations created by the Python parser and the solutions implemented by the DLSODES solver are working properly. See additional information (Fig. A. 1) for a study of the system in which all chemical species are calculated.

Table 1. Reaction list for study case 1. All reaction constants are extracted from the NIST data base (Orkin et al., 2015).

No	Reaction	Rate constants/cm3 s–1
1	OH + O3 → HO2 + O2	$${1.7}^{-12}{e}^{\left(\frac{-940}{T}\right)}$$
2	O3 + HO2 → 2O2 + OH	$${1.0}^{-14}{e}^{\left(\frac{490}{T}\right)}$$
3	OH + HO2 → H2O + O2	$${4.8}^{-11}\left(\frac{250}{T}\right)$$
4	O(1D) + H2O → 2 OH	$${1.63}^{-10}{e}^{\left(\frac{60}{T}\right)}$$
5	H + O3 → OH + O2	$${1.4}^{-10}{e}^{\left(\frac{-470}{T}\right)}$$
6	O + OH → H +O2	$${1.8}^{-11}{e}^{\left(\frac{180}{T}\right)}$$
7	H + O2 → HO2	$$\left(\frac{k_0\left(T\right)\left[M\right]}{1+\frac{k_0\left(T\right)\left[M\right]}{k_{\infty }\left(T\right)}}\right){0.6}^{{\left\{1+{\left[{log}_{10}\left(\frac{k_0\left(T\right)\left[M\right]}{k_{\infty }\left(T\right)}\right)\right]}^2\right\}}^{-1}}$$
8	O + HO2 → OH + O2	$${3.0}^{11}{e}^{\left(\frac{200}{T}\right)}$$
9	H + HO2 → OH + OH	7.2–11
10	H + HO2 → O + H2O	1.6–12

The values of k₀ and k_inf for reaction 7 are 4.4*10^-32 and 7.5*10^-11 respectively.

Table 2. Initial conditions used for case 1

	Number Density/molec. cm–3
Temp (K)	N2	O2	H2O	O3	O(3P)	O(1D)
282	1.8619	5.0018	1.8417	7.3011	3.3503	2.10–3

Fig. 1.

Comparison between benchmark calculation with PATMO (dot lines) and KROME (solid lines). Panels a to e show the time evolution of the chemical species represented in the labels.

Photodissociation of a single chemical species at different altitudes

As a second benchmark calculation, photodissociation of SO₂ was studied by a column of air and SO₂ that is top-down irradiated with UV light. The purpose of this calculation is to check the ability of the code to deal with different light-shielding conditions and how this shielding varies through time due to variation of chemical species. The model was constructed with a 10 ppm SO₂ column in a 100 km atmosphere divided into 100 layers with a total pressure of 1 atm at the bottom of and following hydrostatic equilibrium for the pressure-altitude relation. The rest of the atmospheric composition is assumed to be N₂. The exceedingly large and likely unrealistic condition of 10 ppm throughout a 100 km atmosphere was made by design to observe the behavior of the model when number densities change by several orders of magnitude. This system is exposed to a solar flux (Gueymard, 2004) that permeates into the atmosphere vertically from the top layer. The chemistry consists of a single reaction where SO₂ is photodissociated into SO and O radicals (SO₂ → SO + O). For this calculation the SO₂ a spectra with a 0.1 spectral resolution was used (Endo et al., 2015). Figure 2 shows the photodissociation rate (Eqn. 2) at each wavelength at 8 different altitudes. This figure also shows that most of the spectrum contributes to the total photodissociation value at 100 km (red line). The contribution from shorter wavelengths is reduced at lower altitudes (green, blue, and grey lines) and in the lower part of the atmosphere only long wavelengths contribute to the rate constant. This distribution is caused by the shielding effect described by Eqn. 3 and 4.

Fig. 2.

Photodissociation rate constants at each wavelength calculated at eight different altitudes. The plot shows that at higher altitudes the entire absorption spectra contributes to the total photodissociation rate constant. At lower altitude, low energy photons contribute to the photodissociation rate constant.

The PATMO model was also compared to KROME, and the results are Fig. 3 upper and lower panels where the time evolution of the SO₂ profile with altitude is shown. The photodissociation process indicates that SO₂ on the top layer begins to dissociate as soon as the simulation starts, while the deeper layers remain unaltered. This is because the UV light is completely absorbed in the top layer, and only after the SO₂ concentration decreases (due to photolysis) the lower layers are irradiated. Figure 4 upper and lower panels show the same behavior as Fig. 3 in terms of the evolution of the photodissociation rate. (Eqn. 2). At the beginning of the simulation (t₀) the largest value of J_z is at the top of the atmosphere; this value drops exponentially with the depth of the atmosphere. At t = 2.48 years (21,744 hours) the entire atmosphere has become SO₂ free, therefore the opacity term (Eqn. 4) becomes negligible with a unified reaction rate constant throughout the entire atmosphere (J_z is the same value at all altitudes). The comparison between PATMO and KROME shows that given equal boundary conditions both codes produce the same numerical results. For further analysis, refer to additional information (Fig. A. 2) where time evolution of initial SO₂ concentration (number density) and photodissociation rate constants at 1 and 91 km are monitored through the simulation.

Fig. 3.

SO₂ vertical profiles calculated by PATMO (upper panel) and KROME (lower panel). As the simulation progress SO₂ depletion on the top layers gradually allows for the photochemistry at lower layers. This spatial related property is known as the self-shielding effect.

Fig. 4.

Photodissociation rate constants from the beginning of the simulation (t₀) until 21,768 hours calculated by PATMO (upper panel) and KROME (lower panel). The largest values of photodissociation rate constants are at the top of the atmosphere and they exponentially decrease with atmospheric depth. As the simulation proceeds the SO₂ concentration at the upper layers decreases and consequently the rate constants of layers immediately below increase. By the end of the simulation all the SO₂ is photodissociated and the photodissociation rate

constants are equal throughout atmosphere, due to negligible or 0 self-shielding effect.

Photodissociation rates and isotopic effects, changes in spectral resolution

The following study case is in principle similar to the previous one but extended to photochemically induced isotopic effects and how they are affected by altitude, shielding, opacity, and spectral resolution. The applications of isotopic effects to atmospheric modeling have a wide range of applications such as the effects on oxygen during ozone formation (Cicerone and McCrumb, 1980; Thiemens, 2001), or the sulfur isotopic effect in the Archean atmosphere and its relationship with geological records (Ono, 2017). This model presents the isotopic effect under different spectral resolution during photodissociation of oxygen. Figure 5 presents a simulated ¹⁶O₂, ¹⁶O¹⁸O and ¹⁸O₂ UV absorption spectrum with a resolution of 0.00125 nm. The spectra in this calculation extend from 120 to 250 nm i.e., the Schumann-Runge bands, and the Herzberg continuum. PATMO works with any user-defined spectral resolution. The code then splits the spectrum into linear equi-energetic bins. Table 3 shows the number of bins required to describe the spectra at different spectral resolution. The model was constructed with a 100 km altitude divided into equally spaced 100 layers with a ground atmospheric pressure of 1 atm and 21% of oxygen. Then photodissociation rate constants (J_i,z) are calculated for reaction 1 (R1) isotopologues at all 100 layers.

Fig. 5.

Simulated spectra after measurements of Yoshino et al for ¹⁶O₂, ¹⁶O¹⁸O and ¹⁸O₂ with a resolution of 0.00125 nm. The inbox shows how the isotopic effect changes the structure of the spectra.

Table 3. Oxygen spectra at different resolutions and the number of bins

Spectra Resolution (nm)	Number of bins used to describe the 120–250 nm region.
0.028	5,000
0.014	10,000
0.009	14,000
0.005	23,000
0.004	32,000
0.00125 (Full Resolution)	100,0000

  

$$O_2\stackrel{h\nu }{\to } 2O R1$$

As presented in Eqn. 2, the calculated values of J_i,z at a given altitude (z) is affected by the UV flux depletion from the layers above. This process is usually known as shielding. For this specific case, where oxygen in the upper layers shields (or altitude) absorption in the lower layers; the phenomenon is called self-shielding.

From the photodissociation rate constants isotopic enrichment coefficients were calculated according to the following equation:

  

$${}^{1x}O{}^{1x}O{\epsilon }_z=Ln\left(\frac{{{}^{1x}O{}^{1x}OJ}_z}{{{}^{16}O{}^{16}OJ}_z}\right),$$(6)

where x represents mass 16 or 18 for a single or double substituted oxygen isotopologue. Figure 6 left and right panels show the isotopic enrichment coefficients as function of altitude and spectral resolution. The plot shows that self-shielding dramatically changes the isotopic effect with altitude, and that the variation with spectral resolution is relevant. The specific mechanism that originates the isotopic fractionation constant is the outcome of a complex interaction between the specific absorption and the highly-structured solar irradiance spectra (Gueymard, 2004). For more details on how all these spectra interact among themselves see Figs. A. 3 and A. 4 in the Supplementary Information.

Fig. 6.

¹⁶O¹⁸O and ¹⁸O₂ isotopic enrichment factors produced by photodissociation. Left panel shows enrichment values (ε in ‰ units) for ¹⁶O¹⁸O isotopologues, right panel shows same ε values for ¹⁸O₂ isotopologues. Green, yellow, light blue, red and blue lines represent different spectral resolutions; 0.028 nm, 0.014 nm, 0.009 nm, 0.005 nm, 0.004 nm respectively. Both plots show that isotopic effects for this reaction can be properly calculated with resolution of 0.005 nm.

From this study case it is quite clear that the origin of the isotopic effect is a complex combination of solar flux, absorption spectra, and shielding phenomena.

Photodissociation of SO₂-SO system and their associated isotopic effects

The next benchmark calculation combines the previous two study cases. Here the isotopic fractionation factors (Eqn. 6) are evaluated not only as function of altitude, but also as the change in concentration of the chemical species in the system. Note that, the isotopic studies at this stage only discuss isotopic enrichment factors and not isotopic ratios (δ values as defined at Sharp, 2007). An assessment of isotopic ratios will require mass balance and isotopic effects of all the relevant chemical processes. A more detailed study on the atmospheric isotopic cycle is the aim of a forthcoming study. In this study case, an atmosphere with a constant 1 ppm SO₂ column divided equally spaced 100 layers is exposed to a UV solar flux (Gueymard, 2004) with the following set of reactions:

  

$${SO}_2\stackrel{h\nu }{\to } SO+O R2$$

  

$${SO}_{}\stackrel{h\nu }{\to } S+O R3$$

An initial 1 ppm atmospheric mixing ratio is consistent with a volcanic plume (Hattori et al., 2013). Such settings are slightly excessive for an entire atmospheric column but is optimal for the purpose of testing the ability of the model to account for shielding phenomena throughout the entire atmosphere. The key feature of this study case is that the UV absorption spectrum of both molecules (SO and SO₂) overlaps (Fig. 7). When the SO concentration reaches an equilibrium between production from SO₂ photolysis and depletion from SO photolysis (Reactions R2 and R3), the opacity term (Eqn. 4) has unique values at each altitude and wavelength.

Fig. 7.

Ultraviolet absorption cross sections of ³²SO₂ and ³²SO, these data sets were taken from Endo et al. (2015) and Danielache et al. (2014) respectively. Both spectra significantly overlap and have highly structured waveforms with large differences in absorption intensities at different wavelengths.

Figure 8 upper and lower panels show a fraction of the isotopologues spectra of ^32,33,34,36SO₂ (Endo et al., 2015) and SO (Danielache et al., 2014). The spectral resolution of all isotopologues has been degraded to fit 3000 equally spaced energy bins. It is important to stress that this benchmark calculation focus on the sensitivity of the isotopic effect when the concentrations of UV absorbers change from their initial conditions. These are idealized settings (no other chemistry is involved in the production or depletion of SO₂ and SO nor other UV absorbents are considered) and therefore not suitable for accurate geochemical predictions. The photodissociation rates constants (J_i,z) were calculated with a photodissociation threshold of 219.2 nm for SO₂ and 231.1 nm for SO, however the opacity term for the solar dimming process was calculated with the full range of the absorption spectra of both molecules. At t₀, when SO is not present in the system, calculated isotopic enrichment factors ^33,34,36ε and ^33,36E at the top of the atmosphere are not affected by shielding. Under these conditions ^33,34,36ε and ^33,36E values can be directly compared to those reported by Endo et al. (2015). These results are compared in Figs. 9 (³⁶E vs ³³E) and 10 (³³E vs ³⁴ε) where the modeled enrichment factors accurately reproduce those of Endo et al. (2015). The comparison with the literature data yields a good benchmark which shows that the code is properly calculating Eqn. 2.

Fig. 8.

UV absorption spectra of ^32,33,34,36SO (upper panel) and ^32,33,34,36SO₂ (lower panel). From this figure can be observed that the photodissociation rate constant when calculated at each wavelength (Eqn. 2) is very specific and the overall isotopic effect is a complex combination cross-sections and solar fluxes.

Fig. 9.

Fractionation factors ³⁶E vs ³³E in ‰ units calculated at the top of the atmosphere (100 km) for the SO₂ molecule at the beginning of the simulation. These specific settings represent an unshielded pure SO₂ photodissociation process similar to the one calculated by Endo et al. (2015).

Fig. 10.

Fractionation factors ³³E vs ³⁴ε calculated at the top of the atmosphere (100 km) for the SO₂ molecule at the beginning of the simulation. These specific settings represent an unshielded pure SO₂ photodissociation process similar to the one calculated by Endo et al. (2015). The comparison with the literature data yields a good comparison which means that the code is properly calculating Eqn. 2.

Figure 11 shows the time evolution of SO₂ and SO mixing ratios at 20 km. The model is set to have a constant SO₂ concentration throughout the simulation. SO mixing ratios reach an equilibrium between production from SO₂ photodissociation and depletion from SO photodissociation. Once the model reaches a steady state after 6 years, the mixing ratios remain constant. Figure 12 shows the photolysis rates for SO₂ as function of altitude, the solid green line shows the rates at the initial condition of 1 ppm. The smaller rates at lower altitudes are produced by self-shielding, that is, the reduction of solar flux by absorption in the upper layers. When the system reaches equilibrium, the rates became smaller (Fig. 12 orange line) due to the extra shielding produced by the SO molecule. Table 4 shows the numerical values of the reactions rates at t₀ and after the system reaches the steady state.

Fig. 11.

Time evolution of SO (violet line) and SO₂ (blue line) concentrations, SO₂ number density shows no change throughout the simulation since is set to stay constant. The SO molecule starts with zero number density and increases until it reaches a dynamical balance between production and depletion.

Fig. 12.

SO₂ Photodissociation rate constants from 0 to 25 km at the beginning of the simulation (green line) and after the steady state is achieved (orange line). Photodissociation rate constants are smaller below 25 km due to the shielding effect.

Table 4. Photodissociation rate constants [s^–1] for the dissociation of SO₂ at 20 km at the beginning of the simulation (t₀) and after the system reaches a steady state (t ≈ 6000 hrs).

	32SO2	33SO2	34SO2	36SO2
t0	1.91*10–5	2.06*10–5	2.33*10–5	2.72*10–5
Steady Sate	9.65*10–6	1.00*10–5	1.09*10–5	1.25*10–5

Next, isotopic enrichment factors ^33,34,36ε and ^33,36E were calculated at each layer. The results are shown in Fig. 13 where ³³E and ³⁶E values at the initial conditions (solid lines) are significantly different to those at steady state (dashed lines). These results show the importance of calculating opacities accurately when concentrations change from their initial values, and the role played by shielding from trace species. See additional information (Fig. A. 10) where the variation of the isotopic signal changes over time expressed at the ³³E/³⁶E space. During the preparation of this paper Endo and co-workers (Endo et al., 2022) have reported ^32,33,34,36SO₂ high resolution UV absorption spectra in the 206–220 nm region. A detailed study on the incidence of spectral resolution coupled to spectral errors, spectral range, and variability of the solar spectra deserves a substantial amount of work that would render this paper exceedingly long.

Fig. 13.

SO₂ mass independent isotopic fractionation factors ³³E in ‰ units (upper horizontal axis and red lines) and ³⁶E (lower horizontal axis and black lines) at each altitude for the initial conditions (solid line), and at steady state conditions (dotted line). As observed in both ³³E and ³⁶E, above c.a. 23 km the atmospheres is thin and shielding effects become irrelevant making the ³³E and ³⁶E values fixed throughout the simulation. This values compare quite well with the values presented in Figs. 9 and 10 and also with values reported by Endo et al. (2015) (³³E = 8.18‰, ³³E = –10.05‰).

Present Earth

The next step is to verify the ability of the model to reproduce a specific atmospheric system. The Chapman cycle is an ideal system for a validation study of a 1D model, since an averaged vertical distribution describes the atmosphere accurately. The Chapman cycle is a mechanism that explains the shielding of ultraviolet radiation and production of ozone by photoreactions of oxygen under UV light, and the reactions of atomic oxygen with ozone and oxygen molecules. Table 5 shows the chemical reactions involved in the process while Table 6 shows the UV absorption spectrum and quantum yield of photoactive species oxygen and ozone. The solar spectrum (Gueymard, 2004) extents from 180 to 400 nm, namely the Schuman-Runge band and Herzberg continuum for oxygen and the Hartley, Huggins, and Chappuis bands for ozone. The spectrum was split into 4400 equally and linearly spaced photo-bins (0.05 nm resolution), which is a good balance between accuracy for the calculation of photodissociations constants, and computational cost. The atmosphere was divided into 100 layers of 1 km each and the calculation was executed for 100 years. For this calculation nitrogen and oxygen concentrations were constants throughout the entire simulation. Temperature profiles are assumed to be the US Standard atmosphere of 1976 (United States. National Oceanic and Atmospheric Administration and United States Committee on Extension to the Standard Atmosphere, 1976) and Eddy diffusion profiles are adapted from Massie and Hunten (1981) and presented in the supplementary information section (Fig. A. 5). In order to further validate the consistency of the model, a similar simulation was prepared with KROME.

Table 5. List of reactions used to construct the Chapman system described in Subsection “Present Earth”. The reaction rate constants were taken from (Orkin et al., 2015).

model	Reaction	Rate Constant
Chapman Model	O + O2 + M → O3 + M	6.0*10–34*(T/300)–2.4a
	O + O3 → 2 O2	8.0*10–12*(–2060/T)
	O(1D) + O3 → 2 O2	1.2*10–10
	O(1D) + O3 → O2 + 2 O	1.2*10–10
	O(1D) + N2 → O + N2	2.15*10–11*e(110/T)
	O(1D) + O2 → O + O2	3.3*10–11*e(55/T)
Extended Chapman Model	OH + O3 → HO2 + O2	1.7 × 10–12e(–940/T)
	HO2 + O3 → OH + 2O2	1.0 × 10–12e(–490/T)
	OH + HO2 → H2O + O2	4.8 × 10–11e(250/T)
	O(1D) + H2O → 2OH	1.63 × 10–10e(60/T)
	O(1D) + N2O → N2 + O2	0.39 × 1.19 × 10–10e(20/T)
	O(1D) + N2O → 2NO	0.61 × 1.19 × 10–10e(20/T)
	O + NO2 → NO + O2	5.1 × 10–12e(210/T)
	NO + O3 → NO2 + O2	3.0 × 10–12e(–1500/T)
Extended Chapman Model	NO2 + O3 → NO2 + O2	1.2 × 10–13e(–2450/T)
	NO2 + NO3 + M → N2O5 + M	k0 = 2.4 × 10–30 × (T/300)–3
		k∞ = 4.0 × 10–12 × (T/300)0.1
	NO2 + OH + M → HNO3 + M	k0 = 1.8 × 10–30 × (T/300)–3
		k∞ = 2.8 × 10–11
	HNO3 + OH → NO3 + H2O	k0 = 1.8 × 10–30 × (T/300)–3
		k∞ = 2.8 × 10–11
	HO2 + NO → OH + NO2	3.3 × 10–12e(270/T)
	H + O3 → OH + O2	1.4 × 10–10e(–470/T)
	O + OH → H + O2	1.8 × 10–11e(180/T)
	H + O2 + M → HO2 + M	k0 = 4.4 × 10–32 × (T/300)–1.3c
		k∞ = 7.5 × 10–11 × (T/300)0.2
	O + HO2 → OH + O2	3.0 × 10–11e(200/T)
	H + HO2 → 2OH	7.2 × 10–11
	H + HO2 → O + H2O	1.6 × 10–12
	H + HO2 → H2 + O2	5.0 × 10–12

^a Low Pressure Limit, ^b termolecular reaction has [cm⁶ molecules^–2 s^–1] units.

Table 6. List of photolysis reactions used for in the in the Chapman model in Subsection “Present Earth”

Model	Reaction	Absorption Cross-section (cm2)	Quantum Yield
Chapman Model	O3 + hν → O + O2	180–185 nm[1], 185–232 nm[2] 232–309 nm[3], >309 nm[4]	1.0 – Φ(O(1D))[5]
	O3 + hν → O(1D) + O2
	O2 + hν → 2O	180–203 nm[6], 205–245 nm[5] 246–389 nm[7], 650–788.6 nm[8]	<242 nm: 1.0 >242 nm: 0.0
Extended Chapman Model	H2O + hν → H + OH	180–198 nm[5] >198 nm = 0	1.0
	HO2 + hν → O + OH	180–190 nm Linear Interpolation 190–260 nm [9] >260 nm = 0	1.0
	N2O + hν → N2+ O(1D)	180–240 nm[10] 240–300 nm[11] >300 nm = 0	<230 nm: 1.0 >230 nm: 0.0
	NO2 + hν → O + NO	177.12 nm[12] 185–200 nm[13] 200–238 nm[14] >238 nm[15]	[5]
	NO3 + hν → NO + O2 NO3 + hν → O + NO2	<400 nm = 0 400–691 nm[16] >691 nm[17]	[20]   [5]
	N2O5 + hν → NO2 + NO3 N2O5 + hν → O + NO + NO3	180–198 nm[18] 200–420 nm[5] >420 nm = 0	[21]   [22]
	HNO3 + hν → OH + NO2 HNO3 + hν → H + NO3	180–183 nm[19] 185–350 nm[20] >350 nm: 0	[5]

^[1]Ackerman (1971), ^[2]Molina and Molina (1986), ^[3]Burrows et al. (1999), ^[4]Brion et al. (1998), ^[5]Orkin et al. (2015), ^[6]Yoshino et al. (1987), ^[7]Bogumil et al. (2003), ^[8]Greenblatt et al. (1990), ^[9]Tyndall et al. (2001), ^[10]Selwyn et al. (1977), ^[11]Johnston and Selwyn (1975), ^[12]Au and Brion (1997), ^[13]Bass et al. (1976), ^[14]Olive (2011), ^[15]Vandaele et al. (1998), ^[16]Yokelson et al. (1994), ^[17]Orphal et al. (2003), ^[18]Osborne et al. (2000), ^[19]Suto and Lee (1984), ^[20]Burkholder et al. (1993), ^[21]Harwood et al. (1998), ^[22]Ravishankara et al. (1986).

The results of the simulation were compared to the vertical distributions of O₃ mixing ratios and previous model studies. Figure 14 shows O₂ and O₃ photodissociation rate constants compared to simulations carried out by KROME and the values reported by (Brasseur and Solomon, 2005). Both models show a good level of agreement with the literature data. The difference between KROME and PATMO is, as it was explained in the previous section, the calculation opacity terms at each layer and the approximation for the transport term used for KROME. Figure 15 shows simulated ozone mixing ratios for both codes and the US Standard atmosphere of 1976. Both models show a level of agreement within a factor 2. The differences between the observations and the model are attributed to O₃ depletion by NOx and HOx species.

Fig. 14.

Calculated vertical profiles of photodissociation rate constants for O₃ (dashed lines) and O₂ (solid lines) for a simple Chapman mechanism presented in Subsection “Present Earth”. The calculations presented in this paper (red lines) are compared to simulations with KROME (green lines) and literature values presented in Brasseur and Solomon (2005) (black lines).

Fig. 15.

Altitude profile of O₃ in today’s Earth atmosphere predicted by PATMO (green lines) and KROME (red lines) with a simple Chapman mechanism compared to the US Standard Atmosphere 1976 for mid-latitudes (black lines).

In order to address these inaccuracies, and to test the ability of the model to accurately compute chemistry of short-lived species, the chemical network was expanded to include NOx and HOx chemistry. Table 5 shows an extended version of the Chapman mechanism which includes a chemical reaction network of 25 reactions and 13 photoreactions (Table 6). For this simulation O₃ depletion reactions with NOx and HOx radicals were added to the chemical network. Accurate mixing ratios of such radicals require H₂O, HNO₃, and N₂O mixing ratios which were added as a fixed boundary condition (See Fig. A. 6).

This simulated O₃ profile is compared in Fig. 16 with previous simulations and the US standard atmosphere of 1976 (United States. National Oceanic and Atmospheric Administration and United States Committee on Extension to the Standard Atmosphere, 1976). The additional chemistry adds the catalytic depletion of O₃ reducing the discrepancies with the standard atmosphere. The simulated mixing ratios below 10 km shows a significant discrepancy, coming from tropospheric air pollution, where O₃ is a byproduct of emissions of NOx and various reactive hydrocarbons.

Fig. 16.

Altitude profile of O₃ in today’s Earth atmosphere predicted by the PATMO code with an extended Chapman mechanism (red line) compared to the US Standard Atmosphere 1976 for mid-latitudes (black line) and a simple Chapman mechanism (green line). The extended Chapman mechanism includes NOx and OHx radicals which are key to reproduce field measurements. The large discrepancies below 10 km are expected to be created by anthropogenic emissions of chemical pollutants.

One of the most relevant features of PATMO is the treatment of solar fluxes and how changes in the concentration of trace species affect the solar flux and the photo-kinetics of the middle atmosphere. The ozone/oxygen system presents an ideal case to these the ability of the model to reproduce solar fluxes at different altitudes. Figure 17 shows the actinic flux at 10, 20, 30, 40, 50 km, and at the top of the atmosphere (above 90 km), calculated with the extended Chapman mechanism (solid lines). Modeled fluxes were compared to measured solar spectra at the top pf the atmosphere reported by Gueymard (2004) (solid back line), and modeled fluxes by McLinden et al. (2002) (dotted lines). The high energy section of the spectra shows a clear structure of oxygen’s Schumann-Runge bands. The differences in literature and modeled data can be attributed to lower spectral resolution and the comparison limit of 190 nm. The modeled 210–220 nm solar window reproduces the literature data with a satisfactory level of accuracy. The 230–290 nm range shows the ozone’s Hartley absorption band that absorbs most of the solar flux, while the low energy shoulder (above 300 nm) represents where the solar flux penetrates deeply into the atmosphere. For a detailed analysis of the trace species and solar flux evolutions during the simulation see additional information section D.

Fig. 17.

Actinic flux at 10, 20, 30, 40, 50 km and the top of the atmosphere (extraterrestrial UV irradiance as reported by Gueymard, 2004) calculated with the extended Chapman mechanism (solid lines) and compared to solar spectra reported by Gueymard (2004) (solid back line) and modeled fluxes by McLinden et al. (2002) (dotted lines).

Figure 18 compares modeled photolysis rate constants of O₂ and O₃ to literature data, showing an excellent agreement for both molecules. A slight deviation between the data in this report and the literature can be observed for the O₂ rates at 60–90 km. The source of this difference can be attributed to photo-active chemical species with absorption features in the Schumann-Runge bands that are not accounted for, the computed spectral rang, or spectral resolution of O₂. The photolysis induced isotopic effect on oxygen, as discussed in Subsection “Photodissociation rates and isotopic effects, changes in spectral resolution”, could be added to the reaction network used in this section to draw conclusions about the fate of the isotopic signal within the atmospheric system. To is necessary fractionation factors for all involved reactions. Furthermore, the high reactivity of O(³P) atom and fast atmospheric isotopic exchange of oxygen containing compounds are needed. Although beyond the scope of this paper the PATMO model is capable tool for such analysis. For a more detailed analysis of the calculated photodissociation rate constants, calculated values of trace species NO₂, N₂O, HNO₃, H₂O, NO₃, and N₂O₅ are presented in Fig. 19. The modeled values show a good agreement with the literature data, except for the HNO₃ rates, that present a difference of about a factor of 10 throughout the entire atmosphere. The most likely explanation for this difference is the difficulty to determine accurate UV absorption cross-sections and the molecular-aerosol ratio of HNO₃.

Fig. 18.

Altitude profile of the photodissociation rate constants simulated with the PATMO code (solid lines) implemented with the extended Chapman mechanism and compared to literature data (Brasseur and Solomon, 2005) for O₂ and O₃ photolysis.

Fig. 19.

Altitude profile of the photodissociation rate constants simulated with the PATMO code implemented with the extended Chapman mechanism (solid lines) and compared to literature data (Brasseur and Solomon, 2005) for H₂O, N₂O, HNO₃, N₂O₅, NO₂, and NO₃ photolysis.

Summary

We reported a benchmark study of a chemical solver coupled to a diffusion algorithm for the study of planetary atmospheres. This work extends the work of Ávila et al. (2021) and focuses on different chemical networks and kinetics associated to planetary atmospheres. The study of photodissociation and photo-induced atmospheric chemistry is one of the main objectives of this paper, as well as the calculation of photodissociation rate constants at different altitudes and under different chemical conditions. The code is designed in a way that the spectral resolution used in equations 2 to 4 can be adjusted to have a trade-off between accuracy and computational cost, therefore, the code is suited to study of photo-induced isotopic effects. Subsections “Photodissociation rates and isotopic effects, changes in spectral resolution” and “Photodissociation of SO₂-SO system and their associated isotopic effects” present in detail the effects of spectral resolution on the oxygen photo-induced isotopic effect and the shielding/self-shielding effects at different altitudes for the SO₂/SO system.

The structure of the code allows the addition or removal of chemical reactions as required by specific studies. Subsections “Benchmark calculation”, “Photodissociation of a single chemical species at different altitudes”, “Photodissociation rates and isotopic effects, changes in spectral resolution”, and “Photodissociation of SO₂-SO system and their associated isotopic effects” present study cases with small chemical networks that are not representative of any specific atmosphere but play an important role in calibrating and bench marking the code. Subsection “Present Earth” concentrates on the ability of the code to reproduce the chemical processes of a realistic atmosphere. Starting with a simplified Chapman mechanism the model manages to replicate the basic features of the ozone distribution in today’s Earth atmosphere. With the addition of NOx and HOx species the model successfully reproduces the ozone distribution and the solar actinic flux at different altitudes. This study case exemplifies the main purpose and future applicability of the code, and the study of how chemical reactions affect the chemical composition of planetary atmospheres, with the unique capability of including detailed isotopologues effects and interplay.

Appendix A

Formulation of vertical diffusion flux

The vertical flux operator of the species (Φ_i,j) is derived from the general equation 1 or its explicit form

  

$$\begin{array}{cc}{\Phi }_{i,j}=& -K_j\left[\frac{\partial n_{i,j}}{\partial z}+n_{i,j}\left(\frac{1}{H_0}+\frac{1}{T_j}\frac{dT}{dz}\right)\right]\hfill \\ & -D_j\left[\frac{\partial n_{i,j}}{\partial z}+n_{i,j}\left(\frac{1}{H_{i,j}}+\frac{1+{\alpha }_T}{T_j}\frac{dT}{dz}\right)\right]\hfill \end{array}$$(7)

where K [cm² s^–1] is the turbulent eddy diffusion coefficient along altitude (z-axis), D [cm² s^–1] is the molecular diffusion coefficient, H₀ [cm] is the mean scale height of the atmosphere, H_i [cm] is the molecular scale height for molecule i, T the temperature at a given altitude layer, n_i the number density of a given species, and α_T is the thermal diffusion factor. The physical interpretation of the Eddy diffusion, also known as turbulent diffusion, is that pushes toward an equal mixing ratio curve while the molecular diffusion with its random motion tend to produce a uniform concentration. The right side of Eqn. 2 represents the Eddy diffusion term (first term) and the molecular diffusion term (second term).

The number of moles n_i can be substituted by the molar fraction f_i,j = n_i,j/N_j where N_j is the total number of moles present in the j^th layer, by using the scale height definition (dP/P = –dz/H) the Eddy diffusion term in Eqn. 2 can be written in the following form (i and j are abbreviated)

  

$${\mathrm{\Phi }}_K=-KN\frac{\partial f}{\partial z}$$(8)

The molecular diffusion term takes the following form

  

$${\mathrm{\Phi }}_D=-DN\frac{\partial f}{\partial z}+Dn\left(\frac{1}{H_0}-\frac{1}{H}-\frac{{\alpha }_T}{T}\frac{dT}{dz}\right)$$(9)

From the above the vertical transport term becomes

  

$$\mathrm{\Phi }=-KN\frac{\partial f}{\partial z}-DN\frac{\partial f}{\partial z}+Dn\left(\frac{1}{H_0}-\frac{1}{H}-\frac{{\alpha }_T}{T}\frac{dT}{dz}\right)$$(10)

In order to numerically solve Eqn. 1 the vertical transport term is discretized in the following form

  

$$\frac{\partial n_{i,j}}{\partial t}=P_{i,j}-n_{i,j}{L}_{i,j}-\frac{({\mathrm{\Phi }}_{i,j+1/2}-{\mathrm{\Phi }}_{i,j-1/2})}{∆z}$$(11)

From Eqn. 5 the terms Φ_i,j+1/2, Φ_i,j–1/2 are constructed and constructed:

  

$$\begin{array}{cc}{\mathrm{\Phi }}_{i,j+1/2}=& -\left(K_{j+1/2}+D_{j+1/2}\right)N_{j+1/2}\frac{f_{i,j+1}-f_{i,j}}{\mathrm{\Delta }z}\hfill \\ & +D_{j+1/2}{N}_{j+1/2}\frac{f_{i,j+1}+f_{i,j}}{2}\left(\frac{1}{H_0}-\frac{1}{H_i}-\frac{{\alpha }_T}{T_{j+1/2}}\frac{T_{j+1}-T_j}{\mathrm{\Delta }z}\right)\hfill \end{array}$$(12)

If m_a and m_i are defined as the atmospheric average molecular weight and the molecular weight of molecule i,

  

$$\frac{1}{H_0}-\frac{1}{H_i}=\frac{m_{a}g}{k_{B}T}-\frac{m_{i}g}{k_{B}T}=\frac{(m_a-m_i)g}{k_{B}T}$$(13)

With this definition Eqn. 7 takes the form:

  

$$\begin{array}{cc}{\mathrm{\Phi }}_{i,j+1/2}=\hfill & -\left(K_{j+1/2}+D_{j+1/2}\right)N_{j+1/2}\frac{f_{i,j+1}-f_{i,j}}{\mathrm{\Delta }z}\\ & +D_{j+1/2}{N}_{j+1/2}\frac{f_{i,j+1}+f_{i,j}}{2}\left[\frac{(m_a-m_i)g}{k_B{T}_{j+1/2}}\right.\hfill \\ & \left.-\frac{{\alpha }_T}{T_{j+1/2}}\frac{T_{j+1}-T_j}{\mathrm{\Delta }z}\right]\hfill \end{array}$$(7')

If the same procedure is applied to Φ_i,j–1/2 then:

  

$$\begin{array}{cc}{\mathrm{\Phi }}_{i,j-1/2}=& -\left(K_{j-1/2}+D_{j-1/2}\right)N_{j-1/2}\frac{f_{i,j}-f_{i,j-1}}{\mathrm{\Delta }z}\hfill \\ & +D_{j-1/2}{N}_{j-1/2}\frac{f_{i,j-1}+f_{i,j}}{2}\left[\frac{(m_a-m_i)g}{k_B{T}_{j-1/2}}\right.\hfill \\ & \left.-\frac{{\alpha }_T}{T_{j-1/2}}\frac{T_{j-1}-T_j}{\mathrm{\Delta }z}\right]\hfill \end{array}$$(14)

From the obtained Φ_i,j+1/2, Φ_i,j–1/2 terms are inserted into Eqn. 6 and after re-arrange Eqn. 10 is obtained

  

$$\begin{array}{cc}\frac{\partial n_{i,j}}{\partial t}=& P_{i,j}-n_{i,j}{L}_{i,j}+\left(K_{j+1/2}+D_{j+1/2}\right)N_{j+1/2}\frac{f_{i,j+1}-f_{i,j}}{{\left(\mathrm{\Delta }z\right)}^2}\hfill \\ & {-D}_{j+1/2}\frac{N_{j+1/2}}{2\mathrm{\Delta }z}\left[\frac{\left(m_a-m_i\right)g}{k_B{T}_{j+1/2}}-\frac{{\alpha }_T}{T_{j+1/2}}\frac{T_{j+1}-T_j}{\mathrm{\Delta }z}\right]\left(f_{i,j+1}+f_{i,j}\right)\hfill \\ & -\left(K_{j-1/2}+D_{j-1/2}\right)N_{j-1/2}\frac{f_{i,j}-f_{i,j-1}}{{\left(\mathrm{\Delta }z\right)}^2}+D_{j-1/2}\frac{N_{j-1/2}}{2\mathrm{\Delta }z}\hfill \\ & \left[\frac{(m_a-m_i)g}{k_B{T}_{j-1}}-\frac{{\alpha }_T}{T_{j-1/2}}\frac{T_{j-1}-T_j}{\mathrm{\Delta }z}\right]\left(f_{i,j-1}+f_{i,j}\right)\hfill \end{array}$$(15)

here k_{j ± 1/2} and d_{i, j ± 1/2} are defined as follows

  

$$k_{j\pm 1/2}=\frac{K_{j\pm 1/2}+D_{j\pm 1/2}}{\mathrm{\Delta }{z}^2}$$(16)

  

$$d_{i,j\pm 1/2}=\frac{D_{j\pm 1/2}}{2\mathrm{\Delta }{z}^2}\left[\frac{\left(m_a-m_i\right)g}{k_B{T}_{j\pm 1/2}}∆z-\frac{{\alpha }_T}{T_{j\pm 1/2}}\left(T_{j\pm 1}-T_j\right)\right]$$(17)

By using k_{j ± 1/2} and d_{i, j ± 1/2}, them the 3^rd to 6^th term in Eqn. 10 take the form

3^rd term

  

$$\begin{array}{c}\left(K_{j+1/2}+D_{j+1/2}\right)N_{j+1/2}\frac{f_{i,j+1}-f_{i,j}}{{\left(\mathrm{\Delta }z\right)}^2}\hfill \\ =\frac{K_{j+1/2}+D_{j+1/2}}{{∆z}^2}{N}_{j+1/2}\left(f_{i,j+1}-f_{i,j}\right)\hfill \\ =k_{j+1/2}{N}_{j+1/2}\left(f_{i,j+1}-f_{i,j}\right)\hfill \end{array}$$(18)

4^th term

  

$$\begin{array}{c}{D}_{j+1/2}\frac{N_{j+1/2}}{2\mathrm{\Delta }z}\left[\frac{(m_a-m_i)g}{k_B{T}_{j+1/2}}-\frac{{\alpha }_T}{T_{j+1/2}}\frac{T_{j+1}-T_j}{\mathrm{\Delta }z}\right]\left(f_{i,j+1}+f_{i,j}\right)\hfill \\ =D_{j+1/2}\frac{N_{j+1/2}}{2{\left(∆z\right)}^2}\left[\frac{\left(m_a-m_i\right)g\mathrm{\Delta }z}{k_B{T}_{j+1/2}}-\frac{{\alpha }_T}{T_{j+1/2}}\left(T_{j+1}-T_j\right)\right]\left(f_{i,j+1}+f_{i,j}\right)\hfill \\ =d_{i,j+1/2}{N}_{j+1/2}\left(f_{i,j+1}+f_{i,j}\right)\hfill \end{array}$$(19)

5^th term

  

$$\begin{array}{c}\left(K_{j-1/2}+D_{j-1/2}\right)N_{j-1/2}\frac{f_{i,j}-f_{i,j-1}}{{\left(\mathrm{\Delta }z\right)}^2}\hfill \\ =k_{j-1/2}{N}_{j-1/2}\left(f_{i,j}-f_{i,j-1}\right)\hfill \end{array}$$(20)

6^th term

  

$$\begin{array}{c}{D}_{j-1/2}\frac{N_{j-1/2}}{2\mathrm{\Delta }z}\left[\frac{\left(m_a-m_i\right)g}{k_B{T}_{j-1}}-\frac{{\alpha }_T}{T_{j-1/2}}\frac{T_{j-1}-T_j}{\mathrm{\Delta }z}\right]\left(f_{i,j-1}+f_{i,j}\right)\hfill \\ =d_{i,j-1/2}{N}_{j-1/2}\left(f_{i,j-1}+f_{i,j}\right)\hfill \end{array}$$(21)

By introducing Eqn. 13–16 into Eqn. 10 it can be re-arranged into the form:

  

$$\begin{array}{cc}\frac{\partial n_{i,j}}{\partial t}=& P_{i,j}-n_{i,j}{L}_{i,j}\hfill \\ & +k_{j+1/2}{N}_{j+1/2}\left(f_{i,j+1}-f_{i,j}\right){-d}_{i,j+1/2}{N}_{j+1/2}\left(f_{i,j+1}+f_{i,j}\right)\hfill \\ & -k_{j-1/2}{N}_{j-1/2}\left(f_{i,j}-f_{i,j-1}\right)+d_{i,j-1/2}{N}_{j-1/2}\left(f_{i,j-1}+f_{i,j}\right)\hfill \end{array}$$(22)

Eqn. 17 after rearranging in terms of f_{i, j + 1}, f_{i, j}, f_{i, j – 1} take the form of Eqn. 18

  

$$\begin{array}{cc}\frac{\partial n_{i,j}}{\partial t}=& P_{i,j}-n_{i,j}{L}_{i,j}+f_{i,j+1}\left(k_{j+1/2}{N}_{j+1/2}-d_{i,j+1/2}{N}_{j+1/2}\right)\hfill \\ & -f_{i,j}\left(k_{j+1/2}{N}_{j+1/2}{+d}_{i,j+1/2}{N}_{j+1/2}+k_{j-1/2}{N}_{j-1/2}-d_{i,j-1/2}{N}_{j-1/2}\right)\hfill \\ & +f_{i,j-1}(k_{j-1/2}{N}_{j-1/2}+d_{i,j-1/2}{N}_{j-1/2})\hfill \end{array}$$(23)

If the mixing ratios in Eqn. 18 are converted into number density

  

$$\begin{array}{cc}\frac{\partial n_{i,j}}{\partial t}=& P_{i,j}-n_{i,j}{L}_{i,j}+n_{i,j+1}\left(k_{j+1/2}\frac{N_{j+1/2}}{N_{j+1}}-d_{i,j+1/2}\frac{N_{j+1/2}}{N_{j+1}}\right)\hfill \\ & -n_{i,j}\left(k_{j+1/2}\frac{N_{j+1/2}}{N_j}{+d}_{i,j+1/2}\frac{N_{j+1/2}}{N_j}+k_{j-1/2}\frac{N_{j-1/2}}{N_j}-d_{i,j-1/2}\frac{N_{j-1/2}}{N_j}\right)\hfill \\ & +n_{i,j-1}\left(k_{j-1/2}\frac{N_{j-1/2}}{N_{j-1}}+d_{i,j-1/2}\frac{N_{j-1/2}}{N_{j-1}}\right)\hfill \end{array}$$(24)

As above presented, after discretization and rearranging Eqn. 1 is transformed into Eqn. 19. Equation 19 calculate the chemical concentration of all chemical species at all layers with the transport influence from the upper and lower layers. The boundary conditions for the upper and lower layer are k_{j + 1/2} = 0, d_{i, j + 1/2} = 0, k_{j – 1/2} = 0, d_{i, j – 1/2} = 0 respectively.

Acknowledgments

The authors would like to express their gratitude to the institutions supporting this study and the people that provided valuable comments. This research is founded by Sophia University Special Grant for Academic Research and JSPS KAKENHI Grant Numbers 17H06105, 20H01975, 26108511 and 22H05150.

References

•  Ackerman,  M. (1971) Ultraviolet solar radiation related to mesospheric processes. Mesospheric Models and Related Experiments (Fiocco, G., ed.), 149–159, Springer Netherlands: Dordrecht. https://doi.org/10.1007/978-94-010-3114-1_11
•  Au,  J. W. and  Brion,  C. E. (1997) Absolute oscillator strenghts for the valence-shell photoabsorption (2–200 eV) and the molecular and dissociative photoionization (11–80 eV) of nitrogen dioxide. Chem. Phys. 218, 109–126. https://doi.org/10.1016/S0301-0104(97)00065-7
•  Ávila,  P.,  Grassi,  T.,  Bovino,  S.,  Chiavassa,  A.,  Ercolano,  B.,  Danielache,  S. O. and  Simoncini,  E. (2021) Presence of water on exomoons orbiting free-floating planets: A case study. Int. J. Astrobiology 20, 300–311. https://doi.org/10.1017/S1473550421000173
•  Barge,  L. M.,  Branscomb,  E.,  Brucato,  J. R.,  Cardoso,  S. S. S.,  Cartwright,  J. H. E.,  Danielache,  S. O.,  Galante,  D.,  Kee,  T. P.,  Miguel,  Y.,  Mojzsis,  S.,  Robinson.  K. J.,  Russell,  M. J.,  Simoncini,  E. and  Sobron,  P. (2017) Thermodynamics, disequilibrium, evolution: Far-from-equilibrium geological and chemical considerations for origin-of-life research. Orig. Life Evol. Biosph. 47, 39–56. https://doi.org/10.1007/s11084-016-9508-z
•  Bass,  A. M.,  Ledford,  A. E. Jr. and  Laufer,  A. H. (1976) Extinction Coefficients of NO₂ and N₂O₄. J. Res. Natl. Bur. Stand. A. Phys. Chem. 80A, 143–166. https://doi.org/10.6028/jres.080A.017
•  Bogumil,  K.,  Orphal,  J.,  Homann,  T.,  Voigt,  S.,  Spietz,  P.,  Fleischmann,  O. C.,  Vogel,  A.,  Hartmann,  M.,  Kromminga,  H.,  Bovensmann,  H.,  Frerick,  J. and  Burrows,  J. P. (2003) Measurements of molecular absorption spectra with the SCIAMACHY pre-flight model: Instrument characterization and reference data for atmospheric remote-sensing in the 230–2380 nm region. J. Photochem. Photobiol. Chem. 157, 167–184. https://doi.org/10.1016/S1010-6030(03)00062-5
•  Brasseur,  G. P. and  Solomon,  S. (2005) Aeronomy of the Middle Atmosphere, Springer Dordrecht. 646 pp. https://doi.org/10.1007/1-4020-3824-0
•  Brion,  J.,  Chakir,  A.,  Charbonnier,  J.,  Daumont,  D.,  Parisse,  C. and  Malicet,  J. (1998) Absorption spectra Measurements for the Ozone Molecule in the 350–830 nm Region. J. Atmospheric Chem. 30, 291–299. https://doi.org/10.1023/A:1006036924364
•  Burkholder,  J. B.,  Talukdar,  R. K.,  Ravishankara,  A. R. and  Solomon,  S. (1993) Temperature dependence of the HNO₃ UV absorption cross sections. J. Geophys. Res. Atmospheres 98, 22937–22948. https://doi.org/10.1029/93JD02178
•  Burrows,  J. P.,  Richter,  A.,  Dehn,  A.,  Deters,  B.,  Himmelmann,  S.,  Voigt,  S. and  Orphal,  J. (1999) Atmospheric remote-sencing reference data from GOME-2. Temperature-dependent absorption cross sections of O3 in the 231–794 nm Range. J. Quant. Spectrosc. Radiat. Transf. 61, 509–517. https://doi.org/10.1016/S0022-4073(98)00037-5
•  Catling,  D. C. and  Kasting,  J. F. (2017) Atmospheric Evolution on Inhabited and Lifeless Worlds. Planetary systems and astrobiology, Cambridge University Press, 592 pp.
•  Cicerone,  R. J. and  McCrumb,  J. L. (1980) Photodissociation of Isotopically Heavy O₂ as a source of O₃. Geophys. Res. Lett. 7, 251–254.
•  Claire,  M. W.,  Kasting,  J. F.,  Domagal-Goldman,  S. D.,  Stüeken,  E. E.,  Buick,  R. and  Meadows,  V. S. (2014) Modeling the signature of sulfur mass-independent fractionation produced in the Archean atmosphere. Geochim. Cosmochim. Acta 141, 365–380. https://doi.org/10.1016/j.gca.2014.06.032
•  Danielache,  S. O.,  Tomoya,  S.,  Kondorsky,  A.,  Tokue,  I. and  Nanbu,  S. (2014) Nonadiabatic calculations of ultraviolet absorption cross section of sulfur monoxide: Isotopic effects on the photodissociation reaction. J. Chem. Phys. 140, 044319. https://doi.org/10.1063/1.4862429
•  Endo,  Y.,  Danielache,  S. O.,  Ueno,  Y.,  Hattori,  S.,  Johnson,  M. S.,  Yoshida,  N. and  Kjaergaard,  H. G. (2015) Photoabsorption cross-section measurements of ³²S, ³³S, ³⁴S, and ³⁶S sulfur dioxide from 190 to 220 nm. J. Geophys. Res. Atmospheres 120, 2546–2557. https://doi.org/10.1002/2014JD021671
•  Endo,  Y.,  Danielache,  S. O.,  Ogawa,  M. and  Ueno,  Y. (2022) Absorption spectra measurements at ~1 cm^–1 spectral resolution of 32S, 33S, 34S, and 36S sulfur dioxide for the 206–220 nm region and applications to modeling of the isotopic self-shielding. Geochem. J. 56, 40–56. https://doi.org/10.2343/geochemj.GJ22004
•  Giorgi,  F. and  Chameides,  W. L. (1985) The rainout parameterization in a photochemical model. J. Geophys. Res. Atmospheres 90, 7872–7880. https://doi.org/10.1029/JD090iD05p07872
•  Grassi,  T.,  Bovino,  S.,  Schleicher,  D. R. G.,  Prieto,  J.,  Seifried,  D.,  Simoncini,  E. and  Gianturco,  F. A. (2014) KROME-A Package to Embed Chemistry in Astrophysical Simulations. Mon. Not. R. Astron. Soc. 439, 2386–2419. https://doi.org/10.1093/mnras/stu114.
•  Greenblatt,  G. D.,  Orlando,  J. J.,  Burkholder,  J. B. and  Ravishankara,  A. R. (1990) Absorption measurements of oxygen between 330 and 1140 nm. J. Geophys. Res. Atmospheres 95, 18577–18582. https://doi.org/10.1029/JD095iD11p18577
•  Gueymard,  C. A. (2004) The sun’s total and spectral irradiance for solar energy applications and solar radiation models. Sol. Energy 76, 423–453. https://doi.org/10.1016/j.solener.2003.08.039
•  Harwood,  M. H.,  Burkholder,  J. B. and  Ravishankara,  A. R. (1998) Photodissociation of BrONO₂ and N₂O₅: Quantum yields for NO₃ production at 248, 308, and 352.5 nm. J. Phys. Chem. A 102, 1309–1317. https://doi.org/10.1021/jp9729829
•  Hattori,  S.,  Schmidt,  J. A.,  Johnson,  M. S.,  Danielache,  S. O.,  Yamada,  A.,  Ueno,  Y. and  Yoshida,  N. (2013) SO₂ photoexcitation mechanism links mass-independent sulfur isotopic fractionation in cryospheric sulfate to climate impacting volcanism. Proc. Natl. Acad. Sci. U.S.A. 110, 17656–17661. https://doi.org/10.1073/pnas.1213153110
•  Hindmarsh,  A. C. (1983) A systematized collection of ODE solvers, scientific computing. IMACS Transactions on Scientific Computation 1, 55–56.
•  Hu,  R.,  Seager,  S. and  Bains,  W. (2012) Photochemistry in terrestrial exoplanet atmospheres. I: Photochemistry model and benchmark cases. Astrophys. J. 761, 166. https://doi.org/10.1088/0004-637X/761/2/166
•  Izon,  G.,  Zerkle,  A. L.,  Williford,  K. H.,  Farquhar,  J.,  Poulton,  S. W. and  Claire,  M. W. (2017) Biological regulation of atmospheric chemistry in route to planetary oxygenation. Proc. Natl. Acad. Sci. 114, E2571–E2579. https://doi.org/10.1073/pnas.1618798114.
•  Johnston,  H. S. and  Selwyn,  G. S. (1975) New cross sections for the absorption of near ultraviolet radiation by nitrous oxide (N₂O). Geophys. Res. Lett. 2, 549–551. https://doi.org/10.1029/GL002i012p00549
•  Kasting,  J. F. (1993) Earth’s early atmosphere. Science 259, 920–926. https://doi.org/10.1126/science.11536547
•  Kasting,  J. F.,  Whitmire,  D. P. and  Reynolds,  R. T. (1993) Habitable Zones around Main Sequence Stars. Icarus 101, 108–128. https://doi.org/10.1006/icar.1993.1010
•  Kawamiya,  M.,  Hajima,  T.,  Tachiiri,  K.,  Watanabe,  S. and  Yokohata,  T. (2020) Two decades of earth system modeling with an emphasis on Model for Interdisciplinary Research on Climate (MIROC). Prog. Earth Planet. Sci. 7, 64. https://doi.org/10.1186/s40645-020-00369-5
•  Massie,  S. T. and  Hunten,  D. M. (1981) Stratospheric eddy diffusion coefficients from tracer data. J. Geophys. Res. Oceans 86, 9859–9868. https://doi.org/10.1029/JC086iC10p09859
•  McBride,  B. J.,  Gordon,  S. and  Reno,  M. A. (1993) Coefficients for Calculating Thermodynamic and Transport Properties of Individual Species; NASA-TM-4513.
•  McLinden,  C. A.,  McConnell,  J. C.,  Griffioen,  E. and McElroy, C. T. A (2002) A vector radiative-transfer model for the Odin/OSIRIS project. Can. J. Phys. 80, 375–393. https://doi.org/10.1139/p01-156
•  Molina,  L. T. and  Molina,  M. J. (1986) Absolute absorption cross sections of ozone in the 185- to 350-Nm wavelength range. J. Geophys. Res. Atmospheres 91, 14501–14508. https://doi.org/10.1029/JD091iD13p14501
•  Morgenstern,  O.,  Hegglin,  M.,  Rozanov,  E.,  O’Connor,  F.,  Abraham,  N.,  Akiyoshi,  H.,  Archibald,  A.,  Bekki,  S.,  Butchart,  N.,  Chipperfield,  M.,  Deushi,  M.,  Dhomse,  S.,  Garcia,  R.,  Hardiman,  S.,  Horowitz,  L.,  Jöckel,  P.,  Josse,  B.,  Kinnison,  D.,  Lin,  M., and  Zeng,  G. (2017) Review of the global models used within the Chemistry-Climate Model Initiative (CCMI). Geosci. Model Dev. Discuss. 10, 639–671. https://doi.org/10.5194/gmd-2016-199
•  Nair,  H.,  Allen,  M.,  Anbar,  A. D.,  Yung,  Y. L. and  Clancy,  R. T. (1994) A photochemical model of the martian atmosphere. Icarus 111, 124–150. https://doi.org/10.1006/icar.1994.1137
• Olive, B. S. (2011) Absorption Spectra of Carbon Disulfide, Nitrogen Dioxide, and Anthracene; Results from measurements at the Department of Chemistry and Industrial Hygiene. University of North Alabama, Florence, AL.
•  Ono,  S. (2017) Photochemistry of sulfur dioxide and the origin of mass-independent isotope fractionation in earth’s atmosphere. Annu. Rev. Earth Planet. Sci. 45, 301–329. https://doi.org/10.1146/annurev-earth-060115-012324
• Orkin, V. L., Burkholder, J. B., Sander, S. P., Abbatt, J. P. D. A. D., Barker, J. R., Huie, R. E., Kolb, C. E., Iii, M. J. K., Wilmouth, D. M. and Wine, P. H. (2015) Chemical Kinetics and Photochemical Data for Use in Atmospheric Studies, Evaluation No. 18. JPL Publication 15-10.
•  Orphal,  J.,  Fellows,  C. E. and  Flaud,  P.-M. (2003) The visible absorption spectrum of NO₃ measured by high-resolution Fourier transform spectroscopy. J. Geophys. Res. Atmospheres 108, 4077. https://doi.org/10.1029/2002JD002489
•  Osborne,  B. A.,  Marston,  G.,  Kaminski,  L.,  Jones,  N. C.,  Gingell,  J. M.,  Mason,  N.,  Walker,  I. C.,  Delwiche,  J. and  Hubin-Franskin,  M. J. (2000) Vacuum ultraviolet spectrum of dinitrogen pentoxide. J. Quant. Spectrosc. Radiat. Transf. 64, 67–74. https://doi.org/10.1016/S0022-4073(99)00104-1
•  Pavlov,  A. A.,  Kasting,  J. F.,  Brown,  L. L.,  Rages,  K. A. and  Freedman,  R. (2000) Greenhouse warming by CH₄ in the atmosphere of early Earth. J. Geophys. Res. 105, 11981–11990. https://doi.org/10.1029/1999je001134
•  Pavlov,  A. A. and  Kasting,  J. F. (2002) Mass-independent fractionation of sulfur isotopes in Archean sediments: Strong evidence for an anoxic archean atmosphere. Astrobiology 2, 27–41. https://doi.org/10.1089/153110702753621321
•  Ravishankara,  A. R.,  Wine,  P. H.,  Smith,  C. A.,  Barbone,  P. E. and  Torabi,  A. (1986) N₂O₅ photolysis: Quantum yields for NO₃ and O(³P). J. Geophys. Res. Atmospheres 91, 5355–5360. https://doi.org/10.1029/JD091iD05p05355
•  Rimmer,  P. B. and Helling, C. A (2016) A chemical kinetics network for lightning and life in planetary atmospheres. Astrophys. J. Suppl. Ser. 224, 9. https://doi.org/10.3847/0067-0049/224/1/9
•  Selwyn,  G.,  Podolske,  J. and  Johnston,  H. S. (1977) Nitrous oxide ultraviolet absorption spectrum at stratospheric temperatures. Geophys. Res. Lett. 4, 427–430. https://doi.org/10.1029/GL004i010p00427
•  Sharp,  Z. (2007) Principles of Stable Isotope Geochemistry. Pearson/Prentice Hall, 344 pp.
•  Simoncini,  E.,  Virgo,  N. and  Kleidon,  A. (2013) Quantifying drivers of chemical disequilibrium: Theory and application to methane in the Earth’s atmosphere. Earth Syst. Dynam. 4, 317–331. https://doi.org/10.5194/esd-4-317-2013
•  Suto,  M. and  Lee,  L. C. (1984) Photoabsorption and photodissociation of HONO₂ in the 105–220 nm region. J. Chem. Phys. 81, 1294–1297. https://doi.org/10.1063/1.447816
•  Thiemens,  M. (2001) The mass-independent ozone isotope effect. Science 293, 226. https://doi.org/10.1126/science.1063648
•  Tyndall,  G. S.,  Cox,  R. A.,  Granier,  C.,  Lesclaux,  R.,  Moortgat,  G. K.,  Pilling,  M. J.,  Ravishankara,  A. R. and  Wallington,  T. J. (2001) Atmospheric chemistry of small organic peroxy radicals. J. Geophys. Res. Atmospheres 106, 12157–12182. https://doi.org/10.1029/2000JD900746
• United States. National Oceanic and Atmospheric Administration and United States Committee on Extension to the Standard Atmosphere (1976) U.S. Standard Atmosphere, 1976. National Oceanic and Amospheric [sic] Administration, 227 pp
•  Vandaele,  A. C.,  Hermans,  C.,  Simon,  P. C.,  Carleer,  M.,  Colin,  R.,  Fally,  S.,  Mérienne,  M. F.,  Jenouvrier,  A. and  Coquart,  B. (1998) Measurements of the NO₂ absorption cross-section from 42,000 cm^–1 to 10,000 cm^–1 (238–1000 nm) at 220 K and 294 K. J. Quant. Spectrosc. Radiat. Transf. 59, 171–184. https://doi.org/10.1016/S0022-4073(97)00168-4
•  Yokelson,  R. J.,  Burkholder,  J. B.,  Fox,  R. W.,  Talukdar,  R. K. and  Ravishankara,  A. R. (1994) Temperature dependence of the NO₃ absorption spectrum. J. Phys. Chem. 98, 13144–13150. https://doi.org/10.1021/j100101a009
•  Yoshino,  K.,  Freemen,  D. E.,  Esmond,  J. R. and  Parkinson,  W. H. (1987) High resolution absorption cross-sections and band oscillator strengths of the Schumann-Runge bands of oxygen at 79 K. Planet. Space Sci. 35, 1067–1075. https://doi.org/10.1016/0032-0633(87)90011-0
•  Yung,  Y.,  MA,  A. and  Pinto,  J. (1994) Photochemistry of the atmosphere of Titan - comparison between model and observations. Astrophys. J. Suppl. Ser. 55, 465–506. https://doi.org/10.1086/190963
•  Yung,  Y. L. and  Demore,  W. B. (1982) Photochemistry of the stratosphere of venus: Implications for atmospheric evolution. Icarus 51, 199–247. https://doi.org/10.1016/0019-1035(82)90080-X
• Yung, Y. L. and DeMore, W. B. (1999) Photochemistry of Planetary Atmospheres. Oxford University Press, 456 pp. https://doi.org/10.1093/oso/9780195105018.001.0001