ExoCAM: A 3D Climate Model for Exoplanet Atmospheres

Over the past decade, 3D climate models have been commonly used to study the climate, habitability, and potential observable phenomena of terrestrial extrasolar planets. There are numerous 3D climate models that are currently deployed for this purpose, including the Laboratorie de Météorologie Dynamique Generic (LMD-G) climate model (Wordsworth et al. 2011; Leconte et al. 2013a, 2013b, 2015; Bolmont et al. 2016; Turbet et al. 2016, 2018; Fauchez et al. 2019; Charnay et al. 2021; Lefèvre et al. 2021), the NASA Goddard Institute for Space Studies Resolving Orbital and Climate Keys of Earth and Extraterrestrial Environments with Dynamics (ROCKE-3D) climate model (Way et al. 2016, 2017; Fujii et al. 2017; Way et al. 2018; Aleinov et al. 2019; Colose et al. 2019, 2021; Del Genio et al. 2019; Jansen et al. 2019; Guzewich et al. 2021; Kane et al. 2021), the UK Met Office Unified Model (UM; Boutle et al. 2017; Mayne et al. 2019; Boutle et al. 2020; Sergeev et al. 2020), the National Center for Atmospheric Research (NCAR) Community Earth System Model (CESM; Shields et al. 2013; Yang et al. 2013, 2014, 2019b; Hu & Yang 2014; Kopparapu et al. 2016; Shields & Carns 2018) a growing list of other codes (Abe et al. 2011; Rauscher & Menou 2012; Haqq-Misra & Kopparapu 2014; Popp et al. 2016; Hammond & Pierrehumbert 2018; Ding & Wordsworth 2020; May & Rauscher 2020; Paradise et al. 2021), and ExoCAM, which is a branch of CESM designed specifically for exoplanets, whose description is the focus of this paper. Terrestrial exoplanet climate models typically branch from Earth climate system models that were created and continue to be supported by national modeling centers. While these exoplanet branches can leverage the utilities of their parent Earth model, they often diverge from their parent models in meaningful ways. Such is the relationship between ExoCAM and its parent model, CESM. Thus, it is important to provide definitive descriptions for exoplanet and planetary code branches for community reference (e.g., Way et al. 2017). The purpose of this paper is twofold. First, we provide a description and primary reference for the ExoCAM model, which includes the first presentation of a recent major overhaul to the radiative transfer scheme. Second, we conduct a series of intramodel climate sensitivity tests with ExoCAM. Starting from the TRAPPIST-1 Habitable Atmosphere Intercomparison (THAI) protocol (Fauchez et al. 2020), we test the model's sensitivity to a number of key tunable parameters and settings related to clouds, convection, precipitation, and sea ice to illustrate the robustness of the model's results for simulating habitable terrestrial extrasolar planets.

ExoCAM is an independently supported branch of NCAR CESM version 1.2.1, which facilitates simulations of generic planetary atmospheres. ExoCAM operates as a patch for CESM. The user must first download the core CESM code, which is freely available and provided by NCAR, ⁷ before modifying CESM with ExoCAM-specific configuration scripts and source code, also freely available and provided by the lead author on GitHub. ⁸ Importantly, note that ExoCAM is accompanied by a radiative transfer code for 3D climate models developed by the lead author, ExoRT, ⁹ which must also be downloaded by the user and linked with ExoCAM and CESM. ExoCAM and ExoRT are designed to be used together, with the ExoCAM GitHub repository containing system files, source code files, initial condition files, namelists, and the documentation needed to install and run 3D climate simulations in tandem with CESM, and the ExoRT GitHub repository containing all files needed for flexible radiative transfer calculations. Several years ago, ExoRT was separated from ExoCAM into its own repository in order to facilitate standalone 1D radiative transfer calculations, which significantly eases model development and validation demands for the radiative transfer component. In this paper, any reference to ExoCAM implicitly includes the linked capability of ExoRT.

ExoCAM and ExoRT were first made available to the public on Github in 2018 November. At the time of finalizing this paper, there have been a total of 37 papers published or accessible on arXiv by 20 different lead authors that use ExoCAM, whether in its entirety, by leveraging specifically credited component modules, or by leveraging specifically credited model outputs for extensive analysis as a central part of the work. ExoCAM has been well-received and well-used since being made public. ExoCAM, and its ancestral versions, have been employed for a wide variety of studies, including the Faint Young Sun Paradox for Earth (Wolf & Toon 2013, 2014a), high-CO₂ atmospheres (Wolf et al. 2018; Zhang et al. 2021), runaway and moist greenhouse thresholds for Earth (Wolf & Toon 2014b, 2015), Earth-like extrasolar planets (Wolf et al. 2017; Adams et al. 2019; Kang 2019a, 2019b, 2019c), habitable zone exoplanets around M-dwarf stars (Kopparapu et al. 2017; Wolf 2017; Haqq-Misra et al. 2018; Komacek & Abbot 2019; Komacek et al. 2019, 2020a, 2020b; Wolf et al. 2019; Yang et al. 2019a; Hu et al. 2020; Rushby et al. 2020; Suissa et al. 2020b; Wei et al. 2020; Zhang & Yang 2020; Chen et al. 2021; May et al. 2021; Song & Yang 2021), planets in circumbinary systems (Wolf et al. 2020), and model intercomparison studies (Yang et al. 2016, 2019b; Fauchez et al. 2020, 2021; Sergeev et al. 2021; Turbet et al. 2021). ExoCAM results have also been used off-line in projects to drive photochemistry and observability studies (Afrin Badhan et al. 2019; Suissa et al. 2020a).

ExoCAM is one of the four core 3D climate models used in the THAI project (Fauchez et al. 2020), alongside LMD-G, ROCKE-3D, and UM. The THAI project arose from conversations at the NASA Habitable Worlds Conference in Laramie, WY in 2017 November. It was determined that an important step forward for the growing exoclimate modeling community was to engage in a formal model intercomparison process. This paper describing ExoCAM is part of a collection of papers that have emerged from the THAI project and are featured in this special edition of the Planetary Science Journal (Fauchez et al. 2021; Sergeev et al. 2021; Turbet et al. 2021). The THAI project also led to a community workshop, the proceedings of which are discussed in the THAI workshop report (Fauchez et al. 2021). The THAI project aims to compare and cross-validate 3D climate models that are commonly used for predicting theoretical climate states of habitable zone extrasolar planets while using TRAPPIST-1e as the target object. We refer readers to the other papers in the series for detailed results on the climate states and theoretical observables among the four primary models. This paper provides an addition to the THAI series of papers, for the reader interested in taking a deeper dive into ExoCAM and its intramodel sensitivities.

As introduced in Section 1, ExoCAM is a modeling package designed to be used in tandem with NCAR CESM version 1.2.1 to facilitate flexible modeling of extrasolar planetary atmospheres. In recent years, interest in 3D exoplanet atmosphere modeling has exploded, but the road to proficiency for new users entering into the field is often steep and arduous. Furthermore, the real time and computing time required for full-cycle science (project genesis, model development, simulations, analysis, and publication) with 3D climate models means that current 3D exoclimate modeling groups are often already near or at their carrying capacity, and cannot contribute meaningfully to all of the spontaneous new collaborations that may emerge. The goal of ExoCAM is to lower the barrier for entry into the field for new users, by providing an accessible platform to begin climate modeling of planets and exoplanets while leveraging the NCAR CESM framework. ExoCAM provides accessibility for idealized planet configurations, which are most typical of exoplanet studies, with relatively straightforward handling of modifications of the basic geophysical, orbital, and atmospheric parameters of the modeled planet, and a flexible radiative transfer code. Note that ExoCAM and ExoRT are continually updated on GitHub when new features are developed or bugs are fixed. We encourage users to stay up to date with the code, and to contact the authors for guidance regarding model features, model capabilities, and appropriate use strategies.

2.1. General Features

2.2. Radiative Transfer

At the heart of ExoCAM is its radiative transfer package, ExoRT. Any flexible planetary model must be able to handle radiative transfer for a wide variety of atmospheric compositions. Native Earth-centric radiation schemes tend to perform poorly when pushed outside the bounds of Earth-similar atmospheres. The success of ExoCAM is thus critically dependent on providing a versatile radiative transfer package, and thus the largest share of historical development effort over the years has gone into improving ExoRT. ExoRT is provided via GitHub as a standalone package that is linked to ExoCAM at compile time. ExoRT uses a correlated-k distribution approach, and the two-stream approximation of Toon et al. (1989). The code lineage traces back to Urata & Toon (2013a, 2013b) and Colaprete & Toon (2003), with extensive modifications and reorganization occurring over time. There are now several radiative transfer options available for use with ExoCAM. While the most recent update is now recommended for all use cases moving forward, we describe all available versions below.

The oldest radiation version is named n28archean in the ExoRT repository, and is originally described in Wolf & Toon (2013). This radiative transfer version features 28 spectral intervals across the full spectrum ranging from 10 to 50,000 cm⁻¹ (0.2–1000 μm). Nominally, the longwave calculation is computed over 17 bins between 10 and 4000 cm⁻¹ (2.5–1000 μm), while the shortwave calculation is computed over 23 bins between 980 and 50,000 cm⁻¹ (0.2–10.2 μm). Note that in n28archean, and all other configurations described below, the shortwave bandpass is extended far into the thermal wavelengths in order to accommodate TRAPPIST-1-like spectral energy distributions out of the box. However, the exact bandpasses for the longwave and shortwave streams are modifiable by the user within file exo_init_ref.F90 at compile time, and can be adjusted, in relation to the stellar and planetary temperatures, to squeeze out additional computational efficiency. This is true for all ExoRT configurations. Molecular absorption is included from H₂O, CO₂, and CH₄. Correlated-k distributions were produced using the Atmospheric Environment Research Inc. line-by-line radiative transfer model (LBLRTM; Clough et al. 2005), using the HITRAN 2004 spectroscopic database (Rothman et al. 2005). For all gases, we assumed Voigt line wings cut at 25 cm⁻¹. Water vapor and CO₂ continuum absorption are included using MT_CKD version 2.5, while assuming that the CO₂ self-broadened continuum is 1.3 times the foreign broadening continuum (Kasting et al. 1984; Halevy et al. 2009). Overlapping gas absorption is handled using an amount-weighted scheme, with no more than two overlapping in each band (Mlawer et al. 1997). Collision induced absorption (CIA) is included for N₂–N₂, N₂–H₂, and H₂–H₂ pairs. This version was originally designed for studying the climate of Archean Earth, for CO₂ amounts up to several tenths of a bar and up to several thousand ppm of CH₄, along with Earth-like water vapor amounts. The results were initially compared to line-by-line calculations (see Figure 1 in Wolf & Toon 2013) for Archean-like atmospheres. Note that this version of the radiative transfer is used for the official THAI simulations discussed in Turbet et al. (2021), Sergeev et al. (2021), and Fauchez et al. (2021).

In 2016, an intercomparison project was initiated by Yang et al. (2016) to study the performance of various radiative transfer codes used for simulating warm, moist terrestrial planet atmospheres, up to 360 K. This study found that n28archean (referred to as "CAM4_Wolf" in Yang et al. 2016) performed well in the longwave; however, in the limit of warm moist atmospheres and M-dwarf stellar spectra, it significantly overestimated the strength of near-infrared absorption of the downwelling shortwave radiation stream, resulting in too much radiation being absorbed by the atmosphere and not enough reaching the surface. Motivated by this study, ExoRT was updated to its second iteration, first described in Kopparapu et al. (2017) and referred to in the repository as n42h2o. This version increased the number of spectral intervals to 42 across the full spectrum ranging from 10 to 50,000 cm⁻¹ (0.2–1000 μm). Nominally, the longwave calculation is computed over 19 bins between 10 and 4000 cm⁻¹ (2.5–1000 μm), while the shortwave calculation is computed over 35 bins between 820 and 50,000 cm⁻¹ (0.2–12.2 μm). Spectral resolution was preferentially increased in the near-infrared region in order to address the deficiencies discovered in Yang et al. (2016). This version includes only H₂O gas absorption using the HITRAN 2012 (Rothman et al. 2013) spectroscopic database, along with the BPS water vapor continuum (Paynter & Ramaswamy 2011). Correlated-k distributions for H₂O were produced using HELIOS-K, a GPU-accelerated k-distribution sorting code (Grimm & Heng 2015), which is freely available on GitHub. ¹⁹ Line wings assume a Voigt shape and a 25 cm⁻¹ cut-off. This radiative transfer scheme improved the model performance in the near-infrared around M-dwarf stars, and served the purpose of facilitating studies of water vapor and water cloud interactions on tidally locked planets around M-dwarf stars. CIA is included for N₂–N₂, N₂–H₂, and H₂–H₂ pairs. This scheme was compared favorably against LBLRTM for warm, moist N₂, H₂O atmospheres (see Figure 1 in Kopparapu et al. 2017). Note also that there is an analogous 68 spectral interval version of this radiation scheme, n68h2o, which was used in Wolf et al. (2019) to yield better output spectral resolution for use in phase curve studies. This version has a full spectrum range from 0 to 42,087.00 cm⁻¹ (0.238–infinity μm), with longwave and shortwave calculations computed over the entire spectral range.

In the fall of 2020, ExoRT received a major overhaul, which improved model performance and greatly improved the flexibility of the model to be able to incorporate a wider variety of overlapping gases. This most recent version is called n68equiv. All users of ExoCAM are now recommended to use n68equiv, superseding both n42h2o and n28archean. This new version uses 68 spectral intervals across the full spectrum ranging from 0 to 42,087.00 cm⁻¹ (0.238–infinity μm). Nominally, the longwave calculation is computed over 37 bins between infinity and 4030 cm⁻¹ (2.48–infinity μm), while the shortwave calculation is computed over 53 bins between 875 and 42,087 cm⁻¹ (0.238–11.43 μm). Importantly, n68equiv now uses the equivalent extinction method to handle overlapping gases (Amundsen et al. 2017). In the equivalent extinction method, the primary gas absorber in each spectral interval is treated with an 8 Gauss point k-distribution, while all additional species that have absorption in the band are treated as gray absorbers using their band-averaged absorption value. The major species in each band is dynamically selected at each model grid cell and each model time step by first comparing the band-averaged optical depths of each species before the full radiative calculation is done. The equivalent extinction method for gas overlap offers a significant advantage over other overlap schemes because it permits many overlapping species to be present in each spectral interval without requiring the prior computation of mixed gas k-coefficient tables, which become unwieldy if considering more than two gas species per interval. Presently, n68equiv includes gas absorption from H₂O, CO₂, and CH₄, with k-distributions produced by HELIOS-K (Grimm & Heng 2015) using the HITRAN 2016 spectroscopic database (Gordon et al. 2017). Voigt line shapes with a 25 cm⁻¹ line cut-off were used for H₂O and for CH₄. For H₂O, the plinth, or base-term (Paynter & Ramaswamy 2011), was removed to prevent the double counting of absorption when combined with continuum absorption. Self and foreign water vapor continuum coefficients are included using MT_CKD version 3.2. However, instead of using band-averaged values, we found that improved performance in the near-infrared region is achieved by fitting the water vapor self and foreign continuum coefficients to an 8 point k-distribution, included under the assumption of perfect correlation when water vapor is present in a grid cell. For CO₂, we assume the subLorentzian line shape of Perrin & Hartmann (1989) with a 500 cm⁻¹ line cut-off. In addition to N₂, H₂ CIA pairs, n68equiv includes CIA for CO₂–CO₂ (Wordsworth et al. 2010) and for CO₂–H₂ and CO₂–CH₄ from Turbet et al. (2020). Note that there is also an analogous 84 spectral interval version, n84equiv, which is identical in description to n68equiv, except 16 additional spectral intervals have been added shortward of 0.238 μm. Testing indicates that for stars hotter than about 6500 K, nonnegligible stellar flux is pushed shortward of ∼0.2 μm, and was effectively being lost by n68equiv. If using n84equiv, users are encouraged to examine and tailor the specific shortwave bandpass limits in exo_init_ref.F90 to balance the need for encompassing far-ultraviolet wavelength ranges with computational cost.

Some features are shared by all available versions of ExoRT, and are described here. All versions treat liquid and ice cloud droplets as Mie scattering particles, with refractive indices taken from Segelstein (1981) and Warren & Brandt (2008). 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). Rayleigh scattering is included with the parameterization of Vardavas & Carver (1984), including the effects of N₂, CO₂, H₂O, and H₂. The longwave and shortwave streams share a common spectral grid and k-coefficient tables. The extent of the longwave and shortwave spectral integration streams can be readily adjusted to balance completeness versus computational speed, depending upon the stellar type assumed and planetary temperatures involved (see file exo_init_ref.F90). All k-coefficient sets currently use log₁₀(P₁/P₂) = 0.1 spacing. The n28archean scheme has a pressure range from 100 bars to 0.01 mb, and a temperature range from 80 to 520 K. The n68equiv and n42h2o versions have a pressure range from 10 bars to 0.01 mb, and a temperature range from 100 to 500 K. All versions use the Gauss point binning originated from the Rapid Radiative Transfer Model (RRTM; Mlawer et al. 1997) that is weighted toward 1 and thus favors capturing the peak of the absorption lines in each band. A variety of input stellar spectra are provided for each radiative transfer version, ranging from TRAPPIST-1-like stars (2600 K) up to F-type stars (10,000 K). Additional software is provided, written in IDL, for manipulating input and output data for ExoRT.

Past papers, as noted, illustrate the performance of n42h2o (Kopparapu et al. 2017) and n28archean (Wolf & Toon 2013; Yang et al. 2016). Here, in Figures 1 and 2, we illustrate the basic radiative transfer performance of n68equiv against line-by-line calculations for standard atmospheres. Recall that n68equiv is the recommended version for all use cases; however, n28equiv can still be used for modern Earth and Archean Earth cases, while offering faster computation than n68equiv, and n42h2o remains valid, albeit while lacking greenhouse gases besides H₂O. Figure 1 duplicates the radiative transfer calculations of warm, moist atmospheres originally conducted in in Yang et al. (2016), including comparisons against LBLRTM (Clough et al. 2005) and the Spectral Mapping Atmospheric Radiative Transfer model (SMART; Meadows & Crisp 1996; Goldblatt et al. 2013), considering both the outgoing longwave radiation (OLR) at the top-of-atmosphere (TOA), and the downwelling stellar flux at the surface under irradiation from Solar and M-dwarf (3400 K blackbody) spectra, respectively. The x-axis tracks with the mean surface temperature. These atmospheres are assumed to contain 1 bar N₂, 376 ppm of CO₂, and to be fully saturated with water vapor. Water vapor adds pressure to the atmosphere, thus by 360 K the total surface pressure is ∼1.6 bars. For longwave fluxes, n68equiv performs similarly to LBLRTM, while both diverge from SMART for higher-temperature H₂O-rich atmospheres. This is likely because n68equiv and LBLRTM use the same H₂O continuum treatment (MT_CKD), while SMART uses a χ-factor approach (Goldblatt et al. 2013). The shortwave performance of n68equiv is dramatically improved compared to n28archean (Yang et al. 2016), but differences are still evident for hotter atmospheres, likely due to the complexity of the near-infrared water vapor spectrum, which is difficult to fully resolve at coarse spectral resolutions suitable for 3D climate model operations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 2 shows comparisons of the outgoing longwave radiative flux between n68equiv, SMART, and SOCRATES (a high-resolution band model; Edwards & Slingo 1996) for a 2 bar pure-CO₂ dry atmosphere at Mars gravity with a 250 K surface temperature. The SMART calculations shown were originally published in Kopparapu et al. (2013). Both n68equiv and SOCRATES yield OLR values of ∼93.3 W m⁻², while SMART yields an OLR of ∼88.5 W m⁻². Here, n68equiv notably lacks the spectral resolution to capture the fine detail of CO₂ absorption between 800 and 1200 cm⁻¹. Still, the primary differences between n68equiv and SOCRATES compared to SMART are most noticeable in the wing regions of the primary 667 cm⁻¹ CO₂ absorption band. This may be attributed to differences between the line wing treatments of Perrin & Hartmann (1989), used in n68equiv and SOCRATES, and the stronger line wing formulation used in SMART (Meadows & Crisp 1996). See also Figure 1 in Halevy et al. (2009) for a comparison of relevant CO₂ line wing treatments. Note that several other studies have computed identical CO₂ radiation tests. Identical column tests by Wordsworth et al. (2010) with a correlated-k band model yielded 88.2 W m⁻², in line with SMART. Later, Wordsworth et al. (2017) conducted line-by-line calculations that yielded ∼92 W m⁻². The SOCRATES calculations featured in Guzewich et al. (2021) yielded 96 W m⁻²; however, those calculations used an out-of-date CO₂–CO₂ CIA file that is missing absorption between 250 and 500 cm⁻¹. Note that our goal here is not to anoint any one model as better than the others, but rather to give the reader an adequate benchmarking of performance versus other popular models. Note that the n28archean version of ExoRT, which has been used in prior ExoCAM studies of CO₂-dominated atmospheres (Wolf 2018; Wolf et al. 2018) and in the standard THAI simulations (Sergeev et al. 2021; Turbet et al. 2021), yields an OLR of 81.3 W m⁻², an overestimation of the greenhouse effect for CO₂-dominated atmospheres. Future studies of CO₂-dominated atmospheres with ExoCAM should use n68equiv.

2.3. Additional Features

2.4. Computational Performance

A common criticism of 3D climate models used for exoplanet modeling is that they often contain tunable parameters that were originally calibrated to replicate the climate of the modern Earth, and thus may be a source of uncertainty, leading to untrustworthiness of the modeling results for extrasolar planets. Of particular interest here, tunable parameters are common in the moist physics routines for treating subgrid-scale convection and clouds. Horizontal grid resolutions for 3D climate models are generally large, with horizontal boxes spanning as much as several hundred kilometers on each side. Thus, unresolved small-scale cloud and convection processes, such as localized cloud fractions, cloud overlap, and convective plumes, must be represented statistically through parameterizations. In CESM, the values of the tunable parameters that control these subgrid-scale representations are dependent on the choices of horizontal resolution, dynamical core, and moist physics routines. While the principal THAI simulations used the default model tunings for the selected setup (i.e., 4° × 5° grid, FV dynamical core, and CAM4 moist physics), it is unknown if these tunings would necessarily be appropriate for an alien planet. Thus, an important part of building confidence in our 3D exoplanet climate models is not only in comparing them against other leading 3D models, but also in comparing them against themselves, while exploring the effects of the many important assumptions and choices that are baked into out-of-the-box simulations. In this section, we constrain how large an effect these tunable parameters may have on our subsequent climate predictions of TRAPPIST-1e. While we acknowledge that this section does not constitute a comprehensive analysis of all of the temporal and spatial differences evident due to parameter sensitivities, our first-order analysis supports the robustness of our 3D exoclimate modeling results against plausible changes in important tunable parameters.

We conduct a series of intramodel sensitivity tests with ExoCAM using the THAI protocol (Fauchez et al. 2020) hab1 and hab2 cases as our starting points. We refer the reader to the core series of THAI papers for more details of the THAI protocol and also for detailed analysis of the results (Fauchez et al. 2021; Sergeev et al. 2021; Turbet et al. 2021). To summarize the THAI protocols briefly here, both hab1 and hab2 cases simulate the climate of TRAPPIST-1e, having a radius of 0.91 R_⊕, a surface gravity of 0.93 R_⊕, and a period of 6.1 Earth days in synchronous rotation with a 2600 K host star, while receiving 900 W m⁻² of stellar irradiation. For both the hab1 and hab2 cases, the planet is assumed to be completely ocean-covered with an ocean albedo of 0.06 and a snow and sea ice albedo uniformly of 0.25. The ocean is treated as a 100 m deep slab ocean with no heat transport. The atmospheric compositions used are a 1 bar N₂-dominated atmosphere with 400 ppm of CO₂ for hab1, and a 1 bar CO₂ atmosphere for hab2. In each case, the water vapor, water clouds, and sea ice are allowed to evolve naturally in the system.

Sergeev et al. (2021) discuss the two moist cases (hab1 and hab2), contrasting the model results between ExoCAM, LMD-G, ROCKE-3D, and UM for the standard THAI project intercomparison. Here, including the two moist cases from the standard THAI project discussed in Sergeev et al. (2021), we consider the results from 43 ExoCAM simulations, iterating on hab1 and hab2. Table 1 summarizes the parameters for which sensitivity tests are conducted, including the default values for our model configuration and the ranges explored. The standard THAI simulations used the default values listed. Here, we bracket these default values with end-member values drawn from the CESM parameter selections used for other model configurations that they offer. Recall that ExoCAM is a patch for CESM that enables exoplanet configurations. CESM independently contains numerous different model configurations spanning preindustrial through twenty-first-century Earth, with different combinations of vertical and horizontal resolutions, model physics, and dynamical cores. Each of the many configurations offered requires differing values for the subgrid-scale tuning parameters in question. Here we explore a range parameter values within the context of the ExoCAM configuration, with plausible values informed by other existing CESM configuration requirements.

The parameters explored include the cloud overlap assumption, the liquid cloud droplet radius, cloud fraction relative humidity thresholds for low and high clouds, critical radii for the autoconversion of liquid and ice condensate into precipitation, the autoconversion coefficient for convective clouds, and the dynamical substepping frequency. Because individual clouds cannot be resolved within large grid boxes, their inherent horizontal patchiness and relations to adjacent cloud layers in the vertical dimension are parameterized through cloud fractions linked to cloud overlap assumptions (e.g., Pincus et al. 2003). The cloud fraction represents the fraction of the horizontal grid box that is covered by clouds, with the remaining portion assumed to be clear. Cloud overlap assumptions are used by the radiative transfer model to interpret the relative vertical layerings of clear and cloudy patches of sky. Cloud fraction is generally controlled by simple scalings against the grid-box average relative humidity. When the grid-box average relative humidity exceeds a certain threshold (typically ∼0.8–0.9), clouds begin to form with a cloud fraction that ramps toward unity as the grid-box average relative humidity approaches saturation. Cloud droplet sizes and their subsequent conversion to precipitation are also treated via scaling considerations. Autoconversion thresholds control the maximum size that typical cloud particles can grow to before they become rain or snow, with separate thresholds for liquid, ice, and convectively produced clouds. Autoconversion relies on the available amount of cloud condensate along with assumptions regarding the number density of cloud particles to determine when the size threshold to convert to precipitation is reached. We also test fixing liquid cloud droplet sizes to different values, which affects how fast they fall out of the atmosphere and strongly affects their radiative transfer properties. Note that the default value of 14 μm in the Rasch & Kristjánsson (1998) cloud scheme is taken from Earth observations of liquid water clouds over ocean waters, far from land.

Additionally, we conduct sensitivity tests using the native sea ice and snow albedos. Sea ice and snow albedos are divided into visible and near-infrared bands, and set to 0.67 and 0.3 for sea ice and 0.8 and 0.68 for sea ice and snow, respectively, in typical ExoCAM simulations, following Shields et al. (2013). Note that such two-band sea ice and snow albedos would benefit from recalibration according to the incoming stellar energy distribution of the specific star used (e.g., Shields & Carns 2018); nonetheless, these default albedos provide a meaningful point of comparison against the spectrally uniform albedos specified in the THAI protocol (Fauchez et al. 2020). We also test the effects of using the native CAM5 cloud and convention physics (Morrison & Gettelman 2008; Park & Bretherton 2009; Park et al. 2014; henceforth, MG clouds), instead of the CAM4 cloud and convection physics (Hack 1994; Zhang & McFarlane 1995; Rasch & Kristjánsson 1998) that have been exclusively used in published ExoCAM works to date. As alluded to in Section 2.1, the CAM5 cloud physics uses a two-moment scheme that relies on time- and spatial-dependent aerosols to serve as sites for cloud condensation, and thus can capture the microphysical interactions between clouds and aerosols with greater sophistication than the CAM4 scheme. Granted, this introduces another set of unknowns, as aerosols must be specified in the context of an alien world. Here, we use a default collection of time- and spatial-dependent aerosol fields derived for Earth to drive our CAM5 cloud microphysics. We acknowledge that it is unlikely that TRAPPIST-1e would have an identical aerosol field as the Earth.

Tables 2 and 3 present the global mean climatological results for the hab1 and hab2 simulations, respectively. The model outputs have been averaged over 10 Earth years (∼600 TRAPPIST-1e orbits) after reaching thermal and radiative balance. As mentioned above, for the principal THAI project, the n28archean radiation scheme was used because the development of the n68equiv scheme was not completed in a timely manner. Henceforth, these are referred to as the THAI standard cases. In this paper, we focus on the results using our updated n68equiv radiation scheme, as it provides more trustworthy results, particularly for the hab2 cases, for the reasons discussed in Section 2.2. Henceforth, our nominal simulations with n68equiv are referred to as our control cases. Both the THAI standard and control cases used the default parameter selections from Table 1. The subsequent simulations listed in each of Tables 2 and 3 give the particular tuning parameter that was altered. Included in each table are the global mean surface temperature (T_S ), the planetary albedo, and the integrated column amounts of water vapor, liquid water clouds, and ice water clouds. All cases listed use the CAM4 cloud and convection physics, except for the cases described as "MG clouds." Note that no sensitivity test was conducted using the native sea ice and snow albedos for hab2 because all sea ice and snow are trapped on the night-side of the planet, rendering any ice and snow albedo effects irrelevant due to the assumed synchronous rotation of TRAPPIST-1e.

Figure 3 shows longitude–latitude contour plots for our hab1 and hab2 control simulations. This figure provides a view of several important climatological variables: the surface temperature, OLR, 500 mb vector winds, liquid water cloud column, and ice water cloud column. Note that the white line in the surface temperature plots outlines the approximate sea ice margin. Contour plots of these variables, and numerous other variables, for all hab1 and hab2 sensitivity simulations are included in an accompanying Zenodo repository. ²² Output data of the time mean climates of all simulations are also provided in this repository. For both hab1 and hab2, the overall climate states and the dynamical regimes do not appear sensitive to our suite of tests. A nominal range in global mean surface temperatures of ∼10 K is experienced in our sensitivity tests, with changing the ice albedo and the liquid cloud droplet sizes generally providing the most significant feedback. Still, the climate states of each set are not changed, with hab1 remaining cold and icy, except around the substellar point, and hab1 remaining temperate. A qualitative examination of the wind fields (see the Zenodo supplement) of the hab2 simulations does not reveal any dramatic shifts in the dynamical regimes. While the vector wind speeds can vary in magnitude, 200 and 500 mb vector wind fields appear qualitatively similar, with vortices located on the night-side of the planet near the western terminator, and strong supperrotation in the western hemisphere, which then splits in two across to the eastern hemisphere. A similar statement can be made for the hab1 cases, with the lone exception being the MG case. This case displays the vortices in the eastern hemisphere, as opposed to the western hemisphere, and has weaker equatorial winds in the western hemisphere. This change in dynamics is indicative of the importance of physics–dynamics couplings, and the role of clouds (and their radiative forcings) in influencing the general circulation, and not just in being reactive to the general circulation. Note that permitting the use of the Morrision–Gettlemen two-moment cloud scheme in ExoCAM is a new feature that relies on the specification of a background aerosol field, which can depend on time and space. I encourage others to explore this new model feature within ExoCAM more deeply.

Zoom In Zoom Out Reset image size

To better illustrate the resulting spreads in the climates of our sensitivity tests, we have collected the global mean climatological information in Figures 4–10. Each figure shows a scatter plot comparing two relevant climatological quantities along the x- and y-axes, with the hab1 results in the left panel and the hab2 results in the right panel. The THAI standard simulations, featured in Sergeev et al. (2021), are shown as a black "X"; simulation number 1 in Tables 2 and 3. An orange cross indicates the revised THAI simulations using our updated radiative transfer scheme n68equiv, which we take as our control simulation for further model sensitivity experiments in this work; simulation number 2 in Tables 2 and 3. Blue diamonds show the values for the sensitivity tests exploring different options and tuning parameters in the model. Relative outlier simulations are labeled on each plot to add context; however, generally our sensitivity tests cluster about the control case, implying that significant changes in the overall climate regimes cannot be driven by modifications to subgrid-scale tunable parameters. Note the standard THAI simulations (e.g., Sergeev et al. 2021) show insignificant deviation from our revised control cases for hab1. However, the standard THAI simulations show marked differences to the hab2 cases in most plots. This is a straightforward consequence related to differences between the older radiative transfer version n28archean, which overestimates both the near-infrared absorption by water vapor under M-dwarf spectra and also the greenhouse effect from CO₂, as discussed in Section 2.2. The resulting climate for our THAI standard case is ∼11.5 K warmer than for our control case, which feeds back on most climatologically interesting fields. The results from using n68equiv yield more trustworthy results, which are better in line with the other models discussed in Sergeev et al. (2021). In interpreting these sensitivity simulations, the reader is reminded that tolerances for exoplanet climate science are considerably more forgiving than they would be for Earth climate sciences. For instance, several Kelvin differences in global mean temperature predictions, while enormous in the context of twenty-first-century Earth, do not affect the general conclusions regarding an exoplanet's climate state and potential for habitability.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 4 shows the global mean surface temperature (T_s ) plotted versus the global mean emission temperature (T_emiss). The emission temperature is calculated from a planet's global mean outgoing thermal radiation and inverting the Stefan–Boltzmann Law. Comparing T_s to T_emiss yields information about the greenhouse effect of a planet. Note that the x- and y-axes have different scales for the left and right panels, but each axis spans 15 K. Most parameter choices only have minor effects on both the T_s and T_emiss. There is significant clustering around the control for both the hab1 and hab2 cases, featuring T_s and T_emiss differences of only ∼2 K. Still, several significant outlier cases are noted in Figure 4. In particular, the liquid cloud particle sizes (rliq), which are fixed to constant values in the radiation, have a meaningful impact on both the hab1 and hab2 climates. Reducing the cloud particle size from 14 to 7 μm results in significant global mean cooling, up to ∼8 K, due to increased cloud scattering effects and decreased particle sedimentation, while increasing the cloud particle size to 21 μm can increase the T_s by several K. Changes to cloud properties along with T_s combine to affect the radiating level of the atmosphere, reflected by the values for T_emiss. Naturally, the use of native ice albedos for hab1, which are significantly higher than the 0.25 value used for the control case, cools the hab1 case significantly. Interestingly, despite significant differences in methodology and in the water vapor and condensate partitioning (e.g., Figures 6 and 7), the MG cloud cases have a relatively modest impact on T_s and T_emiss.

In Figure 5, we plot the global mean planetary albedo versus the planetary emissivity. Here, we calculate the planetary emissivity as the ratio of the outgoing longwave radiative flux at the TOA to the outgoing radiative flux at the surface. Using this simple formulation, a value of 1 indicates a planet with no greenhouse effect, while progressively lower planetary emissivities indicate a strengthening greenhouse effect. The greenhouse effect is the totality of CO₂, water vapor, and cloud greenhouse effects. Note that the y-axis (planetary albedo) is identical for both the left and right panels, but the x-axes (planetary emissivity) have different scales. Naturally, the hab2 simulations have a significantly lower emissivity than the hab1 cases, due to their strong CO₂ greenhouse effect. Notable outliers again are driven by the liquid cloud particle size, the MG cloud simulations, and the change in sea ice albedo. While the mean surface temperature of the hab1 MG case does not deviate from the pack, Figure 5 reveals that it has both a weaker greenhouse effect and a lower planetary albedo, which offset each other.

Figures 6 and 7 show the global mean liquid and ice cloud column amounts plotted against the water vapor column amount, respectively. In these figures, note that the cloud and water vapor columns are given in units of kg m⁻², while in Tables 2 and 3 the cloud water columns are given in g m⁻². Note that in Figures 6 and 7, the y-axis scales are identical, but the x-axis scales are different, because there is significantly less water vapor in the cold hab1 atmospheres compared to the temperate hab2 atmospheres that are warmed by 1 bar of CO₂. By plotting the column cloud amounts against the column water vapor amounts, one gets a general idea of the partitioning of water in the atmosphere between thermodynamic phases. In both Figures 6 and 7, for the hab1 cases, we again see significant clustering around the control case. Clustering is also evident for the hab2 sensitivity test simulations; however, the column water vapor amounts of the sensitivity tests are generally larger than is found in the control case. The MG cloud cases are the largest outliers for both the ice and liquid cloud amounts. The MG cloud simulations have significantly less ice and liquid cloud water compared to the control in all cases, and generally have a water vapor column amount toward the lower end of the pack. Additionally, we find that various autoconversion size threshold parameters can have a meaningful effect of the cloud amounts. Autoconversion tuning parameters control the particle size at which condensate is converted to precipitation and thus falls out of the atmosphere. Smaller autoconversion thresholds increase precipitation rates and result in less cloud water remaining in the atmosphere. Interestingly, note also that water vapor and clouds have a nonnegligible sensitivity to the dynamical substepping frequency (nsplit), particularly for the hab2 simulations. Thus, users must be cautious when ramping up the number of dynamical substeps to improve numerical stability midsimulation, as it feeds back on the clouds and climate.

Figure 8 shows the global mean column integrated water amount, including both liquid and ice cloud totals, comparing the substellar hemisphere mean value versus the antistellar hemisphere mean. Interestingly, the hab2 and hab1 simulations contain comparable amounts of substellar cloud water (∼200 g m⁻²). However, the warmer hab2 climates experience more significant day-to-night transport of heat and water vapor, resulting in reduced day-to-night temperature gradients, and greatly increased atmospheric water vapor and cloud condensate on the night-side of the planet. These features can be clearly seen in Figure 3.

We have also computed the global mean longwave and shortwave cloud radiative forcings produced in our runs. In Figure 9, we plot the global mean shortwave cloud forcing versus the global mean longwave cloud forcing. Shortwave cloud forcings are negative, as liquid water clouds on the day-side of the planet predominantly reflect stellar energy back to space, reducing the total amount of absorbed stellar radiation by the planet. Longwave cloud forcings are positive, as high-altitude ice clouds produce a greenhouse effect by lowering the effective emitting temperature of the planet. Perhaps surprisingly, the hab1 cases feature a larger magnitude of shortwave forcing than the hab2 cases. One can observe in Figure 3 that clouds are more concentrated right at the substellar point for hab1, contributing to increased stellar reflection, whereas for hab2, clouds are more evenly spread about the planet. Both the hab1 and hab2 clusters have similar longwave forcings. As expected, changing elements such as the liquid drop effective radii has noticeable effects on the cloud forcing fields. Similar to earlier cloud figures, the MG simulations stand as significant outliers compared to the rest, with reduced magnitudes of both longwave and shortwave cloud forcings. However, their reduction in the magnitudes of cloud forcings approximately balances out, as their net cloud forcing falls among the main cluster, as shown in Figure 10.

Lastly, Figure 10 plots the net cloud radiative forcings (shortwave + longwave) versus the clear-sky greenhouse effect. The clear-sky greenhouse effect is computed as the difference between the upwelling longwave flux at the surface and the upwelling longwave flux at the TOA as calculated in parallel online clear-sky radiative flux calculations. Water vapor and CO₂ contribute to the clear-sky greenhouse effect. Naturally, the magnitude of the clear-sky greenhouse effect is much stronger for the hab2 cases, which have 1 bar CO₂, compared to the hab1 cases, with only have 400 ppm of CO₂. The water vapor greenhouse feedback also plays a role in amplifying the clear-sky greenhouse effect in the hab2 cases. This water vapor amplification is perhaps best seen when comparing the hab2 control and THAI standard simulations, where the standard case has a much stronger clear-sky greenhouse effect than our control simulation due in large part to increased water vapor. Much of the additional warmth of the THAI standard case likely originates from the water vapor greenhouse feedback amplification of the original overestimation of the CO₂ greenhouse effect.

In total, our results show that ExoCAM exhibits some sensitivity to a variety of model parameter selections; however, in the context of TRAPPIST-1e, none of these selections, taken individually, dramatically change the climate outcomes and overall conclusions. Unsurprisingly, bulk differences in the radiative transfer performance between n28archean and n68equiv can yield meaningful differences in climate results for CO₂-dominated atmospheres. Furthermore, a wholesale change in the model cloud schemes, from CAM4 to CAM5, produces significant deviations in cloud fields, but surprisingly little change to the mean surface temperature, as changes to longwave and shortwave cloud forcings remarkably balance out. Among the free parameters studied here, notably the assumptions for the surface albedos and liquid cloud droplet particle sizes play the largest roles in modulating climate. It is also possible that combinations of parameters could drive larger effects than are individually shown here. However, such a multiparameter study would be computationally expensive. Future 3D climate simulations could benefit from more refined treatments of surface albedos and improved cloud and aerosol microphysics.

Note that other published works have also explored intramodel sensitivities for exoplanet modeling using ExoCAM and the NCAR family of models. It is worth briefly reviewing these results here. Wei et al. (2020) tested the climate sensitivity of ExoCAM under different horizontal and vertical resolutions for simulating tidally locked ocean-covered planets around M-dwarf stars. Their simulations show differences in T_s of no more than ∼5 K for changes in horizontal resolution, and negligible differences in T_s from changes in vertical resolution. Presumably, the CESM default resolution-dependent tuning parameters were used, but these are not discussed in their work. Bin et al. (2018) tested standard versions of CAM4 and CAM5 for calculating the inner edge of the habitable zone, and compared their results to ExoCAM (Kopparapu et al. 2017). Bin et al. (2018) found some notable differences between the CAM4, CAM5, and ExoCAM results featured in Kopparapu et al. (2017). Note that the standard version of CAM5 uses the RRMTG radiative transfer package (Mlawer et al. 1997) and the MG cloud scheme tested here, although assumptions regarding aerosol fields are not discussed. In general, Bin et al. (2018) found that standard CAM5 performs similarly to ExoCAM in predicting the inner edge of the habitable zone for synchronously rotating M-dwarf planets. However, the specific physical causes of these model deviations are not identified because there are simultaneous configuration differences in the radiation, clouds, convection, and boundary layer process between standard CAM5 and ExoCAM as used in Kopparapu et al. (2017). Wolf & Toon (2015) explored the model sensitivity of a legacy version of ExoCAM to cloud and convection tuning parameters for moist greenhouse climates possible for future Earth under high solar irradiation. They found that liquid cloud particle sizes and cloud fraction relative humidity thresholds could drive changes of up to ∼10 K, while the tuning of convective timescales generally had little effect. Yang et al. (2013) conducted sensitivity tests for cool, slow-rotating, tidally locked ocean-covered planets around M-dwarf stars, for a variety of NCAR model configurations and parameters, similarly finding that the cloud liquid droplet sizes have the largest effect on climate. Finally, Rushby et al. (2020) used ExoCAM to simulate TRAPPIST-1 planets as completely land-covered planets, and found that the assumed surface albedo can have major ramifications for the global mean surface temperature of the planet.

ExoCAM has proven to be a useful community tool that has enabled numerous different research groups to study the climate and potential for habitability of extrasolar planets. However, to ensure continued success, ExoCAM must continue to grow. Presently, there are funded projects and collaborations to extend the functionality of ExoCAM to new worlds. These will include expanding the current pressure-temperature spaces of the absorption coefficient tables along with adding new absorbing gas species in order to simulate a wider variety of planets, including oxygen-rich planets and subNeptune worlds. Work is also currently ongoing to merge features from Mars-CAM (Urata & Toon 2013a, 2013b) with ExoCAM, including CO₂ condensation and the capability to model low-pressure planets. Other ongoing work will further develop the CARMA functionalities within ExoCAM, including fractal organic hazes (Wolf & Toon 2010), along with dust and ice clouds (Hartwick et al. 2019). Finally, a collaborative effort will begin using ExoCAM with the dynamic ocean model (POP2) to study the long-term evolution of planetary climates in binary star systems. All newly developed features will eventually be ported to the public ExoCAM and ExoRT versions on GitHub for community use, as we believe strongly in the importance of providing useful tools for the community. Where my energy and imagination wane, hopefully others can carry the science forward.

Of further note, one common user critique of ExoCAM is that it is not collected into a single location, but rather the user must download CESM separately (via an NCAR hosted Subversion repository), and then download ExoCAM and ExoRT (via GitHub), and then go through several steps to combine the codes together. However, starting with CESM version 2, CESM is now available on GitHub, and thus can be forked. A long-term goal will be to create a forked version of CESM2 with ExoCAM and ExoRT fully integrated, so that users can obtain the entire code from a single location as a single code package. However, such a large refactoring of the code would be time-intensive, and is unfortunately not on the near horizon.

ExoCAM is a branch of NCAR CESM version 1.2.1 that has been designed specifically for studying planetary and exoplanetary atmospheres. Both CESM and ExoCAM are freely available to the community. CESM is provided by NCAR, while ExoCAM is provided by the authors via GitHub as an easily applied model extension. After its first public release in 2018 November, ExoCAM has been adopted by many groups in the community for simulating exoplanetary atmospheres. At the time of the writing of this paper, ExoCAM has been used as a pillar in 37 peer-reviewed papers, written by 20 different first authors, including numerous student-led projects. ExoCAM is also one of the four core climate models used in the THAI project. In this paper, we have provided the first formal description of ExoCAM, along with details of a recent major upgrade to the radiative transfer module, ExoRT.

In addition, here we have also tested ExoCAM for climate sensitivities to a variety of tunable parameters and model selections. We have tested tunable parameters regarding clouds, convection, time-stepping, and sea ice and snow albedos, along with model configurations regarding the radiative transfer scheme and moist physics packages. Sensitivity tests were conducted in adherence with the THAI protocol (Fauchez et al. 2020) for simulating habitable iterations of TRAPPIST-1e. Changes to primary physical schemes, for instance the radiative transfer module or the cloud microphysics package, can yield meaningful changes to cloud amounts and the global mean surface temperature. Changes to tunable parameters of the snow and sea ice albedos or cloud particle sizes can also drive meaningful changes in the global mean surface temperatures. However, the magnitude of these changes is never more than ∼10 K, indicating that the basic climate states yielded for TRAPPIST-1e are reasonably robust for a variety of intramodel sensitivity tests. Future sensitivity studies may want to explore if the combined sensitivities of multiple parameters changed simultaneously would yield larger perturbations to the climate.

In the coming years, there are several projects planned with ExoCAM that will expand its functionalities. In particular, development will focus on expanding the radiative transfer capabilities to simulate a wider range of atmospheres, adding CO₂ condensation, and leveraging the CARMA cloud and aerosol microphysics package. In continuing to push the model forward, we hope to continue to provide the community with a useful tool for studying planetary and exoplanetary atmospheres.

E.T.W., R.K., and T.J.F. acknowledge support from the GSFC Sellers Exoplanet Environments Collaboration (SEEC), which is supported by NASA's Planetary Science Division's Research Program. E.T.W., R.K., and J.H.M. acknowledge that this work was performed as part of NASA's Virtual Planetary Laboratory, supported by the National Aeronautics and Space Administration through the NASA Astrobiology Institute under solicitation NNH12ZDA002C and Cooperative Agreement Number NNA13AA93A, and by the NASA Astrobiology Program under grant 80NSSC18K0829 as part of the Nexus for Exoplanet System Science (NExSS) research coordination network. E.T.W., J.H.M., and R.K. acknowledge that this material is based upon work supported by the National Aeronautics and Space Administration (NASA) under grant No. 80NSSC21K0905 issued through the Interdisciplinary Consortia for Astrobiology Research (ICAR) program. J.H.M. gratefully acknowledges funding from the NASA Habitable Worlds program under award 80NSSC20K0230. E.T.W. acknowledges that ExoCAM has been developed over many years and over the course of many grants that he has received. Thus, in addition to support from NASA's SEEC, VPL, and ICAR programs noted above, E.T.W. also acknowledges the NASA Earth and Space Science Fellowship award NNX10AR97H, NASA Planetary Atmospheres Program award NNH13ZDA001N-PATM, and NASA Habitable Worlds awards NNX16AB61G, 80NSSC17K0741, 80NSSC17K0257, and 80NSSC20K1421. All projects have provided critical funding for the ongoing development, maintenance, and support of ExoCAM and ExoRT. E.T.W. also acknowledges that this was his first, and perhaps last, lead author paper using Overleaf, and after this experience E.T.W. has doubled down on the thesis that Overleaf and LATEX are overhyped and overrated; however, the coauthors know this thesis to be incorrect.