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

In this part 1 of the THAI trilogy, we present results of the dry atmosphere simulations, named Ben 1 and Ben 2 in the THAI protocol (Fauchez et al. 2020). The Ben 1 and Ben 2 cases were designed to act as benchmarks for the two habitable states Hab 1 and Hab 2, respectively, of TRAPPIST-1e, described in part 2 (Sergeev et al. 2022a). Ben 1 is the colder state with a 1-bar N₂-dominated atmosphere plus 400 ppm of CO₂, while Ben 2 is a warmer state with a 1-bar CO₂-dominated atmosphere.

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

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

• 1.  
  The Exoplanet Community Atmospheric Model (ExoCAM), a branch of the National Center for Atmospheric Research Community Earth System Model version 1.2.1.
• 2.  
  The Laboratoire de Météorologie Dynamique—Generic model (LMD-Generic, aka LMD-G).
• 3.  
  The Resolving Orbital and Climate Keys of Earth and Extraterrestrial Environments with Dynamics (ROCKE-3D Planet_1.0).
• 4.  
  The Met Office Unified Model (UM, science version GA7.0, code version 11.6).

In the following subsections, we provide a brief description of the model components related to the dynamics and dry physics (see part 2 of the THAI intercomparison for details on the moist physics), with a focus on RT, large-scale dynamics (computed using so-called "dynamical cores"), and subgrid-scale dynamics (dry convection and turbulent mixing). Table 2 gives an overview of the dynamical cores, as well as the horizontal and vertical resolutions used in the four 3D models for all THAI simulations, including the Hab 1 and Hab 2 cases analyzed in Sergeev et al. (2022a). Table 3 summarizes the dry physics parameterizations used by the models.

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

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

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

Both the UM and ROCKE-3D share a common RT module SOCRATES, ¹⁶ which is an open-source two-stream correlated-k RT model provided by the Met Office (Edwards & Slingo 1996; Manners et al. 2012). SOCRATES is highly flexible. Spectral files, created for both longwave and shortwave streams separately, are tailored to specific atmospheric compositions and stellar spectra. The spectral files are then input to the GCMs at run-time, containing all information necessary to define the RT problem, including the number of spectral intervals; the absorbing gases, continua, and CIAs; the number of Gauss points used for each gas; the designation of dominant and minor absorbing species in each interval; and cloud and Rayleigh scattering properties. Thus, dramatic changes to the RT can be achieved by careful preprocessing of the spectral files, with no changes to the hard-code being required. When run embedded in a GCM, SOCRATES uses the equivalent extinction method to treat multiple overlapping gas species (Amundsen et al. 2017). In the equivalent extinction method the dominant absorbing gas, as determined from the total transmissivity in a test column, is treated with a full k-distribution, while all minor absorbing gases are treated as gray. The number of Gauss points used for the dominant gas in each spectral interval is determined based on a selected transmission error tolerance (usually 0.001). For the Ben 1 simulations, spectral files were tailored to Earth-like atmospheric compositions, allowing for small amounts of CO₂ in N₂-dominated atmospheres. The spectral files used by ROCKE-3D and the UM are identical (details provided in Table 3). Gas absorption coefficients are derived from the HITRAN 2012 spectroscopic database (Rothman et al. 2013). For Ben 2 simulations, CO₂ absorption is treated, following Wordsworth et al. (2010), with sub-Lorentzian line shape (Perrin & Hartmann 1989) out to 500 cm⁻¹ with CO₂–CO₂ CIA overlaid.

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

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

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

The ROCKE-3D dynamical core utilizes finite differences where atmospheric velocity points are located on the Arakawa B grid (Arakawa 1972; Hansen et al. 1983). The atmospheric vertical layers use a sigma coordinate system from the surface to 150 hPa. From there constant-pressure layers are used to the TOA (≈0.1 hPa). For the THAI experiments the model was run with 72 × 46 grid points, with each grid point being 5° in longitude by 4° in latitude except at the poles, where the latitude coordinate is only 2° in extent. The smaller polar grid boxes are necessary to maintain second-order accuracy of air-mass fluxes, the pressure gradient force, and momentum advection terms (Schmidt et al. 2006). The dynamical time steps used in the Ben 1 and Ben 2 simulations were 225 and 450 s, respectively. Both Hab 1 and Hab 2 simulations use a dynamical time step of 450 s. For all four simulations the parameterized physics uses a time step of ≈30 minutes (1756.8 s to be exact), and less in some submodules, while the SOCRATES radiation was called every two physics time steps (∼1 hr).

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

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

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

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

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

To represent mixing of adiabatically conserved heat, momentum, and tracers throughout the troposphere, including the PBL, the UM uses a first-order turbulence closure of Lock et al. (2000) with modifications described in Lock (2001) and Brown et al. (2008). The scheme treats unstable and stable boundary layers differently. For unstable layers, there are two separate profiles of diffusion coefficients, one for the surface sources of turbulence and one for upper atmospheric sources (e.g., radiative cooling). Initially diagnosed by an adiabatic ascent, the extent of these profiles is adjusted so that the magnitude of the buoyancy consumption of turbulent kinetic energy is limited to a specified fraction of buoyancy production. In moist aquaplanet simulations (see part 2), this permits the cloud layer to decouple from the surface. Additional nonlocal fluxes of heat and momentum can generate more vertically uniform potential temperature and wind profiles in unstable (convective) boundary layers. The scheme is tightly coupled to the convection scheme (see part 2), but the latter is inactive in dry cloudless simulations of the Ben 1 and Ben 2 cases. For stable boundary layers and the free troposphere, diffusion coefficients are calculated using the parameterization of Smith (1990), based on a local Richardson number. The mixing can also be enhanced if the nonlocal coefficients, calculated as described above, are larger than the local ones, which can happen in the neutral boundary layer. The stability dependence is given by the "MES-tail" function within the stable boundary layer and by the "sharp" function at 200 m and above (Brown et al. 2008). Given the resulting coefficient, the diffusion equation is solved implicitly, and the kinetic energy dissipated through the turbulent shear stresses is converted to a local heating term.

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

Zoom In Zoom Out Reset image size

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

The differences between the outputs obtained with our set of GCMs are very small and comparatively similar between the Ben 1 and Ben 2 cases. Similar results are also observed in the shortwave heating rate vertical profiles (see Figures 2(c) and 3(c)), which display similar trends from the surface up to ≈1 hPa. The Ben 2 simulations result in a much stronger shortwave heating rate due to direct absorption of stellar light by CO₂. Both the Ben 1 and Ben 2 cases exhibit shortwave heating rates that increase as altitude increases (or pressure decreases). This is because CO₂ absorption saturates as we go deeper into the atmosphere (because the upper layers absorb all photons in the core of near-infrared CO₂ band lines).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

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

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

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

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

As was mentioned in the previous section, the GCMs agree well on the temperature profile at the substellar point. Subtle differences are revealed by looking at the lapse rate (i.e., the vertical temperature gradient), shown in Figures 2(d) and 3(d). In the troposphere, the lapse rates are about −9 and −11 K km⁻¹ on average for the Ben 1 and Ben 2 cases, respectively. To first order, the difference between the Ben 1 and Ben 2 simulations comes from the variation of the atmospheric heat capacity between the cases (about 1040 and 850 J kg⁻¹ K⁻¹ for Ben 1 and Ben 2, respectively). Note that because there are no changes in the atmospheric composition in our simulations, all four GCMs have the heat capacity set to a fixed parameter in each of the cases.

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

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

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

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

Figure 5 shows the surface temperature maps and the winds at the 250 hPa isobaric surface for the Ben 1 (panel (a)) and Ben 2 (panel (b)) cases. First, we can see that the intermodel spread in the global mean surface temperature is 7 K for the Ben 1 simulations and 6 K for the Ben 2 simulations (see also Table 4). In both cases ExoCAM is on average the warmest and the UM the coldest, similar to what is found for the Hab 1 case (Sergeev et al. 2022a). For the Ben 1 case the day−night thermal contrast is the largest for ROCKE-3D, followed by the UM, LMD-G, and ExoCAM. In the case of ROCKE-3D, two wind gyres form in the western hemisphere at high latitudes, in the coldest regions (see Figure 5(b), bottom left plot). Those gyres are a manifestation of a Rhines rotator dynamical regime, characteristic for slow-rotating tidally locked exoplanet simulations (e.g., Carone et al. 2015; Noda et al. 2017; Haqq-Misra et al. 2018; Hammond et al. 2020; Wang & Yang 2021); this is a regime where a strong upper atmosphere superrotation is combined with a strong upwelling beneath the substellar point. The other three GCM simulations produce a fast-rotator circulation regime: without the gyres, but with two extratropical zonal jets. This is best illustrated by the lower part of Figure 6, which shows the vertical cross sections of the zonal mean zonal wind. Similar to the findings of Sergeev et al. (2020) for TRAPPIST-1e, the thermal contrast between the dayside and the nightside is smaller in this regime than in the slow-rotator regime. This is important not only for TRAPPIST-1e but generally for any tidally locked rocky exoplanet with short enough orbital periods, i.e., short enough for the Rossby deformation radius to be comparable to the planet's radius (Carone et al. 2015, 2016; Haqq-Misra et al. 2018; Penn & Vallis 2018; Pierrehumbert & Hammond 2019; Zhang 2020).

Zoom In Zoom Out Reset image size

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

Zoom In Zoom Out Reset image size

Ben 1 and Ben 2 simulations both exhibit small temporal variability, as illustrated for surface temperatures in Figure 8, due most likely to the manifestation of weather patterns. It has indeed already been shown that on synchronous planets (with a fixed solar forcing) a temporal variability is nevertheless expected (Joshi et al. 1997; Merlis & Schneider 2010; Ding & Wordsworth 2021). The four GCMs show a similar small level of global mean surface temperature high-frequency variability (∼1 K) for Ben 1 and Ben 2 simulations. However, on a period of 10 consecutive orbits, the frequency of temperature oscillations varies from one model to another. The ExoCAM simulations show much greater nightside (Figure 8, right) and dayside (Figure 8, middle) variability than the other GCMs, but these opposite-phase oscillations offset each other on the global mean (Figure 8, left). This behavior, which is also pronounced in the ExoCAM Hab 1 and Hab 2 simulations (part 2; Sergeev et al. 2022a), results from periodic weather systems traveling around the planet. Unfortunately, the design of the Ben 1 and Ben 2 cases of the THAI protocol (Fauchez et al. 2020)—which limited the output files to 10 orbits— does not allow us to probe the low-frequency temporal variability, but there are several indications that it exists, as suggested by rolling averages in Figure 8. This is also indicated by the presence of north–south asymmetries in the different mean values of the models (see, e.g., Figures 5 and 6), somewhat similar to the results of Noda et al. (2017), Section 4.3. There do not appear to be clear correlations between the different model variabilities and the land properties adopted in the models (especially the number and depth of layers; see Table 3, rightmost column) over the time range explored from the THAI simulations. We leave the question of why weather systems seem to be more pronounced on ExoCAM than other GCMs for future studies. We also recall that the level of temporal variability of the different models is far below the detectability threshold of JWST and the astronomical observatories of the coming decades (part 3; Fauchez et al. 2022).

Zoom In Zoom Out Reset image size

The fact that the models reside in different rotational regimes is likely because the radius and rotation period of TRAPPIST-1e place it just at the edge between fast and Rhines (intermediate) rotation regimes (Carone et al. 2015, 2016; Haqq-Misra et al. 2018), as already identified in Sergeev et al. (2020) and Carone et al. (2018) specifically for TRAPPIST-1e. As a matter of fact, Sergeev et al. (2020) showed that just a change in the parameterization of convection can drive a GCM into one or the other dynamical regime. Unlike the simulations in Edson et al. (2011), initial conditions are unlikely to play a role here, as all simulations start from isothermal conditions of 300 K and null winds. We also identified in preliminary THAI simulations (of the UM model) that subtle changes (here, of correlated-k tables) are able to settle the GCM to a different rotation regime. We also note that a similar behavior appears in the aquaplanet THAI simulations, with ROCKE-3D being an outlier and settling in a fast-rotator regime (Sergeev et al. 2022a). Changes in surface friction can also drive GCMs into one or the other climate regime (Carone et al. 2016), but we caution that the surface frictions used in Carone et al. (2016) may be unrealistically high, thus potentially overestimating the importance of this effect (Wicker et al. 2022). More generally, this shows that the tipping point is highly nonlinear and that small parameter changes can lead to a different climate regime.

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

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

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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

Last but not least, we identify that the main source of differences between the models appears in the upper part of the atmosphere. In Figures 2 and 3, temperature differences between models peak in the upper layers of the atmosphere, reaching about 10–20 K in both the Ben 1 and Ben 2 simulations. This is also where variability intrinsic to the models (visible in shading in Figures 2 and 3) is the largest.

In Figure 6, wind distribution is also very different from one model to another above the troposphere. Simulations in the fast-rotator regime (all but ROCKE-3D in the Ben 1 case; ExoCAM in the Ben 2 case) all exhibit two extratropical high-altitude (near 1–10 hPa) prograde jets. In ExoCAM's Ben 2 simulation, these jets extend to the troposphere and merge with the two aforementioned tropospheric jets. Simulations in the Rhines rotator regime (ROCKE-3D in the Ben 1 case; all but ExoCAM in the Ben 2 case) all exhibit an equatorial (near 1–10 hPa) prograde jet and two high-latitude midaltitude (near 10–100 hPa) retrograde jets.

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

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