#  PICASO 3.0: A One-Dimensional Climate Model for Giant Planets and Brown
Dwarfs

Sagnick Mukherjee1 , Natasha E. Batalha2 , Jonathan J. Fortney1 , Mark S.
Marley3 1Department of Astronomy and Astrophysics, University of California,
Santa Cruz, CA 95064, USA
2 NASA Ames Research Center, MS 245-3, Moffett Field, CA 94035, USA
3 Lunar and Planetary Laboratory, The University of Arizona, Tucson, AZ 85721,
USA

###### Abstract

Upcoming James Webb Space Telescope (JWST) observations will allow us to study
exoplanet and brown dwarf atmospheres in great detail. The physical
interpretation of these upcoming high signal–to–noise observations requires
precise atmospheric models of exoplanets and brown dwarfs. While several one-
dimensional and three-dimensional atmospheric models have been developed in
the past three decades, these models have often relied on simplified
assumptions like chemical equilibrium and are also often not open-source,
which limits their usage and development by the wider community. We present a
python-based one-dimensional atmospheric radiative-convective equilibrium
model. This model has heritage from the Fortran-based code (Marley et al.,
1996) which has been widely used to model the atmospheres of Solar System
objects, brown dwarfs and exoplanets. In short, the basic capability of the
original model is to compute the atmospheric state of the object under
radiative–convective equilibrium given its effective or internal temperature,
gravity and host–star properties (if relevant). In the new model, which has
been included within the well-utilized code-base PICASO, we have added these
original features as well as the new capability of self-consistently treating
disequilibrium chemistry. This code is widely applicable to Hydrogen-dominated
atmospheres (e.g. brown dwarfs and giant planets).

###### Subject headings:

Brown Dwarfs, exoplanets

##  1\. Introduction

There are three broad categories of substellar atmosphere modeling frameworks:
one-dimensional (1D) physically self-consistent models, three-dimensional (3D)
General Circulation Models (GCMs) and atmospheric retrieval. Atmospheric
models enable the understanding of physical and chemical processes driving the
behaviour of substellar atmospheres. When combined with observational data,
they can be used to infer the physical state of atmospheres. Additionally, the
adiabat calculated by the atmospheric models provides the upper boundary
conditions of the interior evolutionary models of substellar objects like
giant planets and brown dwarfs (Hubbard, 1977; Burrows et al., 1997; Chabrier
et al., 2000; Saumon & Marley, 2008). Although substellar atmospheres
generally constitute a small fraction ($\sim$ 1%) of the total mass of their
bodies, their radiative properties control the overall cooling of objects
throughout their evolution. Therefore, atmospheric models are crucial for
understanding both physical (e.g. climate and chemistry) and global properties
of substellar objects (e.g. radius/luminosity evolution over time).

A variety of 1D radiative–convective–thermochemical equilibrium models of
exoplanetary and brown dwarf atmospheres have been developed (e.g. Marley et
al., 1996; Barman et al., 2001; Fortney et al., 2005, 2008; Marley et al.,
2021; Phillips et al., 2020; Piskorz et al., 2018; Goyal et al., 2020). These
models originate both from stellar atmosphere modeling codes (e.g. Barman et
al., 2001; Sudarsky et al., 2003; Seager & Sasselov, 1998) and planetary
atmospheric modeling codes (e.g. Marley et al., 1996; Baudino et al., 2015).
These models rely on iterating to a self-consistent radiative-convective-
thermochemical equilibrium (RCTE) solution. Self-consistency here means that
the model iterates on all components of the atmospheric structure
simultaneously such that they are physically consistent with each other (e.g.,
chemical abundances are consistent with the current temperature profile). The
models do not capture 3D dynamics, but instead are able to incorporate the
effects of radiative/convective energy transport, chemistry, and clouds. These
1D models are also computationally efficient, which enhances their exploratory
power significantly over 3D circulation models. For example, (e.g. Zhang et
al., 2021) used the SONORA BOBCAT grid of models, generated using the 1D RCTE
model described in Marley et al. (2021), to infer the physical properties of
55 T-dwarfs in an uniform analysis. The computational speed of these 1D self-
consistent models along with their ability to treat atmospheric physics and
chemistry self-consistently makes them a powerful tool for interpreting
atmospheres of planets and brown dwarfs.

Besides the model used in Marley et al. (2021), several other independent 1D
RCTE models have been widely used in the literature (e.g. Tremblin et al.,
2015; Phillips et al., 2020; Gandhi & Madhusudhan, 2017; Burrows et al., 2008;
Marley et al., 2021; Malik et al., 2017). For example, the GENESIS code from
Gandhi & Madhusudhan (2017) was used for the high significance detection of
various Carbon and Nitrogen bearing species in the atmosphere of HD 209458b by
Giacobbe et al. (2021). The ATMO model (Tremblin et al., 2015; Phillips et
al., 2020) was used by Goyal et al. (2018) to perform a uniform fitting
analysis of the transmission spectra of 10 hot Jupiters which led to important
conclusions about their atmospheric composition. Despite the many insights
gained from these models, their simplicity often makes them insufficient to
explain a variety of atmospheric spectra (e.g. Noll et al., 1997; Oppenheimer
et al., 1998; Saumon et al., 2003; Golimowski et al., 2004; Geballe et al.,
2009; Sorahana & Yamamura, 2012; Leggett et al., 2012; Miles et al., 2020).
For example, Miles et al. (2020) found the presence of disequilibrium
chemistry in a series of late T and early Y-dwarfs by measuring their
significantly enhanced photospheric CO abundances. In another example, Zhang
et al. (2021) concluded that presence of clouds can lead to much better
overall fit to the available infrared spectra of 55 brown dwarfs with a
uniform analysis using the SONORA cloudless grid of models. This motivates
adding physical complexity to 1D RCTE models in order to enable better
interpretation of upcoming data, for example, from JWST.

Models for exoplanetary and brown dwarf atmospheres have been developed in
three-dimensions (3D) (e.g. Showman et al., 2008, 2020; Tan & Showman, 2021;
Menou & Rauscher, 2009; Roman & Rauscher, 2019; Wolf et al., 2022; Lee et al.,
2021) as well. The 3D models work to capture the detailed global dynamics like
horizontal transport/winds in the atmospheres along with other atmospheric
components like radiation and clouds. However, these models are
computationally intensive which limits their power of exploration of the vast
parameter space in question. Moreover, in order to make the models
computationally feasible, they rely on approximations such as gray atmospheres
in the radiative transfer calculations (e.g. Tan & Showman, 2021, 2019; Menou
& Rauscher, 2009; Tan, 2022). As a result, 3D models are unable to resolve
some atmospheric properties predicted by 1D models – like detached convective
zones in brown dwarfs and giant planets.

Atmospheric retrievals are also widely used to interpret observational data.
These studies aim to recover the atmospheric state of any planet/brown dwarf
based on the observed spectral data in a Bayesian framework (e.g. Madhusudhan,
2018; Line et al., 2017; Burningham et al., 2017). The free parameters of
interest typically include the atmospheric temperature-pressure ($T(P)$)
profile, chemical abundances and cloud properties (Line et al., 2017;
Burningham et al., 2017, 2021; Mukherjee et al., 2021; Taylor et al., 2021).
The multi-dimensional parameter space of these free parameters is sampled by
repeatedly calculating model observables (e.g., spectra) with simple and
computationally fast forward models. Comparison of these modeled observables
with the observed data helps in retrieving likelihood distributions for all
the retrieved free parameters. However, these retrieved properties are often
not required to have physical constraints. For example, the retrieval analysis
of hot Jupiter, WASP-18 b, was predicted to have a thermal inversion and
relatively high metallicity (283$\times$Solar), corresponding to a CO
abundance of $>$10% (Sheppard et al., 2017). Later using a self-consistent
grid of models, Arcangeli et al. (2018) showed that this planet is more likely
to have a solar metallicity and hot dayside, which is more consistent with
what is expected from a hot Jupiter with little heat redistribution. This
means that even though retrieval studies are effective in inferring the
various atmospheric properties like abundances and temperature structures,
they are not equipped to provide the physical interpretations of why the
atmosphere is in the retrieved state, and they often result in unphysical
solutions. Therefore, physically motivated self-consistent models are required
to understand the physical and chemical processes that drive the structure of
planetary and substellar atmospheres.

In this work, we focus on 1D physically motivated
radiative–convective–thermochemical equilibrium models for substellar
atmospheres. In addition to include complexities like disequilibrium chemistry
within these models, another important modification we must make to these
codes, is to transform them to open sourced, “FAIR” codes (findable,
accessible, interpretable, and reproducible). Doing so will enable the
community to interpret the upcoming influx of data from missions like JWST
(Pontoppidan et al., 2022). Currently, there are a handful of 1D climate
models that are open source (e.g., HELIOS (Malik et al., 2017, 2019b), and
TLUSTY (Hubeny, 1988; Hubeny et al., 2003; Sudarsky et al., 2003)). HELIOS is
a Python-based GPU dependent model and it has been used for modeling
atmospheres of a variety of substellar objects (e.g. Fossati et al., 2021;
Malik et al., 2019a; Yan et al., 2022; Deline et al., 2022). TLUSTY is a
stellar atmosphere model written in FORTRAN77. This code has been modified to
be applied to substellar atmospheres as well and is called COOLTLUSTY (Hubeny
et al., 2003). COOLTLUSTY also has been widely used for brown dwarfs and
exoplanets (e.g. Burrows et al., 2006; Spiegel et al., 2010; Spiegel &
Burrows, 2012; Lacy et al., 2019; Lacy & Burrows, 2020). In this work, we aim
to add to develop a new open-source Python based one-dimensional radiative-
convective equilibrium model for H-dominated substellar atmospheres – PICASO
3.0 that is both: 1) open source and 2) capable of including disequilibrium
chemistry induced by vertical mixing self-consistently within the 1D
radiative–convective equilibrium framework. PICASO 3.0 has been released
publicly as an extension of the already open-source 1D and 3D radiative
transfer tool PICASO (Batalha et al., 2019). It is also accompanied by
detailed tutorials exploring all of its uses and limitations. Here, we focus
on the description of the numerical techniques used in the code along with
benchmarking to previous studies to demonstrate its functionality.

In §2, we discuss the methodology of PICASO 3.0, the new Python atmospheric
model. We present the benchmarking of our model with other models in §3
followed by recommendations on how to use this model in §4. We briefly discuss
ongoing and future improvements in §6 and conclude in §7.

##  2\. Model Setup of PICASO 3.0

The heritage of our code is the Fortran based EGP sub-stellar atmospheric
model. The EGP model has been used for substellar atmospheres including solar
system planets – Titan and Uranus (McKay et al., 1989; Marley & McKay, 1999),
the L–T–Y brown dwarf sequence (e.g. Marley et al., 1996; Morley et al., 2014,
2012; Marley et al., 2021; Karalidi et al., 2021; Zhang et al., 2021), cloudy
atmospheres (e.g. Cushing et al., 2008; Stephens et al., 2009) and a wide
variety of extrasolar planets (e.g. Fortney et al., 2005, 2007, 2008; Marley
et al., 2012; Fortney et al., 2020; Morley et al., 2015, 2017) for over two
decades. It has also been used as the primary radiative-transfer scheme in the
SPARC GCM (Showman et al., 2009). Briefly, this 1D model solves for the self-
consistent temperature, chemical and cloud structure of H-dominated
atmospheres under the assumption of radiative-convective equilibrium. The code
structure of PICASO 3.0 is composed of Python classes
111https://docs.python.org/3/tutorial/classes.html which generally include
multiple Python functions. We refer to these classes and functions as modules
here. In general, Python is significantly slower than Fortran. However, with
the use of ‘numba‘’s ‘just-in-time‘ framework (Lam et al., 2015) to create on-
the-fly compiled machine code, PICASO 3.0 has comparable run times to the
original Fortran. Figure 1 shows a simplified schematic of PICASO 3.0’s
workflow.

First, the user provides an initial set of physical properties for the object
to be modeled, shown within the red boxes at the top of Figure 1. The model
supports both non-irradiated (e.g., brown dwarf) and irradiated (e.g.,
planets) calculations. Therefore, we specify optional inputs with dashed
outlines. For example, for modeling a field brown dwarf atmosphere, the user
only must specify the $T_{\rm eff}$ ($T_{\rm int}$ in case of a planet),
gravity, atmospheric metallicity and C/O ratio of the brown dwarf. For
modeling an irradiated planet, the user must specify the internal temperature
of the planet $T_{\rm int}$, planet gravity, atmospheric metallicity, C/O
ratio, semi-major axis, and host star properties (e.g., stellar temperature,
metallicity, and gravity). $T_{\rm int}$ represents the temperature obtained
by converting the internal heat flow of the planet via the Stefan-Boltzmann
law. Along with these inputs, the user must also specify an initial pressure-
temperature profile ($T(P)$) guess, which is divided into plane-parallel
logarithmically spaced pressure layers ($\sim$ 60-90). Then, through an
iterative process, which we describe in §2.1, the model iterates through
computing chemistry, opacities, and net upwelling and downwelling fluxes until
a radiative-convective equilibrium threshold has been met. During this
iteration, the model solves for both radiative and convective parts of the
atmosphere and takes into account the possibility of multiple radiative and
convective zones. Ultimately, the model produces the final atmospheric state
of the object ($T(P)$ and associated chemistry). These can then be used to
compute transmission, emission, and/or reflected light spectra used to compare
with observations (see Figure 1). We discuss these physical and chemical
aspects one by one in the following §2.1.

Figure 1.— Schematic diagram of the radiative-convective model in PICASO 3.0.
Solid outlines indicate required modules, while dashed outlines indicate
optional models.

###  2.1. Physics and Chemistry of Substellar Atmospheres

####  2.1.1 Radiative-Convective Equilibrium

The key physical basis of this model is the assumption of radiative-convective
equilibrium in substellar atmospheres. The radiative equilibrium represents
the physical scenario where within each atmospheric layer the energy emitted
must be balanced by the energy absorbed. This means that each atmospheric
layer must allow the transfer of the same amount of radiative energy through
it. Hubeny (2017) provides a detailed derivation of the radiative–convective
equilibrium criteria. For completeness, we present some of the important steps
here. Radiative equilibrium scenario is represented by:

$\int_{0}^{\infty}(\chi_{\nu}J_{\nu}-\eta_{\nu})d{\nu}=0$ (1)

where $\chi_{\nu}$ and ${\eta_{\nu}}$ are absorption and emission
coefficients, respectively and $J_{\nu}$ is the first moment of the intensity
field. Assumption of local thermodynamic equilibrium allows Equation 1 to be
rewritten as:

$\int_{\nu}\kappa_{\nu}(J_{\nu}-B_{\nu})d{\nu}=0$ (2)

where $\kappa_{\nu}$ is the wavelength dependent total absorption opacity of
the layer and $B_{\nu}$ is the local Planck function of the layer. This is the
integral form of the radiative equilibrium condition for each atmospheric
layer. The integral form of the radiative–equilibrium condition (Equation 2)
is applicable and numerically stable throughout the atmosphere. Another method
of computing radiative-equilibrium condition, which is used by our model, is
referred to as the differential form. The differential form is numerically
stable at optically thicker parts of the atmosphere and might become
numerically unstable at parts of the atmosphere which are optically thin
($\tau<<$ 1). In practice, this instability only appears at parts of the
atmosphere with low temperatures and very small optical depths. In our first
release of the climate code, we stick to our original methodology of using the
differential form, despite the instabilities pointed out by Hubeny (2017). In
a later update we will explore and implement the integral form as well within
our model to make our solutions more accurate at the optically thinner parts
of the atmosphere.

The differential form can be derived by using Equation 2 along with the second
moment of the radiative transfer equation,

$\int_{\mu}\mu\dfrac{dI_{\nu}}{d\tau_{\nu}}d\mu=\int_{\mu}\mu(I_{\nu}-S_{\nu})d{\mu}$
(3)

which under local thermodynamic equilibrium leads to the form,

$\dfrac{dH_{\nu}}{dz}=\kappa_{\nu}(J_{\nu}-B_{\nu})$ (4)

where $H_{\nu}$ is the wavelength dependent second moment of the specific
intensity field. Comparing Equation 2 and 4 leads to the condition,

$\int_{\nu}\dfrac{dH_{\nu}}{dz}d{\nu}=0$ (5)

which means that the integral of the second moment of the intensity,
$H_{\nu}$, must be constant at all layers. Mathematically, this can be written
as:

$\int_{\nu}H_{\nu}d{\nu}-\dfrac{\sigma_{sb}T_{\rm eff}^{4}}{4\pi}=0$ (6)

where $\sigma_{sb}$ is the Stefan–Boltzmann constant and $T_{\rm eff}$ is the
effective temperature of the object. This calculation can be found in the
t_start module in the code. The model tries to achieve this radiative
equilibrium in all the atmospheric layers which are stable against convection.

Given Equation 6, which defines the radiative-equilibrium condition, the
general procedure of the model is to start with an initial guess of the $T(P)$
profile of the atmosphere. This can be seen in the required inputs box at the
top of Figure 1. Radiative fluxes through all the atmospheric layers are then
calculated given a chemical state of the atmosphere. These radiative fluxes
are used to compute the wavelength integrated flux carried by each layer
$H(z)$. If the atmospheric layers are not in radiative equilibrium, this
quantity $H(z)$ would not be equal to the target flux given by
$\dfrac{\sigma_{sb}T_{\rm eff}^{4}}{4\pi}$, and the $T(P)$ solution would be
perturbed until the convergence criteria set by Equation 6 is met. The
methodology of the convergence is discussed further in §2.2.

####  2.1.2 Radiative Transfer

In order to compute radiative-convective equilibrium with Equation 6,
radiative fluxes must be calculated for all atmospheric layers. The
fundamental equation describing atmospheric radiative transfer is,

$\displaystyle\mu\dfrac{dI_{\nu}(\mu,\tau_{\nu},\phi)}{d\tau_{\nu}}=$
$\displaystyle I_{\nu}(\mu,\tau_{\nu},\phi)-S_{\nu}(\mu,\tau_{\nu},\phi)$ (7)
$\displaystyle-\dfrac{\omega_{0}}{4\pi}\int_{0}^{2\pi}\int_{-1}^{1}P_{\nu}(\mu,\mu^{{}^{\prime}},\phi,\phi^{{}^{\prime}})$
$\displaystyle\cdot
I_{\nu}(\mu^{{}^{\prime}},\tau_{\nu},\phi^{{}^{\prime}})d\mu^{{}^{\prime}}d\phi^{{}^{\prime}}$

where $\mu=\cos(i)$ where $i$ is the incident angle of the ray, $\phi$ is the
azimuthal angle, $I_{\nu}$ is the specific intensity, $\tau_{\nu}$ is the
frequency dependent optical depth, $S_{\nu}$ is the frequency dependent source
function and $P_{\nu}$ is the scattering phase–function. The first term on the
right–hand–side of Equation 7 describes the attenuation of the intensity with
increasing optical depth. The source function $S_{\nu}$ is used to describe
the emission of the atmosphere itself (e.g. thermal emission) or any external
radiation source (e.g. incident stellar radiation). The third term in Equation
7 represents scattering of the intensity within the atmosphere where the
phase–function $P_{\nu}(\mu,\mu^{{}^{\prime}},\phi,\phi^{{}^{\prime}})$ is the
probability that intensity
$I_{\nu}(\mu^{{}^{\prime}},\tau_{\nu},\phi^{{}^{\prime}})$ will be scattered
from the direction ($\mu^{{}^{\prime}},\phi^{{}^{\prime}}$) to the direction
($\mu,\phi$).

In order to solve Equation 7, we follow the two-stream radiative transfer
methodology described in Toon et al. (1989) for the calculation of radiative
fluxes. The specific intensity $I_{\nu}(\mu,\tau_{\nu},\phi)$ can be
integrated over the azimuthal angle $\phi$ to calculate the azimuthally
integrated intensity $I_{\nu}(\mu,\tau_{\nu})$. The upward and downward fluxes
can also be defined by,

$\displaystyle F_{\nu}^{+}=$
$\displaystyle\int_{0}^{1}{\mu}I^{+}_{\nu}(\mu,\tau_{\nu})d\mu$ (8)
$\displaystyle F_{\nu}^{-}=$
$\displaystyle\int_{0}^{1}{\mu}I^{-}_{\nu}(\mu,\tau_{\nu})d\mu$

where $I^{+}_{\nu}(\mu,\tau_{\nu})$ is the azimuthally integrated intensity
for upward values of $\mu$ and $I^{-}_{\nu}(\mu,\tau_{\nu})$ is the
azimuthally integrated intensity for downward values of $\mu$. Equation 7 can
be integrated to produce two separate coupled equations in terms of the upward
and downward fluxes instead of the $\mu$ and $\phi$ dependent intensities.
These equations for the upward and downward fluxes are222there is a typo in
Eqn. 12 of Toon et al. (1989) which incorrectly swaps in $S^{+}_{\nu}$ for
$S^{-}_{\nu}$ in the negative partial flux,

$\displaystyle\dfrac{{\partial}F_{\nu}^{+}}{\partial\tau_{\nu}}=$
$\displaystyle\gamma_{1}F_{\nu}^{+}-\gamma_{2}F_{\nu}^{-}-S^{+}_{\nu}$ (9)
$\displaystyle\dfrac{{\partial}F_{\nu}^{-}}{\partial\tau_{\nu}}=$
$\displaystyle\gamma_{2}F_{\nu}^{+}-\gamma_{1}F_{\nu}^{-}+S^{-}_{\nu}$

where $\gamma_{1}$ and $\gamma_{2}$ are functions of the scattering properties
of the medium and $S^{+}_{\nu}$ and $S^{-}_{\nu}$ are modified versions of the
source function. The radiative transfer calculation within our model is
divided into two distinct parts: 1) the transfer of thermal emission through
the atmosphere and 2) the transfer of the reflected external (e.g. stellar)
radiation throughout the atmosphere. Equation 9 is applicable to both of these
components but the functions $\gamma_{1}$, $\gamma_{2}$, $S^{+}_{\nu}$, and
$S^{-}_{\nu}$ are defined differently for each. A brief description of the
radiative transfer of thermal radiation in our model is provided below
followed by a discussion of the radiative transfer of external stellar
radiation.

We use the two-stream source function technique described in Toon et al.
(1989) to calculate the thermal upward and downward fluxes in each layer. The
hemispheric mean approximation approach is used in this method where
$\gamma_{1}$ is given by 2-$\omega_{0}(1+g)$ and $\gamma_{2}$ is
$\omega_{0}(1+g)$. $\omega_{0}$ and $g$ are the single scattering albedo and
the scattering asymmetry parameter of the atmospheric layer, respectively.
This hemispheric mean approximation was used to obtain functional forms for
the source functions $S^{+}_{\nu}$ and $S^{-}_{\nu}$ in Toon et al. (1989)
(see Table 3 and Equations 53 & 54 in Toon et al. (1989)). With these source
functions, the upward azimuthally averaged intensities at the top
($I^{+}_{n}(0,\mu)$) and bottom ($I^{+}_{n}(\tau,\mu)$) of an atmospheric
layer with optical depth $\tau$ can be written as,

$\displaystyle I^{+}_{n}(0,\mu)=$ $\displaystyle
I^{+}_{n}(\tau,\mu)e^{-\tau/\mu}+\dfrac{G}{\lambda\mu-1}(e^{\lambda\tau}e^{-\tau/\mu}-1)$
(10)
$\displaystyle+\dfrac{H}{\lambda\mu+1}(1-e^{-\lambda\tau}e^{-\tau/\mu})+\alpha_{1}(1-e^{-\tau/\mu})$
$\displaystyle+\alpha_{2}(\mu-(\tau+\mu)e^{-\tau/\mu})$

This calculation can be found in the get_thermal_1d_gfluxi module. Similarly
the downward azimuthally averaged intensities at the top ($I^{-}_{n}(0,-\mu)$)
and bottom ($I^{-}_{n}(\tau,-\mu)$) of the same atmospheric layer can be
written as,

$\displaystyle I^{-}_{n}(\tau,-\mu)=$ $\displaystyle
I^{-}_{n}(0,-\mu)e^{-\tau/\mu}+\dfrac{K}{\lambda\mu-1}(e^{-\tau/\mu}-e^{-\lambda\tau})$
(11)
$\displaystyle+\dfrac{J}{\lambda\mu+1}(e^{\lambda\tau}-e^{-\tau/\mu})+\sigma_{1}(1-e^{-\tau/\mu})$
$\displaystyle+\sigma_{2}({\mu}e^{-\tau/\mu}+\tau-\mu)$

and can also be found in the get_thermal_1d_gfluxi module. The functions $G$,
$H$, $K$, $I$, $\alpha_{1}$, $\alpha_{2}$, $\sigma_{1}$, and $\sigma_{2}$ have
been computed for the hemispheric-mean approximation in Toon et al. (1989) and
are used in our calculation as well. Solving equations 10 and 11 also require
boundary conditions for the diffuse flux at the top and bottom of the
atmosphere. These boundary conditions are set using the thermal blackbody
intensities at the top and bottom of the atmosphere using,

$\displaystyle B_{\rm top}=$
$\displaystyle(1.-e^{-\tau^{{}^{\prime}}/\mu_{1}})B({T_{\rm top}})$ (12)
$\displaystyle B_{\rm bot}=$ $\displaystyle B({T_{\rm
bot}})+\mu_{1}\dfrac{B({T_{\rm bot}})-B({T_{\rm bot-1}})}{\tau}$

where $B(T)$ represents the blackbody function, $T_{\rm top}$ is the
temperature of the top most atmospheric layer, $T_{\rm bot}$ is the
temperature of the bottom most atmospheric layer, and $\tau$ is the optical
depth of the bottom most atmospheric layer. $\mu_{1}$ is assumed to be 0.5 due
to the hemispheric-mean approximation. Note that the bottom boundary condition
is valid only for gas giant atmospheres where the highest pressure grid point
corresponds to the end of the user-defined grid and does not correspond to a
“surface”. PICASO 3.0 does have the option to swap boundary conditions that
are pertinent to the hard surfaces needed for terrestrial atmospheres, but it
is not relevant for this work on gas giants and as such we do not discuss it
here. $\tau^{{}^{\prime}}$ is given by,

$\tau^{{}^{\prime}}=\tau_{top}\dfrac{{P_{\rm top}}}{({P_{\rm top+1}}-{P_{\rm
top}})}$ (13)

This expression along with the top boundary condition in Equation 12 captures
the downward thermal flux arising from the part of the atmosphere which has
pressure less than the minimum pressure in the used atmospheric pressure grid.
This formulation prevents arbitrary artificial cooling of the top-most
atmospheric layer.

Equations 10 and 11 can be used to calculate the incidence angle ($\mu$)
dependent upward and downward – $I_{n}^{+}(\mu)$ and $I_{n}^{-}(\mu)$
intensity field in each atmospheric layer. But in order to use these radiative
intensities ultimately for the convergence criteria in Equation 6, they need
to be integrated with Equation 8 to calculate the direction independent upward
and downward fluxes – $F^{+}_{\nu}$ and $F^{-}_{\nu}$. In order to compute
these disk-averaged, layer fluxes, the upward and downward thermal intensity
for each atmospheric layer is calculated using Equation 10 and 11 at five
incident angles. The cosine of these incident angles ($\mu$) are determined
using the Gaussian quadrature method with 5 points(Abramowitz & Stegun, 1970).
The choice of 5 Gauss points is the default setting in the code but other
choices are also available in the get_angles_1d module. Table 1 shows the five
default values of $\mu$ and the corresponding Gauss weights used for this
Gaussian quadrature integration. The intensities at different incident angles
are used to compute the integral in Equation 8 with the Gaussian-quadrature
integration technique using the weights in Table 1. The module
get_thermal_1d_gfluxi uses the formulation described above to calculate the
wavelength dependent upward and downward thermal fluxes at the edges of each
atmospheric layer.

$\theta$(deg) | $\mu=cos(\theta)$ | Gauss Weight
---|---|---
84.345 | 0.09853 | 0.015747
72.270 | 0.30453 | 0.073908
55.804 | 0.56202 | 0.146386
36.680 | 0.80198 | 0.167174
16.221 | 0.96019 | 0.096781
Table 1Gauss points and weights used for Gaussian quadrature integration of
thermal flux over different angles.

For the radiative transfer of the reflected stellar light, we use the
quadrature approximation. The calculation is performed for a single incidence
angle of the stellar radiation beam (60${}^{\circ},\mu_{0}$=0.5) unlike the 5
Gauss point method of the thermal counterpart. Under the quadrature
approximation, $\gamma_{1}$ is $0.5\sqrt{3}(2-\omega_{0}(1+g))$, $\gamma_{2}$
is $\omega_{0}\sqrt{3}(1-g)/2$, $\gamma_{3}$ is ($1-\sqrt{3}g\mu_{0}$)/2, and
$\gamma_{4}$ is 1-$\gamma_{3}$. These functions can be found in the
get_reflected_1d_gfluxv module. The source functions $S_{\nu}^{+}$ and
$S_{\nu}^{-}$ for this component are,

$\displaystyle S_{\nu}^{+}=$
$\displaystyle\gamma_{3}{\pi}F_{s}\omega_{0}e^{-(\tau_{c}+\tau)/\mu_{0}}$ (14)
$\displaystyle S_{\nu}^{-}=$
$\displaystyle\gamma_{4}{\pi}F_{s}\omega_{0}e^{-(\tau_{c}+\tau)/\mu_{0}}$

where $\tau_{c}$ is the cumulative optical depth of all the atmospheric layers
above the layer of calculation and $\tau$ is the optical depth of the layer
itself. $F_{s}$ here represents the stellar flux incident on the top of the
atmosphere. $F_{s}$ is interpolated from the PHOENIX grid of models (Allard et
al., 2012) available as a part of the PySynPhot (STScI Development Team, 2013)
package. Using these functions, Equation 9 can be solved for the upward and
downward fluxes for the $n$’th atmospheric layer as has been shown in Toon et
al. (1989),

$\displaystyle F_{\nu,n}^{+}=$ $\displaystyle
k_{1n}e^{\lambda_{n}\tau_{n}}+\Gamma_{n}{k_{2n}}e^{-\lambda_{n}\tau_{n}}+C_{n}^{+}$
(15) $\displaystyle F_{\nu,n}^{-}=$
$\displaystyle\Gamma_{n}k_{1n}e^{\lambda_{n}\tau_{n}}+{k_{2n}}e^{-\lambda_{n}\tau_{n}}+C_{n}^{-}$

where the quantities $\Gamma_{n}$, $k_{1n}$, $k_{2n}$, $\lambda_{n}$,
$C_{n}^{+}$, and $C_{n}^{-}$ are defined in Toon et al. (1989) (Equation 21,
22, 23, and 24). Equation 15 is solved in the get_reflected_1d_gfluxv module.
However, an additional term needs to be added to the downward fluxes solved
from Equation 15 in this formulation. This term is,

$\displaystyle
F_{\nu,n}^{-}=F_{\nu,n}^{-}+\mu_{0}{F_{s}}e^{-\tau_{c}/\mu_{0}}$ (16)

Like the thermal component, complete solutions of Equation 15 also require
boundary conditions on the diffuse flux at the top and bottom of the
atmosphere. We use the boundary conditions,

$\displaystyle F_{\rm top,st}=$ $\displaystyle 0$ (17) $\displaystyle F_{\rm
bot,st}=$ $\displaystyle{R_{s}}\mu_{0}F_{s}e^{-\tau_{\rm c,bot}/\mu_{0}}$

where $F_{\rm top,st}$ and $F_{\rm bot,st}$ are the diffuse flux at the top
and bottom of the atmosphere, $R_{s}$ denotes the reflectivity of the bottom
surface and $\tau_{c,bot}$ is the cumulative optical depth of the deepest
layer. The reflectivity of the bottom surface is typically assumed to be 1% in
the model, but users should check that the bottom surface is not optically
thin. This is further discussed in the §4. This formulation is used by the
get_reflected_1d_gfluxv module to calculate the wavelength dependent upward
and downward reflected stellar light fluxes at the edges of each atmospheric
layer. But the calculation of both the thermal and reflected stellar light
fluxes following the procedure described above requires another important
parameter – the optical depth $\tau$. The implemented $\tau$ calculation
procedure in this model is described below.

Generally, two approaches are used for calculating the optical depth $\tau$ –
1) the line–by–line approach (e.g. Gandhi & Madhusudhan (2017),Burrows et al.
(2008)) or 2) the correlated-k approach (e.g. Marley & McKay (1999), Malik et
al. (2017), Fortney et al. (2008)). This model uses a correlated-k opacity
approach for computationally efficient inclusion of gaseous opacities. In the
correlated-k approach, each wavelength bin is represented by 8 distinct
Gaussian quadrature points ($g_{i}$), each with its associated weight
($w_{i}$). Table 2 shows the Gauss points and the Gauss weights used for the 8
point correlated-k approach in the model. Each of these Gauss-points also have
an optical depth $\tau_{i}$ associated with them. These “Gauss” points should
not be confused with the “Gauss” points used to calculate the disk-integrated
fluxes. This means that the radiative transfer equations for each wavelength
bin must be computed 8 times, once for each value of optical depth $\tau_{i}$
corresponding to the $i$’th Gauss point. This calculation is done by calling
get_thermal_1d_gfluxi and get_reflected_1d_gfluxv modules 8 times by the
climate module. The 8 thermal and radiative upward and downward fluxes are
then added together following,

$F^{\pm}_{n}(\lambda)=\sum_{i=1}^{i=8}w_{i}F^{\pm}_{i,n}(\lambda,\tau_{i,n})$
(18)

where $F^{\pm}_{i,n}$ is the thermal/reflected light flux at the $n$’th
atmospheric layer calculated with the optical depth $\tau_{i,n}$ corresponding
to the $i$’th Gauss point. Our calculation of correlated-k opacities is
described in more details in §2.1.4. This summation with Equation 18 within
the climate module ultimately produces thermal/reflected wavelength dependent
upward and downward fluxes in each atmospheric layer. These fluxes can now be
used to calculate the convergence criteria in Equation 6. The net thermal flux
for the $n$’th layer can be simply obtained with,

$\displaystyle F^{n}_{\rm thermal,net}=$
$\displaystyle\int_{\lambda}(F^{+}_{n}(\lambda)-F^{-}_{n}(\lambda))d\lambda$
(19) $\displaystyle F^{n}_{\rm stellar,net}=$
$\displaystyle\int_{\lambda}(F^{+}_{n,s}(\lambda)-F^{-}_{n,s}(\lambda))d\lambda$

where $F^{+}_{n}(\lambda)$ represents the upward thermal flux from the $n$’th
layer within the wavelength bin between $\lambda$ and $\lambda$+d$\lambda$ and
$F^{-}_{n}(\lambda)$ is the downward thermal flux from the same layer within
the same wavelength bin. $F^{\pm}_{n,s}$ also represent the same quantity for
the reflected light component. The total net radiative flux in the $n$’th
atmospheric layer can be calculated using,

$F^{n}_{\rm net}=r_{\rm th}F^{n}_{\rm thermal,net}+r_{\rm st}F^{n}_{\rm
stellar,net}$ (20)

where rth is the contribution factor of the thermal radiation to the net flux
and is generally fixed at 1 and rst is the contribution factor of the stellar
radiation to the net flux. Non-zero values of rst is only relevant when the
external irradiation on the atmosphere is non-zero. In the scenario when a
user is computing a planet-wide average $T(P)$ profile, the stellar
irradiation is contributing to 50% (one hemisphere) of the planet and as a
result $r_{\rm st}=0.5$. If instead the goal is to compute a night-side
average atmospheric state, $r_{\rm st}$ is set to be 0. On the other extreme,
to compute the day-side atmospheric state of a tidally locked planet $r_{\rm
st}$ should be set at 1.

This full net flux, $F^{n}_{\rm net}$, is same as $4\pi\int_{\nu}H_{\nu}d\nu$
in Equation 6. Therefore, $F^{n}_{\rm net}$ can be used to check if the
convergence criteria (Equation 6) is satisfied at all the radiative layers of
the atmosphere. This check is done in the “Converged?” box of the flow chart
shown in Figure 1. Now that we have described the radiative transfer within
our model, we will discuss other blocks of the model flow chart shown in
Figure 1 starting with atmospheric chemistry.

Gauss Points ($g_{i}$) | Gauss Weights ($w_{i}$)
---|---
0.065960251992824 | 0.165231051440291
0.313509004297193 | 0.309768948559709
0.636490995702807 | 0.309768948559709
0.884039748007175 | 0.165231051440291
0.953471592210149 | 0.008696371128436358
0.966500473910379 | 0.01630362887156367
0.983499526089621 | 0.01630362887156367
0.996528407789851 | 0.008696371128436358
Table 2Gauss points and weights used for calculation of the correlated-k
opacities

####  2.1.3 Equilibrium Chemistry

As shown in Figure 1, in order to compute the radiative fluxes for checking
the convergence criteria, as outlined in the previous §2.1.2, we first need a
method to provide the chemical state of the atmosphere. This is required
because the chemistry dictates the optical depth ($\tau$) required for the
radiative transfer calculations. In the simplest case, we determine the
chemical state of the atmosphere using the $T(P)$ profile of the atmosphere,
atmospheric metallicity, and C/O ratio assuming chemical equilibrium. The
atmospheric metallicity (M/H) is defined as the ratio of the abundances of all
heavy elements to Hydrogen abundance in the atmosphere. Our climate model uses
a pre-calculated grid of molecular abundances on a pressure vs. temperature
vs. [M/H] vs. C/O grid. This chemistry grid is computed using the
thermochemical equilibrium models presented in Gordon & McBride (1994);
Lodders (1999, 2002); Visscher et al. (2006, 2010) and using protosolar
elemental abundances from Lodders (2010). The grid includes 73 temperature
points between 75 K and 4000 K and 20 pressure points logarithmically spaced
between 10-6 and 3000 bars. This corresponds to a pre-computed grid with 1460
grid points. The molecular abundances included in this grid are H2, H, H+, H-,
H${}_{2}{{}^{-}}$, H${}_{2}^{+}$, H${}_{3}^{+}$, He, H2O, CH4, CO, CO2, OCS,
HCN, C2H2, C2H4, C2H6, NH3, N2, PH3, H2S, SiO, TiO, VO, Fe, FeH, MgH, CrH, Na,
K, Rb, Cs, Li, LiOH, LiH, LiCl, and e-. The metallicities included in the grid
are: [M/H] (relative to Solar)= -1, -0.75, -0.5, -0.3, -0.25, 0.0, 0.5, 0.7,
1.0, 1.5, 1.7 and 2. The C/O ratios included in the grid are: C/O (relative to
Solar) = 0.25, 0.5, 1.0, 1.5, 2, and 2.5. In this way, Solar values are
[M/H]=0 and C/O=1. In order to change the C/O ratio for a given atmospheric
metallicity, both the elemental abundances ratios – C/H and O/H are
increased/decreased slightly while maintaining a constant (C+O)/H. The version
of PICASO 3.0, described in this work, does not allow to change other
elemental ratios ( e.g. S/H or N/H ). Future releases will include this
flexibility.

Figure 2 shows the pre-computed volume mixing ratios of four major atmospheric
gases – CH4, CO, H2O, and NH3 for a solar mixture of chemical elements.
Converged brown dwarf $T(P)$ profiles with log(g)=5 and multiple $T_{\rm eff}$
values between 300 K and 2300 K are also overplotted in the pressure-
temperature space to depict parts of the parameter space relevant for objects
with different $T_{\rm eff}$. There are a few important features that strongly
influence the climates of sub-stellar atmospheres that are worth noting.
First, the sharp drop in H2O vapor abundance at $T\lessapprox$200 K which can
be seen in Figure 2 lower left panel. This drop in H2O vapor abundance is
caused by the condensation of H2O into cloud particles. A similar condensation
effect can also be seen in the case of NH3 at $T\leq$100 K. Condensation
induced changes in the gaseous abundances are included in the pre-computed
grid of gaseous abundances. Second, the most abundant carbon bearing gas in
the atmosphere also changes from CO to CH4 between temperatures of $\sim$ 1000
and 1200 K. Figure 2 top right panel shows that CO is the major carbon bearing
gas above $\sim$ 1200 K but at temperatures cooler than 1200 K, CO abundance
decreases rapidly and CH4 abundance rises. This change can be understood by
the net chemical reaction responsible for inter-conversion between CO and CH4
under chemical equilibrium,

${\rm CH_{4}+H_{2}O\leftrightarrow CO+3H_{2}}\\\ $ (21)

This reaction is favored in the forward direction for higher temperatures than
$\sim$ 1200 K and as result CO is more dominant for such temperatures. The
reverse reaction dominates for temperatures lower than $\sim$ 1200 K and CH4
becomes the dominant C-bearing gas there. Many other similar interesting
trends has been explored in the literature for several gases (e.g. Lodders,
2002; Visscher et al., 2006, 2010; Moses et al., 2005, 2011, 2013). These
trends can be directly explored from the available chemical grid within the
model but we move on to another challenge of using a pre-computed chemistry
grid – limited number of grid points.

Due to the finite number of grid points (1460 points), chemical abundances
cannot be calculated exactly at any pressure–temperature point of interest.
Therefore, each time the chemistry routine premix_atmosphere is called, the
climate model uses an interpolation scheme to compute abundances at each
pressure-temperature layer point. This interpolation is performed by the
module chem_interp. The 2D interpolation scheme used relies on finding the
four surrounding grid points namely $T_{\rm low}(P_{\rm low}$), $T_{\rm
low}(P_{\rm high}$), $T_{\rm high}(P_{\rm low}$), and $T_{\rm high}(P_{\rm
high}$). Then the abundances in these surrounding points represented by
$\xi^{i}_{low,low}$, $\xi^{i}_{low,high}$, $\xi^{i}_{high,low}$, and
$\xi^{i}_{high,high}$ are used to interpolate the abundance of each species in
the $T(P)$ point. This interpolation is done using,

$\displaystyle\ln\xi^{i}(T,P)=(1-t_{\rm int})(1-p_{\rm
int})\ln\xi{{}^{i}_{low,low}}+$ (22) $\displaystyle t_{\rm int}(1-p_{\rm
int})\ln\xi{{}^{i}_{high,low}}+t_{\rm int}p_{\rm
int}\ln\xi{{}^{i}_{high,high}}+$ $\displaystyle(1-t_{\rm int})p_{\rm
int}\ln\xi{{}^{i}_{low,high}}$

where $t_{\rm int}$ and $p_{\rm int}$ are given by,

$\displaystyle t_{\rm
int}=\dfrac{1/T-{1}/{T_{low}}}{{1}/{T_{high}}-{1}/{T_{low}}}$ (23)
$\displaystyle p_{\rm int}=\dfrac{lnP-lnP_{low}}{lnP_{high}-lnP_{low}}$ (24)

In the case where the $T(P)$ moves outside the edges of this chemistry grid, a
linear interpolation scheme is instead adopted within the range 50-5200 K.
This allows some flexibility in the iteration scheme. However, beyond these
values, the model is not valid. Note, choices in interpolation routines cause
discrepancies in resulting abundances. For example, only using two nearest-
neighbors causes instabilities in chemical profiles. The interpolation chosen
here was tested in multiple chemical regimes to ensure stability. After
computing the atmospheric chemistry, the next step in Figure 1 is the
computation of the atmospheric opacities, which we describe next.

Figure 2.— Top left panel shows the abundance map of CH4 in the pressure-
temperature chemistry grid used in the model with equilibrium chemistry
described in the paper. Dashed red lines marked on the abundance map depict
different converged $T(P)$ profiles of brown dwarfs from Marley et al. (2021)
with $T_{\rm eff}$ between 300 K and 2300 K and log(g)=5. This grid is used
for interpolation of chemical abundances during iteration of the atmospheric
state in the code. The same for CO, H2O, and NH3 is shown in the top right,
bottom left and bottom right panel. The abundances shown here are calculated
for solar metallicity and C/O ratio.

####  2.1.4 Pre-mixed Opacities

The chemical structure of the atmosphere dictates the optical depth
($\tau(\lambda)$) of each atmospheric layer which is necessary for solving the
radiative transfer equations described in §2.1.2. This calculation of
atmospheric optical depths occurs in the “Opacity” block of the model flow
chart in Figure 1. As previously stated, we use the correlated-k approach to
handle the molecular opacities in this model. An alternative approach is the
line–by–line method which is more precise but comparatively speaking, more
computationally expensive than the correlated-k approach at low resolution,
and over large wavelength ranges. For completeness, we describe the general
methodology of our approach here. However, the data were originally computed
by Lupu et al. (2021), available on Zenodo, and detailed in Marley et al.
(2021). We include a table of the main reference data used to compute the pre-
mixed opacities in Table 3.

C2H2 | Rothman et al. (2013)
---|---
C2H4 | Rothman et al. (2013)
C2H6 | Rothman et al. (2013)
CH4 | Yurchenko et al. (2013); Yurchenko & Tennyson (2014)
CO | Rothman et al. (2010); Gordon et al. (2017); Li et al. (2015)
CO2 | Huang et al. (2014)
CrH | Burrows et al. (2002)
Fe | Ryabchikova et al. (2015); O’Brian et al. (1991); Fuhr et al. (1988); Bard et al. (1991); Bard & Kock (1994)
FeH | Dulick et al. (2003); Hargreaves et al. (2010)
H2 | Gordon et al. (2017)
H${}_{3}^{+}$ | Mizus et al. (2017)
H2–H2 |  Saumon et al. (2012) with added overtone from Lenzuni et al. (1991) Table 8
H2–He | Saumon et al. (2012)
H2–N2 | Saumon et al. (2012)
H2–CH4 | Saumon et al. (2012)
H${}_{2}^{-}$ | Bell (1980)
H- bf | John (1988)
H- ff | Bell & Berrington (1987)
H2O | Polyansky et al. (2018)
H2S | Azzam et al. (2016)
HCN | Harris et al. (2006); Barber et al. (2014); Gordon et al. (2022)
LiCl | Bittner & Bernath (2018)
LiF | Bittner & Bernath (2018)
LiH | Coppola et al. (2011)
MgH |  Yadin et al. (2012); GharibNezhad et al. (2013) assembled in Gharib-Nezhad et al. (2021)
N2 | Rothman et al. (2013)
NH3 | Yurchenko et al. (2011); Wilzewski et al. (2016)
OCS | Gordon et al. (2017)
PH3 | Sousa-Silva et al. (2014)
SiO | Barton et al. (2013)
TiO |  McKemmish et al. (2019) assembled in Gharib-Nezhad et al. (2021)
VO |  McKemmish et al. (2016) assembled in Gharib-Nezhad et al. (2021)
Li,Na,K | Ryabchikova et al. (2015); Allard et al. (2007a, b, 2016, 2019)
Rb,Cs |
Table 3Data used to compute correlated-K opacities. Correlated-K opacities are
available at Lupu et al. (2021) and detailed in Marley et al. (2021).

In the correlated-k approach, the relevant wavelength range for the radiative
transfer is first divided to 196 wavelength bins. These carefully chosen
wavelength bins are shown in Figure 3 and chosen to reflect the approximate
spectral energy distributions of planetary atmospheres. The molecular opacity
$\kappa$ in any one of these wavelength bins has numerous individual molecular
lines, as shown in Figure 4 left panel. The cumulative distribution function
(CDF) of the opacity can be defined as $G(\kappa_{0})=N(\kappa\leq\kappa_{0})$
where $N$ denotes the total number of instances when the condition
$\kappa\leq\kappa_{0}$ is satisfied within the wavelength bin. The CDF
$G(\kappa)$ of the opacity within each wavelength bin is computed. This
function, $G(\kappa)$, is then inverted to obtain the “k-distribution”. An
example k-distribution for the the opacity window shown in the left panel in
Figure 4, is shown in Figure 4 right panel. After computing a k-distribution
for each of the 196 wavelength bins, the 8 Gauss points ($g_{i}$) shown in
Table 2 are used to represent each distribution. Generally, the slope of
k-distributions are shallow and slowly changing at lower Gauss point values
($G$) between 0 and 0.9. Then, the slope steepens rapidly between $G$= 0.9–1,
as can be seen in Figure 4. To capture the complex shape of the k-distribution
with just 8 Gauss points, a double-Gauss method is used for the integration.
The first set of four Gauss points are shown in Figure 4 right panel with
black points. These four Gauss points sparsely sample most of the
k-distribution between $G=$0–0.95 whereas the last four Gauss points, shown
with red points in Figure 4, sample the small, but rapidly changing, part of
the k-distribution between $G=$0.95–1. The first and second set Gauss points
and weights for this double-Gauss point method can be generated from sample
points of the generally used Gauss-Legendre quadrature (Abramowitz & Stegun,
1970) using,

$\displaystyle g_{i,1}=$ $\displaystyle f(g_{i}^{{}^{\prime}}+1)/2$ (25)
$\displaystyle w_{i,1}=$ $\displaystyle fw_{i}^{{}^{\prime}}/2$ $\displaystyle
g_{i,2}=$ $\displaystyle f+(1-f)(g_{i}^{{}^{\prime}}+1)/2$ $\displaystyle
w_{i,2}=$ $\displaystyle(1-f)w_{i}^{{}^{\prime}}/2$

where $g_{i,1}$ and $w_{i,1}$ represent the first set of our Gauss points
(shown with black points in Figure 4) and $g_{i,2}$ and $w_{i,2}$ represents
the second set (shown with red points in Figure 4), $g_{i}^{{}^{\prime}}$ and
$w_{i}^{{}^{\prime}}$ are the sampling points for Gauss-Legendre quadrature of
some order (4 in our case) defined within the interval [-1,1], and $f$ is the
adjustable parameter which sets the division in the values of $G$ which will
be sampled by the first and second set of Gauss points. For example, if
$f=0.8$ then $g_{i,1}$ will sample values between 0 and 0.8 while the rest
will be sampled by $g_{i,2}$. For our purpose, this $f$ is set at 0.95 which
is a reasonable choice based on the typical shape of k-distributions as shown
in Figure 4.

This double-Gauss method for calculating correlated k-coefficients has a large
impact on the time and computational efficiency of the code as it reduces the
number of Gauss points required for radiative transfer calculations but still
maintains sufficient accuracy required for the model. It has been shown that
the radiative transfer with double-Gauss method with 8 Gauss points is as
accurate as a normal set of 20 Gauss points (Michael Line, priv. comm. as
detailed in Marley et al., 2021). The choice of the wavelength-bins for the
correlated-k opacities are also crucial for the model and are described next.

The 196 wavelength bins used span 0.2 $\mu$m – 227 $\mu$m to capture the
general spectral energy distribution of planetary atmospheres. However, these
wavelength bins are not all of equal wavelength width. These wavelength bins
are shown in Figure 3 where three blackbody curves corresponding to
temperatures of 300 K, 1000 K and 2000 K are also shown for comparison. The
bins are narrow between 0.5-10 $\mu$m to capture the large number of molecular
rovibrational bands in this range. At the tail ends of the blackbody
distributions, larger bins help boost computational speed and maintain
sufficient precision in radiative transfer calculations required for the
application cases of this model.

As our model uses a pre-computed chemistry grid during iterations on the
$T(P)$ profile, the k-coefficients are also pre-computed on the same
pressure–temperature–metallicity–C/O ratio grid. For intermediately $T-P$
values, the molecular opacity is interpolated using the same formalism as has
been described for the gas abundances in §2.1.3. Figure 5 shows heat maps of
pre-computed Planck function and abundance weighted molecular cross-sections
in the same pressure-temperature grid as Figure 2. For each $P-T$ point in the
grid, the abundance weighted molecular cross-sections are integrated over all
wavelengths using the Planck function (corresponding to the temperature $T$)
as the integrating kernel (see Equation 2 in Freedman et al. (2014)). Each
panel corresponds to cross-sections at different Gauss points which we use in
our models. Planck mean cross sections at only the first four Gauss points are
shown here in the four panels. Dashed black lines marked on the cross-section
maps depict different converged $T(P)$ profiles of brown dwarfs from (Marley
et al., 2021) with $T_{\rm eff}$ between 300 K and 2300 K and log(g)=5. As the
order of the Gauss points increases in Figure 5, the cross-sections become
higher as higher order Gauss points trace higher opacity parts of the
k-distribution. Achieving atmospheric convergence in regions of the $P-T$
space where the cross-section change rapidly with small changes in temperature
or pressure can be difficult. Such an “opacity cliff” can be seen at $T$
values between$\sim$ 900-1700 K in Figure 5. This cliff appears in the Planck
mean cross-sections due to the overlap of the peak of the Planck function at
900-1700 K with large H2O opacity bands between 1.5-3 $\mu$m. These “opacity
cliffs” were also seen in Freedman et al. (2008, 2014) with Rosseland-mean
opacities. We discuss the effect of these cliffs on model convergence in more
detail in §4.3.

In addition to molecular opacities, collision induced absorption (CIA) of
H2-H2, H2-H, H2-He, H2-N2, H2-CH4, and continuum opacities such as H-bf, and
H-ff are also accounted for. These are included separately from our
correlated-k table. The CIA opacities are pre-calculated between temperatures
of 75 K and 7000 K with 1000 wavelength bins and then interpolated to the
correct $P-T$ combination with a spline. Note that even though the CIA
temperatures are computed up to 7000 K, the temperature valid range of our
model is still limited by the 1460 opacity grid.

Figure 3.— Left panel shows the correlated-k opacity wavelength bins used in
the radiative transfer calculations in our atmospheric model. The blackbody
emission curves of temperatures of 1000 K, 2000 K and 300 K are also shown
here. Right panel shows a zoomed version of the left panel with the
wavelengths between 0.8-13 $\mu$m. The stellar spectrum from a solar
metallicity PHOENIX model for a star with $T_{\rm eff}$= 5000 K is also shown
in both the panels. Figure 4.— Left panel shows the high resolution opacity
$\kappa$ between 1 and 1.3 $\mu$m of a solar mixture of atmospheric gases at a
temperature of 550 K and a pressure of 80 mbar. Right panel shows the
k-distribution computed for the opacity shown in the left panel with the blue
solid line. The black and red points are the Gauss points (from Table 2)
overplotted on the k-distribution for this opacity window. These Gauss points
are applicable for the double Gauss point method. The first set of 4 Gauss
points are shown in black while the second set is shown with red points.
Figure 5.— Top left panel shows the heat map of planck weighted molecular
opacity in the pressure-temperature k-coefficient grid used in the model with
equilibrium chemistry described in the paper. Dashed black lines marked on the
cross-section map depict different converged $T(P)$ profiles of brown dwarfs
from Marley et al. (2021) with $T_{\rm eff}$ between 300 K and 2300 K and
log(g)=5. This grid is used for interpolation of opacities during iteration of
the atmospheric state in the code. The same for the 3rd, 5th, and 7th Gauss
points are also shown in the top right, bottom left and bottom right panel.

We note that the pre-computed opacity grid used with this model is known to be
incomplete with regards to atomic and ionic opacities, which are particularly
important at the high temperature and low pressure parts of the atmosphere
(e.g. Hoeijmakers et al., 2018, 2020). This limits, to some extent, the
current code’s ability to treat the “ultra-hot Jupiters”
($T_{\mathrm{eq}}\gtrapprox$2200 K). In a future update we will add these
opacities to the correlated-K tables.

####  2.1.5 Disequilibrium Chemistry

A significant addition to this model is the capability to treat vertical
mixing induced disequilibrium chemistry self-consistently within the
radiative-convective equilibrium framework. This is an optional part of the
model, as shown in Figure 1. In 1D models, vertical mixing is often
parametrized as a diffusive process which is described by the eddy diffusion
coefficient, $K_{zz}$ (Allen et al., 1981). The mixing timescales of all
relevant gases in each atmospheric layer is given by,

$t_{\rm mix}=\dfrac{H^{2}}{K_{\rm zz}}$ (26)

where $H$ is the local scale height of that atmospheric layer. All
thermochemical reactions within the atmosphere like CO$\leftrightarrow$CH4 and
NH3$\leftrightarrow$N2 also proceed with a characteristic timescale – $t_{\rm
chem}$. Zahnle & Marley (2014) parametrized the $t_{\rm chem}$ of several such
gases using 1D chemical kinetics models. The parametrized $t_{\rm chem}$ from
Zahnle & Marley (2014) and Visscher et al. (2006) are given by,

$\displaystyle t_{\rm CO,CH_{4},H_{2}O}$ $\displaystyle=\dfrac{3\times
10^{-6}}{P_{\rm bar}}exp\left(\dfrac{42000K}{T}\right)$ $\displaystyle t_{\rm
NH_{3},N_{2}}$ $\displaystyle=\dfrac{10^{-7}}{P_{\rm
bar}}exp\left(\dfrac{52000K}{T}\right)$ $\displaystyle t_{\rm CO_{2}}$
$\displaystyle=\dfrac{10^{-10}}{\sqrt{P_{\rm
bar}}}exp\left(\dfrac{38000K}{T}\right)$ $\displaystyle t_{\rm HCN}$
$\displaystyle=\dfrac{1.5\times 10^{-4}}{P_{\rm
bar}m^{0.7}}exp\left(\dfrac{36000K}{T}\right)$ $\displaystyle t_{\rm PH_{3}}$
$\displaystyle=\dfrac{1.9\times
10^{12}}{\rm[OH]}exp\left(\dfrac{6013.6K}{T}\right)$

where $P_{\rm bar}$ is the atmospheric pressure in bars, $T$ is the
temperature and $m$ is the atmospheric metallicity relative to solar
metallicity. $t_{\rm chem}$ is generally short at high pressure-high
temperature regions of the atmosphere. As the atmosphere gets colder at lower
pressures the $t_{\rm chem}$ rises exponentially and becomes large. The volume
mixing ratios of CO, CH4, NH3, CO2, H2O and HCN are expected to follow
equilibrium chemistry at pressures where $t_{\rm mix}$ $>>$ $t_{\rm chem}$
which happens in the deeper atmosphere. But as the atmosphere becomes colder
with lowering pressure, $t_{\rm chem}$ may exceed $t_{\rm mix}$. The pressure
at which this occurs is called the “quench pressure”. At pressures less than
the “quench pressure”, gases are expected to depart from chemical equilibrium
and their mixing ratios become constant.

In this model, we include these disequilibrium chemistry effects for the 8 gas
species shown in Eqn. 2.1.5: H2O, CH4, CO, CO2, NH3, N2, HCN, and PH3. The
quench pressure is determined using the parametrized $t_{\rm chem}$ in Eqn.
2.1.5. The abundance then follows the equilibrium chemistry at pressures more
than the quench level, and is held constant at pressures less than the quench
level.

The uncertain parameter in the disequilibrium modeling framework is the eddy
diffusion parameter, $K_{zz}$ (Phillips et al., 2020; Fortney et al., 2020;
Karalidi et al., 2021). Therefore, this model is flexible in assumptions
regarding $K_{zz}$. Currently, there are two user-options for defining
$K_{zz}$ within a model run: 1) a fully user-defined $K_{zz}$ value which is
either constant or variable throughout the height of the atmosphere, but which
does not change during the iterative climate solution (as was done in Karalidi
et al., 2021; Phillips et al., 2020), and 2) a model-predicted $K_{zz}$, which
is calculated from the $T(P)$ profiles in the convective zones, using mixing
length theory (Gierasch & Conrath, 1985), and in the radiative zones, using
parametrizations (e.g. Moses et al., 2021). In the second case, along with the
$T(P)$ profile, the $K_{zz}$ will also change simultaneously with the
iterations in the model. The $K_{zz}$ in the convective zone can be calculated
with mixing length theory and is given by (Gierasch & Conrath, 1985),

$K_{zz}=\dfrac{H}{3}{\left(\dfrac{L}{H}\right)}^{{4}/{3}}{\left(\dfrac{RF}{\mu{\rho_{a}}c_{p}}\right)}^{{1}/{3}}$
(27)

where $H$ is the local scale height of the atmosphere, $L$ is the turbulent
mixing length, $R$ is the universal gas constant, $\mu$ is the mean molecular
weight of the atmosphere, $\rho_{a}$ is the atmospheric density, $c_{p}$ is
the atmospheric specific heat at constant pressure and $F$ is the convective
heat flux. The convective heat flux can be calculated by the difference
between the net thermal radiative flux and $\sigma{T_{\rm eff}}^{4}$ within
the convective zones of the atmosphere. Therefore, $\sigma{T_{\rm eff}}^{4}$
is the maximum allowed value of $F$ in this framework if the energy transport
within the convective atmosphere is assumed to be completely convective.

While various parametrizations of $K_{zz}$ in the radiative zones of
substellar atmospheres have been discussed in the literature (Wang et al.,
2015; Zhang & Showman, 2018; Parmentier et al., 2013; Tan, 2022), we have
included the parametrization from Moses et al. (2021) as a starting point.
Further parametrizations of $K_{zz}$ can easily be swapped in, in the future.
The Moses et al. (2021) radiative zone $K_{zz}$ is given by,

$K_{zz}=\dfrac{5\times
10^{8}}{\sqrt{P_{\mathrm{bar}}}}\left(\dfrac{H}{{620}{\rm
km}}\right)\left(\dfrac{T_{\rm eff}}{1450{\rm K}}\right)^{4}$ (28)

where $P_{\mathrm{bar}}$ is the pressure of the radiative level in bars and
$H$ is the atmospheric scale height. Both Equation 27 and 28 can be found in
the get_kzz function. As the $T(P)$ profile iterates towards the converged
solution, the quench levels of various gases change as well. This means that
the abundances of quenched species will depart from chemical equilibrium, and
the pre-computed k-coefficient tables described in §2.1.4 are no longer valid.
The methodology to remix the k-coefficients with updated abundances and
recompute resulting optical depth calculation during the convergence process,
is referred to “on–the–fly” mixing, and is described in the following §2.1.6.

####  2.1.6 Mixing Opacities “on–the–fly”

With disequilibrium chemistry, the atmospheric chemistry depends on the quench
pressures, which again depend on the $T(P)$ profile of the atmosphere.
Therefore, because the chemistry of the atmosphere cannot be pre-determined,
the atmospheric opacities also need to be calculated “on–the–fly”.

We mix the correlated-k opacities of individual gases using the methodology of
Amundsen et al. (2017) called the resort-rebin technique. Currently, in this
model we focus on the quenching of CO, CH4, H2O, NH3, CO2, N2, HCN and PH3.
However, the major opacity sources among these gases are mainly CO, CH4, H2O,
and NH3. Meaning, the contribution of N2, HCN and PH3 on the total gas opacity
is negligible for the relatively small departures from the chemical
equillibrium with $\log$(M/H)$\leq$2.5 explored in this analysis. Therefore,
we mix the correlated-k opacities of CO, CH4, H2O, and NH3 with the
correlated-k opacities of all the other sets of background gases which follow
equilibrium chemistry as the volume mixing ratio of these gases evolve due to
quenching with the $T(P)$ profile of the atmosphere. Should the motivation to
include more gases in the “on–the–fly” methodology present themselves in
future observations, they can be included in a future code release.

One of the drawbacks of “on–the–fly mixing” outlined in Amundsen et al. (2017)
is the dependence of the accuracy of the technique on the spectral resolution
of the correlated-k opacities. This effect has been well explored in Karalidi
et al. (2021). Specifically, it was found that the 196 wavelength bins,
traditionally used in the chemical equilibrium version of EGP and this new
python version, is not sufficiently high to use with the resort-rebin
technique. Karalidi et al. (2021) found that 661 wavelength bins in the
correlated-k opacities are required to counteract the inaccuracies in the
resort-rebin technique. Therefore, when disequilibrium calculations are
requested by the user, PICASO automatically switches to 661 wavelength bins.

With the opacities and mixing routines described, the final module, shown in
Figure 1 is the computation of the convective zones, which we outline in the
following §2.1.7.

####  2.1.7 Convective Zones

As the $T(P)$ profile iterates towards radiative equilibrium, parts of the
atmosphere will become unstable against convection, and energy transport will
be expected to occur via convection instead of radiation. These parts of the
atmosphere are forced to follow the local adiabat. Figure 6 shows a heat map
of the adiabatic lapse-rate, ${\rm dlnT}/{\rm dlnP}$, used in PICASO. This
grid is pre-computed assuming a solar H–He mixture with He mass fraction of
Y=0.28. This includes H${}_{2}\leftrightarrow$ 2H dissociation and a detailed
accounting of the molecular vibrational and rotational levels333Raw grid can
be found in GitHub file specific_heat_p_adiabat_grad.json. The calculations
for this grid was done by Didier Saumon as described in Marley et al. (2021).
The pre-computed grid has 53 temperature points between 10 K and 3981 K and 26
pressure points between 10-2 and 103 bars. The lapse-rate of parts of the
profile which become unstable against convection are interpolated from this
grid. For every atmospheric layer, the lapse-rate, ${\rm dlnT}/{\rm dlnP}$, is
first calculated from the $T(P)$ profile. The grid of adiabatic lapse-rates
shown in Figure 6 is then used to interpolate the local adiabatic lapse rate
for each atmospheric layer. This interpolation is done using 2D interpolation
similar to the technique described in §2.1.3.

If the ratio between the lapse-rate obtained from the $T(P)$ profile and the
interpolated, lapse rate, $\nabla$ is greater than $\sim$1 (numerically set to
0.98) in any of the layers, then these layers are considered convective and
are forced to follow the interpolated local adiabatic lapse rate according to
the following equation,

$T^{i+1}=\exp\left({T^{i}+{\nabla}ln\left(\dfrac{P^{i+1}}{P^{i}}\right)}\right)$
(29)

where $T^{i}$ and $P^{i}$ are the temperature and pressure of the i’th
convective layer and $T^{i+1}$ and $P^{i+1}$ are the temperature and pressure
of the $i$+1’th convective layer starting from below the radiative-convective
boundary layer. We should note that this convective adjustment is done for one
layer at a time followed by a iteration of the entire $T(P)$ as this
adjustment for each layer, in principle, would lead to temperature changes in
all the other layers. This approach used in this model is different from
forcing all the layers which are unstable against convection to be convective
at once. The t_start module implements this convective adjustment of the
$T(P)$ profile. This methodology is used to develop and grow convective zones
in the atmosphere with the find_strat module during the iterations of the
model. Note that this module is designed to always form or grow convective
zones and not shrink or remove them. We discuss the important implications of
this for the user in §4.

Figure 6.— A heat map of the adiabatic lapse rate as a function of pressure
and temperature is shown here. This grid is used to determine the lapse rate
in atmospheric layers unstable against convection.

###  2.2. Iteration Scheme And Model Convergence

Each required and optional physical component of the model shown in Figure 1
has been described separately in §2.1. But as our model is iterative in
nature, a description of the iterative scheme of the model along with its
convergence criteria is a necessary component. We provide a brief outline of
the model iterative scheme for both the chemical equilibrium and
disequilibrium runs separately below.

####  2.2.1 Iterative Scheme for Chemical Equilibrium

This model aims to achieve radiative-convective equilibrium using the Newton-
Rhapson iterative scheme. As a first step, the chemistry and the correlated-k
opacities of the atmosphere for the first given guess $T(P)$ profile are
calculated as detailed in §2.1.3 and 2.1.4. Then, the methodology described in
§2.1.2 is used to compute the net radiative fluxes in all atmospheric layers.
The temperature of each layer is then perturbed by a small arbitrary $dT$,
historically set at 0.01% of the current layer temperature, while keeping the
other layer temperatures fixed. The radiative fluxes in all layers are then
recomputed after this perturbation. These two sets of radiative fluxes are
used to compute the Jacobian Aij which quantifies the change in radiative flux
in layer $i$ due to perturbed temperature in layer $j$. The Jacobian is given
by (e.g. Hubeny, 2017),

$A_{ij}=\dfrac{F^{j}_{i,p}-F_{i}}{dT_{j}}$ (30)

where $F_{i}$ is the net radiative flux in layer $i$ before the perturbation
and $F^{j}_{i,p}$ is the perturbed net radiative flux in layer $i$ due to a
change in temperature $dT_{j}$ in layer $j$. The ultimate temperature
correction, $\delta{T}$, needed to ensure radiative equilibrium (in the
radiative zones) is then solved using the equation,

$A\delta{T}=\sigma{T_{\rm eff}}^{4}-F(T(P))$ (31)

where $A$ is the Jacobian and $F(T(P))$ are the net radiative fluxes in each
layer with the current atmospheric state. This process continues iteratively
until a tolerance of maximum allowed radiative flux difference from
$\sigma$$T_{\rm eff}$4 is reached in all the radiative layers of the
atmosphere. This convergence criteria can be expressed as,

$\left|\dfrac{F^{n}_{\rm net}-{\sigma}T_{\rm eff}^{4}}{{\sigma}T_{\rm
eff}^{4}}\right|\leq\epsilon$ (32)

where $\epsilon$ denotes the tolerance parameter, set at $5\times 10^{-3}$ in
our model. This equation is a numerical version of Equation 6. The model is
also considered to be converged if the maximum fractional temperature
correction among all the radiative layers $|\delta{T/T}|$, calculated from
Equation 31, during any iteration is smaller than $\epsilon$. It is important
to note that the model satisfies the convergence criteria in Equation 32 in
the mid-point of each atmospheric layer while the temperature at the edges of
each atmospheric layer (levels) is iterated. For example, if the optical
depths at the bottom and top (edges) of the $i$’th atmospheric layer are
$\tau_{i}$ and $\tau_{i-1}$, then the model aims to satisfy Equation 32 in
this layer for an optical depth of $(\tau_{i}+\tau_{i-1})/2$. But, the
iteration of the temperature profile is done to the temperature at the bottom
and top (or edges) of this atmospheric layer. This is important for the
stability of the iterative process.

As the iteration progresses, the stability of each atmospheric layer against
convective mixing is also checked using the technique described in §2.1.7. If
layers unstable against convection are found, the temperature of these
convectively unstable layers are forced to a H2-He gas mixture adiabat from
the grid shown in Figure 6. Once the convergence criteria of Equation 32 is
met, the model run stops and it produces the outputs outlined by the green
boxes in Figure 1. Outputs such as the converged $T(P)$ profile, chemical
abundances and outgoing radiative fluxes in the 196 wavelength grid are always
returned to the user after model converges. These outputs can further be
optionally used for calculating various observables such as transmission,
emission or reflected light spectra with PICASO.

For a deeper understanding of the iteration scheme of this model, three
atmospheric iteration steps for a brown dwarf with $T_{\rm eff}$ of 1000 K and
log(g)=5 are shown in Figure 7 and iterations for an irradiated planet at a
distance of 0.1 AU from a sun like star with $T_{\rm int}$ of 300 K and
log(g)=3.4 are shown in Figure 8. In Figure 7 and 8, each row corresponds to
the atmospheric state in a certain iteration of the model. The first column
shows the $T(P)$ profile of the atmosphere in each iterative step of this
model with the red solid line whereas SONORA $T(P)$ profile for this case is
shown as a black dashed line. The second column shows the lapse-rate of the
$T(P)$ profile as a function of pressure in each iteration with the solid blue
line and the adiabatic lapse-rate is shown with a black solid line. The third
column shows the volume mixing ratio profiles of four gases – H2O, CH4, NH3,
and CO as a function of pressure and the fourth column shows the emergent
thermal spectrum from the top of the atmosphere. In Figure 7, the first row
shows a simple initial isothermal guess profile of 500 K and the last row
shows the final converged solution. For our initial guess, the bottom four
layers of the atmosphere have been assumed to be convective. The initial
thermal emission spectrum resembles a blackbody because the atmosphere is
isothermal. As the model iterates, the isothermal $T(P)$ profile is perturbed
to reach an atmospheric state such that Equation 32 is satisfied for all
radiative layers of the atmosphere. In the iteration shown in the middle row,
the bottom of the atmosphere becomes unstable against convection and therefore
in the final solution the lapse-rate in that part of the atmosphere is forced
to follow the local adiabatic lapse-rate. As the model iterates towards the
converged solution, the $T(P)$ profile becomes more complex and leads to a
redistribution of flux in each wavelength bin. This leads to molecular
absorption lines becoming more evident in the thermal spectrum with each
iteration. The final converged solution matches well with the SONORA BOBCAT
grid of models (shown as black dashed line). For the irradiated planet
convergence shown in Figure 8, the initial guess is an isothermal profile with
$T=700$ K and the iterations continue until Equation 32 is satisfied.

Figure 7.— Three steps in the iteration towards a converged solution for a
brown dwarf with $T_{\rm eff}$ of 1000 K and log(g) of 5 is shown here in each
row. The first column shows the iteration of the $T(P)$ profile in red with
the SONORA BOBCAT model shown in black dashed. The second columns show the
lapse-rate from the $T(P)$ profiles in blue and the local adiabatic lapse rate
in black. The third column shows the iterations in the volume mixing ratios of
various gases. The fourth column shows the thermal spectrum emergent from the
atmosphere of the object in each iteration. An animated gif of this figure can
be found here Figure 8.— Three steps in the iteration towards a converged
solution for an exoplanet around a sun like star with $T_{\rm int}$ of 300 K
and log(g) of 3.4 is shown here in each row. The left-most column shows the
iteration of the $T(P)$ profile in red with the SONORA BOBCAT model shown in
black. The second left-most columns show the lapse-rate from the $T(P)$
profiles in blue and the local adiabatic lapse rate in black. The third column
shows the iterations in the volume mixing ratios of various gases. The right-
most column shows the thermal spectrum emergent from the atmosphere of the
object in each iteration. An animated gif of this figure can be found here

A second helpful method of visualizing a converged run is to look at the
contribution of radiative and convective fluxes as compared to the target flux
in brown dwarfs given by $\sigma$$T_{\rm eff}$4 or the internal heat flux for
planets given by $\sigma$$T_{\rm int}$4. Figure 9 shows these layer–by–layer
radiative and convective fluxes for converged models of a representative brown
dwarf (left panel), warm Jupiter (middle panel), and hot Jupiter (right
panel). The net radiative flux in each layer is shown with the orange shaded
region. The shaded blue region depicts the convective flux at each layer. The
red dashed line shows the target flux for the brown dwarf – $\sigma$$T_{\rm
eff}$4. This converged brown dwarf model has two convective and radiative
zones as can be seen from its $T(P)$ profile shown in black. As a result, the
convective flux (blue) peaks at the location of these convective zones. The
net radiative flux decays rapidly at the deeper convective zone and convection
carries the majority of the energy in these deep convective layers. The proof
that this is a converged model and Equation 32 is satisfied lies in the fact
that the sum of the net radiative and convective fluxes is equal to the target
flux at all atmospheric layers. For the planet case where there is an
additional energy flux from the host star, we show the additional energy with
black hatched shading (2nd and 3rd column of Figure 9). This incident flux is
downward, compared to the thermal radiative and convective fluxes, which is
upward. Therefore, the upper atmosphere needs larger quantities of upward
thermal radiative fluxes to balance this downward stellar flux. Ultimately,
the goal is to maintain a summed flux of $\sigma$$T_{\rm int}$4 (red dashed
line) at all the atmospheric layers of the planet. Therefore, the upper
atmosphere of irradiated planets, with a given $T_{\rm int}$ , is pushed
toward hotter temperatures when compared to a brown dwarf with a comparable
$T_{\rm eff}$.

Figure 9.— Radiative and convective flux contributions compared to total net
fluxes are shown for three $T(P)$ profiles in order to visualize the concept
of “convergence”. From left to right a representative brown dwarf, warm
Jupiter, and hot Jupiter profile is shown. In each figure, the shaded blue
region depicts the convective flux at each layer. The orange shaded region
depicts the net radiative flux in each layer. The red dashed line shows the
target flux for the brown dwarf in the left panel. For the brown dwarf (non-
irradiated) this target flux is $\sigma$$T_{\rm eff}$4. The red dashed line in
the middle and the right panels shows the internal heat flux of the planets
which is $\sigma$$T_{\rm int}$4. Please note that the scale of the upper and
lower X-axis are different in all the three panels.

####  2.2.2 Iteration Scheme for Chemical Disequilibrium

Figure 10.— Three steps in the iteration towards a converged solution with
disequilibrium chemistry for a brown dwarf with $T_{\rm eff}$ of 700 K and
log(g) of 5.25 is shown here in each row. The first column shows the iteration
of the $T(P)$ profile with the red line and the SONORA BOBCAT chemical
equilibrium is model shown in black dashed. The second column shows the lapse-
rate from the $T(P)$ profiles in blue and the local adiabatic lapse rate in
black. The third column shows the chemical timescales of various gases as a
function of pressure with solid lines of different colors. The pressure
dependant mixing timescale is shown in the third columns with the black dashed
line. The fourth column shows the iterations of the volume mixing ratios of
various gases. An animated gif of this figure can be found here

The primary difference in the iterative scheme of the disequilibrium model is
that it uses a two-step convergence method. With the user specified input
parameters, the code first uses the basic chemical equilibrium model (§2.2.1)
to converge to a chemical equilibrium atmospheric solution. This atmospheric
solution is then used as the initial guess to the disequilibrium chemistry
model in the second step of the iterative process. Of particular importance,
is the need to recalculate the convective zones when disequilibrium chemistry
is turned on. This is important because recent work has shown that there are
major differences in the location, extent and number of convective zones
between the equilibrium chemistry and disequilibrium chemistry atmospheric
solutions (Mukherjee et al. (2022 submitted)). The general procedure is to
partially reset the convective zones before moving to disequilibrium solver.
First, any upper convective zones are removed. Next, the upper boundary of the
deepest convective zone is reset to slightly higher pressures (usually $\sim$5
levels in a 91-level atmosphere from 10-4 to 200 bars). This enables the
convective zones in the presence of vertical mixing to be re-calculated
because we cannot assume there will be minor perturbations from the
equilibrium chemistry solutions.

In order to highlight the convergence of a brown dwarf disequilibrium
chemistry model, a few atmospheric iterations are shown in Figure 10 for a
brown dwarf with $T_{\rm eff}$ of 700 K and log(g)=5.25. The figure is similar
to Figure 7 except here the third column now shows the mixing (black dashed
line) and chemical timescales (solid colored lines) as a function of pressure
and the fourth column shows the volume mixing ratio profiles of various gases.
The first row shows the atmospheric state reached after the first step of the
convergence process, where a converged chemical equilibrium solution is
achieved. The last row shows the final converged solution with disequilibrium
chemistry. The chemical equilibrium solution from the SONORA grid is shown
with the black dashed line in the first column. As the model iterates,
chemical equilibrium $T(P)$ profile is perturbed to reach an atmospheric state
where the $T(P)$ profile is visibly colder than the chemical equilibrium
solution. This is a result of different atmospheric optical depths as a result
of quenching of several gases in the deeper atmosphere. The chemical
equilibrium solution has a single deep convective zone for this brown dwarf.
However, in the iteration shown in the middle row and last row, a second
convective zone also develops when disequilibrium chemistry is treated self-
consistently. This shows that chemical disequilibrium also impacts the
location and number of convective zones in the atmosphere. As a second
convective zone develops in the middle and the third row, a decrease in mixing
timescale can also be seen in that pressure range of the convective zone
because convective zones are assumed to be much more efficient in mixing and
thus have higher $K_{zz}$ and lower $t_{\rm mix}$ than radiative zones.
Comparing the fourth column in the first and the last row shows that mixing
causes several orders of magnitude change in the abundance of CO in the upper
atmosphere when compared to chemical equilibrium solutions. With these
examples, we conclude our detailed discussion on the methodology of the model
and move on to the model benchmarking analysis.

##  3\. Benchmarking Analysis

Figure 11.— Benchmarking tests of the Python based-equilibrium chemistry
developed in this work with the published chemical equilibrium models from
Marley et al. (2021) and Phillips et al. (2020) are shown here. The
benchmarking has been shown in the form of the converged $T(P)$ profiles
predicted by each code for four different $T_{\rm eff}$ values at log(g)=5 in
the left panels while the pressure-derivatives of the $T(P)$ profiles also
have been shown in the right side panels for each case to show that the
convective and radiative zones predicted by both the codes match excellently.
Figure 12.— $T(P)$ profile computed using PICASO 3.0 for WASP-25 b is shown
with the red solid line compared with the ATMO models presented in Goyal et
al. (2020) with dashed black lines. The left panel shows the comparison for
sub-solar metallicity of 0.1$\times$solar whereas the right panel shows
comparisons between solar metallicity models. The same system parameters were
used to compute both the models and they agree with each other with less than
3% differences. Figure 13.— Benchmarking tests of the Python based-
Disequilibrium code used in this work by assuming $K_{zz}$= 0 with the
published chemical equilibrium models from Marley et al. (2021) are shown
here. The benchmarking has been shown in the form of the converged $T(P)$
profiles predicted by each code for four different $T_{\rm eff}$ values at
log(g)=5 in the left panels while the pressure-derivatives of the $T(P)$
profiles also have been shown in the right side panels for each case to show
that the convective and radiative zones predicted by both the codes match
excellently.  Figure 14.— Benchmarking tests of the Python based-
Disequilibrium code used in this work by assuming non-zero $K_{zz}$ values
with the published chemical disequilibrium models from Phillips et al. (2020)
are shown here. The benchmarking has been shown in the form of the converged
$T(P)$ profiles predicted by each code for four different $T_{\rm eff}$ values
and three different log(g) values in the three panels. The $T(P)$ profiles
match well between the two independent models with different sets of opacities
and totally different code methodologies.

We have benchmarked this model in four ways in order to ensure all the
described methodology is functioning as expected: 1) PICASO 3.0 chemical
equilibrium non-irradiated functionality vs. SONORA BOBCAT grid from Marley et
al. (2021), 2) PICASO 3.0 chemical equilibrium non-irradiated functionality
vs. ATMO 2020 code from Phillips et al. (2020), 3) PICASO 3.0 chemical
equilibrium irradiated functionality vs. ATMO code from (Goyal et al., 2020),
and 4) PICASO 3.0 chemical disequilibrium non-irradiated functionality vs.
ATMO 2020 code from Phillips et al. (2020).

###  3.1. Benchmarking chemical equilibrium non-irradiated atmospheres

Figure 11 shows the comparison between the three chemical equilibrium non-
irradiated models. The models described in this paper are shown in red, the
BOBCAT models are shown with a black dashed line, and the blue dot-dashed line
show the models from Phillips et al. (2020). The left columns show a
comparison between the $T(P)$ profiles of the three models. The right column
shows a comparison between the lapse-rate of the profiles. Included in the
benchmarking is four different log(g) = 5 brown dwarf models with $T_{\rm
eff}$ values of 600 K, 700 K, 800 K, and 1000 K, which covers significant
chemical transitions of carbon species. PICASO 3.0 and the models from (Marley
et al., 2021) match well with $T(P)$ profiles differences smaller than 0.1%.
The lapse-rates from this model also match excellently with those from Marley
et al. (2021). Given PIASO’s heritage in the EGP model used to compute the
grid in Marley et al. (2021), this is validation that the update to Python
from Fortran did not introduce numerical issues. The disagreement between
PICASO 3.0 and the chemical equilibrium models from the ATMO 2020 grid are
typically smaller than 1% in the deeper atmosphere below 10 mbars. At
pressures less than 10 mbars, the differences between the two models are
between 1% to 5%. These differences are considered minor given the independent
code setup, opacity calculations, and chemistry routines. For example, models
presented in Phillips et al. (2020) use the “on–the–fly” opacity mixing method
described in §2.1.6 even when modeling atmospheres with chemical equilibrium
whereas PICASO 3.0 uses interpolations of pre-mixed opacities for iterations
of chemical equilibrium models. Differences such as these along with
differences in opacity line lists used by the two models could account for the
1% to 5% differences in $T(P)$ profiles in the upper atmosphere.

###  3.2. Benchmarking chemical equilibrium irradiated atmospheres

We also benchmark the equilibrium chemistry version of the code by comparing
models of irradiated exoplanets calculated using PICASO 3.0 with the published
ATMO models presented in Goyal et al. (2020). We use the RCTE models for
WASP-25 b (Enoch et al., 2011; Southworth et al., 2014) presented in (Goyal et
al., 2020) to perform this benchmarking. WASP-25 b is a hot Saturn with an
estimated equilibrium Temperature of $\sim$ 1210 K (assuming 0 albedo) (Enoch
et al., 2011; Southworth et al., 2014). We use the same system parameters for
WASP-25 b as are used in the models presented in Goyal et al. (2020) and
detailed in the Appendix section of their article. Figure 12 shows the
comparison between the two models at sub-solar and solar metallicities with a
super-solar C/O ratio. The models generally agree for both the metallicities.
For the sub-solar metallicity models, the agreement is better than 1% across
all the pressures. For the solar metallicity models, the maximum disagreement
between the models is $\sim$ 3%. These minor differences were also found
between the PICASO 3.0 and ATMO 2020 models for brown dwarfs. As described
above, these minor differences can be attributed to different opacities and
numerical methodologies used.

###  3.3. Benchmarking non-irradiated atmospheres with diseequilibrium
chemistry

In order to benchmark the version of the model with disequilibrium chemistry
we have done two types of comparisons. First, we set $K_{zz}$ = 0 in our
disequilibrium chemistry model and compare to the results of our equilibrium
chemistry model from the SONORA grid. Second, we benchmark our disequilibrium
chemistry model with results from the ATMO 2020 grid.

In the first test, if $K_{zz}$ is assumed to be zero, $t_{\rm mix}$ becomes
infinitely large according to Equation 26. Therefore, in principle, none of
the gases will be quenched and their volume mixing ratios will follow
equilibrium chemistry throughout the atmosphere. Therefore, a disequilibrium
chemistry model run with $K_{zz}$ = 0 must produce the same result as the
equilibrium chemistry models. Of particular importance, this test checks if
our “on–the–fly” mixing routines, and double-iterative routine are performing
as expected. Figure 13 shows the benchmarking between our disequilibrium
chemistry model with $K_{zz}$ = 0 and the equilibrium chemistry models from
the SONORA BOBCAT grid. The left column shows the comparison of the $T(P)$
profiles between the two models and the right column shows the comparison of
the lapse-rates between the two models. Each row from the top to bottom shows
log(g) = 5.0 brown dwarf models with $T_{\rm eff}$ of 600, 700, 800, and 1000
K. The two models match with differences $\leq$ 1% level. The small $\leq$ 1%
deviations are due to the slight inaccuracies in the “on–the–fly” mixing
technique (Amundsen et al., 2017), which is sensitive to number of wavenumber
bins, versus the pre-weighted technique (see §2.1.6).

In the second test, we directly benchmark our disequilibrium chemistry model
with disequilibrium chemistry models from Phillips et al. (2020).
Disequilibrium chemistry models in Phillips et al. (2020) assume a constant-
with-pressure $K_{zz}$ for the atmospheres. Therefore, we use a constant
$K_{zz}$ in our benchmarking test and show the results in Figure 14. For a
brown dwarf with $T_{\rm eff}$ of 700 K and log(g) of 5.5, the top left panel
in Figure 14 shows: 1) the differences that arise in the $T(P)$ profiles when
disequilibrium is turned on for both the models from Phillips et al. (2020)
and PICASO 3.0, and 2) the differences that arise between Phillips et al.
(2020) and PICASO 3.0. Both the equilibrium chemistry and disequilibrium
chemistry models show excellent matches with among themselves which is an
excellent benchmarking demonstration of our methods and codes. Additionally,
this also shows the impact of disequilibrium chemistry on the $T(P)$ profiles
of this representative brown dwarfs, which results in $\sim$50 K colder $T(P)$
profiles compared to equilibrium chemistry models.

The rest of the panels in Figure 14 show comparisons between various
disequilibrium chemistry $T(P)$ profiles with different $K_{zz}$ and gravity
values produced with our model with the grid of models from Phillips et al.
(2020). Our disequilibrium models generally agree well within 5% levels in the
deeper atmosphere with models from Phillips et al. (2020). However, the
disagreements are higher in the much lower pressure upper atmosphere which was
also the case for the equilibrium chemistry model comparisons. The small
disagreements between the models are mainly due to two reasons – 1) Phillips
et al. (2020) use chemical kinetics models for calculation of quenched
abundances whereas we use the “quench-time” approximation (§2.1.5) to do so
and 2) differences between the accuracy of the “on–the–fly” mixing method used
between the two models due to differences in number of wavelength bins used.
These disagreements show the uncertainty in state–of–the–art models which
perhaps can be explained in future with comparisons with high signal-to-noise
data.

##  4\. Modeling Recommendations

The code is generally well-behaved across the parameter space of interest for
brown dwarfs and exoplanets. However, iterative schemes are sometimes
notoriously tricky to converge. Therefore, here we outline recommendations
regarding the use of code such that users can get meaningful and accurate
results. These recommendations are implemented in the publicly available code
tutorials available via the PICASO documentation page
444https://natashabatalha.github.io/picaso/tutorials.html released with
publication acceptance.

###  4.1. Choosing Model Pressure Grid

We recommend using typically 50-90 atmospheric pressure layers, corresponding
to 51-91 pressure levels. In general, higher number of pressure layers
increases the computational time required for convergence substantially and
lower number of layers makes the atmospheric grid too coarse for an accurate
calculation. The maximum and minimum values of the pressure grid are also
ultimately chosen by the user. While the choice of the minimum pressure can be
made somewhat arbitrarily, we recommend that the model not be run at pressures
lower than 10-6. This is the lowest pressure for which chemistry and opacities
are computed (shown in Figure 2 and 5). Another uncertainty from the double-
Gauss correlated-k approach can arise from running the model at very low
pressures ($\leq$ 10-4 bars). As the molecular lines become very narrow at
such low pressures, most of the opacities resides at very high values of
$G(\kappa)$ (e.g. between 0.995-1.0) in the k-distribution. Therefore, these
opacities are not taken into account even with the double-Gauss method. This
can make the radiative transfer at such low pressures inaccurate.

The maximum pressure of the atmospheric model needs to chosen carefully in the
case of giant planet atmospheres and brown dwarfs. If the maximum pressure of
the atmospheric model is too low, then the model can become transparent to the
deepest atmospheric layer, especially in wavelengths with little gaseous
opacities (e.g. optical wavelengths $\lessapprox 1\mu$m). This would result in
an inaccurate $T(P)$ profile. A good way to check if the profile was run with
low maximum pressure is to plot the wavelength dependent brightness
temperature of the converged model along with the converged temperature of the
deepest atmospheric layer. If the brightness temperature at any wavelength
exceeds the brightness temperature associated with the bottom temperature,
then the pressure grid is not well-suited for the calculation. However, if the
brightness temperature is smaller than this bottom temperature at all
wavelength, then the pressure grid is well-suited for the calculation. In the
first scenario, the pressure grid needs to be extended to higher pressures.
The problem with this approach to choose a pressure-grid is that it can only
be done after running the model once. A more practical alternative to this
approach is to find a comparable model from the SONORA BOBCAT 555Link to
published SONORA BOBCAT models to the temperature, gravity parameter space
wanted and use the maximum pressure corresponding to that model. Also, the
highest pressure for which chemistry and opacities are computed (shown in
Figure 2 and 5) is 3000 bars. Therefore, runs with pressure grids extending to
higher pressures than this will be inaccurate as they would require linear
extrapolation of both gas abundances and opacities.

When computing disequilibrium chemistry runs, the user must take another
aspect into consideration while choosing the maximum pressure of the
atmospheric pressure grid: the quench levels of various gases. If the maximum
pressure of the atmospheric grid is too low such that gases are expected to
quench at deeper parts of the atmosphere than the pressure grid, then the
model run will be incorrect. Therefore, users must check the quench levels of
various gases after performing a disequilibrium chemistry calculation with the
model. If the quench levels of any of the gases are at the last pressure
level, then the user needs to increase the maximum pressure of the run to get
a correct converged solution.

###  4.2. Choosing the Initial Guess T(P) Profile

Climate models that leverage an iterative scheme require a first guess of the
$T(P)$ profile. This guess $T(P)$ profile is then iterated to reach the
converged solution. Even though we used isothermal $T(P)$ profiles as our
first guess in both Figure 7 and 8, it is usually preferable to start with an
initial guess that is close to the expected solution. Using simple profiles
like isothermal profiles can lead to significant increase in run times and in
the worst cases also lead to solutions which do not converge. A far better
alternative to isothermal profiles, are publicly available model grids like
SONORA BOBCAT or the Phillips et al. (2020) models. For exoplanets, we
recommend using parametrized $T(P)$ profiles from, for example Guillot (2010),
as a first guess. Functionalities to browse these profiles directly are
already available in PICASO, via the guillot_pt(), and sonora() functions
666also see various tutorials e.g. Download and Query from Sonora Profile
Grid, and Accessing the pressure-temperature profile parameterizations.

Another first guess required to run the model is the extent of the deepest
convective zone. Specifically, the user needs to specify a guess of the
pressure level of the radiative-convective boundary in the atmosphere which
the code will then modify using methods described in §2.1.7. Even though this
guess is modified within the iteration, we recommend always setting the bottom
3-5 levels to be convective first. If the user sets too many layers to be
convective, the converged solutions may be inaccurate because the model is
only equipped to grow and merge convective zones, not shrink them. Therefore,
the best practice is to start with 3-5 bottom convective levels. If the code
runs into convergence issues after doing so, we recommended increasing this
number slowly to check if it helps converge the solution.

###  4.3. Convergence

In most cloud-free cases, the code will converge. However, there are a few
known cases where extra steps are needed to reach a converged solution. For
equilibrium chemistry brown dwarf models with $T_{\rm eff}$ between 1500
K–1700 K, an oscillating behaviour is seen in the iterations of the $T(P)$
profile. In these cases, the $T(P)$ profile oscillates between two profiles
with a constant temperature offset. This behaviour is not unexpected and is
caused by the sharp “cliff” in gaseous opacity at these effective
temperatures. This can be seen in the cross-section maps in Figure 5. All the
four panels in Figure 5 show a sharp change in gaseous opacity around the 900
K–1700 K temperature range. The overlap of this cliff with the converged brown
dwarf $T(P)$ profile at $T_{\rm eff}$=1700 K can also be seen in Figure 5. As
already mentioned in §2.1.4, this sharp opacity “cliff” makes the convergence
difficult around this temperature range. In order to overcome this, we
recommend performing multiple runs for the same object by recycling the
unconverged final $T(P)$ profile of each run as an initial guess profile for
the next run until a converged solution is reached.

A similar convergence issue also occurs for objects with $T_{\rm eff}$ $<$ 300
K. Figure 2 bottom left panel shows how H2O abundance in the vapor phase falls
drastically due to H2O condensation at such low temperatures. This sharp drop
in the H2O abundance causes the $T(P)$ profiles of these cold objects at the
low pressure regions to oscillate between multiple possible solutions. This
behaviour can be seen in the upper parts of the $T_{\rm eff}$=300 K $T(P)$
profile in Figure 5. We recommend two solutions to this problem: 1) either set
the minimum pressure of the model pressure grid to no less than 1 mbar, or 2)
use a SONORA profile as an initial guess and recycle the unconverged profile
multiple times ($\sim$ 4 times) to reach a converged state. The first
recommendation helps the user exclude the part of the $T(P)$ profile which
causes the instability in the convergence. As these objects are cold, the
excluded pressure region is extremely cold ($<$ 120 K), which means that there
is practically no contribution to the observables (e.g. total emergent flux).

##  5\. Calculating Observables

After using the radiative-convective equilibrium model to calculate the
atmospheric structure of a brown dwarf or an exoplanet, the user may want to
compute various observables. PICASO has been updated to include the capability
of calculating the 1D and 3D thermal emission spectrum, the 1D transmission
spectrum, and the 1D and 3D reflection spectroscopy, and phase curves. With a
1D $T(P)$ profile, users will likely want to produce 1D higher resolution
thermal emission spectra, reflection spectra, and/or transmission spectra of
the planet or brown dwarf. These calculated observables can be directly
compared with observational data.

Figure 15 shows an example of this with an exoplanet. The exoplanet has been
assumed to be a Jupiter mass and size planet with $T_{\rm int}$=300 K at a
distance of 0.1 AU around a Sun like star. The radiative-convective
equilibrium model described in this work has been used to calculate the
planet-wide average atmospheric structure of this hypothetical planet once
with chemical equilibrium and then again with chemical disequilibrium. The top
left panel in Figure 15 shows the $T(P)$ profile of this planet in each case.
It is clear that self-consistent treatment of disequilibrium chemistry also
impact the $T(P)$ profile of irradiated planets compared to equilibrium
chemistry models. We have then used the 1D radiative transfer routines of
PICASO to compute the reflection albedo spectrum of this exoplanet. This
albedo spectrum has been shown in Figure 15 top right panel. The planet is not
very reflective after 0.5 $\mu$m, after which brightness from Rayleigh
scattering tapers off. The transmission spectra of the planet is also modeled
with PICASO and is shown in Figure 15 bottom left panel. This panel shows that
the presence of disequilibrium chemistry can be clearly detected for a hot
Jupiter like this with its transmission spectra between 3-5 $\mu$m. The
thermal spectra are shown in the bottom right panel and the presence of
disequilibrium chemistry in hot Jupiters also can be easily detected with
their thermal spectra in the M-band, or in various NIRSpec and/or NIRCam JWST
modes. Note that for this illustrative example, we have used the planet-wide
average computed $T(P)$ profile for all these viewing geometries. However, a
more self-consistent method to compute these viewing geometries would be to
compute a day-side averaged profile for use in the zero-phase reflection and
thermal spectrum, and a planet-wide average for the transmission spectrum.

PICASO has also been equipped with simple modules to couple the atmospheric
models produced from this code with the evolutionary models from Marley et al.
(2021) to calculate absolute Vega magnitudes of the modeled object in any
filter of choice. The atmospheric model presented in this paper used with
these radiative transfer tool are already available in PICASO and can be
immensely helpful in modeling or planning observations using models.

Figure 15.— Using our atmospheric model to calculate observables. Top left
panel shows the converged $T(P)$ of an exoplanet calculated using our model
with (blue) and without (black) disequilibrium chemistry. Top right, bottom
left and bottom right panel shows the albedo spectrum, transmission spectrum,
and thermal spectrum at a spectral resolution of R=300 of the exoplanet
calculated using radiative transfer tools available in PICASO.

##  6\. Future Improvements

Currently, our publicly released model has the capability to model non-cloudy
brown dwarfs and exoplanets with equilibrium chemistry and disequilibrium
chemistry. However, several future improvements of the model are needed to
enhance its capability to capture a larger parameter space. Here we briefly
discuss these needed future improvements. We note that all PICASO development
is public, and open to community involvement.

###  6.1. Clouds

Our iterative climate model currently does not include the capability to treat
atmospheric clouds self-consistently. The EGP code already includes
condensation clouds with the Ackerman & Marley (2001) model and has been used
to model cloudy L-dwarf atmospheres (e.g. Cushing et al., 2008; Stephens et
al., 2009). Clouds can become important opacity sources in brown dwarf and
exoplanetary atmospheres due to their strong scattering tendencies. This also
has large effects on the $T(P)$ profile of these atmospheres (Morley et al.,
2014, 2012). We plan to couple this model with the Python-based cloud model
VIRGA (Batalha et al., 2020). Recently, VIRGA has been updated to include
variable sedimentation efficiency (fsed) with height (Rooney et al., 2021). We
will couple this updated VIRGA model with our model so that our models can be
used to model hotter L-dwarfs which are generally assumed to be cloudy and
also much colder Y-dwarfs with H2O clouds. This improvement will also enable
us to apply our model to cloudy exoplanets and brown dwarfs.

###  6.2. 1D Chemical Kinetics Model

We use the quench time approximation to model the effects of vertical mixing
on the chemical abundance profiles in this model. This approximation is based
on the parametrized mixing timescales from Zahnle & Marley (2014) outlined in
§2.1.5. However, more robust 1D chemical kinetics models of treating
disequilibrium chemistry caused by vertical mixing are now open-sourced like
the VULCAN model (Tsai et al., 2017, 2021). We plan to couple our Python model
with the VULCAN chemical kinetics model to increase the flexibility of our
model. This will make our disequilibrium chemistry models more robust.
Moreover, VULCAN also includes the capability to treat stellar irradiation
induced photochemistry. This will lead to a significant improvement for
applicability of our model to irradiated exoplanets where photochemistry can
largely impact atmospheric chemistry especially at lower pressures. This
update will allow us to explore the impact of photochemistry on the
atmospheric structure of exoplanets. We have already made progress for this
coupling but still need to test and benchmark several aspects of the coupled
code before this update is publicly available for use.

###  6.3. Time Evolution Version

Recently, Mayorga et al. (2021); Robinson & Marley (2014) enhanced the EGP
model to EGP+ by including a time-stepping version of the code which can model
the dynamic temporal response of the atmosphere to time-varying physical
conditions like stellar irradiation changing with time. This is especially
relevant to planets with highly eccentric orbits as the stellar irradiation
can change by a large amount within one orbital time-period for these planets.
We plan to adapt this improvement within our code in the near future as well.
Studying the atmospheres of eccentric exoplanets is a promising research area
and a subject of future JWST observing campaigns.

##  7\. Conclusions and Summary

We presented a new open-source Python based 1D radiative–convective
equilibrium model as part of the PICASO package. This code derives its
heritage from the EGP code developed by Marley & McKay (1999) based upon a
Titan atmosphere model developed by McKay et al. (1989). The EGP code has been
used to model brown dwarf and exoplanetary atmospheres for almost three
decades now. The model is applicable to H-dominated atmospheres of both brown
dwarfs and irradiated exoplanets. The model includes the capability to do
calculations with both equilibrium chemistry and disequilibrium chemistry due
to vertical mixing. The model includes options to use $K_{zz}$ values constant
with height while performing disequilibrium chemistry runs. We have also
included the capability to use self-consistent prescriptions of $K_{zz}$
within the model where the $K_{zz}$ will also iterate along with the $T(P)$
profile and atmospheric chemistry to ultimately reach a converged atmospheric
state.

We have benchmarked this model with publicly available models from the SONORA
BOBCAT grid and also with results from an independent model used by Phillips
et al. (2020) to produce the ATMO 2020 atmospheric grid. For irradiated
planets, the PICASO 3.0 models were benchmarked against the ATMO grid
presented in Goyal et al. (2020), using the hot Saturn WASP 25-b as a test
case. The chemical equilibrium version of the model was benchmarked both with
models from the SONORA grid, ATMO 2020, and the ATMO models for WASP 25-b. The
chemical disequilibrium models were benchmarked against models from the
disequilibrium chemistry atmospheric models from the ATMO 2020 grid of models.
This benchmarking analysis showed excellent agreement with PICASO.

Our model is open-source and publicly available for the community (we include
several in text links to code throughout this manuscript). Additionally, we
outlined many recommendations in this work for proper usage of the model. This
includes recommendations on choosing the atmospheric pressure grid for a
science object, choosing the initial guess $T(P)$ profile for a particular run
and also the various parts of the parameter space where the model is known to
face convergence issues. We also included ways of resolving these convergence
issues.

We will also release tutorials to apply this model for various science cases
with the code, upon paper acceptance. We plan to actively develop this model
further to include clouds, couple it with 1D chemical kinetics codes for
better robustness and also include time dependent effects like variable
stellar irradiation within our model.

##  8\. Acknowledgments

SM thanks the UC Regents Fellowship award for supporting him for this work.
NEB acknowledges support from NASA Astrophysics Division. NEB and JJF
acknowledge support from NASA’S Interdisciplinary Consortia for Astrobiology
Research (NNH19ZDA001N-ICAR) under award number 19-ICAR19_2-0041. JJF and MSM
acknowledge the support of NASA XRP grant 80NSSC19K0446. We thank the
anonymous referee for helping us to improve the manuscript significantly.

Software: PICASO (Batalha et al., 2019)777upon acceptance we will formally
release PICASO 3.0 and update this with the Zenodo link, pandas (McKinney,
2010), NumPy (Walt et al., 2011), IPython (Pérez & Granger, 2007), Jupyter
(Kluyver et al., 2016), matplotlib (Hunter, 2007)

## References

  * Abramowitz & Stegun (1970) Abramowitz, M., & Stegun, I. A. 1970, Handbook of mathematical functions : with formulas, graphs, and mathematical tables
  * Ackerman & Marley (2001) Ackerman, A. S., & Marley, M. S. 2001, ApJ, 556, 872
  * Allard et al. (2007a) Allard, F., Allard, N. F., Homeier, D., et al. 2007a, A&A, 474, L21
  * Allard et al. (2012) Allard, F., Homeier, D., & Freytag, B. 2012, Philosophical Transactions of the Royal Society of London Series A, 370, 2765
  * Allard et al. (2007b) Allard, N. F., Kielkopf, J. F., & Allard, F. 2007b, European Physical Journal D, 44, 507
  * Allard et al. (2016) Allard, N. F., Spiegelman, F., & Kielkopf, J. F. 2016, A&A, 589, A21
  * Allard et al. (2019) Allard, N. F., Spiegelman, F., Leininger, T., & Molliere, P. 2019, A&A, 628, A120
  * Allen et al. (1981) Allen, M., Yung, Y. L., & Waters, J. W. 1981, Journal of Geophysical Research: Space Physics, 86, 3617
  * Amundsen et al. (2017) Amundsen, D. S., Tremblin, P., Manners, J., Baraffe, I., & Mayne, N. J. 2017, A&A, 598, A97
  * Arcangeli et al. (2018) Arcangeli, J., Désert, J.-M., Line, M. R., et al. 2018, ApJ, 855, L30
  * Azzam et al. (2016) Azzam, A. A. A., Tennyson, J., Yurchenko, S. N., & Naumenko, O. V. 2016, Monthly Notices of the Royal Astronomical Society, 460, 4063
  * Barber et al. (2014) Barber, R. J., Strange, J. K., Hill, C., et al. 2014, MNRAS, 437, 1828
  * Bard et al. (1991) Bard, A., Kock, A., & Kock, M. 1991, A&A, 248, 315
  * Bard & Kock (1994) Bard, A., & Kock, M. 1994, A&A, 282, 1014
  * Barman et al. (2001) Barman, T. S., Hauschildt, P. H., & Allard, F. 2001, ApJ, 556, 885
  * Barton et al. (2013) Barton, E. J., Yurchenko, S. N., & Tennyson, J. 2013, MNRAS, 434, 1469
  * Batalha et al. (2020) Batalha, N., caoimherooney11, & sagnickm. 2020, natashabatalha/virga: Initial Release, doi:10.5281/zenodo.3759888
  * Batalha et al. (2019) Batalha, N. E., Marley, M. S., Lewis, N. K., & Fortney, J. J. 2019, ApJ, 878, 70
  * Baudino et al. (2015) Baudino, J. L., Bézard, B., Boccaletti, A., et al. 2015, A&A, 582, A83
  * Bell (1980) Bell, K. 1980, Journal of Physics B: Atomic and Molecular Physics (1968-1987), 13, 1859
  * Bell & Berrington (1987) Bell, K. L., & Berrington, K. A. 1987, Journal of Physics B: Atomic and Molecular Physics, 20, 801
  * Bittner & Bernath (2018) Bittner, D. M., & Bernath, P. F. 2018, VizieR Online Data Catalog, J/ApJS/235/8
  * Burningham et al. (2017) Burningham, B., Marley, M. S., Line, M. R., et al. 2017, MNRAS, 470, 1177
  * Burningham et al. (2021) Burningham, B., Faherty, J. K., Gonzales, E. C., et al. 2021, MNRAS, 506, 1944
  * Burrows et al. (2008) Burrows, A., Budaj, J., & Hubeny, I. 2008, ApJ, 678, 1436
  * Burrows et al. (2002) Burrows, A., Ram, R. S., Bernath, P., Sharp, C. M., & Milsom, J. A. 2002, The Astrophysical Journal, 577, doi:10.1086/342242
  * Burrows et al. (2006) Burrows, A., Sudarsky, D., & Hubeny, I. 2006, ApJ, 640, 1063
  * Burrows et al. (1997) Burrows, A., Marley, M., Hubbard, W. B., et al. 1997, ApJ, 491, 856
  * Chabrier et al. (2000) Chabrier, G., Baraffe, I., Allard, F., & Hauschildt, P. 2000, ApJ, 542, 464
  * Coppola et al. (2011) Coppola, C. M., Lodi, L., & Tennyson, J. 2011, MNRAS, 415, 487
  * Cushing et al. (2008) Cushing, M. C., Marley, M. S., Saumon, D., et al. 2008, ApJ, 678, 1372
  * Deline et al. (2022) Deline, A., Hooton, M. J., Lendl, M., et al. 2022, A&A, 659, A74
  * Dulick et al. (2003) Dulick, M., Bauschlicher, C. W., J., Burrows, A., et al. 2003, ApJ, 594, 651
  * Enoch et al. (2011) Enoch, B., Cameron, A. C., Anderson, D. R., et al. 2011, MNRAS, 410, 1631
  * Fortney et al. (2008) Fortney, J. J., Lodders, K., Marley, M. S., & Freedman, R. S. 2008, ApJ, 678, 1419
  * Fortney et al. (2007) Fortney, J. J., Marley, M. S., & Barnes, J. W. 2007, ApJ, 659, 1661
  * Fortney et al. (2005) Fortney, J. J., Marley, M. S., Lodders, K., Saumon, D., & Freedman, R. 2005, ApJ, 627, L69
  * Fortney et al. (2020) Fortney, J. J., Visscher, C., Marley, M. S., et al. 2020, AJ, 160, 288
  * Fossati et al. (2021) Fossati, L., Young, M. E., Shulyak, D., et al. 2021, A&A, 653, A52
  * Freedman et al. (2014) Freedman, R. S., Lustig-Yaeger, J., Fortney, J. J., et al. 2014, ApJS, 214, 25
  * Freedman et al. (2008) Freedman, R. S., Marley, M. S., & Lodders, K. 2008, ApJS, 174, 504
  * Fuhr et al. (1988) Fuhr, J. R., Martin, G. A., & Wiese, W. L. 1988, Journal of Physical and Chemical Reference Data, 17
  * Gandhi & Madhusudhan (2017) Gandhi, S., & Madhusudhan, N. 2017, MNRAS, 472, 2334
  * Geballe et al. (2009) Geballe, T. R., Saumon, D., Golimowski, D. A., et al. 2009, ApJ, 695, 844
  * Gharib-Nezhad et al. (2021) Gharib-Nezhad, E., Iyer, A. R., Line, M. R., et al. 2021, ApJS, 254, 34
  * GharibNezhad et al. (2013) GharibNezhad, E., Shayesteh, A., & Bernath, P. F. 2013, MNRAS, 432, 2043
  * Giacobbe et al. (2021) Giacobbe, P., Brogi, M., Gandhi, S., et al. 2021, Nature, 592, 205
  * Gierasch & Conrath (1985) Gierasch, P. J., & Conrath, B. J. 1985, Energy conversion processes in the outer planets., ed. G. E. Hunt, 121–146
  * Golimowski et al. (2004) Golimowski, D. A., Leggett, S. K., Marley, M. S., et al. 2004, AJ, 127, 3516
  * Gordon et al. (2017) Gordon, I. E., Rothman, L. S., Tan, Y., Kochanov, R. V., & Hill, C. 2017, in 72nd International Symposium on Molecular Spectroscopy, TJ08
  * Gordon et al. (2022) Gordon, I. E., Rothman, L. S., Hargreaves, R. J., et al. 2022, J. Quant. Spec. Radiat. Transf., 277, 107949
  * Gordon & McBride (1994) Gordon, S., & McBride, B. J. 1994, Computer program for calculation of complex chemical equilibrium compositions and applications. Part 1: Analysis, Tech. rep.
  * Goyal et al. (2018) Goyal, J. M., Mayne, N., Sing, D. K., et al. 2018, MNRAS, 474, 5158
  * Goyal et al. (2020) Goyal, J. M., Mayne, N., Drummond, B., et al. 2020, MNRAS, 498, 4680
  * Guillot (2010) Guillot, T. 2010, A&A, 520, A27
  * Hargreaves et al. (2010) Hargreaves, R. J., Hinkle, K. H., Bauschlicher, Charles W., J., et al. 2010, AJ, 140, 919
  * Harris et al. (2006) Harris, G. J., Tennyson, J., Kaminsky, B. M., Pavlenko, Y. V., & Jones, H. R. A. 2006, MNRAS, 367, 400
  * Hoeijmakers et al. (2018) Hoeijmakers, H. J., Ehrenreich, D., Heng, K., et al. 2018, Nature, 560, 453
  * Hoeijmakers et al. (2020) Hoeijmakers, H. J., Seidel, J. V., Pino, L., et al. 2020, Astronomy & Astrophysics, 641, A123
  * Huang et al. (2014) Huang, X., Gamache, R. R., Freedman, R. S., Schwenke, D. W., & Lee, T. J. 2014, Journal of Quantitative Spectroscopy and Radiative Transfer, 147, 134
  * Hubbard (1977) Hubbard, W. B. 1977, Icarus, 30, 305
  * Hubeny (1988) Hubeny, I. 1988, Computer Physics Communications, 52, 103
  * Hubeny (2017) Hubeny, I. 2017, Monthly Notices of the Royal Astronomical Society, 469, 841
  * Hubeny et al. (2003) Hubeny, I., Burrows, A., & Sudarsky, D. 2003, ApJ, 594, 1011
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * John (1988) John, T. L. 1988, A&A, 193, 189
  * Karalidi et al. (2021) Karalidi, T., Marley, M., Fortney, J. J., et al. 2021, ApJ, 923, 269
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in ELPUB, 87–90
  * Lacy & Burrows (2020) Lacy, B., & Burrows, A. 2020, ApJ, 892, 151
  * Lacy et al. (2019) Lacy, B., Shlivko, D., & Burrows, A. 2019, AJ, 157, 132
  * Lam et al. (2015) Lam, S. K., Pitrou, A., & Seibert, S. 2015, in Proceedings of the Second Workshop on the LLVM Compiler Infrastructure in HPC, LLVM ’15 (New York, NY, USA: ACM), 7:1–7:6
  * Lee et al. (2021) Lee, E. K. H., Parmentier, V., Hammond, M., et al. 2021, Monthly Notices of the Royal Astronomical Society, 506, 2695
  * Leggett et al. (2012) Leggett, S. K., Saumon, D., Marley, M. S., et al. 2012, ApJ, 748, 74
  * Lenzuni et al. (1991) Lenzuni, P., Chernoff, D. F., & Salpeter, E. E. 1991, ApJS, 76, 759
  * Li et al. (2015) Li, G., Gordon, I. E., Rothman, L. S., et al. 2015, ApJS, 216, 15
  * Line et al. (2017) Line, M. R., Marley, M. S., Liu, M. C., et al. 2017, ApJ, 848, 83
  * Lodders (1999) Lodders, K. 1999, ApJ, 519, 793
  * Lodders (2002) Lodders, K. 2002, The Astrophysical Journal, 577, 974
  * Lodders (2010) Lodders, K. 2010, in Formation and Evolution of Exoplanets, 157
  * Lupu et al. (2021) Lupu, R., Freedman, R., Gharib-Nezhad, E., Visscher, C., & Molliere, P. 2021, Correlated k coefficients for H2-He atmospheres; 196 spectral windows and 1460 pressure-temperature points, doi:10.5281/zenodo.5590989
  * Madhusudhan (2018) Madhusudhan, N. 2018, in Handbook of Exoplanets, ed. H. J. Deeg & J. A. Belmonte, 104
  * Malik et al. (2019a) Malik, M., Kempton, E. M. R., Koll, D. D. B., et al. 2019a, ApJ, 886, 142
  * Malik et al. (2019b) Malik, M., Kitzmann, D., Mendonça, J. M., et al. 2019b, AJ, 157, 170
  * Malik et al. (2017) Malik, M., Grosheintz, L., Mendonça, J. M., et al. 2017, AJ, 153, 56
  * Marley & McKay (1999) Marley, M. S., & McKay, C. P. 1999, Icarus, 138, 268
  * Marley et al. (2012) Marley, M. S., Saumon, D., Cushing, M., et al. 2012, The Astrophysical Journal, 754, 135
  * Marley et al. (1996) Marley, M. S., Saumon, D., Guillot, T., et al. 1996, Science, 272, 1919
  * Marley et al. (2021) Marley, M. S., Saumon, D., Visscher, C., et al. 2021, ApJ, 920, 85
  * Mayorga et al. (2021) Mayorga, L. C., Robinson, T. D., Marley, M. S., May, E. M., & Stevenson, K. B. 2021, ApJ, 915, 41
  * McKay et al. (1989) McKay, C. P., Pollack, J. B., & Courtin, R. 1989, Icarus, 80, 23
  * McKemmish et al. (2019) McKemmish, L. K., Masseron, T., Hoeijmakers, H. J., et al. 2019, MNRAS, 488, 2836
  * McKemmish et al. (2016) McKemmish, L. K., Yurchenko, S. N., & Tennyson, J. 2016, MNRAS, 463, 771
  * McKinney (2010) McKinney, W. 2010, in Proceedings of the 9th Python in Science Conference, ed. S. van der Walt & J. Millman, 51 – 56
  * Menou & Rauscher (2009) Menou, K., & Rauscher, E. 2009, ApJ, 700, 887
  * Miles et al. (2020) Miles, B. E., Skemer, A. J. I., Morley, C. V., et al. 2020, AJ, 160, 63
  * Mizus et al. (2017) Mizus, I. I., Alijah, A., Zobov, N. F., et al. 2017, MNRAS, 468, 1717
  * Morley et al. (2012) Morley, C. V., Fortney, J. J., Marley, M. S., et al. 2012, ApJ, 756, 172
  * Morley et al. (2015) Morley, C. V., Fortney, J. J., Marley, M. S., et al. 2015, The Astrophysical Journal, 815, 110
  * Morley et al. (2017) Morley, C. V., Knutson, H., Line, M., et al. 2017, AJ, 153, 86
  * Morley et al. (2014) Morley, C. V., Marley, M. S., Fortney, J. J., et al. 2014, ApJ, 787, 78
  * Moses et al. (2005) Moses, J. I., Fouchet, T., Bézard, B., et al. 2005, Journal of Geophysical Research (Planets), 110, E08001
  * Moses et al. (2013) Moses, J. I., Madhusudhan, N., Visscher, C., & Freedman, R. S. 2013, ApJ, 763, 25
  * Moses et al. (2021) Moses, J. I., Tremblin, P., Venot, O., & Miguel, Y. 2021, Experimental Astronomy, arXiv:2103.07023
  * Moses et al. (2011) Moses, J. I., Visscher, C., Fortney, J. J., et al. 2011, ApJ, 737, 15
  * Mukherjee et al. (2021) Mukherjee, S., Batalha, N. E., & Marley, M. S. 2021, ApJ, 910, 158
  * Noll et al. (1997) Noll, K. S., Geballe, T. R., & Marley, M. S. 1997, ApJ, 489, L87
  * O’Brian et al. (1991) O’Brian, T. R., Wickliffe, M. E., Lawler, J. E., Whaling, W., & Brault, J. W. 1991, Journal of the Optical Society of America B Optical Physics, 8, 1185
  * Oppenheimer et al. (1998) Oppenheimer, B. R., Kulkarni, S. R., Matthews, K., & van Kerkwijk, M. H. 1998, ApJ, 502, 932
  * Parmentier et al. (2013) Parmentier, V., Showman, A. P., & Lian, Y. 2013, A&A, 558, A91
  * Pérez & Granger (2007) Pérez, F., & Granger, B. E. 2007, Computing in Science & Engineering, 9
  * Phillips et al. (2020) Phillips, M. W., Tremblin, P., Baraffe, I., et al. 2020, A&A, 637, A38
  * Piskorz et al. (2018) Piskorz, D., Buzard, C., Line, M. R., et al. 2018, The Astronomical Journal, 156, 133
  * Polyansky et al. (2018) Polyansky, O. L., Kyuberis, A. A., Zobov, N. F., et al. 2018, MNRAS, 480, 2597
  * Pontoppidan et al. (2022) Pontoppidan, K., Blome, C., Braun, H., et al. 2022, arXiv e-prints, arXiv:2207.13067
  * Robinson & Marley (2014) Robinson, T. D., & Marley, M. S. 2014, ApJ, 785, 158
  * Roman & Rauscher (2019) Roman, M., & Rauscher, E. 2019, ApJ, 872, 1
  * Rooney et al. (2021) Rooney, C. M., Batalha, N. E., Gao, P., & Marley, M. S. 2021, arXiv e-prints, arXiv:2110.05903
  * Rothman et al. (2010) Rothman, L. S., Gordon, I. E., Barber, R. J., et al. 2010, Journal of Quantitative Spectroscopy and Radiative Transfer, 111, doi:10.1016/j.jqsrt.2010.05.001
  * Rothman et al. (2013) Rothman, L. S., Gordon, I. E., Babikov, Y., et al. 2013, Journal of Quantitative Spectroscopy and Radiative Transfer, 130, doi:10.1016/j.jqsrt.2013.07.002
  * Ryabchikova et al. (2015) Ryabchikova, T., Piskunov, N., Kurucz, R. L., et al. 2015, Physica Scripta, 90, 054005
  * Saumon & Marley (2008) Saumon, D., & Marley, M. S. 2008, ApJ, 689, 1327
  * Saumon et al. (2012) Saumon, D., Marley, M. S., Abel, M., Frommhold, L., & Freedman, R. S. 2012, ApJ, 750, 74
  * Saumon et al. (2003) Saumon, D., Marley, M. S., & Lodders, K. 2003, arXiv e-prints, astro
  * Seager & Sasselov (1998) Seager, S., & Sasselov, D. D. 1998, ApJ, 502, L157
  * Sheppard et al. (2017) Sheppard, K. B., Mandell, A. M., Tamburo, P., et al. 2017, ApJ, 850, L32
  * Showman et al. (2008) Showman, A. P., Cooper, C. S., Fortney, J. J., & Marley, M. S. 2008, ApJ, 682, 559
  * Showman et al. (2009) Showman, A. P., Fortney, J. J., Lian, Y., et al. 2009, ApJ, 699, 564
  * Showman et al. (2020) Showman, A. P., Tan, X., & Parmentier, V. 2020, Space Sci. Rev., 216, 139
  * Sorahana & Yamamura (2012) Sorahana, S., & Yamamura, I. 2012, The Astrophysical Journal, 760, 151
  * Sousa-Silva et al. (2014) Sousa-Silva, C., Al-Refaie, A. F., Tennyson, J., & Yurchenko, S. N. 2014, Monthly Notices of the Royal Astronomical Society, 446, 2337
  * Southworth et al. (2014) Southworth, J., Hinse, T. C., Burgdorf, M., et al. 2014, Monthly Notices of the Royal Astronomical Society, 444, 776
  * Spiegel & Burrows (2012) Spiegel, D. S., & Burrows, A. 2012, ApJ, 745, 174
  * Spiegel et al. (2010) Spiegel, D. S., Burrows, A., Ibgui, L., Hubeny, I., & Milsom, J. A. 2010, ApJ, 709, 149
  * Stephens et al. (2009) Stephens, D. C., Leggett, S. K., Cushing, M. C., et al. 2009, ApJ, 702, 154
  * STScI Development Team (2013) STScI Development Team. 2013, pysynphot: Synthetic photometry software package, ascl:1303.023
  * Sudarsky et al. (2003) Sudarsky, D., Burrows, A., & Hubeny, I. 2003, ApJ, 588, 1121
  * Tan (2022) Tan, X. 2022, Monthly Notices of the Royal Astronomical Society, 511, 4861
  * Tan & Showman (2019) Tan, X., & Showman, A. P. 2019, ApJ, 874, 111
  * Tan & Showman (2021) Tan, X., & Showman, A. P. 2021, Monthly Notices of the Royal Astronomical Society, 502, 2198
  * Taylor et al. (2021) Taylor, J., Parmentier, V., Line, M. R., et al. 2021, MNRAS, 506, 1309
  * Toon et al. (1989) Toon, O. B., McKay, C. P., Ackerman, T. P., & Santhanam, K. 1989, J. Geophys. Res., 94, 16287
  * Tremblin et al. (2015) Tremblin, P., Amundsen, D. S., Mourier, P., et al. 2015, ApJ, 804, L17
  * Tsai et al. (2017) Tsai, S.-M., Lyons, J. R., Grosheintz, L., et al. 2017, ApJS, 228, 20
  * Tsai et al. (2021) Tsai, S.-M., Malik, M., Kitzmann, D., et al. 2021, ApJ, 923, 264
  * Visscher et al. (2010) Visscher, C., Lodders, K., & Fegley, B. 2010, The Astrophysical Journal, 716, 1060
  * Visscher et al. (2006) Visscher, C., Lodders, K., & Fegley, Bruce, J. 2006, ApJ, 648, 1181
  * Visscher et al. (2010) —. 2010, ApJ, 716, 1060
  * Walt et al. (2011) Walt, S. v. d., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science & Engineering, 13, 22
  * Wang et al. (2015) Wang, D., Gierasch, P. J., Lunine, J. I., & Mousis, O. 2015, Icarus, 250, 154
  * Wilzewski et al. (2016) Wilzewski, J. S., Gordon, I. E., Kochanov, R. V., Hill, C., & Rothman, L. S. 2016, Journal of Quantitative Spectroscopy and Radiative Transfer, 168, 193
  * Wolf et al. (2022) Wolf, E. T., Kopparapu, R., Haqq-Misra, J., & Fauchez, T. J. 2022, The Planetary Science Journal, 3, 7
  * Yadin et al. (2012) Yadin, B., Veness, T., Conti, P., et al. 2012, MNRAS, 425, 34
  * Yan et al. (2022) Yan, F., Reiners, A., Pallé, E., et al. 2022, A&A, 659, A7
  * Yurchenko et al. (2011) Yurchenko, S. N., Barber, R. J., & Tennyson, J. 2011, Monthly Notices of the Royal Astronomical Society, 413, 1828
  * Yurchenko & Tennyson (2014) Yurchenko, S. N., & Tennyson, J. 2014, Monthly Notices of the Royal Astronomical Society, 440, 1649
  * Yurchenko et al. (2013) Yurchenko, S. N., Tennyson, J., Barber, R. J., & Thiel, W. 2013, Journal of Molecular Spectroscopy, 291, 69
  * Zahnle & Marley (2014) Zahnle, K. J., & Marley, M. S. 2014, ApJ, 797, 41
  * Zhang & Showman (2018) Zhang, X., & Showman, A. P. 2018, ApJ, 866, 1
  * Zhang et al. (2021) Zhang, Z., Liu, M. C., Marley, M. S., Line, M. R., & Best, W. M. J. 2021, ApJ, 921, 95