Photochemical and radiation transport model for extensive use (PROTEUS)

We introduce a new flexible one-dimensional photochemical model named Photochemical and RadiatiOn Transport model for Extensive USe (PROTEUS), which consists of a Python graphical user interface (GUI) program and Fortran 90 modules. PROTEUS is designed for adaptability to many planetary atmospheres, for flexibility to deal with thousands of or more chemical reactions with high efficiency, and for intuitive operation with GUI. Chemical reactions can be easily implemented into the Python GUI program in a simple string format, and users can intuitively select a planet and chemical reactions on GUI. Chemical reactions selected on GUI are automatically analyzed by string parsing functions in the Python GUI program, then applied to the Fortran 90 modules to simulate with the selected chemical reactions on a selected planet. We performed a benchmark test of PROTEUS to confirm its validity, by applying it to the Martian atmosphere and the Jovian ionosphere. PROTEUS can significantly save the time for those who need to develop a new photochemical model; users just need to write chemical reactions in the Python GUI program and just select them on GUI to run a new photochemical model.


Introduction
One-dimensional photochemical and transport models are essential for investigating the typical vertical chemical structure of planetary atmospheres and their evolution throughout the history of the planets.They solve continuity-transport equations considering production and loss of each atmospheric species by numerous chemical reactions including photolysis.So far, plenty of photochemical models have been developed for various planetary atmospheres in the solar system (e.g., Kasting et al. 1979;Nair et al. 1994;Kim and Fox 1994;Fox and Sung 2001;Krasnopolsky 2009Krasnopolsky , 2012;;Chaffin et al. 2017) and terrestrial exoplanets (e.g., Rugheimer et al. 2013;Lustig-Yaeger et al. 2019;Grenfell et al. 2013).As the mass and spectral resolutions of measurements for detecting chemical species in planetary atmospheres increase and the theory of chemical kinetic systems becomes more complex, the need for photochemical models with hundreds or thousands of chemical reactions increases.
There are roughly three approaches to develop a numerical code for solving a lot of chemical reactions (Damian et al. 2002).The first approach is a hard-coding approach, in which the developer analyzes the chemical reactions, derives all the production and loss rate terms for each chemical species, and codes them into a program by hand.This approach is easy to develop when the number of chemical reactions is smaller than a hundred, however, it takes a lot of time to develop a code when the number of reactions becomes larger than hundreds or more.It is also difficult to add new chemical reactions into an already hard-coded program.
The second approach is a totally integrated approach, in which the chemical reactions are listed in a specific file in a certain format, and they are parsed by a program and stored in a memory at a run time.This approach is flexible in adding new chemical reactions after the development of the core program, and easy to deal with hundreds of chemical reactions without developer's manual derivation.This approach was used in the Atmos model in Fortran language for instance, which was originally developed by Kasting et al. (1979), updated by Zahnle et al. (2006) and recently described in Arney et al. (2016).Recently an integrated Martian photochemical model was developed by Chaffin et al. (2017) in Julia language with a more flexible way of describing reactions and rate coefficients.
The third approach is a preprocessing approach, in which the chemical reactions are listed in a specific file like the totally integrated approach but are parsed by a preprocessor to generate a hard-coded program in a high-level language such as Fortran or C language (Damian et al. 2002).This approach is also flexible in implementing hundreds of chemical reactions and is as efficient as the hard-coding approach.This approach was used in the kinetic preprocessor (KPP) originally developed by Damian et al. (2002), which has been widely used for chemical kinetic models for Earth's atmosphere (e.g., Roldin et al. 2019).
In this paper, we present a new integrated photochemical model named Photochemical and RadiatiOn Transport model for Extensive USe (PROTEUS), with a totally integrated approach.PROTEUS couples Python and Fortran modules, which is designed for adaptability to many planetary atmospheres, for flexibility to deal with thousands of or more chemical reactions with high efficiency, and for intuitive operation with a graphical user interface (GUI).A Python GUI program integrates a list of chemical reactions, GUI functions controlling the behavior of GUI operation, and a string parsing functions analyzing chemical reactions that output a Fortran 90 module.Fortran 90 modules solve differential equations numerically.Chemical reactions are written in a simple and flexible string format in the Python GUI program, making it easy to add new chemical reactions into the Python GUI program.The feature of PROTEUS that the Python GUI program outputs a Fortran 90 module is similar to the preprocessing approach, leading to a high efficiency, however, the Fortran modules in PROTEUS are not hard-coded but is rather generic.
PROTEUS has been newly developed and independent of other photochemical models or KPPs that have been developed so far.

Model descriptions
Equations PROTEUS is a one-dimensional photochemical model that solves a system of continuity equations for each species as follows: where n i is the number density of i th species, P i is the production rate of i th species, L i is the loss rate of i th species, z is the altitude and i is the vertical flux of i th species.The vertical flux i for both neutral and ionized species can be expressed as follows: where D i is the binary diffusion coefficient between i th species and the background atmosphere, H i = k B T i /m i g is the scale height of i th species, m i is the mass of i th (1) T e /T i P e species, k B is the Boltzmann constant, g is the gravita- tional acceleration, q i is the charge of i th species, q e is the elementary charge, T e and T i are the temperatures of electrons and i th species, respectively, P e = n e k B T e is the electron pressure, n e is the electron number density, α i is the thermal diffusion coefficient, K is the eddy diffusion coefficient, H = k B T /mg is the mean scale height of the background atmosphere, m is the mean molecular mass of the atmosphere, and T is the neutral temperature.
The temperature profiles are assumed to be stationary in time.The third term in Eq. ( 2) is the ambipolar diffusion term, which is applied only to charged species.The altitude-dependent gravitational acceleration g is calculated by using the mass and radius of the planet.We referred to the numerical method for solving the system of continuity equations for the photochemical model described in Catling and Kasting (2017) with some corrections.The vertical flux in Eq. ( 2) can be rewritten as follows: The spatial derivative of the vertical flux is expressed in a second-order central difference as follows: where superscripts j + 1/2 and j − 1/2 are the upper and lower interfaces of the j th vertical layer.The vertical fluxes at the interfaces j + 1/2 and j − 1/2 are expressed as: where (3.1) T e /T i P e (5.1) T e j+1/2 /T i j+1/2 P e j+1/2 P e j+1 − P e j Then the spatial derivative of the vertical flux can be described as follows: Equation (1) can then be written as follows: where A i j , B i j , C i j , and D i j are expressed as follows except for lower and upper boundaries: and D i 1 are expressed as follows: (6.2) T e j−1/2 /T i j−1/2 P e j−1/2 P e j − P e j−1 where v i 1/2 and � i 1/2 are the fixed velocity and flux of i th species at the lower boundary.Either v i 1/2 or � i 1/2 is given for each i th species as a lower boundary condition, or they are set to zero if not needed.If the lower boundary condition of the i th species is the constant density, , and D i 1 are set to zero.Note that the velocity and flux with positive values at the lower boundary conditions direct vertically upward.At the upper boundary, A i J , B i J , C i J , and D i J are expressed as follows: where J is the number of vertical layers, and v i J +1/2 and � i J +1/2 are the fixed velocity and flux of i th species at the upper boundary escaping the atmosphere.Either v i J +1/2 or � i J +1/2 is given for each i th species as an upper boundary condition, or they are set to zero if not needed.Note that the velocity and flux with positive values at upper boundary conditions also direct vertically upward.
The basic equations of the photochemical model are stiff equations in which some of the variables such as the number densities of short-lived species change more quickly than others.Therefore, PROTEUS has applied an implicit scheme described by Catling and Kasting (2017), in which the density and flux are solved simultaneously.Note that the production rate P i j and the loss rate L i j are functions of the number density.Jacobian components related to the production and loss rates are analyzed by the Python GUI program.The timestep can gradually increase from an initial value 10 -8 s to ideally more than 10 14 s, allowing us to investigate from a sporadic event response in several minutes to the evolution of planetary atmospheres in a billion-year time scale.
PROTEUS is basically a one-dimensional photochemical model and currently does not solve horizontal (11.1) transportation, but considers the rotation of a planet as options to obtain a simplified global distribution.PRO-TEUS has four options for the simulation geometry: (1) one-dimensional simulation at a given latitude at noon, (2) two-dimensional simulation at noon at each latitude from the north pole to the south pole, (3) two-dimensional simulation at a given latitude with rotation, and (4) three-dimensional simulation at all latitudes with rotation.The three-dimensional simulation has already been applied by Nakamura et al. (2022) for the Jovian ionosphere.Note that our model is fundamentally a onedimensional model that can only deal with the temporal variation of the solar zenith angle and is not able to solve horizontal transport, which is different from the state-ofthe-art three-dimensional photochemical and transport models (e.g., Braam et al. 2022;Cooke et al. 2023;Chen et al. 2018).

Radiative transfer
PROTEUS uses the solar EUV irradiance model for aeronomic calculations (EUVAC) (Richards et al. 1994) for the reference irradiance spectrum of solar extreme ultraviolet (EUV) flux to calculate the photoionization rates of atmospheric species.EUVAC model provides the solar EUV flux in 37 wavelength bins ranging from 5 to 105 nm.EUVAC model requires the input of F10.7 value and its 81 days running average value.High resolution solar EUV reference flux model such as the highresolution version of EUVAC (HEUVAC) (Richards et al. 2006) and the Flare Irradiance Spectral Model-Version 2 (FISM2) (Chamberlin et al. 2020) will be implemented into PROTEUS in the future.For calculating the photodissociation rates of atmospheric species, we used the reference irradiance spectrum of the solar flux in the wavelength range 0.05-2499.5nm taken from Woods et al. (2009).Adopting the solar flux taken from EUVAC model and Woods et al. (2009), the radiative transfer is solved by considering the absorption of the solar irradiation by atmospheric species.In PROTEUS, users can flexibly change the wavelength bin size for the solar irradiance of Woods et al. (2009) and absorption/dissociation cross sections of chemical species at each wavelength.The solar flux and cross section data can be provided in any wavelength bins, which are automatically interpolated linearly and binned to the wavelength bin given by the user.The automatically binning algorithm is especially useful when users need high resolution wavelength bins in a limited wavelength range.For example, if users need to resolve the Schumann-Runge bands of the oxygen molecule, the user can set a 0.01 nm resolution at 176-192.6 nm and 1 nm at other wavelength range, which could reduce the computational cost in solving radiation transfer and dissociation rate of atmospheric species even fully resolving the structured Schumann-Runge bands of the oxygen molecule.This algorithm is also useful for resolving a slight difference in absorption cross sections between isotopes (Yoshida, et al. 2023).
The reference solar irradiance spectra are then divided by the square of the distance r in the unit AU between the planet and the Sun.The distance r between the planet and the Sun at a given solar longitude L s is given by where r m is the mean distance in unit AU between the planet and the Sun, e is the eccentricity of the planetary orbit, and L s,P is the solar longitude at perihelion.The solar longitude is defined as the angle along the planet's orbit, starting from zero at the northern spring equinox, dividing the season from the northern spring equinox to the northern summer solstice into 90 degrees.The solar zenith angle at latitude θ and an hour angle η is given by where δ is solar declination.δ and η are given by where ε is the tilt angle of the rotational axis, t L is the time in second measured from the local noon, and T p is the rotational period of the planet.The optical depth τ (z, , χ) and the photon flux I(z, , χ) at the altitude z and the solar zenith angle χ as a function of wavelength can be calculated as follows: where n s is the number density of the s th species, σ a s is the absorption cross section of the s th species, I ∞ is the solar flux at the top of the atmosphere, and Ch(z, χ) is the Chapman's grazing incidence integral described in Smith III and Smith (1972) with the improved approximation of complementary error functions from Ren and ManKenzie (2007).

Structure of PROTEUS
PROTEUS consists of a Python GUI program and of Fortran 90 modules.The Python GUI program contains a list of chemical reactions for each planet, string parsing functions that parse the reactions and reaction rate coefficients selected in GUI, and GUI functions that control the behavior of GUI.Python language was adopted because of its flexibility in parsing strings and its capability in operating GUI.Fortran language was adopted to solve differential equations numerically.The structure of PROTEUS is illustrated in Fig. 1.Chemical reactions listed in the Python file are first read by GUI functions (arrow 1 in Fig. 1).Then, users select reactions on GUI (arrows 2 in Fig. 1), which will be analyzed by string parsing functions (arrow 3 in Fig. 1) to output a Fortran 90 module named "v__in.f90"and data files (in directories "PLJ_list" and "settings") including information of the selected chemical reactions (arrow 4 in Fig. 1).The Fortran 90 module "v__in.f90" is the only module that includes information of chemical reactions, and other Fortran 90 modules are independent of chemical reactions to be used in the simulation.Datafiles of initial density profiles and temperature profiles are read by "v__ in.f90" (arrow 5 in Fig. 1).The selected reactions, settings such as temperature profile and initial density profiles are applied to the Fortran 90 model when the main Fortran 90 routine "e__main.f90"call subroutines in "v__in.f90"(arrow 6 in Fig. 1).Users can then run the Fortran model by compiling all the Fortran 90 modules (indicated as "7" in Fig. 1).
The structure of the Fortran codes is illustrated in Fig. 2. It should be noted that each Fortran 90 file contains several modules with distinct functions, but only the Fortran 90 file is indicated for simplicity.The main routine "e__main.f90"consists of three parts, (1) initialization, (2) calculation, and (3) finalization.The description of each parts and each Fortran 90 files are as follows.
(1) All derived data types are defined in "v__tdec.f90".Information of the chemical reactions, boundary conditions, and calculation settings are defined in "v__in.f90".Physical constants and parameters of the planetary orbit are given in "c__prm.f90".Production rates calculated by other models (e.g., ionization rate calculated by a meteoroid model (Nakamura et al. 2022)) can be input in "p__Mars.f90"and "p__Jupiter.f90″.The solar EUV flux is calculated by the EUVAC model and the absorption and ionization cross sections are defined in "p__EUVAC.f90″.The solar flux of Woods et al. (2009) is defined and absorption and dissociation cross sections are calculated in"p__UV.f90".
(2) The radiative transfer is solved and the optical depth is calculated in "p__photochem_opticaldepth.f90".Ionization and dissociation rates, reaction rate coefficients and production and loss rates of each species are calculated in "p__photochem_rate.f90".The vertical diffusion flux is calculated in "p__photochem_transport.f90".The eddy and binary diffusion coefficients are defined in "p__eddy_diffusion.f90"and "p__molecular_diffusion. f90", respectively."p__photochem_scheme.f90" calculates the Jacobian matrix and advances the timestep using the implicit method.

Graphical user interface
The Python GUI program uses the tkinter package, a standard library of Python to use the Tcl/Tk GUI toolkit (https:// docs.python.org/3/ libra ry/ tkint er.html).
The Python GUI allows users to easily and intuitively select a planet, chemical reactions of interest, and run the simulation.An example of GUI is shown in Fig. 3, and the operation of GUI is as follows.Once the user runs the Python GUI program, one can select a planet as indicated ("1" in the upper panel of Fig. 3).After selecting a planet, one can create a new project directory, or select or rename a pre-existing project directory as indicated ("2" in the upper panel of Fig. 3).Then the chemical reaction list for the selected planet appears in the window (the lower panel of Fig. 3).Chemical reactions and their rate coefficients written in the Python GUI program are automatically converted into Unicode and displayed, making them easy to read.One can select or clear reactions by clicking on the checkbox at each reaction ("3" in the lower panel of Fig. 3).By inserting chemical species, reference or label into a search box, only related reactions appear in the window.One can set upper and lower boundary conditions, initial density profiles, vertical grid size, and other calculation settings such as dimension of the simulation, season, latitude, integration time, maximum timestep ("4" in the lower panel of Fig. 3).After all the settings are done, one can press "Output f90 module", then the following files are generated in the selected project directory: a Fortran 90 module named "v__in.f90",setting files stored in the directory "settings", and information about production and loss reactions of each chemical species, rate coefficient labels, and Jacobian matrix analyzed by the string parsing function stored in the directory "PLJ_list".The list of chemical species, the number of chemical species, indices, mass, and charge of each chemical species are automatically determined and written in "v__in.f90"at this time.Those text files are read by the Fortran 90 module "v__in.f90".One can also output those files, compile and run the Fortran codes by pressing "Output f90 module & Run model", which requires the installation of an open source software CMake into user's computer ("5" in the lower panel of Fig. 3).All the settings and selected reactions are saved, and users can use the same settings and selected reactions the next time they run GUI.At the end of the simulation, users can quickly plot the density Fig. 3 Overview of GUI and instruction of the operation profiles by pressing "Plot setting" button ("6" in the lower panel of Fig. 3).

Format of chemical reaction list
The main feature of PROTEUS is the simple format of chemical reaction list in the Python GUI program and on the GUI.Format of chemical reaction list in the Python GUI program and some examples are illustrated in Fig. 4.
Any reactions and their rate coefficients are described in the following string format in the Python GUI program.Reaction and rate coefficient are separated by a colon ":", and left-and right-hand sides of the reaction are separated by an arrow "-> ".Chemical species can be written simply as string.For instance, ionized species "N 2 + ", "CO 2 + ", and "H + (H 2 O) 4 " are simply described as "N2+", "CO2+", and "H+(H2O)4", respectively, and electron is described as "e-".Isotope species such as " 13 CO 2 " can be written as "^13CO2".Each chemical species and an addition operator " " or an arrow "->" should be separated by at least one space.PROTEUS also deals with three-body reactions with the expression "M" describing the total atmospheric number density.Temperature-dependent rate coefficient equation can be simply described as string in infix notation.Addition operator " + ", subtraction operator "−", multiplication operator "*", division operator "/", exponentiation operator "^" or "**", exponential function "exp()" square root function"sqrt()", neutral, ion and electron temperatures "Tn", "Ti", and "Te", respectively, altitude in km "h", decimal fraction values such as "1.16", integer values such as "300", and values in E notation such as "4.9e−11" can be used in the rate coefficient equation.If there is a temperature range T 1 -T 2 [K] in which the rate coefficient is valid, one can describe the temperature range by "for T = T 1 ~ T 2 [K]".
Reactions and their rate coefficients selected on GUI are parsed by the string parsing functions in the Python GUI program.Index for each chemical species is automatically determined by the string parsing function, and mass and charge of each chemical species are also automatically identified by the string parsing function.String of each species are automatically divided into constituent elements and mass is calculated by the sum of mass of all the elements, and charge is calculated by counting the number of " + " and "−" in the string of each species.The string parsing function analyzes which reaction produces or lose each chemical species.The rate coefficient expressions written in infix notation are first separated into tokens.Then the order of tokens in infix notation are converted into reverse Polish notation (i.e., postfix notation) and automatically labeled.The Fortran 90 modules calculates the reaction rate coefficient using the labeled tokens arranged in reverse Polish notation.This method allows PROTEUS to process a variety of expression of temperature-and altitude-dependent rate coefficients (as seen in Fig. 4) at high computational speed.All the information needed to calculate production rate, loss rate and Jacobian matrix are output as text files, which will be read by Fortran 90 module to apply the information about chemical reactions selected.The number of chemical species and reactions, mass and charge of chemical species, rate coefficient of each reaction and contribution of each reaction to production/loss of each species are automatically applied to Fortran 90 modules by reading those text files.Those features make PROTEUS a flexible  Nair et al. (1994), Fox and Sung (2001), Scott et al. (1999), Dotan et al. (1997), Chaffin et al. (2017), Troe (2005), Pavlov (2014), andVerronen et al. (2016) photochemical model that can be applied to many planetary atmospheres with different set of chemical reactions.

Software requirements
PROTEUS is compatible with Python 3 and does not support Python 2. The Python libraries, Matplotlib and Numpy, need to be installed on the user's computer.
Additionally, the open-source software CMake needs to be installed on the user's system.

Mars
For the application to the Martian atmosphere, parameters of Mars and its orbit are implemented into PRO-TEUS; The mean distance between Mars and the sun is r m = 1.524AU, the eccentricity is e = 0.0934, the solar longi- tude at perihelion is L s,P = 250°, tilt angle of the rotational axis is ε = 25.2°, the rotational period is T p = 88,775 s, the mass of Mars is 6.417 × 10 23 kg, and the mean radius of Mars is 3389.5 km (Patel et al. 2002;Williams 2021).
The cross sections implemented into PROTEUS for the application to Mars are as follows.Ionization cross sections of CO 2 , CO, O 2 , N 2 , and O are taken from Schunk and Nagy (2009).Absorption cross sections and quantum yields for calculating dissociation rates of atmospheric molecules are listed in Table 1.In order to validate PROTEUS, we compared with one-dimensional Martian photochemical model by Chaffin et al. (2017) (hereafter called as C17 model), using the same chemical reactions and their rate coefficients, boundary conditions, temperature and water vapor profiles, and binary and eddy diffusion coefficient profiles.The neutral density profiles simulated by PROTEUS and C17 model are shown in Fig. 5. PROTEUS and C17 model are in good agreement except for small differences for O 3 , OH, HO 2 and H 2 O 2 .Those differences could be due to the difference in the photo-absorption and dissociation cross sections used in the two models.

Jupiter
For the application to the Jovian parameters of Jupiter and its orbit are implemented into PROTEUS; The mean distance between Jupiter and the sun is r m =5.2 AU, the eccentricity e , the solar longitude at perihelion L s,P , and tilt angle of the rotational axis ε are set to zero for simplicity, the rotational period is assumed to be the System III period related to the period of radio burst T p = 35,729.71s, the mass of Jupiter is 1.898 × 10 27 kg, and the equatorial radius of Jupiter is 71,492 km (Williams 2021; Russell et al. 2001).
Ionization cross sections of hydrogen molecule and atom, helium atom, hydrocarbon molecules (CH 4 , C 2 H 2 , C 2 H 4 , and C 2 H 6 ) and metallic atoms (Fe, Mg, Si, and Na) implemented into PROTEUS for the application to the Jovian ionosphere are found in Appendix of Nakamura et al. (2022) and references therein.
PROTEUS has recently been applied to the Jovian ionosphere by Nakamura et al. (2022).Chemical reactions regarding hydrocarbon ion chemistry used in the simulation are described in Nakamura et al. (2022).Ion density profiles calculated by PROTEUS with 218 reactions are shown in Fig. 6.Simulated ion density profiles are in good agreement with Kim and Fox (1994), as discussed in Nakamura et al. (2022).Slight differences seen in the shape of profiles of hydrocarbon ions could result from the difference in the input density profiles of hydrocarbon molecules, which are not indicated in Kim and Fox (1994).

Conclusions
We have newly developed a flexible one-dimensional photochemical model named PROTEUS, which consists of a Python GUI program and of Fortran 90 modules.Chemical reactions can be easily implemented into Python code as a simple string format, and users can intuitively select a planet and chemical reactions to be considered in their calculation on GUI.Chemical reactions selected on GUI are automatically analyzed by a string parsing code written in Python, which will be applied to Fortran 90 modules to simulate with selected chemical reactions on a selected planet.This paper presents examples of PROTEUS application to the Martian atmosphere and the Jovian ionosphere, which are in good agreement with previous numerical models.PROTEUS can significantly save time for those who need to develop a new photochemical model; they just need to add chemical reactions in the Python code and just select them on GUI to run a new photochemical model.PROTEUS can be easily extended to other planets and satellites, e.g., Venus, Earth, Titan, and exoplanets in the future.

Appendix
See Tables 1 and 2. σ a : Absorption cross section, σ d : dissociation cross section, φ : quantum yield, a : data files are taken from The MPI-Mainz UV/VIS Spectral Atlas (Keller-Rudek et al. 2013),   b   : data files are taken from PHIDRATES (Huebner and Mukherjee 2015), c : quantum yields for each photolysis reaction of HNO 3 were estimated by quantum yield of each product (OH, O, and O( 1 D)) obtained by Johnston et al. (1974), Turnipseed et al. (1992), and Margitan and Watson (1982).

Fig. 1
Fig. 1 Schematic illustration of the structure of PROTEUS

Fig. 5
Fig. 5 Application of PROTEUS to the Martian atmosphere.Vertical profiles of neutral density simulated by PROTEUS (solid) and one-dimensional Martian photochemical model by Chaffin et al. (2017) (dashed).The same boundary conditions, chemical reactions and their rate coefficient, binary and eddy diffusion coefficients, and temperature profile, were used in both simulations for validation

Table 1
List of cross sections and quantum yields implemented into PROTEUS

Table 1
(continued)Photoabsorption, dissociation cross sections and quantum yield implemented into PROTEUS are listed in Table1.Photoionization cross sections implemented into PROTEUS are listed in Table2.