# The TRAPPIST-1 Habitable Atmosphere Intercomparison (THAI). I. Dry Cases-The Fellowship of the GCMs

Martin Turbet, Thomas J. Fauchez, Denis E. Sergeev, Ian A. Boutle, Kostas

Tsigaridis, Michael J. Way, Eric T. Wolf, Shawn D. Domagal-Goldman,

François Forget, Jacob Haqq-Misra, et al.

► To cite this version:

Martin Turbet, Thomas J. Fauchez, Denis E. Sergeev, Ian A. Boutle, Kostas Tsigaridis, et al. The TRAPPIST-1 Habitable Atmosphere Intercomparison (THAI). I. Dry Cases-The Fellowship of the GCMs. The Planetary Science Journal, 2022, 3, 10.3847/PSJ/ac6cf0. insu-03847094

HAL Id: insu-03847094

https://insu.hal.science/insu-03847094

Submitted on 10 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The TRAPPIST-1 Habitable Atmosphere Intercomparison (THAI). I. Dry Cases—The Fellowship of the GCMs

Martin Turbet<fnref n="1" />,<fnref n="2" /> ID, Thomas J. Fauchez<fnref n="3" />,<fnref n="4" />,⁵ ID, Denis E. Sergeev<fnref n="6" /> ID, Ian A. Boutle⁷,<fnref n="8" /> ID, Kostas Tsigaridis⁹,<fnref n="10" /> ID, Michael J. Way⁵,<fnref n="10" />,<fnref n="11" /> ID, Eric T. Wolf⁵,<fnref n="12" />,¹³ ID, Shawn D. Domagal-Goldman<fnref n="3" />,⁵,¹³ ID, François Forget<fnref n="1" /> ID, Jacob Haqq-Misra¹³,<fnref n="14" /> ID, Ravi K. Kopparapu<fnref n="3" />,⁵,¹³ ID, F. Hugo Lambert<fnref n="6" /> ID, James Manners⁷ ID, Nathan J. Mayne<fnref n="8" /> ID, and Linda Sohl⁹,<fnref n="10" /> ID

Received 2021 September 24; revised 2022 April 9; accepted 2022 May 3; published 2022 September 15

## Abstract

With the commissioning of powerful, new-generation telescopes such as the James Webb Space Telescope (JWST) and the ground-based Extremely Large Telescopes, the first characterization of a high molecular weight atmosphere around a temperate rocky exoplanet is imminent. Atmospheric simulations and synthetic observables of target exoplanets are essential to prepare and interpret these observations. Here we report the results of the first part of the TRAPPIST-1 Habitable Atmosphere Intercomparison (THAI) project, which compares 3D numerical simulations performed with four state-of-the-art global climate models (ExoCAM, LMD-Generic, ROCKE-3D, Unified Model) for the potentially habitable target TRAPPIST-1e. In this first part, we present the results of dry atmospheric simulations. These simulations serve as a benchmark to test how radiative transfer, subgrid-scale mixing (dry turbulence and convection), and large-scale dynamics impact the climate of TRAPPIST-1e and consequently the transit spectroscopy signature as seen by JWST. To first order, the four models give results in good agreement. The intermodel spread in the global mean surface temperature amounts to 7 K (6 K) for the N₂-dominated (CO₂-dominated) atmosphere. The radiative fluxes are also remarkably similar (intermodel variations less than 5%), from the surface (1 bar) up to atmospheric pressures ~5 mbar. Moderate differences between the models appear in the atmospheric circulation pattern (winds) and the (stratospheric) thermal structure. These differences arise between the models from (1) large-scale dynamics, because TRAPPIST-1e lies at the tipping point between two different circulation regimes (fast and Rhines rotators) in which the models can be alternatively trapped, and (2) parameterizations used in the upper atmosphere such as numerical damping.

*Unified Astronomy Thesaurus concepts: Planetary climates (2184); Exoplanet atmospheres (487); Habitable planets (695); Atmospheric circulation (112); Extrasolar rocky planets (511); Hydrodynamical simulations (767); M dwarf stars (982); Radiative transfer (1335)*

## 1. Introduction

"Don't adventures ever have an end? I suppose not. Someone else always has to carry on the story." —J. R. R. Tolkien, *The Fellowship of the Ring* (1954).

With the commissioning of next-generation telescopes and instruments, the astronomical community is for the first time in the process of detecting and characterizing the chemical

composition of the atmospheres of rocky exoplanets receiving a moderate instellation (Greene et al. 2016; Morley et al. 2017; Wunderlich et al. 2019; Gillon et al. 2020; Turbet et al. 2020). Specifically, observations of the TRAPPIST-1 planets (Gillon et al. 2017) with the James Webb Space Telescope (JWST) are likely to be our first real opportunity to characterize the surface and atmosphere of planets with a similar mass, radius, and bolometric irradiation to Earth but orbiting other stars than the Sun (Morley et al. 2017; Fauchez et al. 2019; Lustig-Yaeger et al. 2019; Wunderlich et al. 2019, 2020; Gillon et al. 2020; Turbet et al. 2020).

Preparing for these observations with sophisticated numerical climate models has recently become an active topic of research aimed at optimizing the scientific return from observational campaigns (see Fauchez et al. 2021 and references therein). First, model predictions guide observation proposals, and second, they

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ Laboratoire de Météorologie Dynamique/IPSL, CNRS, Sorbonne Université, École Normale Supérieure, PSL Research University, École Polytechnique, F-75005 Paris, France; martin.turbet@lmd.ipsl.fr</footnote>

<footnote n="2">² Observatoire Astronomique de l'Université de Genève, Université de Genève, Chemin des Maillettes 51, 1290 Versoix, Switzerland</footnote>

<footnote n="3">³ NASA Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20771, USA</footnote>

<footnote n="4">⁴ Goddard Earth Sciences Technology and Research (GESTAR), Universities Space Research Association (USRA), Columbia, MD 7178, USA
⁵ NASA GSFC Sellers Exoplanet Environments Collaboration, USA</footnote>

<footnote n="6">⁶ Department of Mathematics, College of Engineering, Mathematics, and Physical Sciences, University of Exeter, Exeter, EX4 4QF, UK
⁷ Met Office, FitzRoy Road, Exeter, EX1 3PB, UK</footnote>

<footnote n="8">⁸ Department of Astrophysics, College of Engineering, Mathematics, and Physical Sciences, University of Exeter, Exeter, EX4 4QL, UK
⁹ Center for Climate Systems Research, Columbia University, New York, NY 10025, USA</footnote>

<footnote n="10">¹⁰ NASA Goddard Institute for Space Studies, 2880 Broadway, New York, NY 10025, USA</footnote>

<footnote n="11">¹¹ Theoretical Astrophysics, Department of Physics and Astronomy, Uppsala University, Uppsala, Sweden</footnote>

<footnote n="12">¹² Laboratory for Atmospheric and Space Physics, University of Colorado Boulder, Boulder, CO 80303, USA
¹³ NASA NExSS Virtual Planetary Laboratory, Seattle, WA, 98195, USA</footnote>

<footnote n="14">¹⁴ Blue Marble Space Institute of Science, Seattle, WA 98104, USA</footnote>
are used to help in the analysis and interpretation of the data. However, little work has been done so far to identify sources of uncertainty between atmospheric models employed in the exoplanet community, as well as how this may affect the feasibility and the interpretation of observations of exoplanet atmospheres. To this end, global climate model intercomparison projects, which have been pioneered and successfully used to study the past, present, and future evolution of the Earth's climate (e.g., Eyring et al. 2016), have also been used more recently in the context of exoplanets (e.g., Yang et al. 2019).

The TRAPPIST-1 Habitable Atmosphere Intercomparison (THAI) project (Fauchez et al. 2020) brings together four state-of-the-art 3D global climate models to simulate the climate and observability with JWST of TRAPPIST-1e under various atmospheric composition scenarios. TRAPPIST-1e is to date the best potentially habitable exoplanet target we have (Wolf et al. 2017; Turbet et al. 2018; Fauchez et al. 2019; Sergeev et al. 2020). As of 2021 September, four transit observations of TRAPPIST-1e are already planned in the first observation cycle of JWST using the NIRSpec instrument (GTO Program 1331; PI: Nikole Lewis).

The goal of the THAI intercomparison project is twofold. The first aim is to identify the extent to which, for given atmospheric composition scenarios, the different models agree or disagree on the representation of the climate of TRAPPIST-1e. This first investigation is key to identifying the processes responsible for most of the disparities between models, and therefore where the community must deploy significant effort. The second aim is to assess how uncertainties in climate modeling affect observables, which provides a degree of confidence in our interpretation of observations of the atmospheric composition and climate of TRAPPIST-1e.

The results of the THAI intercomparison project are presented in three parts detailed in three distinct manuscripts. In the first part (this manuscript), we present and compare the results of dry (i.e., without water cycle and without clouds) 3D numerical climate simulations of TRAPPIST-1e performed with four GCMs (ExoCAM, LMD-Generic, ROCKE-3D, UM). These simulations are used as benchmark cases to better understand the differences in the water cycle (moist convection, clouds), described in Sergeev et al. (2022a, referred to as part 2 of THAI), and their effects on the synthetic observations, described in Fauchez et al. (2022, referred to as part 3 of THAI). The benchmark simulations presented here were designed to test the impact of the radiative transfer (RT), the dynamical core, the boundary-layer parameterization, and the subgrid-scale dry dynamics on the mean climate state and ultimately on the synthetic observables of TRAPPIST-1e (Fauchez et al. 2022).

The manuscript is structured as follows: In Section 2 the methods and tools used in this study are presented, with particular emphasis on the description of the parameterizations (RT, dynamical core, subgrid dynamics) used. In Section 3, we present the results of this first intercomparison phase on the climate states of TRAPPIST-1e and discuss the possible source of differences between the models. Finally, conclusions are given in Section 4.

## 2. Methods

### 2.1. Planetary Configuration

In this part 1 of the THAI trilogy, we present results of the dry atmosphere simulations, named Ben 1 and Ben 2 in the THAI protocol (Fauchez et al. 2020). The Ben 1 and Ben 2

<table><caption>**Table 1**
TRAPPIST-1 Stellar Spectrum and Planetary Parameters of TRAPPIST-1e (Grimm et al. 2018) Used for the THAI Simulations</caption><thead><tr><td>Parameter</td><td>Units</td><td>Value</td></tr></thead><tbody><tr><td>Star and spectrum</td><td></td><td>2600 K BT-Settl with Fe/H = 0</td></tr><tr><td>Semimajor axis</td><td>au</td><td>0.029 28</td></tr><tr><td>Orbital period</td><td>Earth days</td><td>6.1</td></tr><tr><td>Rotation period</td><td>Earth days</td><td>6.1</td></tr><tr><td>Obliquity</td><td></td><td>0</td></tr><tr><td>Eccentricity</td><td></td><td>0</td></tr><tr><td>Instellation</td><td>W m<sup>-2</sup></td><td>900.0</td></tr><tr><td>Planet radius</td><td>km</td><td>5798</td></tr><tr><td>Gravity</td><td>m s<sup>-2</sup></td><td>9.12</td></tr></tbody></table>

cases were designed to act as benchmarks for the two habitable states Hab 1 and Hab 2, respectively, of TRAPPIST-1e, described in part 2 (Sergeev et al. 2022a). Ben 1 is the colder state with a 1-bar N₂-dominated atmosphere plus 400 ppm of CO₂, while Ben 2 is a warmer state with a 1-bar CO₂-dominated atmosphere.

In these simulations, the planet is assumed to be tidally locked to its host star TRAPPIST-1. Both simulations are started from an isothermal (300 K) dry atmosphere at rest. The models are integrated until top-of-atmosphere (TOA) radiative balance is reached. Model outputs are then provided for 10 orbits (i.e., 61 Earth days) and every 6 hr. Readers can refer to Fauchez et al. (2020) for more details on the experimental setup of the Ben 1 and Ben 2 simulations discussed in this manuscript. All Ben 1 and Ben 2 simulations use fixed surface albedo (0.3), according to the THAI protocol (Fauchez et al. 2020). More information on the planet's parameters used in THAI simulations is also provided in Table 1. Note that CO₂ condensation—although it is expected for Ben 2 simulations (Turbet et al. 2018)—has been disabled because not all four GCMs possess this parameterization. Note also that gravity wave parameterizations have been turned off.

### 2.2. Models

As specified in the THAI protocol (Fauchez et al. 2020), four GCMs were used in this study:

1. The Exoplanet Community Atmospheric Model (ExoCAM), a branch of the National Center for Atmospheric Research Community Earth System Model version 1.2.1.

2. The Laboratoire de Météorologie Dynamique—Generic model (LMD-Generic, aka LMD-G).

3. The Resolving Orbital and Climate Keys of Earth and Extraterrestrial Environments with Dynamics (ROCKE-3D Planet_1.0).

4. The Met Office Unified Model (UM, science version GA7.0, code version 11.6).

In the following subsections, we provide a brief description of the model components related to the dynamics and dry physics (see part 2 of the THAI intercomparison for details on the moist physics), with a focus on RT, large-scale dynamics (computed using so-called “dynamical cores”), and subgrid-scale dynamics (dry convection and turbulent mixing). Table 2 gives an overview of the dynamical cores, as well as the horizontal and vertical resolutions used in the four 3D models for all THAI simulations, including the Hab 1 and Hab 2 cases analyzed in Sergeev et al. (2022a). Table 3 summarizes the dry physics parameterizations used by the models.
<table><caption>**Table 2**
Resolution Used in the Models</caption><thead><tr><th>GCM</th><th>Number of Grid Points in the Horizontal</th><th>Number of Vertical Levels (Top Layer)</th></tr></thead><tbody><tr><td>ExoCAM</td><td>72 × 46</td><td>51 (1 Pa)</td></tr><tr><td>LMD-G</td><td>72 × 46</td><td>40 (4 Pa)</td></tr><tr><td>ROCKE-3D</td><td>72 × 46</td><td>40 (10 Pa)</td></tr><tr><td>UM</td><td>144 × 90</td><td>41 and 38 (4 and 13 Pa)</td></tr></tbody></table>

## 2.3. Radiative Transfer

To take into account how the atmosphere, clouds, and surface interact with the light emitted by the star and the planet, GCMs employ RT schemes.

ExoCAM includes a two-stream correlated-k distribution RT scheme with a code lineage that traces back to Urata & Toon (2013a, 2013b) and Colaprete & Toon (2003) with extensive modifications and reorganization occurring over time (Wolf & Toon 2013; Kopparapu et al. 2017; Wolf et al. 2022). This RT scheme, known as ExoRT, has been separated from the 3D model and is provided on GitHub<fnref n="15" /> in order to provide an accessible 1D column model version that facilitates more rapid development and testing. There are now several different RT versions available with ExoRT for use as part of ExoCAM that permit the simulation of different atmospheres. For the THAI simulations, we utilized the oldest and most published RT version, named *n28archean* in the ExoRT repository. This RT version was originally designed and tested for Archean Earth-like conditions, CO₂ amounts up to no more than 30%, Earth-like water vapor column amounts, and Sun-like stellar inputs (Wolf & Toon 2013). Later studies have found this version of the RT to overestimate the greenhouse effect from pure CO₂ atmospheres (Wolf et al. 2022). Note that later iterations of ExoRT have corrected the issue with RT as is discussed in Kopparapu et al. (2017) and Wolf et al. (2022). The *n28archean* version features 28 spectral intervals that extend from 0.2 to 1000 µm, with eight Gauss points per bin. The Gauss point binning used originates from RRTMG (Mlawer et al. 1997) and is weighted toward 1, thus favoring capturing the peak of the absorption lines in each band. While the longwave and shortwave stream bandpasses are customizable to optimize computational expense versus completeness, here we used a brute-force approach and calculated the longwave and shortwave streams over the full spectrum in each case. (Note that, for planets in the TRAPPIST-1 system, about 1% of the stellar energy received lies longward of 9 µm.) Molecular absorption by CO₂ is included, with correlated-k distributions produced using the Atmospheric Environment Research Inc. line-by-line RT model (LBLRTM; Clough et al. 2005) using the HITRAN 2004 spectroscopic database (Rothman et al. 2005) and Voigt line wings cut at 25 cm⁻¹ for both H₂O and CO₂. Water vapor and CO₂ continuum absorption are included using MT_CKD version 2.5, while assuming that CO₂ self-broadened continuum is 1.3 times the foreign broadening continuum (Kasting et al. 1984; Halevy et al. 2009). Recent improvements in the RT module of ExoCAM (Wolf et al. 2022) subsequent to this THAI intercomparison project have shown that this approach overestimates the greenhouse effect

of CO₂ compared to using a sub-Lorentzian CO₂ profile (Perrin & Hartmann 1989; Tran et al. 2011) with a 500 cm⁻¹ cutoff, in agreement with the calculations of Wordsworth et al. (2010). To treat overlapping gas absorption, the amount-weighted scheme of Mlawer et al. (1997) is used. Collision-induced absorption (CIA) is included for N₂-N₂ pairs. Rayleigh scattering is treated using the parameterization of Vardavas & Carver (1984) including N₂ and CO₂.

The LMD-G GCM includes a flexible RT historically originating from the NASA Ames RT scheme, as described in Wordsworth et al. (2011). The RT is performed here for variable gaseous atmospheric compositions made of varying concentrations of CO₂ and N₂, using the correlated-k method. The CO₂ line list was taken from HITRAN 2012 (Rothman et al. 2013) and used to compute high-resolution absorption spectra (for multiple temperatures, pressures, and molecular mixing ratios) based on Voigt profiles (truncation at 500 cm⁻¹, following Wordsworth et al. 2010) adjusted using experimentally based χ-factors (Perrin & Hartmann 1989). These high-resolution spectra were then converted into correlated-k coefficients, using 16 nonregularly spaced grid points for the g-space integration (aka Gauss points), where g is the cumulative distribution function of the absorption data for each band. The correlated-k coefficients were also built using between 32 and 38 spectral bands in the thermal infrared (from 2.3 to 1000 µm) and between 32 and 36 spectral bands in the visible domain (from 0.3 to 5.1 µm), depending on the atmospheric composition considered. Besides, several continuum absorptions (CIA, dimer, far wings) were added for N₂ (N₂-N₂ CIA from Richard et al. 2012) and CO₂ (CO₂-CO₂ CIA, and dimer absorption from Gruszka & Borysow 1997 and Baranov et al. 2004). LMD-G RT uses a two-stream scheme (Toon et al. 1989) to take into account the scattering effects of the atmosphere (Rayleigh scattering) and the clouds (Mie scattering; but note that this is only relevant for Hab 1 and Hab 2 simulations), using the method of Hansen & Travis (1974).

Both the UM and ROCKE-3D share a common RT module SOCRATES,<fnref n="16" /> which is an open-source two-stream correlated-k RT model provided by the Met Office (Edwards & Slingo 1996; Manners et al. 2012). SOCRATES is highly flexible. *Spectral files*, created for both longwave and short-wave streams separately, are tailored to specific atmospheric compositions and stellar spectra. The spectral files are then input to the GCMs at run-time, containing all information necessary to define the RT problem, including the number of spectral intervals; the absorbing gases, continua, and CIAs; the number of Gauss points used for each gas; the designation of dominant and minor absorbing species in each interval; and cloud and Rayleigh scattering properties. Thus, dramatic changes to the RT can be achieved by careful preprocessing of the spectral files, with no changes to the hard-code being required. When run embedded in a GCM, SOCRATES uses the equivalent extinction method to treat multiple overlapping gas species (Amundsen et al. 2017). In the equivalent extinction method the dominant absorbing gas, as determined from the total transmissivity in a test column, is treated with a full k-distribution, while all minor absorbing gases are treated as gray. The number of Gauss points used for the dominant gas in each spectral interval is determined based on a selected transmission error tolerance (usually 0.001). For the Ben 1 simulations, spectral files were tailored to Earth-like

<footnote n="15">¹⁵ https://github.com/storyofthewolf/ExoRT</footnote>

<footnote n="16">¹⁶ https://code.metoffice.gov.uk/trac/socrates</footnote>
<table><caption>**Table 3**
Dry Physics Parameterizations in the THAI GCMs</caption><thead><tr><th>GCM</th><th>Radiative Transfer</th><th>Dynamical Core</th><th>Subgrid Dynamics</th><th>Land Module</th></tr></thead><tbody><tr><td>ExoCAM</td><td>Correlated-k<br>RT time step of 1 hr<br>0.2–1000 µm<br>17 × 23 bands (SW × LW)<br>HITRAN 2004<br>CO<sub>2</sub> self-broadening = 1.3 × that of N<sub>2</sub><br>MT_ CKD v2.5 CO<sub>2</sub> continuum; N<sub>2</sub>–N<sub>2</sub> CIA</td><td>Finite volume<br>Dynamical time step of 56 s<br>Arakawa C grid<br>Fourth-order divergence damping<br>With Laplacian sponge<br>Second-order velocity component<br>Damping for top layers</td><td>First-order closure<br>With explicit nonlocal flux of heat<br>Roughness = 10<sup>-2</sup> m</td><td>15 ground layers<br>total depth = 35 m<br>Heat capacity = 2.1 × 10<sup>6</sup> J m<sup>-3</sup> K<sup>-1</sup><br>Thermal inertia = 4 × 10<sup>3</sup> J m<sup>-2</sup> K<sup>-1</sup> s<sup>-1/2</sup></td></tr><tr><td>LMD-G</td><td>Correlated-k<br>RT time step of 15 minutes/7.5 minutes (Ben1/Ben2)<br>0.3–1000 µm<br>36 × 38 bands (SW × LW) for Ben 1<br>36 × 32 bands (SW × LW) for Ben 2<br>HITRAN 2012<br>CO<sub>2</sub> self-broadening = that of N<sub>2</sub><br>CO<sub>2</sub> sub-Lorentzian far wings<br>CO<sub>2</sub>–CO<sub>2</sub> dimer + CIA; N<sub>2</sub>–N<sub>2</sub>CIA</td><td>Finite difference<br>Dynamical time step of 90 s/45 s<br>(for Ben1/Ben2)<br>Arakawa C grid<br>Scale-selective horizontal dissipation<br>Based on an n time iterated Laplacian</td><td>1.5-order (TKE) closure<br>With dry convective adjustment<br>Roughness = 10<sup>-2</sup> m</td><td>18 ground layers<br>Total depth = 20 m<br>Heat capacity = 2 × 10<sup>6</sup> J m<sup>-3</sup> K<sup>-1</sup><br>Thermal inertia = 2 × 10<sup>3</sup> J m<sup>-2</sup> K<sup>-1</sup> s<sup>-1/2</sup></td></tr><tr><td>ROCKE-3D</td><td>Correlated-k<br>RT time step of 1 hr<br>0.2–10,000 µm<br>21 × 12 bands (SW × LW) for Ben 1<br>42 × 17 bands (SW × LW) for Ben 2<br>HITRAN 2012<br>CO<sub>2</sub> self-broadening = that of N<sub>2</sub><br>CO<sub>2</sub> sub-Lorentzian far wings<br>CO<sub>2</sub>–CO<sub>2</sub> CIA; no N<sub>2</sub>–N<sub>2</sub> CIA</td><td>Finite difference<br>Dynamical time step of 225 s/450 s<br>(for Ben1/Ben2)<br>Arakawa B grid<br>Numerical sponge in the top six layers</td><td>1.5-order (TKE) closure<br>With explicit nonlocal flux of heat<br>Roughness = 4.1 × 10<sup>-2</sup> m</td><td>Six ground layers<br>Total depth = 3.5 m<br>Heat capacity = 0.85 × 10<sup>6</sup> J m<sup>-3</sup> K<sup>-1</sup><br>Thermal inertia = 1.2 × 10<sup>3</sup> J m<sup>-2</sup> K<sup>-1</sup> s<sup>-1/2</sup></td></tr><tr><td>UM</td><td>Correlated-k<br>RT time step of 1h<br>0.2–10,000 µm<br>21 × 12 bands (SW × LW) for Ben 1<br>42 × 17 bands (SW × LW) for Ben 2<br>HITRAN 2012<br>CO<sub>2</sub> self-broadening = that of N<sub>2</sub><br>CO<sub>2</sub> sub-Lorentzian far wings<br>CO<sub>2</sub>–CO<sub>2</sub> CIA; no N<sub>2</sub>–N<sub>2</sub> CIA</td><td>Finite difference<br>Dynamical time step of 1200 s/600 s<br>(for Ben1/Ben2)<br>Arakawa C grid<br>Numerical damping of vertical wind<br>With altitude and latitude</td><td>First-order closure<br>with explicit nonlocal<br>Fluxes of heat and momentum<br>Roughness = 10<sup>-2</sup> m</td><td>1 ground layer<br>total depth = 1 m<br>Heat capacity = 2 × 10<sup>6</sup> J m<sup>-3</sup> K<sup>-1</sup><br>Thermal inertia is implicitly infinite</td></tr></tbody></table>
atmospheric compositions, allowing for small amounts of CO₂ in N₂-dominated atmospheres. The spectral files used by ROCKE-3D and the UM are identical (details provided in Table 3). Gas absorption coefficients are derived from the HITRAN 2012 spectroscopic database (Rothman et al. 2013). For Ben 2 simulations, CO₂ absorption is treated, following Wordsworth et al. (2010), with sub-Lorentzian line shape (Perrin & Hartmann 1989) out to 500 cm⁻¹ with CO₂–CO₂ CIA overlaid.

## 2.4. Dynamical Cores (and Numerical Damping)

A key part of a GCM is its dynamical core, which solves a system of Navier–Stokes equations on a spherical grid to simulate movement of mass and energy in the atmosphere. To solve these equations, three THAI GCMs (LMD-G, ROCKE-3D, the UM) use a finite-difference approach, while ExoCAM is based on a finite-volume approach; there are also important differences in the implementation and approximations used.

ExoCAM uses a finite-volume scheme (Lin & Rood 1996). Modifications have been brought to its dynamical core to improve numerical stability in more strongly forced atmospheres. They consist in incrementally applying physics tendencies for temperature and wind speed evenly throughout the dynamical time step, rather than only at the beginning of it (Bardeen et al. 2017). The number of points in the latitudinal and longitudinal directions is 72 and 46, respectively. This gives the horizontal grid spacing of 5° × 3°.75. ExoCAM uses 51 levels in the vertical, stretching from the surface to 1 hPa in the atmosphere. Regarding numerical damping, a fourth-order divergence damping with Laplacian sponge is applied. A second-order velocity-component damping is also used for top layers. There is therefore an implicit numerical sponge layer produced by the degradation of the order of the numerical scheme.

The LMD-G's GCM dynamical core solves the primitive equations of meteorology using an entrophy and total angular momentum conserving finite-difference dynamical core on an Arakawa C grid (for more details, see Forget et al. 1999; Hourdin et al. 2006). A filter is applied at high latitudes to satisfy the Courant–Friedrichs–Lewy numerical stability criterion. For the THAI simulations, LMD-G uses a uniform resolution of 5° in longitude and 3°.75 in latitude. The dynamical time step is 45 s. The model uses hybrid vertical coordinates (terrain-following coordinate system in the lower atmosphere; fixed pressure levels in the upper atmosphere). For the THAI simulations, we used 40 vertical layers, with the lowest midlayer level at 5 m and the top midlayer level at ≈4 Pa. Nonlinear interactions between explicitly resolved scales and subgrid-scale processes (which are necessary to ensure numerical stability) are parameterized by applying a scale-selective horizontal dissipation operator based on an *n* time iterated Laplacian as described in Forget et al. (1999), Section 3.

The ROCKE-3D dynamical core utilizes finite differences where atmospheric velocity points are located on the Arakawa B grid (Arakawa 1972; Hansen et al. 1983). The atmospheric vertical layers use a sigma coordinate system from the surface to 150 hPa. From there constant-pressure layers are used to the TOA (≈0.1 hPa). For the THAI experiments the model was run with 72 × 46 grid points, with each grid point being 5° in longitude by 4° in latitude except at the poles, where the latitude coordinate is only 2° in extent. The smaller polar grid

boxes are necessary to maintain second-order accuracy of air-mass fluxes, the pressure gradient force, and momentum advection terms (Schmidt et al. 2006). The dynamical time steps used in the Ben 1 and Ben 2 simulations were 225 and 450 s, respectively. Both Hab 1 and Hab 2 simulations use a dynamical time step of 450 s. For all four simulations the parameterized physics uses a time step of ≈30 minutes (1756.8 s to be exact), and less in some submodules, while the SOCRATES radiation was called every two physics time steps (~1 hr).

The UM employs the ENDGAME dynamical core with a semi-implicit semi-Lagrangian formulation to solve the non-hydrostatic fully compressible deep-atmosphere equations of motion (Wood et al. 2014). Model fields are discretized horizontally onto a regular longitude–latitude Arakawa C grid, while in the vertical a terrain-following hybrid height coordinate with a Charney–Phillips staggering is used. THAI experiments are run with 144 × 90 grid points in the horizontal, giving a grid spacing of 2°.5 in longitude and 2° in latitude. In the vertical, 41 levels between the surface and a fixed model lid, located at ≈63 km height, are used for the Ben 1 experiment. For Ben 2, the number of levels is 38 and the model lid is at ≈40 km to improve the model stability. The primary dry prognostic variables in the UM are the three wind components, virtual dry potential temperature, Exner pressure,<fnref n="17" /> and density. Increments to these variables are obtained on each dynamical time step (1200 s in Ben 1, 600 s in Ben 2), within which processes are split into an outer loop. The semi-Lagrangian departure points are obtained within the outer loop using the latest estimates for the wind components. Fields are then interpolated to the updated departure points. Within the inner loop, a linear Helmholtz problem is solved to calculate the pressure increment. Estimates for all variables are then obtained from the pressure increment via a back-substitution process (for a schematic of the time-stepping algorithm and more details, see Walters et al. 2017). There is no explicit diffusion or dissipation in the model. At the TOA, a sponge layer is used to suppress numerical instabilities by explicitly damping vertical velocity using a sin² function. The damping is applied in the upper half of the model domain above at the equator, dropping in altitude with latitude, eventually to the surface at the poles.

## 2.5. Subgrid-scale Mixing

To take into account how the atmosphere acts on scales too small to be resolved by the dynamical cores, GCMs employ subgrid-scale parameterizations. For a dry simulation, they include parameterizations of turbulence and dry convection, which are usually handled by the boundary-layer or convection schemes, or a combination of both.

The ExoCAM parameterization package (CAM4; Neale et al. 2010) can be described by a sequence of four components: the moist precipitation processes, the cloud and radiation, the surface model, and the turbulent mixing. In the Ben 1 and Ben 2 cases, only the shallow part of convection is active. This is done by ensuring that the lapse rate is stable and stratified. The surface model includes surface fluxes obtained from land, ocean, and sea ice models, or computed from

<footnote n="17">¹⁷ Exner pressure is defined as $\Pi = (p/p_0)^{R_d/c_p}$, where *p* is the air pressure, *p₀* is the reference pressure, *R_d* is the specific gas constant, and *c_p* is the specific heat capacity at constant pressure.</footnote>
specific conditions. These surface fluxes provide lower flux boundary conditions for the turbulent mixing, which is composed of the planetary boundary-layer (PBL) parameterization and vertical diffusion.

In LMD-G, subgrid-scale dynamical processes including turbulent mixing and convection are parameterized as in Forget et al. (1999). In practice, the boundary-layer dynamics are accounted for by the Mellor & Yamada (1982) unstationary 2.5-level closure scheme (using the implementation of Galperin et al. 1988) plus a convective adjustment that rapidly mixes the atmosphere in the case of unstable temperature profiles. In that case, the thermal profile is brought back to the dry or wet adiabat (in the case where water vapor condenses, which is only relevant for Hab 1 and Hab 2 simulations described in THAI part 2). In the Ben 1 and Ben 2 (dry) simulations, only the former case applies. Turbulence and convection mix energy (potential temperature), momentum (wind), and tracers (gases and aerosols) as described in Forget et al. (1999), Section 6.1. For the THAI simulations, a standard roughness coefficient $z_0 = 1 \times 10^{-2}$ m is used for both dry (Ben 1 and Ben 2) and ocean (Hab 1 and Hab 2; see THAI part 2) surfaces for simplicity.

In ROCKE-3D, level 2 of the Cheng et al. (2002) model is used in the free troposphere (above the PBL—which is determined dynamically). This scheme has the ability to produce fluxes when turbulence is weak, at large Richardson numbers (large buoyancy to shear forcing ratios). From the midpoint of the first atmospheric model layer to the surface there are eight sublayers used to compute surface turbulent fluxes with the level 2.5 model from Cheng et al. (2002). In between the two it depends on the stability of the PBL. In the stable case, we use level 2 of the Cheng et al. (2002) model, which increases the critical Richardson number (from the Mellor & Yamada 1982 model value 0.2–1), in agreement with several data sets, and with the potential to increase the PBL height (the latter is Richardson number related). In the unstable case, for diffusivities, we use the nonlocal model of Holtslag & Boville (1993) with countergradient fluxes, and for the turbulent kinetic energy (TKE) we use the parameterization based on the large eddy simulation (LES) by Moeng & Sullivan (1994). For details see Schmidt et al. (2006) and references therein.

To represent mixing of adiabatically conserved heat, momentum, and tracers throughout the troposphere, including the PBL, the UM uses a first-order turbulence closure of Lock et al. (2000) with modifications described in Lock (2001) and Brown et al. (2008). The scheme treats unstable and stable boundary layers differently. For unstable layers, there are two separate profiles of diffusion coefficients, one for the surface sources of turbulence and one for upper atmospheric sources (e.g., radiative cooling). Initially diagnosed by an adiabatic ascent, the extent of these profiles is adjusted so that the magnitude of the buoyancy consumption of turbulent kinetic energy is limited to a specified fraction of buoyancy production. In moist aquaplanet simulations (see part 2), this permits the cloud layer to decouple from the surface. Additional nonlocal fluxes of heat and momentum can generate more vertically uniform potential temperature and wind profiles in unstable (convective) boundary layers. The scheme is tightly coupled to the convection scheme (see part 2), but the latter is inactive in dry cloudless simulations of the Ben 1 and Ben 2 cases. For stable boundary layers and the free troposphere,

diffusion coefficients are calculated using the parameterization of Smith (1990), based on a local Richardson number. The mixing can also be enhanced if the nonlocal coefficients, calculated as described above, are larger than the local ones, which can happen in the neutral boundary layer. The stability dependence is given by the “MES-tail” function within the stable boundary layer and by the “sharp” function at 200 m and above (Brown et al. 2008). Given the resulting coefficient, the diffusion equation is solved implicitly, and the kinetic energy dissipated through the turbulent shear stresses is converted to a local heating term.

### 3. Results

We present and compare results from the four GCMs on board the THAI intercomparison, for the dry Ben 1 and Ben 2 simulations of TRAPPIST-1e. We first discuss radiation fluxes, then convection and turbulence, and eventually large-scale dynamics.

#### 3.1. Radiative Budget and Fluxes

In order to compare the radiative budgets of the GCMs, we first start with the shortwave fluxes (which should not be very sensitive to the thermal structure of the atmosphere in the Ben 1 and Ben 2 dry cases) and then compare the thermal, longwave fluxes. Figure 1 shows the absorbed stellar radiation (ASR) horizontal map for the Ben 1 (left) and Ben 2 (right) cases. ASR is defined here as the total flux from the star (TRAPPIST-1) that is absorbed by the planetary surface and atmosphere. In both the Ben 1 and Ben 2 cases, the mean and maximum values agree between the models within a few $W m^{-1}$. The minimum values are equal to $0 w m^{-2}$ for all GCMs, as the nightside is not illuminated. The highest ASR values are for ROCKE-3D, followed by the UM, ExoCAM, and LMD-G GCMs.

Since the surface albedo is fixed at 0.3 in the THAI Ben 1 and Ben 2 simulations (as defined in the THAI protocol of Fauchez et al. 2020) and the stellar spectral flux is the same for the four GCMs, differences in absorbed flux may arise from Rayleigh scattering but more likely from molecular absorption (from $CO_2$ and marginally from $N_2-N_2$ CIA, in the Ben 1 and Ben 2 simulations) because the spectrum of the star TRAPPIST-1 is very strongly shifted toward infrared wavelengths (see, e.g., Turbet & Selsis 2021, Figure 1), where $CO_2$ can absorb (and where Rayleigh scattering has a very limited role). For the Ben 1 simulations, the mean Bond albedo is equal to 0.27/0.28/0.29/0.3 for ROCKE-3D/UM/ExoCAM/LMD-G, respectively (see Table 4). These values are very close to the surface albedo, indicating little influence from Rayleigh scattering and $CO_2$ absorption (only ≈400 ppm of $CO_2$ in the Ben 1 simulations). For the Ben 2 simulations, the mean Bond albedo is equal to 0.18/0.19/0.21/0.22 for ROCKE-3D/UM/ExoCAM/LMD-G, respectively (see Table 4). These values are significantly lower than that of the Ben 1 simulations because here the atmosphere is essentially composed of 1 bar of $CO_2$ (about 2500× the amount in Ben 1) that can efficiently absorb the incoming stellar radiation in the near-infrared.

The differences between the outputs obtained with our set of GCMs are very small and comparatively similar between the Ben 1 and Ben 2 cases. Similar results are also observed in the shortwave heating rate vertical profiles (see Figures 2(c) and 3(c)), which display similar trends from the surface up to ≈1 hPa. The Ben 2 simulations result in a much stronger
## THAI Ben1

<figure><img src="image_4.png" /><figcaption>(a) Results for the Ben 1 simulations</figcaption></figure>

## THAI Ben2

<figure><img src="image_5.png" /><figcaption>(b) Results for the Ben 2 simulations</figcaption></figure>

**Figure 1.** ASR maps for the Ben 1 and Ben 2 THAI simulations, for the four GCMs. The temporal average minimum, mean, and maximum values are also shown below each map.
<caption>Table 4</caption>

Global Mean Surface Temperature ($T_s$, K), Global Mean Top-of-atmosphere Outgoing Longwave Radiation ($F_{IR}$, W m⁻²), and Planetary Albedo ($\alpha_p$) in Ben 1 and Ben 2 Simulations

<table><thead><tr><th rowspan="2">GCM</th><th colspan="3">Ben 1</th></tr><tr><th>T<sub>s</sub></th><th>F<sub>IR</sub></th><th>α<sub>p</sub></th></tr></thead><tbody><tr><td>ExoCAM</td><td>220</td><td>160</td><td>0.29</td></tr><tr><td>LMD-G</td><td>217</td><td>160</td><td>0.30</td></tr><tr><td>ROCKE-3D</td><td>216</td><td>164</td><td>0.27</td></tr><tr><td>UM</td><td>213</td><td>162</td><td>0.28</td></tr></tbody><tbody><tr><th rowspan="4">ExoCAM</th><th colspan="3">Ben 2</th></tr><tr><td>245</td><td>179</td><td>0.21</td></tr><tr><td>242</td><td>174</td><td>0.22</td></tr><tr><td>242</td><td>184</td><td>0.18</td></tr><tr><td>UM</td><td>239</td><td>182</td><td>0.19</td></tr></tbody></table>

shortwave heating rate due to direct absorption of stellar light by CO₂. Both the Ben 1 and Ben 2 cases exhibit shortwave heating rates that increase as altitude increases (or pressure decreases). This is because CO₂ absorption saturates as we go deeper into the atmosphere (because the upper layers absorb all photons in the core of near-infrared CO₂ band lines).

Most of the differences in the RT likely come from differences in the parameterization of CO₂ absorption (CO₂ is the main absorbing species in Ben 1 and Ben 2 simulations). This is supported in particular by the fact that the two models—ROCKE-3D and UM—that share the same RT scheme, namely, SOCRATES (as well as the same spectral files), display very similar shortwave heating rates from the surface up to ≈1 hPa (see orange and green lines in Figures 2(c) and 3(c)). The two models diverge mainly for the two uppermost levels of UM (especially the last altitude level), presumably due to different RT boundary conditions between UM and ROCKE-3D. At the pressures where the differences between shortwave heating rates are the most significant between models, CO₂ absorption is in principle dominated by the cores of CO₂ band lines (and not the far wings or CIA/dimer absorptions caused by molecular collisions). Note that the shortwave heating rate variability within individual models (identified by shades in Figures 2(c) and 3(c)) is small, indicating that the near-infrared CO₂ absorption weakly depends on atmospheric temperature. A more detailed comparison of the correlated-k tables used by the models is required to identify the origins of the variations in the high-altitude shortwave heating rates.

Figure 4 is similar to Figure 1, but for the TOA outgoing longwave radiation (OLR). TOA OLR is the total amount of infrared radiation that is emitted at the top of the model. For the Ben 1 case, the intermodel agreement is quite good, with LMD-G and ExoCAM showing exactly the same mean value of 160 W m⁻², while the UM and ROCKE-3D have higher values by only 2 and 4 W m⁻², respectively (see Table 4). The agreement is also good for the Ben 2 case (see Table 4), but the difference slightly increases up to 10 W m⁻² between the mean values of LMD-G and ROCKE-3D, and up to 35 W m⁻² between the maximum values of LMD-G and the UM (slightly eastward of the substellar point). The global mean values of OLR match those for ASR, confirming that the models are close to the radiative balance at the TOA. The OLR maps are controlled by the temperature of the surface (Figure 5) and the atmosphere, both of which are

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Vertical profiles at the substellar point for the Ben 1 simulations of (a) the temperature, (b) the net longwave heating rate (i.e., the heating rate of each layer resulting from a balance between thermal infrared emission of the layer and the absorption of thermal infrared emission of the other layers and the surface), (c) the shortwave heating rate (i.e., the heating rate of each layer produced by the absorption of the stellar flux), and (d) the lapse rate (i.e., the vertical temperature gradient). Time-averaged values are represented by thick lines, while the 1σ deviations (over the 10 orbits) are represented by shades.</figcaption></figure>

heated and cooled by the global heat redistribution. Broadly, the OLR peak is located at or a few degrees to the east of the substellar point, which is where the surface absorbs most of the instellation and thus is the warmest. Weakest OLR regions are also colocated with coldest surface temperature regions.

At the substellar point, the tropospheric (from the surface up to 100 hPa) temperature structure is very similar across all models (in both Ben 1 and Ben 2 simulations), as seen in Figures 2(a) and 3(a), respectively. Under these conditions, the substellar longwave heating rate profiles in the troposphere (see Figures 2(b) and 3(b)) agree well with each other, in both the Ben 1 and Ben 2 simulations. The longwave heating rates are net quantities: in each atmospheric layer they result from its thermal emission, as well as the absorption of thermal radiation of the other layers and the surface. The longwave heating is relatively small near the surface for all models and for both the Ben 1 and Ben 2 simulations, while in the upper atmosphere strong longwave cooling dominates owing to the effective infrared emission by CO₂, particularly prominent for the Ben 2 simulations (Figure 3(b)). In the upper atmosphere (i.e., the radiative part), net longwave heating rate profiles are strongly
<figure><img src="image_7.png" /><figcaption>Figure 3. Same as Figure 2, but for the Ben 2 simulations.</figcaption></figure>

anticorrelated with temperature profiles (in other words, longwave cooling is strongly correlated with temperature) because the temperature controls the emission of the layers.

In summary, the RT calculations, shown here for the substellar point, across all GCMs produce very similar results (with variations of radiative heating rates between the four models of always less than 5%, from the surface up to ≈5 hPa), except in the uppermost part of the atmosphere. Vertical variations of the net longwave heating rates are strongly linked to the temperature structure. The longwave heating varies horizontally across the planet too, which is a function of the temperature structure discussed in the next sections.

## 3.2. Turbulence and Convection

In order to compare the behavior of our GCMs for the convective boundary layer, we focus on the substellar point, where direct stellar heating of the surface drives intense convection (evident, e.g., in the vertical wind patterns in Figures 9 and 10).

As was mentioned in the previous section, the GCMs agree well on the temperature profile at the substellar point. Subtle differences are revealed by looking at the lapse rate (i.e., the vertical temperature gradient), shown in Figures 2(d) and 3(d). In the troposphere, the lapse rates are about −9 and −11 K km$^{-1}$ on average for the Ben 1 and Ben 2 cases, respectively. To first order, the difference between the Ben 1

and Ben 2 simulations comes from the variation of the atmospheric heat capacity between the cases (about 1040 and 850 J kg$^{-1}$ K$^{-1}$ for Ben 1 and Ben 2, respectively). Note that because there are no changes in the atmospheric composition in our simulations, all four GCMs have the heat capacity set to a fixed parameter in each of the cases.

ExoCAM, ROCKE-3D, and the UM show superadiabatic behavior near the surface, while LMD-G—whose temperature profile is driven by a dry convective adjustment scheme— follows closely the $-g/c_p$ lapse rate. Dry convection is triggered by the instability induced by the presence of a superadiabatic layer near the surface. In the case of LMD-G, the convective adjustment instantly removes the superadiabatic layer, forcing it to follow the $-g/c_p$ lapse rate. In the other three models, the superadiabatic layer persists owing to the (nonzero) timescale needed for the hot near-surface air parcels to rise.

The models all have a troposphere that reaches about 100 hPa, for both Ben 1 and Ben 2 simulations, in line with the theoretical prediction of Robinson & Catling (2014). However, there are notable differences between models, in particular regarding the temperature of the tropopause. The intermodel variability is about 15 K for both Ben 1 and Ben 2 simulations (see Figures 2(a) and 3(a)). LMD-G consistently produces the coldest tropopause, as in the Hab 1 and Hab 2 aquaplanet simulations presented in THAI part 2 (Sergeev et al. 2022a). This is most likely due to the fact that LMD-G also has the weakest stratospheric shortwave heating rates in both Ben 1 and Ben 2 simulations (see Figures 2(c) and 3(c)).

Despite these differences, we note that the convective layer of the atmosphere is the part that produces the smallest temperature variations among the models (as well as the smallest internal variability in the models, represented by the shades in Figures 2(a) and 3(a)). This is a good marker that stellar energy input is consistent across the four GCMs.

## 3.3. Large-scale Dynamics

To understand and interpret the horizontal variations in OLR and temperature between the GCMs, and more generally the climate and circulation regime as a whole, it is necessary to compare how the models represent the large-scale dynamics.

Figure 5 shows the surface temperature maps and the winds at the 250 hPa isobaric surface for the Ben 1 (panel (a)) and Ben 2 (panel (b)) cases. First, we can see that the intermodel spread in the global mean surface temperature is 7 K for the Ben 1 simulations and 6 K for the Ben 2 simulations (see also Table 4). In both cases ExoCAM is on average the warmest and the UM the coldest, similar to what is found for the Hab 1 case (Sergeev et al. 2022a). For the Ben 1 case the day-night thermal contrast is the largest for ROCKE-3D, followed by the UM, LMD-G, and ExoCAM. In the case of ROCKE-3D, two wind gyres form in the western hemisphere at high latitudes, in the coldest regions (see Figure 5(b), bottom left plot). Those gyres are a manifestation of a Rhines rotator dynamical regime, characteristic for slow-rotating tidally locked exoplanet simulations (e.g., Carone et al. 2015; Noda et al. 2017; Haqq-Misra et al. 2018; Hammond et al. 2020; Wang & Yang 2021); this is a regime where a strong upper atmosphere superrotation is combined with a strong upwelling beneath the substellar point. The other three GCM simulations produce a fast-rotator circulation regime: without the gyres, but with two extra-tropical zonal jets. This is best illustrated by the lower part of
## THAI Ben1

<figure><img src="image_8.png" /><figcaption>(a) Results for the Ben 1 simulations</figcaption></figure>

## THAI Ben2

<figure><img src="image_9.png" /><figcaption>(b) Results for the Ben 2 simulations</figcaption></figure>

**Figure 4.** Same as Figure 1, but for the outgoing longwave radiation (OLR).
THAI Ben1

<figure><img src="image_10.png" /><figcaption>(a) Results for the Ben 1 simulations</figcaption></figure>

THAI Ben2

<figure><img src="image_11.png" /><figcaption>(b) Results for the Ben 2 simulations</figcaption></figure>

**Figure 5.** Surface temperature (shading, K) maps on which we superimposed horizontal wind vectors in the upper troposphere (≈250 hPa).
<figure><img src="image_12.png" /><figcaption>**Figure 6.** Vertical cross sections of the zonal mean zonal wind (m s⁻²). The gray dashed horizontal line marks the 250 hPa level used in Figures 9 and 10.</figcaption></figure>

Figure 6, which shows the vertical cross sections of the zonal mean zonal wind. Similar to the findings of Sergeev et al. (2020) for TRAPPIST-1e, the thermal contrast between the dayside and the nightside is smaller in this regime than in the slow-rotator regime. This is important not only for TRAPPIST-1e but generally for any tidally locked rocky exoplanet with short enough orbital periods, i.e., short enough for the Rossby deformation radius to be comparable to the planet's radius (Carone et al. 2015, 2016; Haqq-Misra et al. 2018; Penn & Vallis 2018; Pierrehumbert & Hammond 2019; Zhang 2020).

In the Ben 2 case, the circulation regime for two of the models (LMD-G and the UM) changes to the Rhines rotator (Figure 5(b)). The day-night contrast stays the largest in ROCKE-3D's simulation, followed then by that in ExoCAM. Despite the circulation regime change, LMD-G and the UM still have a lower day-night temperature contrast than that in ExoCAM. Note, however, that the magnitude of this contrast is overall smaller than in the Ben 1 case owing to overall more efficient heat redistribution (see Figure 7).

Ben 1 and Ben 2 simulations both exhibit small temporal variability, as illustrated for surface temperatures in Figure 8, due most likely to the manifestation of weather patterns. It has indeed already been shown that on synchronous planets (with a fixed solar forcing) a temporal variability is nevertheless expected (Joshi et al. 1997; Merlis & Schneider 2010; Ding & Wordsworth 2021). The four GCMs show a similar small level of global mean surface temperature high-frequency variability (~1 K) for Ben 1 and Ben 2 simulations. However, on a period of 10 consecutive orbits, the frequency of temperature oscillations varies from one model to another. The ExoCAM simulations show much greater nightside (Figure 8, right) and dayside (Figure 8, middle) variability than the other GCMs, but these opposite-phase oscillations offset each other on the global mean (Figure 8, left). This behavior, which is also pronounced in the ExoCAM Hab 1 and Hab 2 simulations (part 2; Sergeev et al. 2022a), results from periodic weather systems traveling

around the planet. Unfortunately, the design of the Ben 1 and Ben 2 cases of the THAI protocol (Fauchez et al. 2020)—which limited the output files to 10 orbits—does not allow us to probe the low-frequency temporal variability, but there are several indications that it exists, as suggested by rolling averages in Figure 8. This is also indicated by the presence of north-south asymmetries in the different mean values of the models (see, e.g., Figures 5 and 6), somewhat similar to the results of Noda et al. (2017), Section 4.3. There do not appear to be clear correlations between the different model variabilities and the land properties adopted in the models (especially the number and depth of layers; see Table 3, rightmost column) over the time range explored from the THAI simulations. We leave the question of why weather systems seem to be more pronounced on ExoCAM than other GCMs for future studies. We also recall that the level of temporal variability of the different models is far below the detectability threshold of JWST and the astronomical observatories of the coming decades (part 3; Fauchez et al. 2022).

The fact that the models reside in different rotational regimes is likely because the radius and rotation period of TRAPPIST-1e place it just at the edge between fast and Rhines (intermediate) rotation regimes (Carone et al. 2015, 2016; Haqq-Misra et al. 2018), as already identified in Sergeev et al. (2020) and Carone et al. (2018) specifically for TRAPPIST-1e. As a matter of fact, Sergeev et al. (2020) showed that just a change in the parameterization of convection can drive a GCM into one or the other dynamical regime. Unlike the simulations in Edson et al. (2011), initial conditions are unlikely to play a role here, as all simulations start from isothermal conditions of 300 K and null winds. We also identified in preliminary THAI simulations (of the UM model) that subtle changes (here, of correlated-k tables) are able to settle the GCM to a different rotation regime. We also note that a similar behavior appears in the aquaplanet THAI simulations, with ROCKE-3D being an outlier and settling in a fast-rotator regime (Sergeev et al. 2022a). Changes
<figure><img src="image_13.png" /><figcaption>**Figure 7.** Overturning circulation shown by the tidally locked mass stream function $\Psi^{\text{TL}}$ (10$^{11}$ kg s$^{-1}$). The mass flux is clockwise where $\Psi^{\text{TL}}$ is positive and counterclockwise where $\Psi^{\text{TL}}$ is negative. The maximum of $\Psi^{\text{TL}}$ in each GCM is shown in the upper left corner. The substellar point (SS) is at 90° latitude, and the antistellar point (AS) is at −90° latitude.</figcaption></figure>

in surface friction can also drive GCMs into one or the other climate regime (Carone et al. 2016), but we caution that the surface frictions used in Carone et al. (2016) may be unrealistically high, thus potentially overestimating the importance of this effect (Wicker et al. 2022). More generally, this shows that the tipping point is highly nonlinear and that small parameter changes can lead to a different climate regime.

Similar regime changes are seen between N$_2$- and CO$_2$-rich atmospheres in the cases including moisture (Sergeev et al. 2022a), and it is apparent that the simulations are sensitive to several parameters for this regime. Model parameterizations affect the eddy generation in our GCMs in different ways, which may then tip the balance one way or another. Further work has been performed to explore the potential bistability of this climate state (Sergeev et al. 2022b).

To get a more detailed view of these dynamical regimes, we used the Helmholtz decomposition methodology of Hammond & Lewis (2021). The total wind flow at 250 hPa is decomposed into its zonal mean rotational, eddy rotational, and divergent flows, as can be seen in Figures 9 and 10 for the Ben 1 and Ben 2 simulations, respectively. The zonal means of the rotational wind clearly show the prograde jets induced by the fast rotation of TRAPPIST-1e (Figures 9(e)–(h) for the Ben 1 case and Figures 10(e)–(h) for the Ben 2 case). The eddy component (as the deviation from the zonal mean) reveals planetary-scale stationary waves straddling the substellar longitude (Figures 9(i)–(l) for Ben 1 and Figures 10(i)–(l) for Ben 2). For the Ben 1 simulations, only ROCKE-3D displays stationary gyres at midlatitudes. As ExoCAM, LMD-G, and the UM do not show those gyres, this is likely a sign of a faster atmospheric superrotation placing them in the fast-rotator regime. For the Ben 2 simulations, LMD-G, ROCKE-3D, and the UM display

stationary gyres at midlatitudes. Only ExoCAM produces the fast-rotating regime, with no midlatitude gyres. The asymmetry visible in the UM simulation of the Ben 2 case is likely an artifact of a relatively short simulation length.

These differences in circulation regime are also visible in the tidally locked mass stream function (Figure 7), which describes the overturning circulation as a function of the substellar latitude (between 90° at the substellar and −90° at the anti-substellar points). In all eight simulations (Ben 1 and Ben 2, and for the four GCMs), the mass stream function $\Psi^{\text{TL}}$ shows a single circulation cell between the dayside and the nightside (with air parcels rising at the substellar point and subsiding on the nightside). There is some variability in the strength of the circulation across models, but the circulation is consistently stronger in the Ben 2 simulations than in the Ben 1 simulations. This is likely due to the fact that the intensity of convection is stronger at the substellar point in Ben 2 (due mostly to the fact that longwave CO$_2$ absorption by the dayside atmosphere leads to surface warming that further increases convection), which is further strengthened by the strong radiative nightside cooling (due to stronger longwave emission by CO$_2$). This reflects directly on Figure 4, which shows that the OLR is weaker on the dayside and stronger on the nightside in the Ben2 simulations than in the Ben1 simulations. We observe that the circulation cells extend to higher substellar latitudes for simulations trapped in the Rhines rotation regime (ROCKE-3D for Ben 1; LMD-G, ROCKE-3D, and the UM for Ben 2) than those trapped in the fast rotation regime, as previously identified in Haqq-Misra et al. (2018).

Last but not least, we identify that the main source of differences between the models appears in the upper part of the atmosphere. In Figures 2 and 3, temperature differences between models peak in the upper layers of the atmosphere,
<figure><img src="image_14.png" /><figcaption>**Figure 8.** Temporal variability of the surface temperatures (K) for the Ben 1 (top row) and Ben 2 (bottom row) simulations. Each column shows global mean, dayside mean, and nightside mean. Thin lines show the raw data; thick lines show a 2.5-orbit rolling mean.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 9.** Helmholtz decomposition of the horizontal wind at 250 hPa in the Ben 1 cases (quivers): (a)–(d) total wind, (e)–(h) zonal mean rotational component, (i–l) eddy rotational component, (m)–(p) divergent component. Note the different scaling of the eddy rotational and divergent components. Also shown is the upward wind velocity (shading, m s⁻¹), with the 0.05 m s⁻¹ highlighted by a black contour. Note that for ExoCAM only the pressure velocity (ω, Pa s⁻¹) is available in the output, so the vertical velocity (w, m s⁻¹) is approximated as w = −ω/ρg, where ρ is air density and g is the acceleration due to gravity.</figcaption></figure>

reaching about 10–20 K in both the Ben 1 and Ben 2 simulations. This is also where variability intrinsic to the models (visible in shading in Figures 2 and 3) is the largest.

In Figure 6, wind distribution is also very different from one model to another above the troposphere. Simulations in the fast-rotator regime (all but ROCKE-3D in the Ben 1 case;
<figure><img src="image_16.png" /><figcaption>**Figure 10.** Same as Figure 10, but for the Ben 2 simulations.</figcaption></figure>

ExoCAM in the Ben 2 case) all exhibit two extratropical high-altitude (near 1–10 hPa) prograde jets. In ExoCAM’s Ben 2 simulation, these jets extend to the troposphere and merge with the two aforementioned tropospheric jets. Simulations in the Rhines rotator regime (ROCKE-3D in the Ben 1 case; all but ExoCAM in the Ben 2 case) all exhibit an equatorial (near 1–10 hPa) prograde jet and two high-latitude midaltitude (near 10–100 hPa) retrograde jets.

While we do not show wind patterns above 1 hPa in the atmosphere here, because some of the models (the UM specifically) have a model top that stops just above this pressure, we do observe strong discrepancies near 1 hPa and above. These discrepancies can originate from differences in radiative forcings (see, e.g., Figures 2(c) and 3(c)) and differences in the location and spacing of vertical levels (this choice was left free in the THAI protocol), but most likely from differences in numerical damping parameterization (e.g., by the presence and formulation of a sponge layer). As detailed in Section 2, each model uses its own recipe to simulate wind braking (as well as dissipation to smaller scales) with varying degrees of amplitude and levels of sophistication. For example, winds above 1 hPa are close to zero in ROCKE-3D simulations (for both Ben 1 and Ben 2 simulations), which comes from the presence of enhanced Rayleigh damping in the topmost layers of the model, while LMD-G exhibits strong, superrotating winds above 1 hPa (for both Ben 1 and Ben 2 simulations), which is made possible by the absence of a sponge layer. Wind braking has indeed already been shown to have a significant impact on the stratospheric circulation (Carone et al. 2018). Additional sensitivity studies performed with ROCKE-3D using moderate Rayleigh damping (the default ROCKE-3D sponge layer configuration) reestablish a high-altitude stratospheric superrotating jet. These sensitivity studies also reveal that the impact of the upper atmosphere damping parameterization on the climate state of the lower atmosphere is very weak (e.g., mean surface temperature only differs by 0.1 K).

Although the direct impact on the lower atmosphere and the surface of the presence (or absence) of a superrotating jet appears to be limited, this could have an indirect impact by affecting the photochemistry and the distribution of aerosols in the stratosphere, which could also impact the observations of the planet (Chen et al. 2019; Boutle et al. 2020). There is also a causality in the opposite direction: the radiative heating due to different gases drives the atmospheric circulation in the upper atmosphere (and may be one of the causes of intermodel differences). Given the proximity of these dynamical features to the top of our model domains and the additional processes that become important at low pressures, we reserve a full analysis of this element to a future work, where the additional model treatments can be included and maximum altitude raised.

## 4. Conclusions

In this manuscript, we reported the results of the first part of the THAI project, which compares simulation results of four state-of-the-art 3D global climate models (ExoCAM, LMD-Generic, ROCKE-3D, the Unified Model) for the exoplanet TRAPPIST-1e. This first part consisted of comparing and analyzing the results of the so-called “Ben 1” and “Ben 2” cases, which assume that the planet has a dry surface and an N₂-dominated and CO₂-dominated atmosphere, respectively.

These simulations were designed as benchmarks to test how radiative processes, small-scale dynamics (dry turbulence and convection), and large-scale dynamics in 3D GCMs impact the climate of TRAPPIST-1e, whether the GCMs agree with each other, and what implications might the intermodel differences have for observable features of this planet.

To first order, the four models give results in fairly good agreement. The intermodel spread in the global mean surface temperature amounts to 7 K for the N₂-dominated atmospheric composition and 6 K for the CO₂-dominated one. The radiative
fluxes (both stellar and thermal) are also remarkably similar (variations of radiative heating rates between the four models of always less than 5%), from the surface (≈1000 hPa) up to atmospheric pressures ≈5 hPa.

Moderate differences between the models appear in the atmospheric circulation pattern (winds) and the (stratospheric) thermal structure. These differences arise between the models from (1) large-scale dynamics, because TRAPPIST-1e lies at the tipping point between two different circulation regimes (fast and Rhines rotators) in which the models can be alternatively trapped, and (2) parameterizations used in the upper atmosphere such as numerical damping (e.g., the presence and formulation of a sponge layer).

Transit spectra (computed in THAI part 3; see Fauchez et al. 2022) are shown to be sensitive to atmospheric pressures down to 1 and $1 \times 10^{-2}$ Pa for the Ben 1 and Ben 2 cases, respectively, in the core of the 4.3 µm CO$_2$ band (Fauchez et al. 2022), and considering the typical spectral resolution of the instruments ($R \sim 10^2-10^3$) on board JWST. At the typical spectral resolution of ground-based spectrographs ($R \sim 10^5$), transit spectra are found to be sensitive to atmospheric pressures down to $1 \times 10^{-3}$ Pa and $1 \times 10^{-7}$ Pa for the Ben 1 and Ben 2 cases, respectively, in the core of CO$_2$ lines (of the 4.3 µm CO$_2$ band). Efforts should thus be made to improve the description of the processes affecting the thermal structures of planets above 1 hPa or so in GCMs because temperature-induced variations of the stratospheric scale height can significantly alter the amplitude of the strongest molecular absorption bands (Fauchez et al. 2022).

Overall, the main objective of these dry, benchmark simulations has been achieved. The simulations demonstrate that the main differences between the models presented for the aquaplanet simulations (described in part 2 of the THAI trilogy; see Sergeev et al. 2022a) come mostly from water vapor and clouds (although differences due to convection and RT are amplified when water vapor and clouds are present). Similarly, most of the differences in the simulated transit spectra of TRAPPIST-1e (described in part 3 of the THAI trilogy; see Fauchez et al. 2022) are dominated by the variations in the 3D distribution of water vapor and clouds.

This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 832738/ESCAPE. M.T. thanks the Gruber Foundation for its generous support to this research. M.T. acknowledges support from the PORTAL BRAIN-be 2.0 BELSPO project. M.T. thanks Ray Pierrehumbert and Jeremy Leconte for useful discussions on this work. M.T. thanks Shuang Wang for helping to track a mistake in the original version of the manuscript. M.T. was granted access to the High-Performance Computing (HPC) resources of Centre Informatique National de l'Enseignement Supérieur (CINES) under the allocations N° A0020101167 and A0040110391 made by Grand Équipement National de Calcul Intensif (GENCI). This work has been carried out within the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation. M.T. acknowledges the financial support of the SNSF.

M.T. and F.F. thank the LMD Generic Global Climate team for the teamwork development and improvement of the model. T.F. and M.J.W. acknowledge support from the GSFC Sellers Exoplanet Environments Collaboration (SEEC), which is

funded in part by the NASA Planetary Science Divisions Internal Scientist Funding Model. M.J.W. was supported by NASA's Nexus for Exoplanet System Science (NExSS). Resources supporting the ROCKE-3D simulations were provided by the NASA High-End Computing (HEC) Program through the NASA Center for Climate Simulation (NCCS) at Goddard Space Flight Center. D.E.S., I.A.B., F.H.L, J.M., and N.J.M. acknowledge use of the Monsoon system, a collaborative facility supplied under the Joint Weather and Climate Research Programme, a strategic partnership between the Met Office and the Natural Environment Research Council. We acknowledge support of the Met Office Academic Partnership secondment program. This work was partly supported by a Science and Technology Facilities Council Consolidated grant (ST/R000395/1), UKRI Future Leaders Fellowship (MR/T040866/1), and the Leverhulme Trust (RPG-2020-82). J.H.M. acknowledges funding from the NASA Habitable Worlds program under award 80NSSC20K0230.

The authors thank the two anonymous reviewers whose comments helped to improve the quality and clarity of this manuscript.

The THAI GCM intercomparison team is grateful to the Anong's THAI Cuisine restaurant in Laramie for hosting its first meeting on 2017 November 15.

Numerical experiments performed for this study required the use of supercomputers, which are energy-intensive facilities and thus have nonnegligible greenhouse gas emissions associated with them. We estimate that the final versions of Ben 1 and Ben 2 experiments resulted in ≈840 kg CO$_2$e, including ≈750 kg from ExoCAM runs, ≈20 kg from LMD-G, ≈70 kg from ROCKE-3D, and ≈3 kg from the UM.

**Software:** Scripts to process and visualize THAI data are available on GitHub, https://github.com/projectcuisines/thai_triology_code, and are dependent on the following Python libraries: aeolus (Sergeev & Zamyatina 2021), iris (Met Office 2021), jupyter (Kluyver et al. 2016), matplotlib (Hunter 2007), numpy (Harris et al. 2020), windspharm (Dawson 2016), xarray (Hoyer & Hamman 2017). ExoCAM (Wolf & Toon 2015) is available on GitHub: https://github.com/storyofthewolf/ExoCAM. The radiative transfer component of ExoCAM is available on GitHub: https://github.com/storyofthewolf/ExoRT. The National Center for Atmospheric Research Community Earth System Model version 1.2.1, a prerequisite for using ExoCAM, is available publicly at https://www.cesm.ucar.edu/models/cesm1.2/tags/. The Met Office Unified Model is available for use under licence; see https://www.metoffice.gov.uk/research/approach/modelling-systems/unified-model. ROCKE-3D is public domain software and available for download for free from https://simplex.giss.nasa.gov/gcm/ROCKE-3D/. Annual tutorials for new users take place annually, whose recordings are freely available online at https://www.youtube.com/user/NASAGISStv/playlists?view=50&sort=dd&shelf_id=15. LMD-G is available upon request from Martin Turbet (martin.turbet@lmd.ipsl.fr) and François Forget (francois.forget@lmd.ipsl.fr).

## Appendix

### Data Accessibility

All our GCM THAI data are permanently available for download at https://ckan.emac.gsfc.nasa.gov/organization/thai, with variables described for each data set. If you use these data, please cite the current paper and add the following
statement: "THAI data have been obtained from https://ckan.
emac.gsfc.nasa.gov/organization/thai, a data repository of the
Sellers Exoplanet Environments Collaboration (SEEC), which
is funded in part by the NASA Planetary Science Divisions
Internal Scientist Funding Model."

Scripts to process the THAI data are available on GitHub:
https://github.com/projectcuisines

ORCID iDs

Martin Turbet https://orcid.org/0000-0003-2260-9856

Thomas J. Fauchez https://orcid.org/0000-0002-5967-9631

Denis E. Sergeev https://orcid.org/0000-0001-8832-5288

Ian A. Boutle https://orcid.org/0000-0002-1485-4475

Kostas Tsigaridis https://orcid.org/0000-0001-5328-819X

Michael J. Way https://orcid.org/0000-0003-3728-0475

Eric T. Wolf https://orcid.org/0000-0002-7188-1648

Shawn D. Domagal-Goldman https://orcid.org/0000-0003-0354-9325

François Forget https://orcid.org/0000-0002-3262-4366

Jacob Haqq-Misra https://orcid.org/0000-0003-4346-2611

Ravi K. Kopparapu https://orcid.org/0000-0002-5893-2471

F. Hugo Lambert https://orcid.org/0000-0002-4664-1327

James Manners https://orcid.org/0000-0003-4402-6811

Nathan J. Mayne https://orcid.org/0000-0001-6707-4563

Linda Sohl https://orcid.org/0000-0002-6673-2007

References

Amundsen, D. S., Tremblin, P., Manners, J., Baraffe, I., & Mayne, N. J. 2017, A&A, 598, A97

Arakawa, A. 1972, Design of the UCLA General Circulation Model, Numerical Simulation of Weather and Climate, Tech. Rep. 7, NASA, https://ntrs.nasa.gov/citations/19730012781

Baranov, Y. I., Lafferty, W. J., & Fraser, G. T. 2004, JMoSp, 228, 432

Bardeen, C. G., Garcia, R. R., Toon, O. B., & Conley, A. J. 2017, PNAS, 114, E7415

Boutle, I. A., Joshi, M., Lambert, F. H., et al. 2020, NatCo, 11, 2731

Brown, A. R., Beare, R. J., Edwards, J. M., et al. 2008, BoLMe, 128, 117

Carone, L., Keppens, R., & Decin, L. 2015, MNRAS, 453, 2412

Carone, L., Keppens, R., & Decin, L. 2016, MNRAS, 461, 1981

Carone, L., Keppens, R., Decin, L., & Henning, T. 2018, MNRAS, 473, 4672

Chen, H., Wolf, E. T., Zhan, Z., & Horton, D. E. 2019, ApJ, 886, 16

Cheng, Y., Canuto, V. M., & Howard, A. M. 2002, JAtS, 59, 1550

Clough, S. A., Shephard, M. W., Mlawer, E. J., et al. 2005, JQSRT, 91, 233

Colaprete, A., & Toon, O. B. 2003, JGRE, 108, 5025

Dawson, A. 2016, JORS, 4, 31

Ding, F., & Wordsworth, R. D. 2021, PSJ, 2, 201

Edson, A., Lee, S., Bannon, P., Kasting, J. F., & Pollard, D. 2011, Icar, 212, 1

Edwards, J. M., & Slingo, A. 1996, QJRMS, 122, 689

Eyring, V., Bony, S., Meehl, G. A., et al. 2016, GMD, 9, 1937

Fauchez, T. J., Turbet, M., Sergeev, D. E., et al. 2021, PSJ, 2, 106

Fauchez, T. J., Turbet, M., Villanueva, G. L., et al. 2019, ApJ, 887, 194

Fauchez, T. J., Turbet, M., Wolf, E. T., et al. 2020, GMD, 13, 707

Fauchez, T. J., Villanueva, G. L., Sergeev, D. E., et al. 2022, PSJ, 3, 213

Forget, F., Hourdin, F., Fournier, R., et al. 1999, JGR, 104, 24155

Galperin, B., Kantha, L. H., Hassid, S., & Rosati, A. 1988, JAtS, 45, 55

Gillon, M., Meadows, V., Agol, E., et al. 2020, BAAS, 52, 0208

Gillon, M., Triaud, A. H. M. J., Demory, B.-O., et al. 2017, Natur, 542, 456

Greene, T. P., Line, M. R., Montero, C., et al. 2016, ApJ, 817, 17

Grimm, S. L., Demory, B.-O., Gillon, M., et al. 2018, A&A, 613, A68

Gruszka, M., & Borysow, A. 1997, Icar, 129, 172

Halevy, I., Pierrehumbert, R. T., & Schrag, D. P. 2009, JGRD, 114, D18112

Hammond, M., & Lewis, N. T. 2021, PNAS, 118, e2022705118

Hammond, M., Tsai, S.-M., & Pierrehumbert, R. T. 2020, ApJ, 901, 78

Hansen, J., Russell, G., Rind, D., et al. 1983, MWRv, 111, 609

Hansen, J. E., & Travis, L. D. 1974, SSRv, 16, 527

Haqq-Misra, J., Wolf, E. T., Joshi, M., Zhang, X., & Kopparapu, R. K. 2018, ApJ, 852, 67

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Holtslag, A. A. M., & Boville, B. A. 1993, JCI, 6, 1825

Hourdin, F., Musat, I., Bony, S., et al. 2006, CIDy, 27, 787

Hoyer, S., & Hamman, J. J. 2017, JORS, 5, 1

Hunter, J. D. 2007, CSE, 9, 90

Joshi, M. M., Haberle, R. M., & Reynolds, R. T. 1997, Icar, 129, 450

Kasting, J. F., Pollack, J. B., & Crisp, D. 1984, JAtC, 1, 403

Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and Power in Academic Publishing: Players, Agents and Agendas, ed. F. Loizides & B. Scmidt (Netherlands: IOS Press), 87, https://eprints.soton.ac.uk/403913/

Kopparapu, R. k., Wolf, E. T., Arney, G., et al. 2017, ApJ, 845, 16

Lin, S.-J., & Rood, R. 1996, MWRv, 124, 2046

Lock, A. P. 2001, MWRv, 129, 1148

Lock, A. P., Brown, A. R., Bush, M. R., Martin, G. M., & Smith, R. N. B. 2000, MWRv, 128, 3187

Lustig-Yaeger, J., Meadows, V. S., & Lincowski, A. P. 2019, AJ, 158, 27

Manners, J., Vosper, S. B., & Roberts, N. 2012, QJRMS, 138, 720

Mellor, G. L., & Yamada, T. 1982, RvGSP, 20, 851

Merlis, T. M., & Schneider, T. 2010, JAMES, 2, 13

Met Office 2021, Iris: A Python Library for Analysing and Visualising Meteorological and Oceanographic Data Sets, http://scitools.org.uk/

Mlawer, E. J., Taubman, S. J., Brown, P. D., Iacono, M. J., & Clough, S. A. 1997, JGR, 102, 16,663

Moeng, C.-H., & Sullivan, P. P. 1994, JAtS, 51, 999

Morley, C. V., Kreidberg, L., Rustamkulov, Z., Robinson, T., & Fortney, J. J. 2017, ApJ, 850, 121

Neale, R. B., Gettelman, A., Park, S., et al. 2010, Description of the NCAR Community Atmosphere Model (CAM 4.0), Tech. Note NCAR/TN-486 +STR, NCAR, https://www.cesm.ucar.edu/models/ccsm4.0/cam/docs/description/cam4_desc.pdf

Noda, S., Ishiwatari, M., Nakajima, K., et al. 2017, Icar, 282, 1

Penn, J., & Vallis, G. K. 2018, ApJ, 868, 147

Perrin, M. Y., & Hartmann, J. M. 1989, J. Quant. Spec. Radiat. Transf., 42, 311

Pierrehumbert, R. T., & Hammond, M. 2019, AnRFM, 51, 275

Richard, C., Gordon, I. E., Rothman, L. S., et al. 2012, JQSRT, 113, 1276

Robinson, T. D., & Catling, D. C. 2014, NatGe, 7, 12

Rothman, L. S., Gordon, I. E., Babikov, Y., et al. 2013, JQSRT, 130, 4

Rothman, L. S., Jacquemart, D., Barbe, A., et al. 2005, JQSRT, 96, 139

Schmidt, G. A., Ruedy, R., Hansen, J. E., et al. 2006, JCI, 19, 153

Sergeev, D. E., Fauchez, T. J., Turbet, M., et al. 2022a, PSJ, 3, 212

Sergeev, D. E., Lewis, N. T., & Lambert, F. H. 2022b, PSJ, 3, 214

Sergeev, D. E., Lambert, F. H., Mayne, N. J., et al. 2020, ApJ, 894, 84

Sergeev, D. E., & Zamyatina, M. 2021, Aeolus—a Python Library for the Analysis and Visualisation of Climate Model Output, Zenodo, doi:10.5281/zenodo.5145604

Smith, R. N. B. 1990, QJRMS, 116, 435

Tolkien, J. R. R. 1954, The Fellowship of the Ring (London: George Allen & Unwin)

Toon, O. B., McKay, C. P., Ackerman, T. P., & Santhanam, K. 1989, JGR, 94, 16287

Tran, H., Boulet, C., Stefani, S., Snels, M., & Piccioni, G. 2011, JQSRT, 112, 925

Turbet, M., Bolmont, E., Bourrier, V., et al. 2020, SSRv, 216, 100

Turbet, M., Bolmont, E., Leconte, J., et al. 2018, A&A, 612, A86

Turbet, M., & Selsis, F. 2021, arXiv:2104.06391

Urata, R. A., & Toon, O. B. 2013a, Icar, 226, 229

Urata, R. A., & Toon, O. B. 2013b, Icar, 226, 336

Vardavas, I. M., & Carver, J. H. 1984, P&SS, 32, 1307

Walters, D., Boutle, I., Brooks, M., et al. 2017, GMD, 10, 1487

Wang, S., & Yang, J. 2021, ApJ, 907, 28

Wicker, W., Greatbatch, R. J., & Claus, M. 2022, QJRMS, 148, 199

Wolf, E. T., Kopparapu, R., Haqq-Misra, J., & Fauchez, T. J. 2022, PSJ, 3, 7

Wolf, E. T., Shields, A. L., Kopparapu, R. K., Haqq-Misra, J., & Toon, O. B. 2017, ApJ, 837, 107

Wolf, E. T., & Toon, O. B. 2013, AsBio, 13, 1

Wolf, E. T., & Toon, O. B. 2015, JGRD, 120, 5775

Wood, N., Staniforth, A., White, A., et al. 2014, QJRMS, 140, 1505

Wordsworth, R., Forget, F., & Eymet, V. 2010, Icar, 210, 992

Wordsworth, R. D., Forget, F., Selsis, F., et al. 2011, ApJL, 733, L48

Wunderlich, F., Godolt, M., Grenfell, J. L., et al. 2019, A&A, 624, A49

Wunderlich, F., Scheucher, M., Godolt, M., et al. 2020, ApJ, 901, 126

Yang, J., Leconte, J., Wolf, E. T., et al. 2019, ApJ, 875, 46

Zhang, X. 2020, RAA, 20, 099