# The TRAPPIST-1 Habitable Atmosphere Intercomparison (THAI). II. Moist Cases-The Two Waterworlds

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

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

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

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

► To cite this version:

Denis E. Sergeev, Thomas J. Fauchez, Martin Turbet, Ian A. Boutle, Kostas Tsigaridis, et al. The TRAPPIST-1 Habitable Atmosphere Intercomparison (THAI). II. Moist Cases-The Two Waterworlds. The Planetary Science Journal, 2022, 3, 10.3847/PSJ/ac6cf2. insu-03847093

HAL Id: insu-03847093

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

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_2.png" /></figure>

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

# The TRAPPIST-1 Habitable Atmosphere Intercomparison (THAI). II. Moist Cases—The Two Waterworlds

Denis E. Sergeev¹, Thomas J. Fauchez²,³,⁴, Martin Turbet⁵,⁶, Ian A. Boutle⁷,⁸, Kostas Tsigaridis⁹,¹⁰, Michael J. Way⁴,¹⁰,¹¹, Eric T. Wolf⁴,¹²,¹³, Shawn D. Domagal-Goldman²,⁴,¹³, François Forget⁶, Jacob Haqq-Misra¹³,¹⁴, Ravi K. Kopparapu²,⁴,¹³, F. Hugo Lambert¹, James Manners⁷, and Nathan J. Mayne⁸

¹Department of Mathematics, College of Engineering, Mathematics, and Physical Sciences, University of Exeter, Exeter, EX4 4QF, UK; d.sergeev@exeter.ac.uk

²NASA Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20771, USA

³Goddard Earth Sciences Technology and Research (GESTAR), Universities Space Research Association (USRA), Columbia, MD 7178, USA

⁴NASA GSFC Sellers Exoplanet Environments Collaboration

⁵Observatoire Astronomique de l'Université de Genève, Université de Genève, Chemin des Maillettes 51, 1290 Versoix, Switzerland

⁶Laboratoire de Météorologie Dynamique/IPSL, CNRS, Sorbonne Université, École Normale Supérieure, PSL Research University, École Polytechnique, F-75005 Paris, France

⁷Met Office, FitzRoy Road, Exeter, EX1 3PB, UK

⁸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, USA

¹⁰NASA Goddard Institute for Space Studies, 2880 Broadway, New York, NY 10025, USA

¹¹Theoretical Astrophysics, Department of Physics and Astronomy, Uppsala University, Uppsala, Sweden

¹²Laboratory for Atmospheric and Space Physics, University of Colorado Boulder, Boulder, CO, USA

¹³NAŞA NExSS Virtual Planetary Laboratory, Seattle, WA, 98195, USA

¹⁴Blue Marble Space Institute of Science, Seattle, WA, USA

Received 2021 September 24; revised 2022 March 30; accepted 2022 April 25; published 2022 September 15

## Abstract

To identify promising exoplanets for atmospheric characterization and to make the best use of observational data, a thorough understanding of their atmospheres is needed. Three-dimensional general circulation models (GCMs) are one of the most comprehensive tools available for this task and will be used to interpret observations of temperate rocky exoplanets. Due to parameterization choices made in GCMs, they can produce different results, even for the same planet. Employing four widely used exoplanetary GCMs—ExoCAM, LMD-G, ROCKE-3D, and the UM— we continue the TRAPPIST-1 Habitable Atmosphere Intercomparison by modeling aquaplanet climates of TRAPPIST-1e with a moist atmosphere dominated by either nitrogen or carbon dioxide. Although the GCMs disagree on the details of the simulated regimes, they all predict a temperate climate with neither of the two cases pushed out of the habitable state. Nevertheless, the intermodel spread in the global mean surface temperature is nonnegligible: 14 K and 24 K in the nitrogen- and carbon dioxide-dominated case, respectively. We find substantial intermodel differences in moist variables, with the smallest amount of clouds in LMD-Generic and the largest in ROCKE-3D. ExoCAM predicts the warmest climate for both cases and thus has the highest water vapor content and the largest amount and variability of cloud condensate. The UM tends to produce colder conditions, especially in the nitrogen-dominated case due to a strong negative cloud radiative effect on the day side of TRAPPIST-1e. Our study highlights various biases of GCMs and emphasizes the importance of not relying solely on one model to understand exoplanet climates.

*Unified Astronomy Thesaurus concepts: Exoplanet atmospheres (487); Planetary atmospheres (1244); Habitable planets (695); Habitable zone (696); Water vapor (1791); Atmospheric circulation (112); Atmospheric clouds (2180); Atmospheric composition (2120); Exoplanet atmospheric variability (2020); Planetary climates (2184)*

## 1. Introduction

"You must understand, young Hobbit, it takes a long time to say anything in Old Entish. And we never say anything unless it is worth taking a long time to say." —J.R.R. Tolkien, *The Two Towers* (1954)

With the advent of powerful new telescopes, such as the James Webb Space Telescope (JWST) and the Extremely Large Telescope, the characterization of a temperate rocky planet

orbiting an M-dwarf star becomes highly likely (Snellen et al. 2015; Anglada-Escudé et al. 2016; Lovis et al. 2017; Gillon et al. 2017, 2020; Fauchez et al. 2019; Turbet et al. 2020). In order to make the most use of the vast amount of data these telescopes will collect, a solid understanding of extraterrestrial atmospheres is needed. Three-dimensional general circulation models (or global climate models, GCMs) represent the complexity of planetary atmospheric dynamics to the best of our present knowledge. Being complex pieces of code, different GCMs typically produce different results even for the well-constrained climate of Earth (e.g., Webb et al. 2015; Eyring et al. 2016; Ceppi et al. 2017). Almost all previous exoplanetary 3D GCM studies have employed only a single GCM due to the challenges of simulating exotic atmospheres on one hand and the exciting variety of possible targets on the other. It is therefore unclear how accurate predictions from a single GCM will be when used for interpreting observations of a particular exoplanet.

<figure><img src="image_4.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.
A concerted effort is required to compare predictions of exoplanet GCMs against each other. Such model intercomparison will provide a more robust reference for explaining the observational data and will serve as a rigorous benchmark for future GCM studies. Along with theory and abundant observations, model intercomparisons are the bedrock of our understanding of modern and future Earth climate (e.g., Eyring et al. 2016). For exoplanets, observations are extremely scarce, making a model intercomparison an especially valuable tool for our scientific understanding of planetary atmospheres. The only 3D model intercomparison conducted for a generic exoplanet setup was carried out by Yang et al. (2019), who reported significant intermodel differences for slowly rotating planets irradiated by M-dwarf stars.

Yang et al. (2019) found that the water vapor radiation feedback was the major cause of the intermodel discrepancy, in agreement with previous 1D intercomparison (Yang et al. 2016); radiative properties of clouds in GCMs were also found to be important. More recently, participants in the 2020 THAI workshop noted that intermodel differences in subgrid parameterizations continue to be a key source of disagreement between simulated climates (Fauchez et al. 2021). The intermodel disagreement due to the treatment of clouds is well known in Earth climate science. For paleoclimate reconstructions, it continues to fuel the debate around the Faint Young Sun paradox (e.g., Goldblatt et al. 2021) and Snowball deglaciations (e.g., Abbot et al. 2012). For future climate projections, it is thought to be the largest source of uncertainty in equilibrium climate sensitivity<fnref n="15" /> (Ceppi et al. 2017). Interestingly, the equilibrium climate sensitivity did not reduce in the latest climate model intercomparison—Coupled Model Intercomparison Project phase 6 (CMIP6)—mainly because of stronger positive cloud feedbacks from decreasing low cloud amount and, as a result, albedo of the extratropics (Zelinka et al. 2020). Clouds are tightly coupled with moist convection and so can be affected by convective parameterizations in GCMs (Sherwood et al. 2014), though the direct effect of convection on Earth is thought to be of secondary importance (Webb et al. 2015). However, the direct effect of convection on exoplanetary climate predictions from GCMs is still unclear (Sergeev et al. 2020).

The purpose of this study is to analyze intermodel differences for two scenarios of an aquaplanet climate of TRAPPIST-1e. In particular, we aim to revisit the moist climate on a tidally locked exoplanet residing within the habitable zone (HZ) of TRAPPIST-1, an ultracool M-dwarf star, and discuss how moist physics, parameterized differently in different GCMs, affects the temperature regime of the planet, its atmospheric circulation, and cloud variability. Our results highlight various biases of the THAI GCMs, providing guidance for terrestrial exoclimatology studies, the majority of which so far are done with these four models. This part of THAI is a logical next step after the dry cases (Part I; see Turbet et al. 2021), increasing the model complexity by including the effects of water vapor and clouds. Moreover, understanding synthetic observations (Part III; see Fauchez et al. 2022) cannot be complete without a thorough understanding of the simulated climates.

This manuscript is structured as follows. In Section 2, we describe the experimental setup and key differences in the representation of moist physics in the four THAI GCMs. In Section 3, we present the details of the global climate in each

<table><caption>**Table 1**
TRAPPIST-1 Stellar Spectrum and Planetary Parameters of TRAPPIST-1e
(Grimm et al. 2018)</caption><thead><tr><th>Parameter</th><th>Units</th><th>Value</th></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 day</td><td>6.1</td></tr><tr><td>Rotation period</td><td>Earth day</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>5797</td></tr><tr><td>Gravity</td><td>m s<sup>-2</sup></td><td>9.12</td></tr></tbody></table>

group of simulations: N₂-dominated (Section 3.1) and CO₂-dominated (Section 3.2), starting with the radiative fluxes and thermodynamic profiles and continuing to the cloud cover and atmospheric circulation, before briefly discussing the time variability. In Section 4, we discuss overall commonalities and discrepancies between the GCMs. Section 5 concludes the paper.

## 2. Methodology

### 2.1. Planetary Configuration

In this Part II of the THAI trilogy, we conduct moist atmosphere aquaplanet simulations, referred to as Hab 1 & Hab 2 here and in the THAI protocol (Fauchez et al. 2020). These simulations are configured to represent two potentially habitable states of TRAPPIST-1e, assuming the planet is tidally locked to its host star (Table 1). Several modeling studies involving GCMs (Wolf et al. 2017; Turbet et al. 2018; Fauchez et al. 2019) show that it was possible for TRAPPIST-1e to retain water in its atmosphere or on its surface once it achieved temperate conditions in the HZ (Kopparapu et al. 2013). Hab 1 is a colder state with an N₂-dominated atmosphere, 400 ppm of CO₂, while Hab 2 is a warmer state with a CO₂-dominated atmosphere. In both cases, the H₂O is a condensible species, and the mean atmospheric pressure at the surface is 1 bar. CO₂ condensation has been disabled, as not all four GCMs possess this parameterization. Hab 1 represents a habitable atmosphere broadly similar to that of modern Earth. Hab 2 could be representative of the early Hadean period on Earth or Hesperian period on Mars (e.g., Wordsworth 2016). The planet's surface in both Hab 1 and Hab 2 is prescribed to be a slab ocean with the heat capacity of water (≈4 × 10⁶ J m⁻³ K⁻¹) and a constant bolometric albedo, 0.06 and 0.25, above and below freezing, respectively. The slab ocean setup in our GCMs is slightly different, though it does not affect the time averaged steady state of the atmosphere. ExoCAM and the UM use a single-layer slab ocean, while in LMD-G and ROCKE-3D, the slab ocean is split into layers and therefore depends not only on the total heat capacity but also on the thermal inertia of the layers. The slab ocean depths are as follows: 100 m in ExoCAM, 20 m in LMD-G, 90 m in ROCKE-3D, and 1 m in the UM (Table 2). See Fauchez et al. (2020) for more details on the experimental setup and its motivation.

Both sets of simulations are started from an isothermal (300 K) dry atmosphere at rest. The models are integrated until they reach a steady state<fnref n="16" /> and then for a further 100 orbits,

<footnote n="15">¹⁵ The global mean surface temperature increase after the climate system fully adjusts to a sustained doubling of CO₂ relative to pre-industrial conditions.</footnote>

<footnote n="16">¹⁶ Whether a GCM has reached a steady state is determined by the absence of a long-term trend in variables such as the global mean surface temperature and the top of the atmosphere (TOA) net energy flux.</footnote>
<table><caption>**Table 2**
Grid Resolution and Slab Ocean Setup in THAI GCMs</caption><thead><tr><th>GCM</th><th>Number of Grid Points in the Horizontal<br>(Longitude × Latitude)</th><th>Number of Vertical Levels (Top Layer<br>Pressure)</th><th>Slab Ocean Depth (Number of<br>Layers)</th></tr></thead><tbody><tr><td>ExoCAM</td><td>72 × 46</td><td>51 (1 Pa)</td><td>100 m (1)</td></tr><tr><td>LMD-G</td><td>72 × 46</td><td>40 (4 Pa)</td><td>20 m (18)</td></tr><tr><td>ROCKE-3D</td><td>72 × 46</td><td>40 (10 Pa)</td><td>90 m (13)</td></tr><tr><td>UM</td><td>144 × 90</td><td>Hab 1: 41 (4 Pa) Hab 2: 38 (13 Pa)</td><td>1 m (1)</td></tr></tbody></table>

with the analysis for this study performed on the data from these 100 orbits (610 Earth days). The output is provided every 6 hr in order to sufficiently resolve atmospheric variability.

## 2.2. Models

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

1. the Exoplanet Community Atmospheric Model (Exo-CAM, 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-G),

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

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

The horizontal and vertical resolution of each model are given in Table 2. We refer to Part I (Turbet et al. 2021) of this trilogy for a further description of dynamical cores and model grids. Here we provide a brief description of the model components related to the moist physics. Note that no model parameterizations were specifically adapted for the THAI intercomparison and thus the model configurations are close to what has been used in their recent single-model studies. A comparison of moist physics parameterizations in THAI GCMs is given in Table 3.

### 2.2.1. Radiative Transfer

While the main description of radiative transfer (RT) parameterizations is given in THAI Part I (Turbet et al. 2021), we provide additional details in this section on how these parameterizations deal with water vapor and cloud particles. The RT schemes used by THAI GCMs share a number of characteristics: they are based on the two-stream correlated-k approach; the HITRAN database is used for spectroscopic data, and MT_CKD formalism is used for continuum absorption.<fnref n="17" /> They take into account the scattering effects of the atmosphere (Rayleigh scattering) and H₂O clouds (Mie scattering). There are however, several key differences among the schemes, as outlined below.

ExoCAM simulations are performed with ExoRT, whose oldest and most published version is used here; it is referred to as n28archean (Wolf & Toon 2013; Yang et al. 2016). This version has 28 bands covering the spectral range between 0.2

and 1000 µm. Note that subsequent intercomparison studies have found this version of the RT to overestimate the greenhouse effect from pure CO₂ atmospheres (Wolf et al. 2022), and overestimate the near-infrared absorption of H₂O for warm moist atmospheres irradiated by M-dwarf stars (Yang et al. 2016). Thus, particularly for Hab 2, we would expect this version of RT to have a bias toward producing a climate that is too warm based on clear-sky RT considerations only. Later iterations of ExoRT have corrected these issues, primarily by increasing the spectral resolution in the near-infrared and by using a newer version of MT_CKD continuum, as is discussed in Kopparapu et al. (2017) and Wolf et al. (2022). Liquid and ice cloud particles are treated as Mie scattering particles with refractive indices taken from Segelstein (1981) and Warren & Brandt (2008). Liquid cloud particle effective radii are fixed at 14 µm, while ice crystal sizes are allowed to vary depending on the air temperature (see Fauchez et al. 2022, for more detail on the particle sizes and their implications for observations). Note that the cloud condensation nuclei (CCN) number used for the effective radii of cloud particles in ExoRT does not affect particle sizes in the large-scale cloud parameterization (see below). The radiative effect of cloud overlap is treated using the Monte Carlo Independent Column Approximation while assuming maximum-random overlap (Pincus et al. 2003; Barker et al. 2008).

In LMD-G's RT module (Wordsworth et al. 2011), the correlated-k coefficients were 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. Mie scattering by cloud particles is parameterized following Hansen & Travis (1974), and their refractive indices cloud particles are taken from Warren (1984). The mean cloud particle radius is calculated from (i) the CCN number and (ii) the amount of condensed H₂O (Equation (3) in Turbet et al. 2020). As described in Charnay et al. (2013), RT is computed twice: in a clear-sky column, and in a cloudy sky column whose area is equal to the total cloud cover.

Both ROCKE-3D and the UM share a common RT module SOCRATES,<fnref n="18" /> which is an open-source model developed by the Met Office (Edwards & Slingo 1996; Manners et al. 2022). SOCRATES is highly flexible. Spectral files, created for both shortwave (0.2–20 µm) and longwave (3.33–10,000 µm) streams separately, contain all necessary information to define the RT problem and are tailored to specific atmospheric compositions and stellar spectra, then input to the GCMs at run-time. For both Hab 1 and Hab 2 simulations, the UM and ROCKE-3D use identical pairs of spectral files.<fnref n="19" /> For Hab 1,

<footnote n="17">¹⁷ ExoCAM and LMD-G use the MT_CKD v2.5, while ROCKE-3D and the UM use MT_CKD v3.0. The newer version leads to a better treatment of H₂O shortwave absorption in the near-infrared, which can be particularly important to H₂O-rich climates (Yang et al. 2016; Wolf et al. 2022).</footnote>

<footnote n="18">¹⁸ https://code.metoffice.gov.uk/trac/socrates</footnote>

<footnote n="19">¹⁹ Available at https://portal.nccs.nasa.gov/GISS_modelE/ROCKE-3D/spectral_files/.</footnote>
<table><caption>**Table 3**
Summary of Moist Physics Parameterizations in THAI GCMs</caption><thead><tr><th></th><th>Scheme Components</th><th>ExoCAM</th><th>LMD-G</th><th>ROCKE-3D</th><th>UM</th></tr></thead><tbody><tr><td rowspan="6">Radiative transfer</td><td>Type of scheme</td><td>Two-stream, correlated-k</td><td>Two-stream, correlated-k</td><td>Two-stream, correlated-k</td><td>Two-stream, correlated-k</td></tr><tr><td>SW spectral range &amp; number of bands</td><td>0.2–1000 µm; 28 bands</td><td>Hab 1: 0.3–5.1 µm; 36 bands<br>Hab 2: 0.3–5.1 µm; 36 bands</td><td>Hab 1: 0.2–20 µm; 21 bands<br>Hab 2: 0.2–20 µm; 42 bands</td><td>Hab 1: 0.2–20 µm; 21 bands<br>Hab 2: 0.2–20 µm; 42 bands</td></tr><tr><td>LW spectral range &amp; number of bands</td><td></td><td>Hab 1: 2.3–1000 µm; 38 bands<br>Hab 2: 2.3–1000 µm; 32 bands</td><td>Hab 1: 3.3–10,000 µm; 12 bands<br>Hab 2: 3.3–10,000 µm; 17 bands</td><td>Hab 1: 3.3–10,000 µm; 12 bands<br>Hab 2: 3.3–10,000 µm; 17 bands</td></tr><tr><td>Spectroscopy database</td><td>HITRAN2004</td><td>HITRAN2012</td><td>HITRAN2012</td><td>HITRAN2012</td></tr><tr><td>Continuum</td><td>MT_CKD v2.5</td><td>MT_CKD v2.5</td><td>MT_CKD v3.0</td><td>MT_CKD v3.0</td></tr><tr><td>Liquid particles effective radius</td><td>14 µm</td><td>Set by the CCN number and cloud water mixing ratio</td><td>Set by the CCN number and liquid water content; parameterization of thick averaging with Pade fits</td><td>Set by the CCN number and liquid water content; parameterization of thick averaging with Pade fits<br>See Edwards et al. (2007)</td></tr><tr><td rowspan="2">Convection</td><td>Ice particles effective radius</td><td>Temperature-dependent</td><td></td><td></td><td></td></tr><tr><td>Subgrid-scale cloud variability</td><td>MCICA</td><td>None</td><td></td><td>Scaling factors</td></tr><tr><td rowspan="6">Large-scale clouds &amp; precipitation</td><td>Type of scheme</td><td>Mass-flux</td><td>Adjustment</td><td>Mass-flux</td><td>Mass-flux</td></tr><tr><td>Subtypes</td><td>Deep, shallow</td><td>None</td><td>None</td><td>Deep, mid-level, shallow</td></tr><tr><td>Trigger</td><td>Conditional instability diagnosed by entraining plume ascent</td><td>Conditional instability</td><td>Conditional instability diagnosed by one undilute and one entraining plume</td><td>Conditional instability diagnosed by undilute parcel ascent</td></tr><tr><td>Closure</td><td>CAPE-based</td><td>None</td><td>Mass flux reaching neutral buoyancy</td><td>CAPE-based and dependent on the vertical velocity for deep and mid-level convection, cloud-based for shallow convection</td></tr><tr><td>Downdrafts</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>Convective momentum transport</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td rowspan="5">Large-scale clouds &amp; precipitation</td><td>Treatment of cloud condensate</td><td>Prognostic</td><td>Prognostic</td><td>Prognostic</td><td>Prognostic</td></tr><tr><td>Treatment of cloud fraction</td><td>Diagnostic</td><td>Diagnostic</td><td>Diagnostic</td><td>Prognostic</td></tr><tr><td>Liquid/ice cloud particles separation</td><td>Temperature-dependent</td><td>Temperature-dependent</td><td>Temperature-dependent</td><td>Process-based</td></tr><tr><td>Supercooled water limit</td><td>233.15 K</td><td>255.15 K</td><td>233.15 K</td><td>233.15 K</td></tr><tr><td>Column cloud fraction</td><td>Maximum-random overlap</td><td>Maximum overlap</td><td>Maximum-random overlap</td><td>Exponential random overlap</td></tr><tr><td></td><td>CCN number (m<sup>-3</sup>)</td><td>Fixed (see the text for details)</td><td>10<sup>7</sup> for liquid; 10<sup>4</sup> for ice particles</td><td>Fixed (see the text for details)</td><td>10<sup>8</sup> for liquid; variable for ice particles</td></tr><tr><td></td><td>Precipitation</td><td>Diagnostic</td><td>Diagnostic</td><td>Diagnostic</td><td>Prognostic</td></tr></tbody></table>
both models use a spectral file tailored to present-day Earth-like conditions having 21 shortwave bins (*sp_sw_21_dsa*) and 12 longwave bins (*sp_lw_12_dsa*). For Hab 2, both models use a spectral file with 42 shortwave (*sp_sw_42_dsa_mars*) and 17 longwave (*sp_lw_17_dsa_mars*) bands tailored to a CO₂-dominated atmosphere. The parameterization of cloud droplets uses the method of “thick averaging” (Edwards & Slingo 1996) and Padé fits for the variation with effective radius. Ice crystal treatment uses the parameterization of Edwards et al. (2007). The CCN number is treated consistently between the RT and cloud modules. The subgrid-scale cloud inhomogeneity is taken into account using scaling factors, which depend on whether the cloud is stratiform or convective and consists of liquid or frozen condensate.

### 2.2.2. Moist Convection

To account for the subgrid-scale transport of heat and moisture by cumulus clouds, GCMs employ convection parameterizations.

In ExoCAM, deep penetrating moist convection is parameterized following Zhang & McFarlane (1995) with modifications to the numerical solver to improve numerical stability for warm climates (Wolf & Toon 2015). This mass-flux scheme represents deep convection using a plume ensemble approach, triggering updrafts and downdrafts by conditional instability, which is diagnosed by entraining plume ascent. The parameterization is closed, i.e., the mass fluxes are determined as a function of the rate at which convective available potential energy (CAPE) is consumed. A separate parameterization for shallow convection (i.e., between three or fewer layers) is treated following Hack (1994).

In contrast to the other three GCMs, LMD-G employs a convective adjustment scheme, which is a conceptually simpler parameterization, but does not represent the full complexity of subgrid-scale processes associated with moist convection. This parameterization is based on the concept of radiative-convective equilibrium and lets convection rapidly mix heat and moisture in the vertical, effectively relaxing an unstable temperature profile to a reference state. In other words, this convective adjustment scheme tries to maintain a specified moist adiabatic lapse rate in the atmospheric column. Note that the Manabe-Wetherald scheme (Manabe & Wetherald 1967) is used instead of the more common Betts-Miller, because it is more robust for a wide range of pressures, though at the cost of giving enhanced precipitation near the substellar point (Charnay et al. 2013).

Like ExoCAM, ROCKE-3D accounts for cumulus convection via the mass-flux approach (Schmidt et al. 2014). The scheme is triggered by conditional instability calculated using two plumes: undilute and dilute (entraining). This allows for the representation of different types of cumulus clouds without invoking a separate shallow convection parameterization. The scheme's closure is based on the assumption of neutral buoyancy achieved by the mass flux at the cloud base. Relative to its parent code, the Earth GCM ModelE 2, ROCKE-3D includes a few modifications to the convection scheme: e.g., a stronger entrainment rate, more vigorous vertical transport of condensed water, and a relaxed limit for the convection top pressure (Way et al. 2017).

The UM also uses a mass-flux convection scheme, which is based on Gregory & Rowntree (1990) and developed further to improve the representation of downdrafts and convective

momentum transport (Walters et al. 2019). The convection is triggered by the surface buoyancy flux diagnosed by the boundary-layer scheme and an undilute parcel ascent to the level of neutral buoyancy. The type of convection is then classified as shallow, mid-level or deep, depending on the height of the level of neutral buoyancy relative to the freezing level and the level of strong vertical updrafts. For deep and mid-level convection, the UM uses a CAPE-based closure with a dependency on the grid-scale or vertical velocity; for shallow convection, a cloud-based closure is used.

### 2.2.3. Clouds and Precipitation

While cumulus physics in GCMs are handled by the convection scheme, grid-scale cloud processes are the responsibility of the large-scale cloud and precipitation schemes. Depending on the complexity of the model, cloud-related processes directly impact the state of the atmosphere by changing its radiative properties and by releasing latent heat of condensation or freezing. There is a spectrum of cloud parameterizations, from purely diagnostic (e.g., Smith 1990; Liu et al. 2021), to mixed prognostic-diagnostic (e.g., Sundqvist 1978), to fully prognostic (e.g., Wilson et al. 2008), as well as statistical schemes based on various moments of subgrid-scale variability. While diagnostic schemes offer transparency and simplicity of representing clouds in a GCM, the advantage of the prognostic approach is more physically based evolution of clouds, including the web of complex interactions between vapor and condensate phases, and better treatment of cloud optical properties by the RT code. As an example, a prognostic cloud scheme makes it possible for a cloud to be transported away from where it formed, as in the case of anvilsetrained from convection that can persist after the convection itself has ceased. Three THAI GCMs rely on Sundqvist-type prognostic-condensate/diagnostic-cloud schemes, while the UM uses a fully prognostic approach.

ExoCAM relies on the CAM4 moist physics package (Neale et al. 2010) to control water vapor, liquid cloud, and ice cloud condensate fields via prognostic bulk microphysical parameterizations. Cloud fraction is diagnosed in three categories: marine stratus clouds dependent on the temperature profiles; convective clouds defined by the convective mass-flux; and layered clouds forming at a sufficient relative humidity threshold. The total cloud fraction within each grid box is diagnosed via the maximum overlap assumption. ExoCAM cloud scheme uses a fixed cloud droplet number concentration, equal to $1.5 \times 10^8 \text{ m}^{-3}$, which is used to compute the conversion of liquid water to rain (Neale et al. 2010). Precipitation is treated diagnostically and can be in the form of rain or snow, or a mixture of both.

In LMD-G, the large-scale condensation is based on the prognostic equations of Le Treut & Li (1991) that track the evolution of the total cloud condensate, which is split into the liquid or ice phase according to Equation (2) of Charnay et al. (2013). Liquid cloud droplets are allowed to exist in the supercooled state down to the limit of 255.15 K, below which clouds are only composed of ice crystals. Cloud fraction is diagnosed from relative humidity. The total cloud cover of each vertical atmospheric column is taken equal to the partial cover of the optically thickest cloud of the column. LMD-G then assumes that each individual cloud cover is equal to this total cloud cover with a maximum overlap (Charnay et al. 2013)—equivalent to the “LMDG_max” simulation in Yang et al. (2019). The CCN number is set to $10^7$ and $10^4 \text{ kg}^{-1}$ for liquid and ice cloud particles, respectively, following Turbet et al. (2020). Water precipitation is divided into rainfall (Boucher et al. 1995) and
<table><caption>**Table 4**
Global Mean Surface Temperature ($T_s$, K) Global Mean Top-of-atmosphere Shortwave (CRE<sub>SW</sub>), Longwave (CRE<sub>LW</sub>), and Net Cloud Radiative Effect (CRE, W m<sup>-2</sup>), Day Side Net CRE (CRE<sub>day</sub>), Night Side Net CRE (CRE<sub>night</sub>), Planetary Albedo ($\alpha_p$), All-sky (G<sub>as</sub>) and Clear-sky (G<sub>cs</sub>) Greenhouse Effect (K), Water Vapor Volume Mixing Ratio at the 1 hPa Level (VMR<sub>1 hPa</sub>, mol mol<sup>-1</sup>) in Hab 1 and Hab 2 Simulations</caption><thead><tr><th rowspan="2">GCM</th><th colspan="10">Hab 1</th></tr><tr><th>T<sub>s</sub></th><th>CRE<sub>SW</sub></th><th>CRE<sub>LW</sub></th><th>CRE</th><th>CRE<sub>day</sub></th><th>CRE<sub>night</sub></th><th>α<sub>p</sub></th><th>G<sub>as</sub></th><th>G<sub>cs</sub></th><th>VMR<sub>1 hPa</sub></th></tr></thead><tbody><tr><td>ExoCAM</td><td>245.6</td><td>-37.9</td><td>15.5</td><td>-22.4</td><td>-46.0</td><td>-0.3</td><td>0.24</td><td>10.3</td><td>5.8</td><td>4.6 × 10<sup>-7</sup></td></tr><tr><td>LMD-G</td><td>242.2</td><td>-21.8</td><td>5.1</td><td>-16.8</td><td>-33.2</td><td>-1.2</td><td>0.22</td><td>6.0</td><td>4.8</td><td>5.8 × 10<sup>-8</sup></td></tr><tr><td>ROCKE-3D</td><td>244.0</td><td>-44.6</td><td>9.3</td><td>-35.3</td><td>-70.9</td><td>0.3</td><td>0.27</td><td>9.3</td><td>6.7</td><td>5.1 × 10<sup>-7</sup></td></tr><tr><td>UM</td><td>231.6</td><td>-48.1</td><td>10.1</td><td>-38.0</td><td>-70.9</td><td>-5.0</td><td>0.28</td><td>2.6</td><td>0.0</td><td>5.1 × 10<sup>-7</sup></td></tr><tr><th colspan="10">Hab 2</th></tr><tr><td>ExoCAM</td><td>295.2</td><td>-22.0</td><td>13.2</td><td>-8.8</td><td>-26.9</td><td>8.0</td><td>0.15</td><td>53.0</td><td>48.9</td><td>1.5 × 10<sup>-5</sup></td></tr><tr><td>LMD-G</td><td>270.9</td><td>-32.3</td><td>3.2</td><td>-29.1</td><td>-60.9</td><td>1.0</td><td>0.21</td><td>34.3</td><td>33.2</td><td>3.8 × 10<sup>-6</sup></td></tr><tr><td>ROCKE-3D</td><td>284.0</td><td>-35.0</td><td>8.6</td><td>-26.4</td><td>-58.7</td><td>5.9</td><td>0.19</td><td>44.7</td><td>41.9</td><td>4.5 × 10<sup>-5</sup></td></tr><tr><td>UM</td><td>280.4</td><td>-27.3</td><td>8.3</td><td>-19.0</td><td>-41.7</td><td>3.6</td><td>0.16</td><td>39.1</td><td>36.5</td><td>9.2 × 10<sup>-6</sup></td></tr></tbody></table>

snowfall, and is computed as in Charnay et al. (2013). For both cases, the LMD-G assumes that precipitation is instantaneous (i.e., it goes directly to the surface), though it can evaporate while falling through subsaturated layers, thus producing water vapor and latent heat cooling. The precipitation re-evaporation scheme is based on Gregory (1995) and is discussed in more detail in Charnay et al. (2013).

ROCKE-3D uses a prognostic cloud water approach with diagnostic cloud fraction (Schmidt et al. 2014; Way et al. 2017). Like in ExoCAM and LMD-G, a single cloud water variable is predicted; its phase is then diagnosed from temperature and corrected for precipitation from upper model layers. Note, however, that the cloud condensate is not advected by resolved winds or turbulence. Cloud volume fraction is a function of relative humidity, with a different threshold for the boundary layer and free troposphere. The column cloud fraction is computed using maximum-random overlap. The cloud particle number is fixed over the ocean surface (and thus for the whole atmosphere in Hab 1 and Hab 2 simulations) and is equal to $\approx 6 \times 10^7 \text{ m}^{-3}$ for liquid clouds and to $\approx 6 \times 10^4 \text{ m}^{-3}$ for ice clouds (Genio et al. 1996). Precipitation is treated diagnostically.

The UM represents large-scale clouds via the Prognostic Cloud fraction and Prognostic Condensate (PC2) scheme (Wilson et al. 2008) along with additional parameterizations of cloud erosion and critical relative humidity. The scheme includes three prognostic mixing ratios (water, liquid, and ice) and three prognostic cloud fractions (liquid, ice, and mixed-phase). The prognostic variables are updated by process-driven increments, and all cloud types are handled by single scheme, with the assumption of minimum overlap of phases. The column cloud fraction is a product of exponential random vertical overlap. The CCN number is fixed at $10^8 \text{ m}^{-3}$ for liquid droplets and is variable for ice crystals. Precipitation is formed according to a single-moment scheme based on Wilson & Ballard (1999). The warm rain is based on Boutle et al. (2014) and includes prognostic rain, allowing for 3D transport and subgrid variability of precipitation. Rain-rate-dependent particle size and fall velocities allow the UM to represent better sedimentation and evaporation of small droplets (Walters et al. 2019).

### 3. Results

For each case, Hab 1 (Section 3.1) and Hab 2 (Section 3.2), we present results from the four GCMs, comparing and contrasting them against each other and previous studies. We

report the following key aspects of the climate in our numerical simulations. We first discuss radiation fluxes at the top of the atmosphere (TOA) and at the planet's surface, following up with the analysis of temperature and humidity profiles in each of the GCMs. We then discuss the cloud amount and its spatial distribution. We then move on to comparing the GCMs with respect to the mean large-scale atmospheric circulation. We conclude the analysis of each of the cases by analyzing the time variability of cloudiness at the planet's limb.

#### 3.1. Global Climate in Hab 1 Simulations

##### 3.1.1. Radiation Fluxes and Thermodynamic Profiles

At a planetary scale, the climate of TRAPPIST-1e is determined by the amount of the shortwave stellar radiation absorbed at the TOA on the day side and the amount of longwave radiation lost to space. In this study, as in many other GCM studies, the shortwave (SW) and longwave (LW) streams of radiation refer to the radiation emitted by the host star and by the planet, respectively. In all four THAI GCMs, there is an overlap in wavelengths of the SW and LW streams, as noted in Table 3. The clear-sky absorbed radiation peaks at the substellar point, as in the dry experiments (see Part I; Turbet et al. 2021), but the cloud cover reflects a large part of shortwave radiation back to space. The impact of clouds on radiation is evident in all four GCMs and can be quantified as the cloud radiative effect (CRE), $\text{CRE} = F_{\text{clear-sky}} - F_{\text{all-sky}}$, where $F_{\text{clear-sky}}$ is the upward TOA radiative flux assuming clouds are absent<sup>20</sup> and $F_{\text{all-sky}}$ is the actual upward TOA radiative flux (e.g., Ceppi et al. 2017). Here the CRE is averaged over the whole planet. The shortwave CRE, $\text{CRE}_{\text{SW}}$, is negative in the THAI simulations (Table 4), because clouds have a higher albedo than the planet's surface and thus reduce the amount of absorbed energy. In Hab 1, the strongest $\text{CRE}_{\text{SW}}$ is exhibited by the UM ($-48.1 \text{ W m}^{-2}$) and is close to the average value for the modern Earth's climate, $-45 \text{ W m}^{-2}$ (e.g., Henderson et al. 2013). Even though the UM does not have the highest cloud amount among the GCMs on average, the cloud cover is concentrated and is the thickest on the day side of the planet (see Section 3.1.2, Figures 8 and 9). Hab 1 climate $\text{CRE}_{\text{SW}}$ estimates (multiplied by 2 because there is no $\text{CRE}_{\text{SW}}$ on the night side) are also close to the typical values

<footnote n="20"><sup>20</sup> The clear-sky radiative fluxes are calculated within the RT schemes in each GCM and provided as standard output in the THAI simulations.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Key radiation fluxes in Hab 1 simulations. Panels (a)–(d): TOA outgoing longwave radiation (TOA OLR, W m⁻²); panels (e)–(h): downward shortwave radiation flux at the surface (W m⁻²); and panels (i)–(l): net upward longwave radiation flux at the surface (W m⁻²).</figcaption></figure>

observed over the warm pools of the tropical Pacific and Indian Oceans on Earth, where $CRE_{SW}$ is 60–80 W m⁻² (Wall et al. 2019). The weakest $CRE_{SW}$ is exhibited by LMD-G (−21.8 W m⁻²) and is smaller by more than a half of that found for simulations using the other GCMs. This is primarily due to a reduced cloud fraction and small cloud water content produced by LMD-G on the day side (Section 3.1.2). The second reason is the size of cloud particles, which is set to a constant or parameterized depending on variables such as temperature and CCN number (Turbet et al. 2020; see also Section 2.2.3 for details). The effective radius of cloud droplets and crystals tends to be the largest in LMD-G compared to that in the other three models (not shown). Larger cloud particles result in a decrease of the shortwave optical depth of clouds, which leads to a decrease in the reflected radiation flux at TOA, thus making the $CRE_{SW}$ less negative.

The cloud cover, produced by the interaction between the convective (Section 2.2.2) and large-scale cloud (Section 2.2.3) parameterizations and shaped by the atmospheric circulation, leaves a distinct imprint in the TOA radiation fluxes and shifts the shortwave absorption maximum away from the substellar point in all GCMs (not shown, but the pattern is similar to that shown in Figure 1(e)–(h)). The resulting planetary albedo repeats the trend in $CRE_{SW}$ (Table 4), with LMD-G having an albedo of about 0.22, followed by ExoCAM (0.24), then ROCKE-3D (0.27), and the UM (0.28).

While the strongest TOA longwave emission is on the day side (exceeding 240 W m⁻²), a significant portion of it is on the night side too, especially in low latitudes of the eastern hemisphere (Figures 1(a)–(d)). Due to high-altitude clouds, the longwave emission has a pronounced minimum at and to the east of the substellar point, especially in ExoCAM, ROCKE-3D and the UM. The global longwave CRE, $CRE_{LW}$, is positive, because high cloud cover reduces the amount of outgoing thermal radiation; reaches 15.5 W m⁻² in the ExoCAM simulation, dropping by a third for ROCKE-3D and the UM and then even more for LMD-G. Compared to

27 W m⁻² for Earth (Henderson et al. 2013), $CRE_{LW}$ in the Hab 1 simulations is smaller, which is explained chiefly by the low $CRE_{LW}$ on the night side, where the low altitude of clouds coincides with the highest temperature due to the boundary-layer thermal inversion (Sergeev et al. 2020). The effective radii of cloud particles offers a small contribution to the intermodel spread of $CRE_{LW}$, due to its smaller effect on the longwave fluxes (e.g., Abbot 2014).

The overall effect of clouds on TRAPPIST-1e is cooling, confirming the earlier theory of clouds stabilizing the climate of tidally locked terrestrial exoplanets and expanding the inner edge of the habitable zone (Joshi 2003; Yang et al. 2013, 2014; Kopparapu et al. 2016). The magnitude of the net CRE is the largest in the UM (−38.0 W m⁻²) and ROCKE-3D (−35.3 W m⁻²), diminishing by almost a half in ExoCAM (−22.4 W m⁻²) and lower still in LMD-G (−16.8 W m⁻²; Table 4). These numbers summarize that the radiative influence of clouds on the Hab 1 climate is the strongest in the UM and the weakest in the LMD-G. The tendency of the LMD-G to produce weaker CRE than that in other GCMs, mostly due to its shortwave part, was also reported by Yang et al. (2019), although the amplitude of CRE is significantly lower in our experiments due to lower stellar flux and due to faster rotation of the planet (Kopparapu et al. 2016). To sum up, the models have different CREs because of the distribution of clouds, especially on the day side of the planet, and because of the size of cloud particles, which affects primarily the shortwave radiation fluxes.

Due to the reflection and scattering by the day-side clouds, the stellar radiation flux reaching the planet's surface is unevenly distributed around the substellar point (Figures 1(e)–(h)). Its maximum, located broadly to the west of the substellar longitude, reaches values as high as 438 W m⁻² in the ROCKE-3D case, while it is only 296 W m⁻² in the UM case. The arrow-shaped distribution of clouds in the ROCKE-3D simulation (Figure 8(c)) shields the equatorial region from the stellar radiation, instead allowing the flux to reach its maximum
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Hab 1 results: time variability of (a)–(c) surface temperature (K), (d) ice-free ocean area fraction, (e)–(h) cloud water path (kg m⁻²), and (i)–(l) the altitude of the cloud mass mixing ratio (MMR) maximum (km). Each column shows ((a), (d), (e), (i)) global mean, ((b), (f), (j)) day-side mean, ((c), (g), (k)) night-side mean, and ((h), (l)) terminator mean. Thin lines show the raw data, while thick lines show a five-orbit rolling mean.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Hab 1 results: time-mean surface temperature (K) in (a) ExoCAM, (b) LMD-G, (c) ROCKE-3D, and (d) the UM.</figcaption></figure>

at ≈45° latitude. This cloud pattern, and the radiation flux pattern as a consequence, is due to the dynamical regime in ROCKE-3D being distinctly different from that in the other 3 GCMs (see Section 3.1.3). Namely, the ROCKE-3D simulation is in a rapid-rotator regime, in which the meridional circulation and convective activity sustaining the cloud cover is confined to a narrow equatorial band, similar to the simulations of a planet orbiting a late M-dwarf in Komacek & Abbot (2019). In terms of the global average of the surface stellar flux, GCMs fall into two groups: those with a flux of about 75 W m⁻² (ExoCAM and UM) and those with a flux of about 50% higher (LMD-G and ROCKE-3D). This dichotomy is explained by the overall low amount of clouds in the LMD-G case and the concentrated distribution of clouds in the equatorial band in the ROCKE-3D case (see Section 3.1.2). The ice-free ocean surface, occupying almost a half of the day-side area (Figure 2(d)), absorbs 94% of the available stellar radiation, while the “ice-covered” periphery of the day side absorbs 75% (due to the fixed albedo; see Section 2.1).

Our models predict mean global surface temperatures of ≈232–246 K (Table 4, Figures 2(a) and 3), but the intermodel differences in temperature do not simply follow the differences between model planetary albedos. While ExoCAM, LMD-G, and ROCKE-3D have similar estimates of the mean temperature as well as the temperature extrema, the UM produces a consistently colder climate (Figures 3 and 4(a)). Because this difference is mostly due to the night-side temperatures of the surface below the stationary cyclonic mid-latitude gyres (Figure 3(d)), the day–night contrast in the UM is larger by ≈20 K than that in the other three models. The surface temperature in ExoCAM is somewhat higher than that in LMD-G, which contradicts previous studies in which LMD-G was much warmer than CAM-family models (Yang et al. 2019). The opposite outcome found here is not likely due to model heritage issues, because the model configurations changed minimally since that intercomparison was conducted. Rather, the opposite outcome is likely a function of natural model behaviors arising from two different circulation regimes. Yang et al.’s (2019) intercomparison studied slow rotating planets (60 day periods)
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Surface area fraction with corresponding time-mean temperature binned every 10 K in (a) Hab 1 and (b) Hab 2.</figcaption></figure>

receiving the same stellar insolation as modern Earth but from a 3400 K star with a blackbody radiation spectrum. Here, TRAPPIST-1e receives less flux but from a much redder star, and perhaps most importantly, its dramatically faster rotation rate shifts the atmospheric circulation into a distinct dynamical state (e.g., Haqq-Misra et al. 2018; Sergeev et al. 2020) compared to the planets studied by Yang et al. (2019). Therefore it is perhaps unsurprising that the complex coupling between atmospheric circulation, cloud, and radiative processes yields nonlinear model sensitivities when making changes to fundamental planetary configuration. This finding is another incentive to perform model intercomparisons for other planetary configurations, as will be done within the CUISINES framework (Fauchez et al. 2021).

The day-side surface temperature is, on the other hand, relatively similar in all models (Figure 3). This is perhaps unsurprising given that the day side is under a constant forcing from the stellar radiation flux, while the night-side temperature is susceptible to the changes in global circulation and intermodel differences in subgrid parameterizations. The day-side mean surface temperature ranges from 261 K in the UM to 268 K in ExoCAM, the latter being the warmest simulation overall (Figure 2(b)). ExoCAM-simulated climate is in fact the warmest at almost every height at the substellar point (Figure 5(a)) and in the troposphere on average (Figure 6(a)). As Figure 5(a) shows, LMD-G appears to have an even warmer boundary layer at the substellar point than that in ExoCAM and thus a steep lapse rate of about ≈8 K km⁻¹ below 750 hPa. The steepest lapse rate, exceeding ≈25 K km⁻¹ at the substellar point is produced by the UM, whose boundary-layer parameterization is based on a nonlocal closure and is able to maintain this temperature gradient very close to the surface. Additionally, the UM is the only THAI GCM that has a nonhydrostatic dynamical core (Turbet et al. 2021), which

allows for explicit convection and is another likely reason for the steep lapse rate in this model. Above the surface layer, all four models have a relatively similar profile of temperature, with the lapse rate fluctuating close to the moist adiabatic lapse rate—between 7 K km⁻¹ and 4 K km⁻¹—throughout most of the troposphere (Figure 5(b)).

The key distinction between the results from the GCMs is seen at the level of the tropopause (≈60 hPa), where the lapse rate becomes positive (Figure 5(b)). The lowest temperature minimum of <160 K is predicted by LMD-G, while in the UM it is ≈170 K; in ExoCAM and ROCKE-3D, the tropopause has the least pronounced cold trap with a temperature minimum of ≈180 K (Figure 5(a)). Above ≈50 hPa, the temperature begins to rise—gradually in the ExoCAM, ROCKE-3D, and LMD-G simulations and more sharply in the UM—outlining a stratosphere of TRAPPIST-1e. The temperature increase in the upper layers is driven by the shortwave absorption, which heats the atmosphere at a rate of up to 15 K d⁻¹ at the substellar point (Figure 5(c)). This extreme heating rate is present mostly at the topmost level of the GCMs, and so the upper level spacing and the model top height are what partially explain the intermodel differences. Compared to the dry benchmark case Ben 1 analyzed in Turbet et al. (2021),<fnref n="21" /> shortwave heating rates here are “bottom-heavy,” meaning that their magnitude is comparable, if not higher, than that at the TOA. This is due to intense absorption of stellar radiation by water vapor concentrated in the lowest half of the atmosphere (≈100 hPa and below).

The excess of energy received by the day side is redistributed vertically in the troposphere by convection and laterally to the night side by the global atmospheric circulation (see Section 3.1.3). Temperature profiles at the substellar point (Figure 5(a)) and at its antipode on the night side (not shown) reveal the low spatial variation of the free tropospheric temperature, confirming the weak temperature gradient assumption (similar to earlier studies, e.g., Turbet et al. 2016; Boutle et al. 2017). This can be seen in a relatively low gradient of the tropospheric temperature, whose meridional-mean varies from ≈230 K on the day side to ≈260 K on the night side, with ROCKE-3D having the smallest variation (Figure 6(a)). The weak variation of tropospheric temperature in ROCKE-3D with longitude is a consequence of a more zonally symmetric global circulation, which is characteristic of synchronously rotating planets in a “fast rotation” dynamical regime (Haqq-Misra et al. 2018; Sergeev et al. 2020). This regime is dominated by two extratropical zonal jets (see Section 3.1.3) and a stronger equator-to-pole temperature gradient relative to the substellar-antistellar gradient.

Vertically integrated specific humidity, i.e., the water vapor path, is the lowest in ROCKE-3D (up to 12 kg m⁻²), while the highest amount of water vapor (up to almost 25 kg m⁻²) is simulated by ExoCAM and LMD-G (Figure 6(b)). In the case of LMD-G, however, moisture is mostly concentrated in the lowest layers close to the substellar point and decreases rapidly with height, making the stratosphere more than an order of magnitude drier compared to the other three models (see last column in Table 4). As a result, the absorption of shortwave radiation by the water vapor in the stratosphere is the lowest in the LMD-G

<footnote n="21">²¹ Ben 1 and Ben 2 are dry benchmark cases equivalent to the Hab 1 and Hab 2 cases, but have no moisture. Ben 1 has a 1 bar N₂-dominated atmosphere with 400 ppm of CO₂, while Ben 2 has a 1 bar CO₂-dominated atmosphere.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Hab 1 results: substellar vertical profiles of (a) air temperature (K), (b) lapse rate (K km⁻¹), (c) SW heating rate (K d⁻¹), and (d) LW heating rate (K d⁻¹). Time-mean values are shown by solid lines, the 1σ deviation is shown by shading.</figcaption></figure>

(Figure 5(c)). The low humidity in LMD-G can be explained by the temperature minimum at the tropopause (Figure 5(a)), which controls the amount of water entering the stratosphere in the convective region. In addition, the amount of water vapor that reaches the upper atmosphere is strongly affected by the type of a convection parameterization in the GCM. As mentioned in Section 2.2.2, LMD-G uses the convective adjustment scheme, which works by adjusting the temperature and humidity on each level toward a pre-defined adiabat. This scheme does not allow for convective overshooting or entrainment, and so does not

affect neutral or stable layers. ExoCAM, ROCKE-3D, and the UM, on the other hand, use the mass-flux parameterizations. These schemes carry information in the vertical from their original perturbation, modifying the entraining and detraining rates along the way. If the convective plume is particularly warm and moist when it reaches the tropopause, it can travel a few model levels above it before becoming neutrally buoyant. During this overshooting, warm and humid air from the plume is heavily detrained into the lower stratosphere. As a result, the mass-flux schemes tend to deposit more heat and moisture in the
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Hab 1 results: meridional-mean profiles of (a) the mass-weighted tropospheric mean of air temperature (K), (b) water vapor path (kg m⁻²), (c) cloud water path (kg m⁻²): total (solid lines), liquid (dashed lines), and ice (dotted lines), and (d) cloud area fraction (%). Here, the troposphere's boundary is defined by the lapse rate decreasing below 2 K km⁻¹ in the upper atmosphere.</figcaption></figure>

upper levels of the model where the adjustment schemes would stop.

### 3.1.2. Clouds

Confirming previous studies of cloudy tidally locked planets (e.g., Joshi 2003; Kopparapu et al. 2016; Turbet et al. 2016; Boutle et al. 2017; Wolf et al. 2017; Turbet et al. 2018; Fauchez et al. 2019; Komacek & Abbot 2019; Yang et al. 2019; Eager et al. 2020), all of our Hab 1 simulations have a deep cloud layer covering the substellar point, demonstrated by the cloud area fraction maximum in Figure 6(d). The cloud fraction varies significantly between the four GCMs (Figures 7(i)–(l)). The LMD-G is the least cloudy, having a mean cloud fraction barely above 10% throughout most of the atmosphere (Figure 7(j)). Its global mean cloud water path, i.e., the vertical integral of cloud condensate, is also the lowest among the models at a value of 3.3 × 10⁻² kg m⁻² (Figure 8(b)), which is supported by previous intercomparisons where this model was likewise less cloudy than CAM and other GCMs (e.g., Abbot et al. 2012; Yang et al. 2019). ROCKE-3D, on the other hand, has the highest total cloud fraction (77%), presenting TRAPPIST-1e as a planet completely enshrouded in clouds (Figure 8(g)). Note that the total cloud fraction in ROCKE-3D appears much higher than in the meridional-mean cross section in Figures 6(d) and 7(c), which may be an artifact of how the model integrates the cloud field vertically.

On the day side, the lowest part of the cloud deck consists of liquid water droplets, while above ≈600 hPa and farther away from the substellar point, the clouds are dominated by ice crystals (Figures 7(a)–(h)). The total cloud water path is

partitioned differently in the four GCMs: ExoCAM's and LMD-G's clouds almost entirely comprise liquid, while ROCKE-3D's clouds are mostly icy, and the UM is roughly half-way between these two extremes. Despite the high cloud ice content, the longwave CRE in ROCKE-3D is not as high as in ExoCAM, because in the former, the cloud water path is concentrated at the substellar point (Figure 8(c)), exposing the rest of the day side to the loss of thermal radiation. The eastward tilt of the cloudy area in the meridional-mean cross sections is the consequence of a superrotating equatorial jet and leads to the difference in cloud content at the western and eastern terminators (for more details, see Paper III; Fauchez et al. 2022). Illustrating the process of detrainment of convective cloud condensate, an anvil cloud is visible at the top of the day-side atmosphere in all four GCMs, and is especially pronounced in ExoCAM and the UM (Figures 7(i), (l)). This anvil is composed of ice crystals (Figures 7(a)–(d)) and is the key contributor to the longwave CRE warming the day side, which is indeed the largest in ExoCAM and the UM (Table 4). LMD-G exhibits a weak cloud anvil and a low-altitude maximum in the icy clouds, which explains the small shortwave and longwave CRE in this model mentioned above.

The largest discrepancy in clouds between our models is on the night side of TRAPPIST-1e. In the meridional average, the total cloud fraction ranges from ≈15% (LMD-G) to ≈70% (ROCKE-3D). Clouds cover the night side differently: they occur mostly either in the high latitudes (in LMD-G), at the equator (in the UM), or in both regions (in ExoCAM and ROCKE-3D). They occupy the lowest portion of the troposphere, having the maximum cloud content at 1–3 km above the surface, although the ROCKE-3D case is an outlier that has
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Hab 1 results: meridional-mean cross sections of (a)–(d) mass mixing ratio (MMR) of ice cloud particles (log[kg kg⁻¹]), (e)–(h) MMR of liquid cloud particles (log[kg kg⁻¹]), and (i)–(l) total cloud fraction (%).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Hab 1 results: (a)–(d) cloud water path (vertically integrated cloud condensate content, kg m⁻²), and (e)–(h) cloud fraction (%).</figcaption></figure>

cloud maximum at ≈8 km above the surface (Figure 2(k)). Note that even though in ROCKE-3D the cloud content peaks at ≈8 km, there are plenty of low-level clouds on the night side too (Figure 7(c)). The night side is generally characterized by low-level cloudiness, because convection (especially deep convection) is suppressed, and most night-side clouds are of stratiform type, formed due to the advection and condensation of water vapor from the day side of the planet. Their formation is largely controlled by the relative humidity, and thus temperature. Due to a different circulation regime (see Section 3.1.3), the temperature distribution in ROCKE-3D has a local minimum in mid-troposphere in high latitudes, just below the eastward jet cores (Figure 10(c)). This results in a

local maximum of the relative humidity (not shown), leading to more cloud formation.

From the perspective of the total cloud content and its composition, the ExoCAM case is a distinct outlier, because its cloud water path on the night side is consistently several times higher than that in the other three GCMs (Figures 2(g) and (c)) and because its night-side clouds are dominated by liquid water (Figure 7(e)). Due to the absence of the shortwave radiation on the night side, the water phase of clouds (and thus their optical properties) is less important than their altitude for the net night-side CRE: for ExoCAM, LMD-G, and the UM, it is negative, while for ROCKE-3D, it is positive but close to zero (Table 4). The altitude and optical properties of clouds at the terminators
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Helmholtz decomposition of the horizontal wind at 250 hPa in the Hab 1 case (quivers): (a)–(d) total wind, (e)–(h) zonal mean rotational component, (i)–(l) eddy rotational component, and (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>

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Hab 1 results: vertical cross sections of the zonal mean zonal wind (m s⁻¹). The gray dashed horizontal line marks the 250 hPa level shown in Figure 9.</figcaption></figure>

affect the transmission spectrum of the atmosphere, and so it is
crucial to simulate correctly in GCMs, because transmission
spectroscopy is the primary mode of atmospheric characteriza-
tion of TRAPPIST-1e. This is discussed in more detail in
Part III (Fauchez et al. 2022).

3.1.3. Atmospheric Circulation

The stationary pattern of the global circulation in the
troposphere of tidally locked planets is a combination of
several key components. This is elucidated by Figure 9, where
the Hab 1 wind field at 250 hPa is split into its rotational and
divergent components, using the Helmholtz decomposition
methodology of Hammond & Lewis (2021). First, on a
relatively fast-rotating planet such as TRAPPIST-1e, the zonal
(eastward) component of the wind forms prograde jets, evident
in the zonal mean of the rotational wind field (Figures 9(e)–(h)
and 10). Second, its eddy component (here, the deviation from
the zonal mean) reveals planetary-scale stationary waves

straddling the substellar longitude (Figures 9(i)–(l)). Third, the divergent wind is relatively weak at this altitude, but clearly shows the outflow branch of the overturning circulation. At the substellar point, strong divergence aloft is compensated by the convergence in the boundary layer (not shown) and is associated with intense convection (positive upward wind velocity contours in Figure 9). Note that even though the magnitude of the divergent wind is small, it is one of the dominant branches of the moist static energy transport between the day and night side on tidally locked planets (Hammond & Lewis 2021).

The substellar-antistellar overturning circulation is driven by
the divergent wind and, due to its isotropic structure, is best
presented in a tidally locked coordinate system (Koll &
Abbot 2015). In this coordinate system, the tidally locked
latitude is −90 degrees at the antistellar point, 0 degrees at the
north and south Poles, and 90 degrees at the substellar point.
The tidally locked longitude is 0 on the line connecting the
substellar point, North Pole, and antistellar point, increasing
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Hab 1 results: overturning circulation shown by the tidally locked mass streamfunction $\Psi^{TL}$ (kg s$^{-1}$). The mass flux is clockwise where $\Psi^{TL}$ is positive and counterclockwise where $\Psi^{TL}$ is negative. The maximum of $\Psi^{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>

toward the (original) eastern terminator. The mass streamfunction $\Psi^{TL}$ in the tidally locked latitude-pressure plane shows a single circulation cell between the day side and the night side: the air rises at the substellar point, moves toward the antistellar point aloft, subsides on the night side or at the terminators and finally returns to the day side in the lower troposphere (Figure 11).

ExoCAM, LMD-G, and the UM show a classic wind field pattern in the free atmosphere, consisting of a strong superrotating jet at the equator, stationary gyres in mid-latitudes, and a moderate divergence at the substellar point (Matsuno 1966; Gill 1980). Such a pattern emerges in many simulations of moderate climates of tidally locked exoplanets in the other 3D GCMs and idealized models (e.g., Carone et al. 2014, 2015, 2016; Haqq-Misra et al. 2018; Hammond et al. 2020). ROCKE-3D, however, produces a different pattern: strong zonal jets in mid-latitudes and a weak superrotation at the equator, reminiscent of the regime of rapidly rotating Earth-like planets (Figure 9(g)). This pattern is associated with a weaker stationary wave response (Figure 9(k)), but a stronger substellar divergence than in other GCMs (Figure 9(o)). Correspondingly, the tidally locked mass streamfunction $\Psi^{TL}$ reaches $7.5 \times 10^{11}$ kg s$^{-1}$ in ROCKE-3D (Figure 11(c)). In the other three models, $\Psi^{TL} = 3.9-5.8 \times 10^{11}$ kg s$^{-1}$, while the overturning cell extends beyond the terminator to the night side (Figures 11(a), (b), and (d)).

The outlier case of ROCKE-3D provides more evidence of the delicate balance in atmospheric circulation of tidally locked but relatively fast-rotating planets such as TRAPPIST-1e, as hinted at by earlier studies (e.g., Edson et al. 2011; Kopparapu et al. 2017). The global circulation on TRAPPIST-1e can settle on two different regimes: with a dominant eastward jet at the equator (Figures 10(a), (b), and (d)) or a pair of mid-latitude eastward jets (Figure 10(c)). This confirms the findings of Sergeev et al. (2020), who demonstrated that this transition is susceptible to change in the parameterization of convection in the UM. Given the differences in parameterizations between the THAI GCMs, it is unsurprising to capture both circulation regimes in this intercomparison, though it is reassuring to do it with GCMs other than the UM. A similar dichotomy between circulation regimes can exist for fast-rotating tidally locked planets, assuming both a dry and wet atmosphere, as was first shown by Edson et al. (2011) using the GENESIS GCM (a fork of the NCAR CESM). The authors found that for an Earth-like case, an abrupt transition occurs at rotation periods, especially in their wet case (3–4 days), lower than for TRAPPIST-1e (6.1 days; see Table 1). It is thus possible that due to intermodel

differences, only ROCKE-3D is susceptible to the circulation change at the rotation period of 6.1 days, while other models would change their circulation if the rotation period is lowered even more. Another reason for the discrepancy between the studies is that Edson et al. (2011) used the stellar spectrum and irradiation of the Sun. However, the study by Carone et al. (2014, 2015, 2016) used the MIT GCM to explore a broad parameter space of synchronously rotating Earth-like atmospheres; this model configuration included a present-Earth Newtonian relaxation scheme as radiative forcing, which showed this dynamical regime transition to occur at a period of 6–7 days. Carone et al. (2014, 2015, 2016) attributed this discrepancy with Edson et al. (2011) to differences in each model’s treatment of the night-side forcing. Just as in the sensitivity experiments of Edson et al. (2011), Carone et al. (2014, 2015, 2016), and Sergeev et al. (2020), the two circulation regimes are associated with different cloud patterns, which is indeed seen in ROCKE-3D’s cloud fraction map (Figures 8(c) and (g)). The nature of the two distinct global circulation regimes for planets like TRAPPIST-1e epitomizes the importance of multimodel studies and is an important avenue of research (Sergeev et al. 2022).

Even more intermodel differences are visible in the upper layers of the atmosphere, as shown in Figure 10. At $\approx 50$ hPa, the flow weakens and in two models (ExoCAM and the UM) becomes retrograde, while farther above the prograde flow, it increases again. This eastward flow takes the shape of two jets near the poles in the UM, while in other models, it is concentrated at the equator and is especially strong in ExoCAM ($>100$ m s$^{-1}$). The reason why the four GCMs exhibit large differences in the upper atmosphere is the different altitude of the model top, as well as different radiative heating rates (Figure 5). Another possible factor is the numerical damping at the top of the domain, which is implemented differently in each model (Turbet et al. 2021).

### 3.1.4. Time Variability

Fluctuations around the time mean seen in Figure 2 are the manifestation of weather patterns growing and moving in the simulated atmosphere of TRAPPIST-1e. It is clear that different GCMs have a different amplitude and frequency of the weather variability within the Hab 1 climate. The largest amplitude of fluctuations is exhibited by ExoCAM, and its principal oscillation period is 12.5 orbits for the global mean cloud content (Figure 12(a)). A secondary peak of variability amplitude can be discerned at about three orbits, similar to ROCKE-3D. LMD-G’s variability is overall similar to
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Time spectrum of the global mean cloud water path variability. This is the time series shown in Figures 2(e) and 14(e) transformed to a time frequency domain in log-log axes for (a) Hab 1 and (b) Hab 2, respectively. The time frequency units are TRAPPIST-1e orbits (see Table 1). The dashed vertical lines mark the maximum amplitude for each of the models. Periods corresponding to the maximum amplitude are shown in the upper-right corner.</figcaption></figure>

ExoCAM, albeit at a lower amplitude, while ROCKE-3D and the UM exhibit the lowest period of cloud content fluctuations and a moderate amplitude. In the UM, the most significant mode of variability is high-frequency fluctuations at 1.1 orbits, though secondary peaks are seen at lower frequency.

The temperature variability in the Hab 1 simulations is also the largest in ExoCAM (Figure 2(a)). Its amplitude is the largest on the night side of the planet (Figure 2(c)) and is partially canceled by the day-side variability (Figure 2(b)). Preliminary analysis of this variability reveals that it is driven by periodic shifts of the planetary-scale atmospheric waves, resulting in variable heat and water vapor transport to the night side of the planet. The surface heat budget for the night side is controlled only by the advection of warm and moist air from the day side (e.g., Lewis et al. 2018), because the ocean heat transport is absent from our simulations. As the night-side atmosphere becomes warmer and moister, downward fluxes of sensible heat and net LW radiation on the night side increase, raising the surface temperature to ≈235 K (Figure 2(c)). In the opposite phase of the oscillation, the air becomes drier and colder, allowing the night-side surface temperature to drop to ≈215 K.

The period of this global-scale oscillation is about 12.5 TRAPPIST-1e orbits, or approximately 76 Earth days. This rather long period as well as a preference toward positive zonal wavenumbers (not shown) makes the ExoCAM variability somewhat similar to the Madden-Julian Oscillations (MJO) observed on Earth (see, e.g., Madden & Julian 1971; Kiladis et al. 2009; Matthews 2021). Faithful representation of the MJO remains a challenge even for sophisticated climate and weather models (Vallis & Penn 2020). In particular, components of the convective parameterization such as entrainment rate and convection triggering affect the initiation and propagation of the MJO (e.g., Holloway et al. 2013). It is thus possible that a combination of moist physics parameterizations and the overall warmer climate in ExoCAM's simulation is able to excite such internal variability, while in other models (e.g., ROCKE-3D and the UM) it is not. Precisely why this happens in ExoCAM is left for the future, because the amount of additional experiments and outputs required to tackle this question is beyond the scope of the THAI project.

Our time-variability analysis is limited by the 610 Earth-day sample size of the THAI project. Therefore, the presence of any

significant variability on Earth-yearly scales, such as the Quasi-Biennial Oscillation, are missed from our analysis. On the other hand, the spectra in all four GCMs peak at 12.5 orbits or less, making the THAI sample able to capture at least some of the important modes of variability. Time oscillations with shorter periods are also more likely to be observed because of the limitations in telescope time.

## 3.2. Global Climate in Hab 2 Simulations

### 3.2.1. Radiation Fluxes and Thermodynamic Profiles

At the top of the CO₂-dominated atmosphere of Hab 2, the absorbed shortwave radiation is larger than that for Hab 1 by 2–32 W m⁻², but the overall pattern is the same (not shown). This increase is mostly due to stronger absorption by CO₂ and H₂O and, in all models except LMD-G, due to a weaker CRE<sub>SW</sub> (Table 4). In LMD-G, CRE<sub>SW</sub> reaches −31.8 W m⁻² and is the second largest after that in ROCKE-3D (−35.0 W m⁻²). Accordingly, the Bond albedo α<sub>p</sub> of TRAPPIST-1e is the highest in these two GCMs compared to ExoCAM and the UM, though it is still lower than that in Hab 1 (Table 4).

The longwave part of CRE in Hab 2 is somewhat smaller compared to that in Hab 1, but is likewise the lowest in LMD-G and the highest in ExoCAM (Table 4). The spatial pattern of the TOA OLR is broadly similar to that in Hab 1 (Figures 13(a)–(d)). However, the contrast between the day side and the night side is reduced dramatically. What ExoCAM illustrates in particular, is that the OLR has two sharp peaks in its longitudinal distribution, thus having an imprint in the emission spectra. This would also be apparent in broadband spectra for which the minimum would be shifted eastward and the maximum would be closer to 90° (not shown).

The net effect of clouds on the TOA radiation fluxes is still negative, though much weaker relative to that in Hab 1, especially in the cases of ExoCAM and the UM (Table 4). The exception here is LMD-G, which, due to its stronger shortwave part, has a net CRE of −28.6 W m⁻². Thus, in three out of four models, clouds play a lesser role in controlling energy fluxes in the climate of a CO₂-dominated atmosphere (Hab 2) than that in an N₂-dominated atmosphere (Hab 1); in LMD-G though, the trend is the opposite.

The distribution of shortwave radiation at the surface of the planet reaches its maximum to the west of the substellar point
<figure><img src="image_17.png" /><figcaption>**Figure 13.** Key radiation fluxes in Hab 2 simulations. Panels (a)–(d): TOA OLR (W m⁻²); panels (e)–(h): downward shortwave radiation flux at the surface (W m⁻²); and panels (i)–(l): net upward longwave radiation flux at the surface (W m⁻²).</figcaption></figure>

(Figures 13(e)–(h)). This flux is the lowest in ROCKE-3D and the highest in the UM, with their respective distributions nevertheless being similar to those in Hab 1. The average incident shortwave radiation is similar in all four models, between 55 and 60 W m⁻² and almost twice as low as in Hab 1 due to larger absorption of radiation by the atmosphere. Note that the bulk of the shortwave radiation is absorbed by the oceanic surface with a relatively low albedo of 0.06, because most of the day side is above the freezing point of water, which is the threshold for the albedo change in the THAI protocol (Figure 14(b)). While the net downward surface shortwave radiation flux is positive, its longwave counterpart is positive only in the Hab 2 cases of ROCKE-3D and UM. In ExoCAM and especially in LMD-G, the net upward longwave flux at the surface is negative at almost all grid points (Figure 13(j)), implying that the atmosphere is radiating heat to the surface—in stark contrast to the Hab 1 simulations (Figures 1(i)–(l)) and the Hab 2 simulations of ROCKE-3D and the UM (Figures 13(k) and (l)). This is likely due to the differences in the gaseous absorption in the RT schemes. ROCKE-3D and the UM employ the same scheme and the same set of spectral files (see Section 2.2.1), thus simulating overall similar (weaker) atmospheric absorption (Figure 15(c)), while ExoCAM and LMD-G use schemes that tend to overestimate shortwave absorption, especially by H₂O in the case of ExoCAM.

At the substellar point, the surface temperature exceeds 300 K in all four models, but the range of its maximum values is wider than that in the Hab 1 case: from 301 K in LMD-G to 319 K in ExoCAM (Figure 16). The ExoCAM simulation is the warmest throughout all of the troposphere at the substellar point,<fnref n="22" /> before being overtaken above by ROCKE-3D and the

UM (Figure 15(a)). The transition between the troposphere and stratosphere is located slightly higher and is less pronounced in Hab 2 than in Hab 1. In the upper atmosphere, the heat balance is primarily between the absorption of stellar radiation and emission of planetary radiation, both of which are an order of magnitude larger than those in Hab 1 (Figures 15(c) and (d)). Due to strong absorption of stellar flux by CO₂ in the Hab 2 case, which is overwhelming the absorption by water vapor in the troposphere, the radiative heating profile is larger at higher altitudes, just as in the benchmark simulations (Turbet et al. 2021). Compared to the benchmark simulations, both Hab 1 and Hab 2 cases see larger differences in the radiative heating profiles between the GCMs (even those employing the same RT code, as is the case for ROCKE-3D and the UM), which is due to the radiative effects of water vapor and clouds.

Due to heavily enhanced greenhouse effect of the CO₂ atmosphere (the *G* parameter in Table 4), the global surface temperatures of the Hab 2 cases are 40–50 K higher than those of Hab 1 (Figure 16). This is mostly due to the significant warming of the night side of the planet (Figure 14(c)). The intermodel spread in predicted surface temperature is larger than that in the Hab 1 simulations, slightly more on the day side and slightly less on the night side (see Figures 2(b) and (c) and Figures 14(b) and (c)). This spread reduces somewhat when the whole troposphere is concerned, which is substantially warmer than that for the Hab 1 case (Figure 17(a)). ExoCAM’s troposphere is warmer than those of the other models and thus contains a markedly larger amount of water vapor (Figure 17(b)). The water vapor path is an order of magnitude larger in Hab 2 than in Hab 1 due to the Clausius-Clapeyron law.

### 3.2.2. Clouds

Despite the climate being much warmer, the Hab 2 simulations bear many similarities to those of Hab 1 in terms of cloud patterns. The bulk of cloud water is on the day side of the planet in all four GCMs (Figure 17(c)), having an average

<footnote n="22">²² The ExoCAM THAI simulations use an old version of the RT code, which overestimates the shortwave absorption by H₂O and the longwave absorption by CO₂ for dense atmospheres. The new THAI simulations were not ready to be included in the intercomparison in time, but preliminary results show that with the updated CO₂ RT, ExoCAM gets a mean surface temperature of ≈282 K and as such is no longer an outlier in Hab 2.</footnote>
<figure><img src="image_18.png" /><figcaption>**Figure 14.** Hab 2 results: time variability of (a)–(c) surface temperature (K), (d) ice-free ocean area fraction, (e)–(h) cloud water path (kg m⁻²), and (i)–(l) the altitude of the cloud content maximum (km). Each column shows ((a), (d), (e), (i)) global mean, ((b), (f), (j)) day-side mean, ((c), (g), (k)) night-side mean, and ((h), (l)) terminator mean. Thin lines show the raw data, and thick lines show a five-orbit rolling mean.</figcaption></figure>

day-side cloud water path of 0.1–0.3 kg m⁻² (Figures 14(f) and (g)). However, vertically integrated cloud diagnostics, such as the cloud water path or the cloud area fraction, show that TRAPPIST-1e is cloudier in the Hab 2 case than in the Hab 1 case, especially on the night side of the planet (Figures 17(c) and (d) and 14(g)).

For Hab 2 simulations, the clouds are more extended in the vertical, which is visible in the meridional-mean cross sections (Figure 18). In fact, nonnegligible cloud condensate extends to the top of the model domain in two GCMs: ExoCAM and LMD-G. From an observational perspective, simulated trans-mission spectra are affected differently when the upper atmosphere is cloudy, so using ExoCAM or LMD-G may contradict the interpretations obtained from ROCKE-3D or the UM (see Paper III; Fauchez et al. 2022). The global mean altitude of the cloud content maximum is also higher in the Hab 2 case relative to that in Hab 1 in three out of four models: ExoCAM (most pronounced), LMD-G, and the UM (Figure 14(i)). While ROCKE-3D predicts the most amount of clouds among the models (Figures 14(e) and 19(c) and (g)), its clouds are closer to the surface, especially liquid clouds on the night side (Figure 18(g)). The second cloudiest simulation is by ExoCAM, where the cloud content (mostly liquid droplets) reaches its maximum at the highest altitude (11 km) among the four GCMs (Figure 14(i)). Interestingly, the minimum of cloudiness in ExoCAM is not on the night side as in other models, but in the western half of the day hemisphere (e.g., Figures 18(a), (e), and (i)), which may have a detectable imprint in the reflected starlight (Wolf et al. 2019).

### 3.2.3. Atmospheric Circulation

The tropospheric wind structure in all four Hab 2 simulations is similar to that in their respective Hab 1 simulations (cf. Figures 9(a)–(d) and 20(a)–(d)). ExoCAM, LMD-G, and the UM develop the equatorial jet pattern, while ROCKE-3D develops two extratropical jets (Figure 21). Accordingly, the day-side vertical velocity maximum is more pronounced at the equator in the ROCKE-3D simulation and spread out across latitudes in the other three simulations (color shading in Figure 20).

Figure 22 shows that the strength of the day-night circulation does not change significantly compared to that of the Hab 1 case. The mass streamfunction $Ψ^{TL}$ in each GCM has broadly the same shape in both the Hab 1 and Hab 2 experiments. Namely, $Ψ^{TL}$ in LMD-G is bottom-heavy, showing a strong flow toward the substellar point in the boundary layer, while in the UM, it is top-heavy, showing strong divergence in the upper atmosphere. As is the case for Hab 1, the ROCKE-3D’s overturning cell in the Hab 2 simulations is confined to the day side, while in the other three GCMs, it extends to the night side (as in, e.g., Zhang & Yang 2020; Hammond & Lewis 2021).

In the troposphere, the differences in the mean eastward wind speed between the Hab 1 and Hab 2 cases are rather small, with the maximum value staying within 50 to 70 m s⁻¹ (Figure 21). The stratospheric circulation, however, differs more noticeably: in all four models, the uppermost levels do not develop an additional region of superrotation (not shown). Instead, they are dominated by a westward flow, which is especially strong in the UM (Figure 21(d)). As discussed in
<figure><img src="image_19.png" /></figure>

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

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

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

<figure><img src="image_23.png" /><figcaption>**Figure 15.** Hab 2 results: substellar vertical profiles of (a) air temperature (K), (b) lapse rate (K km⁻¹), (c) SW heating rate (K d⁻¹), and (d) LW heating rate (K d⁻¹). Time-mean values are shown by solid lines, the 1σ deviation is shown by shading.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 16.** Hab 2 results: time-mean surface temperature (K) in (a) ExoCAM, (b) LMD-G, (c) ROCKE-3D, and (d) the UM.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>**Figure 17.** Hab 2 results: meridional-mean profiles of (a) the mass-weighted tropospheric mean of air temperature (K), (b) water vapor path (kg m⁻²), (c) cloud water path (kg m⁻²): total (solid lines), liquid (dashed lines), and ice (dotted lines), and (d) cloud area fraction (%). Here, the troposphere's boundary is defined by the lapse rate decreasing below 2 K km⁻¹ in the upper atmosphere. Note the different y-axis limits compared to Figure 6.</figcaption></figure>

more detail in Part I (Turbet et al. 2021), the source of intermodel differences near the top lid is due to the differences in clear-sky RT, location and spacing of vertical levels, and numerical damping (sponge layers). For example, ROCKE-3D has relatively high numerical damping coefficients to increase model stability, which results in a low wind speed near the top boundary compared to that in the other three GCMs (Figure 21(c)). The impact on the stratospheric circulation by the numerical drag has been discussed in Carone et al. (2018), but this question would further benefit from a dedicated GCM intercomparison. Thus, the overall characteristic of the time-mean zonal flow for the Hab 2 simulations is that a single structure of the equatorial (ExoCAM, LMD-G, the UM) or extratropical (ROCKE-3D) jet extends throughout most of the lower atmosphere.

### 3.2.4. Time Variability

Compared to Hab 1 case, the warmer climate of the Hab 2 case sees an increase of the vapor and condensed water at the terminators (Figure 14), which consequently leads to an increased amplitude of the cloud content variability (Figure 12(b)). The periodicity of these fluctuations, however, stays roughly the same for each of the GCMs, with relatively longer periods in ExoCAM and LMD-G and shorter periods in ROCKE-3D and the UM. This is also true for the global surface temperature fluctuations (spectra not shown, but compare Figure 2(a) to Figure 14(a)), though for LMD-G and the UM, they exhibit longer periods than those for the cloud content. Unlike the cloud variability, the amplitude of

temperature variability in the ExoCAM simulation of Hab 2 drops by more than an order of magnitude compared to that of Hab 1, which was noted in the preliminary THAI results (Fauchez et al. 2020). However, it is not the case for the other GCMs: the amplitude of temperature fluctuations increases in the UM and LMD-G and stays the same in ROCKE-3D (not shown).

## 4. Synthesis and Discussion

### 4.1. Inter-model Similarities and Differences

To a large degree, all four THAI GCMs predict temperatures within 20–25 K of each other in both the Hab 1 and Hab 2 cases. On the largest scale, the atmosphere stays habitable and sustains a stable water cycle without entering a moist greenhouse state or an atmospheric collapse due to the CO₂ condensation (see Section 4.2). While the intermodel spread in the global mean surface temperature is nonnegligible, it is smaller in the Hab 1 case (14 K) and larger in Hab 2 (24 K, see Table 4). These differences are two and four times larger than those found in the Ben 1 and Ben 2 simulations, respectively (Turbet et al. 2021), demonstrating the amplification of intermodel spread by moist physics in GCMs. Our results suggest that while in the dry climates the intermodel spread is similar (Turbet et al. 2021), in moist climates it is larger for a warmer case—at least for the four THAI GCMs. The intermodel spread is smaller on the day side (especially in the Hab 1 case) and larger on the night side of the planet. Interestingly, the largest temperature difference in the Hab 2
<figure><img src="image_26.png" /><figcaption>**Figure 18.** Hab 2 results: meridional-mean cross sections of (a)–(d) mass mixing ratio (MMR) of ice cloud particles (log[kg kg⁻¹]), (e)–(h) MMR of liquid cloud particles (log[kg kg⁻¹]), and (i)–(l) total cloud fraction (%).</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 19.** Hab 2 results: (a)–(d) cloud water path (vertically integrated cloud condensate content, kg m⁻²), and (e)–(h) cloud fraction (%).</figcaption></figure>

simulations is between ExoCAM and LMD-G. These two models were at the extremes in a previous exoplanet GCM intercomparison (Yang et al. 2019), although LMD-G was warmer than ExoCAM, whereas in our study it is the opposite.

Global mean energy fluxes at the TOA generally agree well (for example, the model spread in the TOA flux is below 7%). Radiative heating rates in the atmosphere are also rather similar between the models in their respective cases; the models disagree substantially only at the very top of the atmosphere. The tropospheric circulation in all models except for ROCKE-3D is in a Rhines-rotator regime with a strong equatorial superrotation. The overall strength of the circulation is comparable between the models and the two cases, with the

mean zonal winds reaching 50–70 m s⁻¹ and the day–night overturning streamfunction $\Psi_{\text{max}}^{\text{TL}}$ staying within 4.2–7.4 × 10¹¹ kg s⁻¹. Thus, the agreement between the models is relatively better in the dry thermodynamics and radiation, as well as in the atmospheric circulation than it is for moist physics. While there are certainly important intermodel differences in radiation and dynamics, as we discuss in more detail in Part I of the THAI trilogy (Turbet et al. 2021), we find a much larger discrepancy in the moist aspects of the simulated Hab 1 and Hab 2 climates.

Our experiments show that the largest intermodel differences appear in the water vapor content, cloud amount, and distribution. The wettest troposphere is simulated by ExoCAM
<figure><img src="image_28.png" /><figcaption>**Figure 20.** The same as in Figure 9, but for Hab 2.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Figure 21.** Hab 2 results: vertical cross sections of the zonal mean zonal wind (m s⁻¹). The gray dashed horizontal line marks the 250 hPa level shown in Figure 20.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Figure 22.** The same as in Figure 11, but for Hab 2.</figcaption></figure>

both in the Hab 1 and Hab 2 cases, which can be attributed to its overall warmer climate. In the stratosphere, the clear outlier is LMD-G, predicting specific humidity at least an order of magnitude lower than that in the other three models (VMR₁ hPa in Table 2). This outcome contradicts previous studies, in which ExoCAM and LMD-G were used to simulate fast-rotating planets orbiting a G-dwarf star (Wolf & Toon 2015) and tidally locked planets orbiting an M-dwarf (Yang et al. 2019). As mentioned in Section 3.1.1, the main causes of

LMD-G's dry stratosphere are the cold tropopause and the inability of the convective adjustment scheme to account for moisture detrainment beyond the level of neutral buoyancy at the tropopause. Additional evidence for the latter is the scant amount of clouds at high altitudes (in a shape of an “anvil”) above the substellar point (see Section 3.1.2). As a result, the water absorption band at ≈6 µm is much weaker in LMD-G than that in the simulations by other models (see Paper III; Fauchez et al. 2022).
Not only does LMD-G have the driest upper atmosphere, its cloud amount is consistently the smallest among THAI models. This contradicts the Earth-like simulations of Wolf & Toon (2015), but generally supports the more recent simulations of Yang et al. (2019), when compared to ExoCAM (then called CAM4_Wolf). The cloudiest simulation, in terms of the global mean cloud fraction and integrated cloud condensate, is produced by ROCKE-3D. Similar to the UM simulations of TRAPPIST-1e in Sergeev et al. (2020), increased cloud fraction as well as the cloud condensate concentrated at the substellar point, are congruent with a change in the tropospheric circulation regime. Unlike the other three models, ROCKE-3D favors the accumulation of eastward momentum in two zonal jets in the extratropics instead of the equator. Such a regime of the two extratropical jets was shown by Turbet et al. (2018) for a case of an N₂ atmosphere with 10% CH₄ in their 3D simulations of TRAPPIST-1e, f, h, and g. This suggests that a change in the gaseous absorption and a resulting change in the static stability of the atmosphere may induce a change in circulation, which, in turn, modulates the cloud cover.

Cloud metrics shown in this paper also demonstrate that one must be cautious when interpreting GCM output, as has been also noted by Kopparapu et al. (2017). Cloud area fraction can give a general impression of how “cloudy” a planet is and thus how reflective its atmosphere is, i.e., the albedo. At the same time, this parameter is highly model-dependent. Moreover, it is usually insensitive to the optical thickness or the altitude of the cloud and can give a false impression of its radiative importance. One solution is to filter optically thin clouds (Boutle et al. 2017). Another is to show vertically integrated cloud content, i.e., the cloud water path (e.g., Kopparapu et al. 2017; Wolf et al. 2017, 2019). The difference between these two diagnostics can be seen in Figures 8 and 19. For example, in the ROCKE-3D simulation of the Hab 1 case, almost all grid boxes have a cloud fraction >50% (Figure 8(g)), but the cloud water path is quite low in the tropical band on the night side of the planet (Figure 8(c)), indicating that that region is covered by optically thin high clouds (Figures 7(c), (g), and (k)). This demonstrates how looking at only one cloud diagnostic in a GCM can be misleading and should be taken into account when connecting GCM output to atmospheric retrieval tools.

Clouds are a serious obstacle to atmospheric characterization, because they raise the altitude of the continuum level of a transmission spectrum (e.g., Fauchez et al. 2019), therefore reducing the amplitude of the absorption lines. Even though ROCKE-3D is cloudier overall than other GCMs in both the Hab 1 and Hab 2 cases, its cloud condensate is concentrated at the substellar point (Figures 8(c) and 19(c)). At the terminators, the largest cloud water path is simulated by ExoCAM. Furthermore, the cloud mixing ratio maximum in ExoCAM reaches its maximum at the highest altitude among the THAI models, especially at the eastern terminator, where the cloud content is higher. Consequently, its continuum level is the highest and, compared to other models, requires a higher number of transits to detect CO₂ absorption peaks using JWST (see Paper III; Fauchez et al. 2022). The UM, on the other hand, tends to produce lower cloud decks at the terminators with smaller temporal variability than that in the ExoCAM, which means it has the lowest continuum level and requires the smallest number of orbits for CO₂ detection. Accurate representation of cloud opacity and cloud top heights, including their spatial variability, has a substantial impact on the

detectability estimates, making 3D GCMs a more appropriate tool than 1D models with a fixed cloud layer (e.g., Lin et al. 2021). The intermodel differences in cloudiness of the terminator stem from the differences in how GCMs distinguish between convective and stratiform clouds. For example, if one GCM produces a lot of convective clouds, which would be concentrated at the substellar point, and few stratiform clouds, the terminator region in its simulation would be less cloudy. While this separation between cloud types is artificial, it is common for terrestrial GCMs and thus should be investigated further in future intercomparisons.

A key potentially observable feature of cloudy dynamic atmospheres is the time variability of opacity at the terminators. This variability will most likely stay below the noise level of JWST's NIRSpec PRISM instrument and not affect retrieved abundances (May et al. 2021). While our simulations support May et al.'s (2021) finding, they provide a multimodel estimate and may prove useful if a similar planet is found closer to Earth or if TRAPPIST-1e is probed by better instruments (Fauchez et al. 2022). As discussed in Section 3.1.4 and 3.2.4, ExoCAM has the highest amplitude of cloud variability among the THAI models, and this corresponds to the highest variability in the atmospheric transit depth relative to other GCMs, as shown for different wavelengths of the transmission spectrum in Paper III (Fauchez et al. 2022). The periodicity of ExoCAM's cloud variability is approximately 12.5 orbits, which is comparable to the minimum number of transits required to detect CO₂ features and thus may need to be taken into account if TRAPPIST-1e is observed during successive transits. Differences between the simulations of the Hab 1 and Hab 2 climate lead to an increased amplitude of the cloud variability, which affects ExoCAM and the UM the most (Fauchez et al. 2022). Note, however, that for other models, the differences between the variability in the spectra between the Hab 1 and Hab 2 cases are muted, likely because the highest clouds have a bigger impact than the total cloud condensate in the column.

## 4.2. Habitability and Climate Stability

Differences in the mean climate shown in Section 3 imply that the prediction of habitability of TRAPPIST-1e is model-dependent. With respect to the inner edge of the habitable zone, none of the THAI simulations enter the water-loss stage or a runaway greenhouse state. Stratospheric water content in the Hab 1 case, taken here as the H₂O volume mixing ratio at the 1 hPa, is about 5–50 × 10⁻⁸ mol mol⁻¹ (Table 4) and stays comfortably below the classical water-loss limit of 3 × 10⁻³ (Kasting et al. 1993), especially in the case of LMD-G. The intermodel difference in the stratospheric water content is an order of magnitude (due to LMD-G having a much lower value), which for a warmer climate, e.g., a planet orbiting the star closely such as TRAPPIST-1d, may result in the eventual loss of water on geological timescales. While the atmosphere is substantially warmer and thus moister in the Hab 2 simulations relative to Hab 1, the stratospheric volume mixing ratio for the Hab 2 case is still one to two orders of magnitude lower than the limit—similar to the fast-rotating simulations for a 2600 K star demonstrated by Kopparapu et al. (2017). Note, however, that the Kasting et al. (1993) limit is a conservative estimate derived from a cloudless 1D Earth simulation.

The Hab 1 and Hab 2 simulations also avoid the other extreme, namely the CO₂ condensation limit, which can trigger a runaway atmospheric collapse (Turbet et al. 2018). In the
coldest simulation—the UM’s simulation of the Hab 1 case—
the night side exhibits a strong temperature inversion in the
boundary layer, exceeding 30 K in the lowest 100 m. Even in
the coldest places on Earth, the inversion rarely reaches this
strength, counteracted by turbulent mixing by topography-
induced mesoscale eddies (Joshi et al. 2020). Our simulations
have flat topography and therefore potentially overestimate the
surface cooling in the night-side cold traps. Even with this
overestimation, the lowest temperature stays above the
condensation point of CO₂, both for the Hab 1 (≈125 K) and
in Hab 2 cases (≈195 K). Note that the Hab 1 surface
temperature minimum in the UM is 20–30 K lower than that
in the other three GCMs (see Figures 3(d) and 4(a)); so in a
colder state, it would be closer to a potential atmospheric
collapse, all other factors being equal.

Within these extremes demarcating the habitable zone boundaries, habitability of the temperate climate of TRAPPIST-1e can be defined by the abundance of stable surface liquid water. We can use the planet-average open ocean fraction as a metric of this abundance, with the obvious caveat of the lack of a dynamic ocean in our simulations (Del Genio et al. 2019). Figure 2(d) demonstrates a fairly close estimate of the ice-free ocean area at around 0.20–0.23 in the GCMs. Despite the fact that the UM predicts a lower surface temperature both globally and on the day side, its open ocean fraction is essentially the same as that in ExoCAM (0.23). LMD-G and ROCKE-3D have a similar day-side average surface temperature, and their ocean fraction fluctuates around 0.20.

While being relatively close to entering the water-loss limit, Hab 2 simulations predict a larger fraction of deglaciated habitable area at the surface: the surface is warmer and the day- night contrast is roughly twice as small as in Hab 1. The ice- free fraction varies from 0.36 in LMD-G to ≈0.75 in the UM and ROCKE-3D to 1 in ExoCAM (Figure 14(d)). Thus, in line with a larger intermodel spread in the global mean surface temperature, the ice-free fraction disagreement in the Hab 2 case is larger than that in Hab 1. This is a result of the surface temperature distribution with respect to the freezing point of water being different in the two experiments (Figure 4). In the Hab 1 simulations, the freezing point is closer to the right tail of this distribution, where the agreement between the models is higher; for Hab 2, it is closer to the left tail of the distribution (i.e., the night-side temperatures), where models have the largest discrepancy.

4.3. *Caveats and Future Work*

We have used four state-of-the-art 3D GCMs to re-examine
the climate of TRAPPIST-1e in a more structured way than
previous studies. Nevertheless, our study has some caveats.
First, we chose a simplified gas composition for both the Hab 1
and Hab 2 cases: an N₂- and CO₂-dominated atmosphere with
H₂O as the condensible species. Furthermore, we chose a total
atmospheric pressure of 1 bar. While these setups cover two
key types of atmospheric composition, not least for the sake of
easier comparison with the climate dynamics on Earth,
TRAPPIST-1e may possess a completely different atmosphere
(Turbet et al. 2020) and more 3D GCM experiments are needed
to widen the parameter space (e.g., Fauchez et al. 2019; Zhang
& Yang 2020).

Second, the Hab 1 and Hab 2 experiments presented here are
predicated on the assumption of abundant water at the surface

of the planet. It is possible that TRAPPIST-1e has lost all water
throughout its lifetime, but some studies suggest that the
hydrodynamic escape may have stopped once it entered the
habitable zone (Bolmont et al. 2017; Bourrier et al. 2017).
Additionally, late-stage outgassing of the mantle may have
replenished water on the surface (Moore & Cowan 2020).
TRAPPIST-1e may never enter a water-loss phase and retain its
oceans over a timescale of 100 Gyr (Bolmont et al. 2017;
Kopparapu et al. 2017). Thus, together with the dry planet
scenarios explored in Paper I (Turbet et al. 2021), scenarios
with abundant water shown in the present paper cover the two
main pathways for rocky planets orbiting M-dwarfs (Tian &
Ida 2015).

Third, the surface is assumed to be completely covered with a static slab ocean without any continents. The presence of continents and orography may significantly change the planet's hydrological cycle (e.g., Lewis et al. 2018) and boundary-layer mixing on the night side (Joshi et al. 2020). In some cases and in the absence of large continents, the oceanic heat transport may also change the simulated climate (e.g., Del Genio et al. 2019), but the oceanic component of exoplanet models requires a separate intercomparison (Fauchez et al. 2021).

Fourth, it is unknown which model is closer to the truth,
since there are no observations of the atmosphere of
TRAPPIST-1e available yet. However, the present study and
similar future multimodel intercomparisons will help to identify
biases of one GCM relative to others, whether in a specific
scenario or systematically across the parameter space. Before
observations of temperate rocky exoplanets become available,
high-resolution cloud-resolving models or large eddy simula-
tions could be used to validate moist physics parameterizations
of exoplanet GCMs (e.g., Sergeev et al. 2020; Lefèvre et al.
2021). More detailed intercomparisons of individual parame-
terizations of convection, boundary layer, and microphysics are
also needed.

Future studies of temperate climates from a 3D perspective
are encouraged to build upon the experiments analyzed in this
paper. It is pertinent to address the obvious caveats mentioned
above, as well as expand this effort both in the parameter space
and by comparing more GCMs. More habitable scenarios
should be explored by varying the atmospheric composition
and pressure of TRAPPIST-1e. This has been started in a few
previous works (e.g., Wolf et al. 2017; Turbet et al. 2018;
Fauchez et al. 2019), but requires a coordinated multimodel
approach, several benefits of which have been demonstrated by
the present study. This type of study would bring the best value
if performed periodically (e.g., every 5 yr), similar to the CMIP
framework (e.g., Eyring et al. 2016). We note, however, that
future intercomparisons of exoplanet GCMs should find a
balance between increasing the model complexity and reducing
the intermodel spread.

5. Conclusions

In this Part II of the THAI study, we analyzed a moist
climate of TRAPPIST-1e in an aquaplanet configuration
simulated by the four 3D GCMs—ExoCAM, LMD-G,
ROCKE-3D, and the UM. We hope this study will propel
future model intercomparisons of potentially habitable extra-
solar worlds. Our conclusions are as follows.

1. Despite many intermodel differences, Hab 1 and Hab 2 climates in all four models stay within the habitability
limits. They sustain an active hydrological cycle without entering the water-loss stage: the stratospheric water content is at most $O(10^{-7} \text{mol mol}^{-1})$ in the Hab 1 case and $O(10^{-5} \text{mol mol}^{-1})$ in the Hab 2 case—below the classical water-loss limit. The lowest surface temperature is also above the CO$_2$ condensation limit, and thus the condition for night-side atmospheric collapse is not triggered.

2. The global mean surface temperature is 232–246 K in the Hab 1 simulations and 271–295 K in the Hab 2 simulations, with a relatively larger intermodel spread on the night side than on the day side. The models agree more in their estimate of the ice-free area fraction for Hab 1 (0.20–0.23) and less for Hab 2 (0.36–1.00).

3. The largest intermodel differences exist in the variables directly affected by moist physics parameterizations. LMD-G produces the driest stratosphere, especially in the Hab 1 case due to a “cold trap” at the tropopause. ExoCAM has the wettest troposphere, due to it being substantially warmer, which is a consequence of strong stellar flux absorption by CO$_2$. Globally averaged, LMD-G is also consistently less cloudy, both in terms of mean cloud area fraction and column-integrated cloud content. Note that this is not the case at the terminators (for further discussion, see Part III; Fauchez et al. 2022). The cloudiest Hab 1 and Hab 2 climates are predicted by ROCKE-3D, which simulates TRAPPIST-1e enveloped by clouds with a high concentration of cloud water at the substellar point.

4. Tropospheric circulation in three models (ExoCAM, LMD-G, and the UM) is characterized by a Rhines-rotator regime with a strong equatorial superrotation both in the Hab 1 and Hab 2 cases. ROCKE-3D is an outlier and predicts a zonal circulation dominated by two extratropical jets. A different circulation regime in ROCKE-3D partly explains the differences in cloud patterns relative to other GCMs.

5. The models exhibit different amplitude and periodicity of atmospheric variability. At the terminators, ExoCAM has a distinctly higher amplitude of cloud variability than other models, and a period of 12.5 orbits. In other models, the variability is generally of higher frequency (≈2.5–11 orbits). The amplitude grows in Hab 2 relative to Hab 1, especially in the ExoCAM and UM cases.

The authors are grateful to the two anonymous reviewers whose comments helped to improve this paper. D.E.S., I.A.B., F.H.L., J.M., and N.J.M. acknowledge use of the Monsoon2 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).

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. 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. was granted access to the High-Performance Computing (HPC) resources of Centre Informatique National de l’Enseignement Supérieur (CINES) under the allocations Nos. 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. acknowledges support from the PORTAL BRAIN-be 2.0 BELSPO project. M.T. and F.F. thank the LMD Generic Global Climate team for the teamwork development and improvement of the 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. J.H.M. acknowledges funding from the NASA Habitable Worlds program under award 80NSSC20K0230.

The authors acknowledge the help of Andrew Ackerman to set up the cloud diagnostics in ROCKE-3D.

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 the Hab 1 and Hab 2 experiments resulted in roughly ≈1280 kg kgCO$_2$e, including ≈520 kg from ExoCAM runs, ≈80 kg from LMD-G, ≈650 kg from ROCKE-3D, and ≈30 kg from the UM. These numbers are the product of the following three factors: (i) CPU-hours used by each of the model to reach the steady state, (ii) kWh per CPU-hour, which depends on the efficiency of the supercomputer used, and (iii) kgCO$_2$e per kWh of energy, which varies a lot from country to country and so is the biggest source of uncertainty in our estimates.

**Software:** Scripts to process and visualize THAI data are available on GitHub: https://github.com/projectcuisines/thai_trilogy_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), and xarray (Hoyer & Hamman 2017). ExoCAM (Wolf & Toon 2015) is available on GitHub: https://github.com/storyofthewolf/ExoCAM. LMD-G is available upon request from Martin Turbet (martin.turbet@lmd.jussieu.fr) and François Forget (francois.forget@lmd.jussieu.fr). ROCKE-3D is public domain software and available for download from https://simplex.giss.nasa.gov/gcm/ROCKE-3D/. Annual tutorials for new users take place annually, whose recordings are freely available on line at https://www.youtube.com/user/NASAGISStv/playlists?view=50&sort=dd&shelf_id=15. The Met Office Unified Model is available for use under licence; see http://www.metoffice.gov.uk/research/modeling-systems/unified-model.
## Data accessibility

All of our GCM THAI data are permanently available for download here: https://ckan.emac.gsfc.nasa.gov/organization/thai, with variables described for each data set. If you use those 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.”

## ORCID iDs

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

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

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

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

## References

Abbot, D. S. 2014, JCLI, 27, 4391

Abbot, D. S., Voigt, A., Branson, M., et al. 2012, GeoRL, 39, L20711

Anglada-Escudé, G., Amado, P. J., Barnes, J., et al. 2016, Natur, 536, 437

Barker, H. W., Cole, J. N. S., Morcrette, J. H., et al. 2008, QJRMS, 134, 1463

Bolmont, E., Selsis, F., Owen, J. E., et al. 2017, MNRAS, 464, 3728

Boucher, O., Le Treut, H., & Baker, M. B. 1995, JGR, 100, 16,395

Bourrier, V., de Wit, J., Bolmont, E., et al. 2017, AJ, 154, 121

Boutle, I. A., Eyre, J. E. J., & Lock, A. P. 2014, MWRv, 142, 1655

Boutle, I. A., Mayne, N. J., Drummond, B., et al. 2017, A&A, 601, A120

Carone, L., Keppens, R., & Decin, L. 2014, MNRAS, 445, 930

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

Ceppi, P., Brient, F., Zelinka, M. D., & Hartmann, D. L. 2017, Wiley Interdisciplinary Reviews: Climate Change, 8, e465

Charnay, B., Forget, F., Wordsworth, R., et al. 2013, JGRD, 118, 414

Dawson, A. 2016, JORS, 4, 31

Del Genio, A. D., Way, M. J., Amundsen, D. S., et al. 2019, AsBio, 19, 99

Eager, J. K., Reichelt, D. J., Mayne, N. J., et al. 2020, A&A, 639, A99

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

Edwards, J. M., Havemann, S., Thelen, J.-C., & Baran, A. J. 2007, AtmRe, 83, 19

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

Genio, A. D. D., Yao, M.-S., Kovari, W., & Lo, K. K.-W. 1996, JCLI, 9, 270

Gill, A. E. 1980, QJRMS, 106, 447

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

Goldblatt, C., McDonald, V. L., & McCusker, K. E. 2021, NatGe, 14, 143

Gregory, D. 1995, MWRv, 123, 2716

Gregory, D., & Rowntree, P. R. 1990, MWRv, 118, 1483

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

Hack, J. J. 1994, JGR, 99, 5551

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

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

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

Henderson, D. S., L'Ecuyer, T., Stephens, G., Partain, P., & Sekiguchi, M. 2013, JAPMC, 52, 853

Holloway, C. E., Woolnough, S. J., & Lister, G. M. S. 2013, JAtS, 70, 1342

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

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

Joshi, M. 2003, AsBio, 3, 415

Joshi, M. M., Elvidge, A. D., Wordsworth, R., & Sergeev, D. 2020, ApJL, 892, L33

Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, Icar, 101, 108

Kiladis, G. N., Wheeler, M. C., Haertel, P. T., Straub, K. H., & Roundy, P. E. 2009, RvGeo, 47, RG2003

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. Schmid (Netherlands: IOS Press), 87

Koll, D. D. B., & Abbot, D. S. 2015, ApJ, 802, 21

Komacek, T. D., & Abbot, D. S. 2019, ApJ, 871, 20

Kopparapu, R. K., Ramirez, R., Kasting, J. F., et al. 2013, ApJ, 765, 131

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

Kopparapu, R. K., Wolf, E. T., Haqq-Misra, J., et al. 2016, ApJ, 819, 84

Le Trent, H., & Li, Z.-X. 1991, CIDy, 5, 175

Lefèvre, M., Turbet, M., & Pierrehumbert, R. 2021, ApJ, 913, 101

Lewis, N. T., Lambert, F. H., Boutle, I. A., et al. 2018, ApJ, 854, 171

Lin, Z., MacDonald, R. J., Kaltengegger, L., & Wilson, D. J. 2021, MNRAS, 505, 3562

Liu, Q., Collins, M., Maher, P., Thomson, S. I., & Vallis, G. K. 2021, GMD, 14, 2801

Lovis, C., Snellen, I., Mouillet, D., et al. 2017, A&A, 599, A16

Madden, R. A., & Julian, P. R. 1971, JAtS, 28, 702

Manabe, S., & Wetherald, R. T. 1967, JAtS, 24, 241

Manners, J., Edwards, J. M., Hill, P., & Thelen, J.-C. 2022, SOCRATES (Suite Of Community Radiative Transfer codes based on Edwards and Slingo) Technical Guide, Met Office, UK, https://code.metoffice.gov.uk/trac/socrates

Matsuno, T. 1966, JMeSJ, 44, 25

Matthews, A. J. 2021, QJRMS, 147, 4310

May, E. M., Taylor, J., Komacek, T. D., Line, M. R., & Parmentier, V. 2021, ApJL, 911, L30

Met Office 2021, Iris: A Python library for analysing and visualising meteorological and oceanographic data sets, http://scitools.org.uk/

Moore, K., & Cowan, N. B. 2020, MNRAS, 496, 3786

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, Natl. Cent. for Atmos, 2009

Pincus, R., Barker, H. W., & Morcrette, J. J. 2003, JGRD, 108, 4376

Schmidt, G. A., Kelley, M., Nazarenko, L., et al. 2014, JAMES, 6, 141

Segelstein, D. R. 1981, M.S. Thesis, Univ. of Missouri

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

Sergeev, D. E., Lewis, N. T., Lambert, F. H., et al. 2022, PSJ, 3, 214

Sergeev, D. E., & Zamyatina, M. 2021, Aeolus—a Python library for the analysis and visualisation of climate model output, doi:10.5281/zenodo.5145604

Sherwood, S. C., Bony, S., & Dufresne, J.-L. 2014, Natur, 505, 37

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

Snellen, I., de Kok, R., Birkby, J. L., et al. 2015, A&A, 576, A59

Sundqvist, H. 1978, QJRMS, 104, 677

Tian, F., & Ida, S. 2015, NatGe, 8, 177

Tolkien, J. R. R. 1954, The Two Towers (London: George Allen & Unwin)

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., Fauchez, T. J., Sergeev, D. E., et al. 2022, PSJ, 3, 211

Turbet, M., Gillmann, C., Forget, F., et al. 2020, Icar, 335, 113419

Turbet, M., Leconte, J., Selsis, F., et al. 2016, A&A, 596, A112

Vallis, G. K., & Penn, J. 2020, QJRMS, 146, 2297

Wall, C. J., Hartmann, D. L., & Norris, J. R. 2019, GeoRL, 46, 12495

Walters, D., Baran, A. J., Boutle, I., et al. 2019, GMD, 12, 1909

Warren, S. G. 1984, ApOpt, 23, 1206

Warren, S. G., & Brandt, R. E. 2008, JGR, 113, D14220

Way, M. J., Aleinov, I., Amundsen, D. S., et al. 2017, ApJS, 231, 12

Webb, M. J., Lock, A. P., Bretherton, C. S., et al. 2015, RSPTA, 373, 20140414

Wilson, D. R., & Ballard, S. P. 1999, QJRMS, 125, 1607
Wilson, D. R., Bushell, A. C., Kerr-Munslow, A. M., Price, J. D., & Morcrette, C. J. 2008, *QJRMS*, **134**, 2093

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

Wolf, E. T., Kopparapu, R. K., & Haqq-Misra, J. 2019, *ApJ*, **877**, 35

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*, **7**, 656

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

Wordsworth, R. D. 2016, *AREPS*, **44**, 381

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

Yang, J., Cowan, N. B., & Abbot, D. S. 2013, *ApJ*, **771**, L45

Yang, J., Leconte, J., Wolf, E. T., et al. 2016, *ApJ*, **826**, 222

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

Yang, J., Liu, Y., Hu, Y., & Abbot, D. S. 2014, *ApJL*, **796**, L22

Zelinka, M. D., Myers, T. A., McCoy, D. T., et al. 2020, *GeoRL*, **47**, e85782

Zhang, G., & McFarlane, N. A. 1995, *AtO*, **33**, 407

Zhang, Y., & Yang, J. 2020, *ApJ*, **901**, L36