TRAPPIST Habitable Atmosphere Intercomparison (THAI) Workshop Report

Preliminary results of the THAI intercomparison (Fauchez et al. 2020) have been published. Four atmospheric compositions have been simulated by the four GCMs (see Appendix A.1). The simulated atmospheres include two dry (no surface liquid water) benchmark cases, "Ben1" and "Ben2," presented by Martin Turbet (https://www.youtube.com/watch?v=B8a2-G8NmmA), with the atmospheric compositions of 1 bar of N₂, 400 ppm of CO₂, and a pure 1 bar of CO₂, respectively, and two moist habitable cases "Hab1" and "Hab2," with a global ocean and the same respective atmospheric compositions, presented by Thomas J. Fauchez (https://www.youtube.com/watch?v=kYLbp_BrJFs). The overall outcome of the Ben1 and Ben2 cases is a good agreement between the four GCMs in terms of surface and atmospheric fields. However, we note some differences in the circulation regime, which manifest a sensitivity of TRAPPIST-1e simulations to GCM setup due to a combination of planetary parameters, noted, e.g., in Sergeev et al. (2020). Synthetic transmission spectra have been produced by the Planetary Spectrum Generator (PSG; Villanueva et al. 2018), and they are in good agreement between the models as long as the top of the atmosphere is extended (assuming an isothermal atmosphere and fixed gas mixing ratio) up to about 100 km (10⁻⁷ and 10⁻¹⁰ bar for Ben1 and Ben2, respectively). Without this extension, the strongest absorption features of CO₂ are truncated. Concerning the Hab1 and Hab2 cases, clouds are the largest source of discrepancies between the models, as expected, due to differences in convection (Section 4.3), bulk condensation, cloud microphysics, boundary layer, and other parameterizations, and their coupling with atmospheric dynamics and radiation. The altitude and thickness of the cloud deck at the terminators impact the simulated transmission spectra, leading to about 20% differences between the models on the number of transits required to detect those atmospheres with the James Webb Space Telescope (JWST) at a 5σ confidence level. More details on the Ben1 and Ben2 simulations, Hab1 and Hab2 simulations, and the impact on observable transmission spectra and thermal phase curves will be presented in three follow-up papers. We also welcome other modeling groups to join THAI at any time. Interest has been shown from the THOR (Deitrick et al. 2020) and Isca (Vallis et al. 2018) GCM groups that we hope to host and compare results soon, while the Exeter Exoplanet Theory Group (EETG) will contribute results from the Unified Model's (UM) replacement, LFRic (Adams et al. 2019), in the future.

Once JWST data for TRAPPIST-1e are available, GCM output will be compared to observational data. We expect such comparison will lead to a new set of simulations and to further validation of model performance against terrestrial exoplanet data. It is important to maintain and improve the level of collaboration between the exoplanet GCM community (including THAI) and the observational community. The "TRAPPIST-1 JWST Community Initiative" (Gillon et al. 2020) is particularly relevant as it aims to develop a coordinated framework to study TRAPPIST-1 planets with JWST from both the observational and theoretical/modeling levels.

Yang et al. (2016) compared the differences in 1D radiative transfer calculations among two line-by-line codes (SMART and LBLRTM), a moderate resolution code (SBART), and four low-resolution codes that are used in GCMs (CAM3, CAM4_Wolf, LMD-G, and AM2). Note that CAM4_Wolf would eventually become ExoCAM; see Appendix A.1.1. The atmospheric composition was set to 1 bar of N₂, 376 ppmv CO₂, and variable H₂O. They showed that there are small differences between the models when the surface temperature is lower than about 300 K. At higher temperatures, such as 320–360 K, the differences between the models could be tens of watts per square meter. The source of the differences is mainly from water vapor radiative transfer calculations in both shortwave and longwave. The differences are larger for shortwave than longwave and also for an M-dwarf spectrum than a solar spectrum. These results suggest that radiative transfer codes should be verified first (such as the absorption and continuum behavior of water vapor) before being used in an exoplanet GCM especially when targeting planets with hot climates or to estimate the inner edge of the habitable zone. Notably, an important lesson learned from this study is that the adequate performance of shortwave radiative transfer for warm moist atmospheres is contingent upon sufficiently resolving the near-IR H₂O spectral absorption bands and windows, particularly when considering irradiation from M-dwarf stars.

Yang et al. (2019c) compared the differences of 3D GCM simulation results on a rapidly rotating aqua planet receiving a G-star spectral energy distribution and on a tidally locked aqua planet receiving an M-star spectral energy distribution. Several GCMs were considered, including various versions of CAM, the LMD-G, and the AM2 GCM (Anderson et al. 2004) They found a relatively small difference (<8 K) in global mean surface temperature predicted by various GCM for cloudy planets orbiting a G star but rather large differences (20–30 K) for cloudy planets orbiting M stars. These discrepancies are due to differences in the atmospheric dynamic, clouds, and radiative transfer. Clouds are the largest difference between the models. The interactions between radiative transfer (such as shortwave absorption by water vapor) and atmospheric circulation can influence the atmospheric relative humidity and therefore affect the surface temperature.

Mars—Several GCM intercomparison efforts were arranged in the Mars atmosphere modeling community as soon as enough teams could contribute to such a project. They were organized in advance of two workshops that helped structure the community just like the 2020 THAI workshop did, although no official reports were published in the literature. The first meeting was the "Mars GCM intercomparison workshop" organized at Oxford University, United Kingdom, 1996 July 22–24. It was later followed by the first "Mars atmosphere modeling and observations," which took place in Granada (Spain) on 2003 January 13–15. In both cases, instructions were sent to the different teams to prepare comparable simulations in advance, and volunteers analyzed the simulations by comparing the zonal mean fields and the predicted planetary waves. Radiative transfer models were also compared for specific test cases (as announced by Harri et al. 2003). These efforts showed that, most of the time, significant differences could always be attributed to different settings in the models and that it was difficult to reach profound scientific conclusions from such organized intercomparisons. Mars GCM teams later focused on comparing their models with the numerous observational data sets that became available in the 2000s. Nevertheless, specific intercomparison studies have continued being conducted to study phenomena or altitude ranges for which not much observational data are available (e.g., González-Galindo et al. 2010; Wilson et al. 2014) or to compare model predictions at a mission landing site (e.g., Kass et al. 2003; Newman et al. 2020).

Venus and Titan—Similarly, intercomparison campaigns have been organized in the Titan and Venus GCM communities, which reached a sufficient stage of development a few years after the Martian case. In practice, because on Venus and Titan the problem of superrotation is so striking and challenging for the GCMs, the most interesting comparisons actually focused on the behavior of the various dynamical cores and their ability to simulate superrotation and conserve angular momentum (Lebonnois et al. 2012, 2013; see also Section 4.2.1). These studies revealed that various dynamical cores, which would give very similar results in Earth or Mars conditions, can predict very different circulation patterns in Venus-like conditions. Recently, a detailed Titan climate model intercomparison has been organized, with the motivation of preparing the planned Dragonfly Titan lander mission (Lora et al. 2019), like intercomparison studies done for Mars to prepare for a mission landing. On the basis of the acceptable agreement between the different models, the authors conclude that the "low-latitude environment on Titan at this season is now fairly well constrained," which is reassuring when preparing an ambitious mission like Dragonfly.

Snowball Earth—A noticeable GCM intercomparison has been presented in Abbot et al. (2012) concerning the impact of clouds in the snowball Earth deglaciation. Six different GCMs were on board the intercomparison. They found that clouds could warm a snowball Earth enough to reduce the amount of CO₂ required for deglaciation. But because the amount of clouds varies from one model to another, the amount of CO₂ required differs by one order of magnitude depending on the model. This intercomparison highlights clouds as an important source of discrepancies between the GCMs.

In planning future community-based exoplanet model intercomparisons ("ExoMIPs") like the THAI project, it is useful to consider the 25 yr history of the future Earth-focused Coupled Model Intercomparison Project (CMIP). This discussion has been presented at the THAI workshop by Linda Sohl.

CMIP is perhaps best known now for its contributions to the periodic assessments issued by the Intergovernmental Panel on Climate Change (see, e.g., IPCC report 2013, https://www.ipcc.ch/report/ar5/wg1/), but it began in 1995 as an independent project of the World Climate Research Program (WCRP). Over the years, CMIP has grown from an effort to use simple global coupled ocean–atmosphere GCM experiments, with interactive sea ice and land models, to separate natural climate variability from human-induced climate change (Covey et al. 2003), to the exploration of a variety of sophisticated climate change scenarios, using GCMs with ever-more complex capabilities that include chemistry and higher-resolution dynamical interactions (Eyring et al. 2016), as well as specialized ancillary investigations (e.g., the various projects under the Paleoclimate Modelling Intercomparison Project, or PMIP ³¹ ). The history of the CMIP experience highlights four key considerations that would benefit any future exoMIP effort:

Context/Rationale: Establish how the MIP will advance the state of knowledge and/or the state of the art. CMIP's overall experiment designs are developed as an outgrowth of the WCRP Grand Challenges, ³² which are updated periodically via community input. These Grand Challenges, which encompass observational, theoretical, and modeling-based research, are meant to:

• 1.  
  Identify the key research questions needing to be addressed in order to move the field forward in a substantive way, as well as the barriers to progress (what do we need to learn next, and what stands in the way?);
• 2.  
  Define effective and measurable performance metrics (how will we know we have been successful in achieving our goals?);
• 3.  
  Provide storylines that engage a broad interested audience, from the media and general public to scientists from other disciplines (how can we attract future talent and improving interdisciplinary connections?).

Planetary science and astrobiology do not have an internationally defined set of grand challenges as such. However, documents such as the NASA Astrobiology Strategy ³³ and AstRoMap European Astrobiology Roadmap ³⁴ outline research topics of interest for advancing the field, and not surprisingly, there is overlap that can help narrow the context and rationale of an exoplanet MIP. An ExoMIP should strive to make connections with as many of these topics as is plausible: linking the MIP rationale to broad themes of community-wide interest is to our advantage in connecting to our fellow researchers whose focus is on theory, field work, or observations.

Experiment Design: Encourage broad participation by planning core experiments with low entry barriers for most groups, with specialized subprojects as needed. One of the benefits of MIPs, from a model development standpoint, is that comparisons across multiple models can illustrate which model design/parameterization approaches provide the most robust results. Thus, the more models we can encourage to participate in a given MIP, the better for the community. The participation of any one modeling group in an MIP is going to be limited by three factors: the technical requirements of the MIP (how intensive is the setup process for the experiments?), the available resources (how much computing time is available?), and expertise on hand (are there enough people with the necessary model knowledge and time to run all the experiments?).

CMIP has taken the approach of defining a limited set of required core experiments that are easy to implement (e.g., Eyring et al. 2016). This lowers the entry barrier to groups interested in joining and improves the chances of useful outcomes beneficial to the community. More complex scenarios, specialized topics, and extended parameter space are addressed through related MIPs, ³⁵ which attract groups that have the additional resources, interest, and expertise available. For an ExoMIP, it is conceivable that a similar design could focus on global-scale core experiments with simple forcing changes (stellar insolation, relatively thin atmospheres/oceans), and some of the related MIPs could engage with 1D/EBM models on details of atmospheric composition/radiative forcing.

MIP Logistics: Plan realistic schedules for experiment completion and group analyses/manuscript preparation. Keeping model groups focused and reaching project milestones in a timely fashion are important for the overall success of an MIP. CMIP and the related MIPs typically establish these schedules via community planning workshops, where ideas for experiments and additional "rules of engagement" regarding MIP participation are also defined in advance. These MIP protocols should then be published as close to the official start of a MIP project as possible so that additional groups not involved in the planning workshops can also make a timely decision to participate.

Data Sharing: When/how to release experiment results for broadest impact? On this topic, CMIP and some of the related MIPs—mainly the PMIP projects—handle data sharing differently. Because the groups contributing to CMIP are frequently not entirely the same as researchers conducting multimodel evaluations of the experiment results, data are released immediately to the community. In contrast, the specialized PMIP project experiments are often run and evaluated by the MIP participants themselves. In the latter case, a data embargo is often declared until the first group papers are published, as part of the agreed-upon project schedule. An ExoMIP might want to consider a similar data embargo as part of a community agreement not to publish each other's work prematurely.

Data-sharing logistics are a more complicated issue. No ExoMIP will have the vast resources currently available to CMIP for data sharing, so at present any data sharing is likely to happen on an ad hoc basis. However, it is possible to develop community standards for what and how data should be shared. While raw model output and some accompanying post-processing scripts might provide maximum flexibility to fellow modelers, model results that have already been postprocessed into file formats such as netCDF for map views and plain text for line plots allow the broadest possible audience—from fellow scientists who are not modelers, to educators and students—to work with and learn from the output with the help of free apps.

The beginning of the first network of exoplanet model intercomparison is a great time to address organizational challenges, and it should be a top priority for the community.

In order to advance intermodel comparisons for exoplanet study, a first-order requirement is a collaboration workshop for the community to discuss key issues, especially regarding data to be shared (and how to share it), which would be very important for establishing standards and best practices for intercomparisons going forward. A formal intercomparison workshop could be organized for roughly one year after the THAI workshop (fall 2021), which would allow time for planning the workshop and getting funding support to encourage participation.

We can expect that an intercomparison workshop would produce documentation about best practices for model intercomparisons and, most importantly, produce a community consensus on how to share data so that a common repository (to be identified) would host a comprehensive set of common diagnostic outputs that are most important for addressing the science questions asked (e.g., which diagnostic attributes contribute most to the synthetic spectral signatures of interest for observational programs?) and yet do not overwhelm potential users or the repository itself with ancillary diagnostics/large data files.

To compare GCMs, even to one another, we must necessarily address a range of cases (by bracketing single-point cases with increasingly complex physics included sequentially, etc.). This requires a clearly defined question (e.g., "is this climate sensitive predominantly to clouds or surface water reservoir effects?") with the goal of producing concrete test cases (that is, with an exhaustive list of clearly stated parameters) and testable predictions so that it is possible to differentiate possible states. Having a clear, well-defined question, in conjunction with a concise list of simulations and deliverables, is necessary to ensure that requirements for participation can be met within realistically allocatable work efforts for such projects. Transparency with respect to set parameters also ensures reproducibility and provides diagnostic access for future tests and observations. Future observational modes further constrain important test deliverables (wavelength coverage and spectral resolution if spectroscopic, or dynamic/geologic/phenomenological diagnostics; spatial scales; arrival times; duration of observation or mission lifetime, etc.). If the models are too far away from one another, then how do each of the GCMs motivate those disparate results, and does this lead to additional testable predictions?

Given the abundance of "hidden" parameterizations in models, it can be difficult to assess whether a given model in a given part of parameter space matches observations for the "right" reason (e.g., the same physical driver as the primary control in both the model and the planetary environment) versus a confluence of other effects. Further intercomparison work can elucidate some of these factors, but it is unlikely that they will all become explicit dependencies, even with additional documentation. Note that potential new observables outside current capabilities and ones where additional precision would refine parameter ranges also help to steer future mission development, which then feeds back into how close the models are to ground truth. Continued intercomparisons require continued funding, and given the dependence of adequate documentation, this suggests the need for clear funding lines for the development of testing frameworks (either for a single model or as part of a new or ongoing collaboration), validation, documentation, etc. GCM model development work is costly, and national laboratories like NCAR or GISS typically hire software engineers to support the scientists. The exoclimate community is young, and the scientists generally lack such support, hampering progress. Lastly, reducing the model output to a single metric (e.g., generic climate state or spectral signature produced) is an additional constraint (metasensitivity?), suggesting that disagreeing models may "agree" in some diagnostic sense. This would help to identify additional directions to explore, such as broadening the parameter space identified initially, or through secondary observables (diurnal/seasonal variability, etc.).

Finally, we would also recommend that potential intercomparison contributors think beyond the goal of "what does this planet simulation look like from the observation perspective?" Three-dimensional GCMs in particular produce a wide variety of diagnostic outputs that are interesting and relevant to understanding the potential habitability of a particular world configuration as well as model performance, but these do not necessarily produce directly observable results. This is especially true knowing that some observables that we have currently identified may be in fact unobservable and that new ones will eventually be found later. Current modeling should therefore not be only constrained by the current set of possible observables.

In the upcoming era of JWST, it becomes essential to focus community effort on benchmarking/comparing/validating the performance of exoplanet climate models, both with respect to other models and to observations (when available). As noted in Section 2.2.1, model intercomparisons have been widely used for decades by the Earth science community in this way, as a very valuable means to improve model reliability, mitigate model dependencies, track down bugs, and provide benchmarks for new models. While individual intercomparison projects should have their own clearly defined protocols, the exoplanet community would benefit also from a metaframework—essentially, a framework for designing model intercomparison projects. This metaframework is what we propose with CUISINES. This framework would be open from 0D to 3D models as well as radiative transfer models and not limited only to rocky exoplanets. One of the first steps in establishing this metaframework will be to create a CUISINES committee and then to prepare a workshop on best practices for intercomparisons.

At the end of the THAI workshop, two future intercomparisons were already discussed: one between EBMs for ice belts and one between GCMs for cloud-free mini Neptunes (see Section 2.2.4 below). In the era of JWST, mini Neptunes and hot Jupiters in particular will require focused modeling efforts from the community. Other ideas for intercomparisons under the CUISINES metaframework are welcome.

A future THAI-equivalent GCM intercomparison project for mini Neptunes has been proposed during a breakout discussion. Mini Neptunes are the most abundant category of exoplanets that have been discovered so far and thanks to their larger size they will be more easily characterized through transmission spectroscopy with JWST. For now, there is a wide range of approaches considering mini Neptunes modeling so there is a significant need for an intercomparison to see what the differences are when everyone makes the same assumptions (more important for now than thinking about the more complicated physics we need to implement). For instance, aerosols are very challenging to include in mini Neptunes (Charnay et al. 2015a, 2015b, 2021); therefore, it has been suggested clear-sky simulations are considered as a start. Cloudiness has been shown to decrease with decreasing equilibrium temperature (Crossfield & Kreidberg 2017), which motivates the case for relatively cold temperatures (K2-18b; Benneke et al. 2019; Tsiaras et al. 2019; and colder), where less photochemical haze is expected. Gliese-436b could be a very good candidate with respect to the amount of data potentially available and the quality of constraints on the planetary parameters (Demory et al. 2007; Gillon et al. 2007; Lanotte et al. 2014; Ehrenreich et al. 2015; Bourrier et al. 2016; dos Santos et al. 2019, and references therein). Atmospheric compositions should be limited to common gases expected for these planets such as hydrogen, helium, water, methane, and carbon dioxide and surface pressures could range from few millibars to tens of bars. However, in order to take into account deep atmosphere effects on the upper atmosphere, which have been shown to be important for Titan (Lebonnois et al. 2012), it may require inclusion of pressures up to 10–100 bar which can subsequently increase the computational time (Wang & Wordsworth 2020). In a next step, the atmospheric compositions can be refined to include more processes to match with observations that JWST would have provided. Also, other models such as EBM and 1D radiative-convective models could be engaged, too.

Most GCMs require numerical diffusion or a filter that is applied in addition to the existing terms of the Euler or primitive equations (Lauritzen et al. 2011, Chapter 13). This mechanism typically serves two practical purposes, which are intimately related: providing numerical stability and achieving a kinetic energy spectrum that is consistent with our understanding of turbulent cascades.

In the case of numerical diffusion, diffusivity is generally a tunable parameter. For Earth GCMs, the diffusivity can be tuned to achieve a spectral slope that matches empirically measured values (Nastrom & Gage 1985; Lauritzen et al. 2011). For exoplanets, there is little hope of measuring the kinetic energy spectrum, but we can use the expectation of a −3 power law (Charney 1971) or −5/3 power law (Pope 2000; Cabanes et al. 2020) as guidance.

More specifically, the turbulent cascade causes energy to build up at the grid scale, well above the level at which molecular viscosity would act to convert this energy to heat (see Lauritzen et al. 2011, Figure 13.7). A numerical diffusion term is thus usually included and selected to have a form that preferentially diffuses the fields at the smallest scales in the model by using, for example, iterated Laplacian operators (see, e.g., Spiga et al. 2020, Appendix A.2). However, selecting the strength and form of this term is somewhat of an art, as it will depend on the resolution, time-step size, solver, grid, and numerous other factors (Lauritzen et al. 2011; Thrastarson & Cho 2011). Fortunately, most exoplanet observables are relatively insensitive to the strength of numerical diffusion (Heng et al. 2011a; Deitrick et al. 2020). Nonetheless, the exoplanet GCM community should bear in mind that other properties may be sensitive to ad hoc numerical settings. For instance, differences have been noted between the GCM prediction of Titan superrotation possibly due to numerical diffusion (Newman et al. 2011).

One numerical issue deserves further attention: the need for so-called "sponge layers," that is, enhanced diffusion or drag near the model top and/or bottom. This need arises because GCMs typically use reflecting boundary conditions, allowing waves (usually gravity waves) to be reflected back into the model domain. These reflected waves are unphysical and can amplify and trigger numerical instabilities (Lauritzen et al. 2011). Thus, an additional drag mechanism is often used to eliminate these reflections. Various types of sponge layer exist, for example, Rayleigh friction, which directly damps wind speeds toward zero or another value such as the zonal mean (see, for example, Mayne et al. 2014b; Mendonça et al. 2018). This type of sponge is easy to implement but is nonconservative (for instance, terrestrial studies by Shaw & Shepherd 2007 show that sponge layers adversely impact the angular momentum balance, thus the simulated circulations). It is instructive to note that some GCM simulations of solar system gas giants do not employ a sponge layer so as to avoid altering the angular momentum balance of the atmosphere (Schneider & Liu 2009; Liu & Schneider 2010; Spiga et al. 2020). Another commonly used sponge layer, which is conservative in finite-volume formulations, reduces the order of the numerical diffusion, from (for example) fourth order to second order (Lauritzen et al. 2011). In either case, it is not always clear how to tune the strength and size of the sponge layer, which needs to damp waves without strongly affecting the general circulation or inducing additional reflections. While sponge layers have been carefully calibrated for Earth simulations, the effects of sponge layers and reflected waves arguably deserve more attention in exoplanet atmospheres. Indeed, the exact settings required for these various damping mechanisms are currently unknown, as there is little constraint from observations, although in cases physically motivated (e.g., capturing dissipation from subgrid eddies or emulating the propagation of waves into space) their main use is for numerical stability (see Heng et al. 2011b for an example of the dissipation and maximum wind speed).

There has been some debate in the literature on the sensitivity of hot-Jupiter simulations to initial conditions (Thrastarson & Cho 2010; Liu & Showman 2013). Other recent work has hinted at the possibility that zonal wind speeds on these planets may be sensitive to the initial temperature–pressure profile used in the deep atmosphere (Sainsbury-Martinez et al. 2019). More investigation should be done on initial conditions in exoplanet GCMs, although it is a challenge to explore many possibilities with such computationally expensive models.

The atmospheres of exoplanets present new territories in which physical processes may be unfamiliar and poorly constrained, compared to Earth and other solar system bodies. Unlike bodies in our solar system, for exoplanets we have little, if any, spatial information on the atmospheric structure. Thus in modeling these atmospheres we must utilize any and all available criteria to ensure physical realism.

The dynamical cores of GCMs are formulated using conservation laws (e.g., the Euler equations). As such, the global conservation of properties such as mass, energy, and angular momentum provides a diagnostic of the model's performance and accuracy. Thuburn (2008) provides some guidance on which properties may be conserved and the desirable degree of conservation. We reiterate a few of those concepts here but more details on the dynamical cores will be given in Section 4.2.

First, as pointed out in Thuburn (2008), while the continuous forms of the equations of fluid dynamics can be formulated to conserve all physical properties, the discrete forms of the equations do not. Choices must be made regarding which properties to conserve (to numerical precision). One example of this is the thermodynamic equation (i.e., the first law of thermodynamics), which can be written in terms of different variables, such as potential temperature, pressure, internal or total energy, etc. Many GCMs use potential temperature, which is convenient for modeling convective processes. This leads to a conservation law for entropy, whereas sometimes a conservation law for energy (total or internal) may be preferred (Satoh 2002). A similar choice must be made in the momentum equations, as these may be written in terms of linear momenta, angular momenta, vorticity, etc.

Conservation of mass is particularly important, as noted by Thuburn (2008), because it affects all other conservation laws, and should be robust in the absence of significant sources and sinks (e.g., escape to space or volatile freezeout). In other words, errors in mass conservation will lead to a cascade of errors elsewhere.

As we develop GCMs further and explore novel territory, conservation also provides a critical way to identify coding errors. In his talk, Russell Deitrick briefly discussed using mass and angular momentum conservation to identify bugs in the THOR GCM. Finite-volume models (such as THOR) should conserve properties naturally to roughly machine precision because the equations are discretized in flux form—fluxes flow across boundaries such that control volumes on either side experience the exact same flux. Model discretized in other ways (e.g., spectral models) may ensure conservation by use of "fixers" (see Lauritzen et al. 2011, Chapter 13).

Traditionally, GCMs have used a latitude–longitude (lat-lon) spherical grid, which is easy to construct and has operators that are well known and intuitive. It does, however, suffer from singularities and resolution clustering due to the convergence of meridional lines at the poles (Staniforth & Thuburn 2012). Many models have solved this issue using a combination of semi-implicit time integration and numerical filters, for example, the UM.

Another solution is to use an alternative horizontal grid structure. A comprehensive review of grid types, and their advantages and disadvantages, is provided in Staniforth & Thuburn (2012). The most commonly used of these seem to be the cubed sphere, used in versions of the FMS, for example, Lin (2004), and the icosahedral grid, used in NICAM (Tomita & Satoh 2004), THOR (Mendonça et al. 2016), and DYNAMICO (Dubos et al. 2015). These quasi-uniform grids succeed in making the resolution more uniform across the sphere, avoiding numerical complications at polar regions in the lat-lon grid. These grids also scale very well with a high number of processors, which is usually not the case for the lat-lon grid (Staniforth & Thuburn 2012). However, these grids are a true challenge to work with—the divergence, gradient, and curl operators must be written using Gaussian integrals on the icosahedral grid, for example (Tomita & Satoh 2004). Further, they are not completely uniform and thus still admit the possibility of grid imprinting, wherein errors build up to a level that makes the underlying grid visible to the eye (Staniforth & Thuburn 2012). Also, the core utilization of icosahedral grid is far below the one of the lat-lon grid. Unfortunately, there is no known "perfect" grid, so an understanding of the shortcomings of a particular grid is essential, as is a comparison between models utilizing different grids.

A further word of caution is warranted here: some models, as in Dobbs-Dixon & Lin (2008), have avoided numerical issues with polar regions by omitting them entirely. While some features of the resulting simulations may be qualitatively reasonable, it is likely that the polar regions are particularly crucial for the circulation of tidally locked exoplanets.

Dissipation and accuracy—When GCMs are used to model more "extreme" planets, such as hot Jupiters, changes in the physical conditions can lead to reductions in the stability and potential accuracy of the simulation results. As discussed, GCMs rely on several forms of dissipation in the model (Jablonowski & Williamson 2011) to control subgrid "noise," which can lead to model instability. These can take several forms, from a diffusion of the winds themselves, "filtering" over polar regions in lat-lon grid GCMs and so-called "sponge" layers (see Section 4.1.2, paragraph on sponge layers)

Rayleigh drag and closing the angular momentum budget—The question of the conservation of axial angular momentum is paramount in atmospheric modeling, as is illustrated for instance in studies of slow-rotating bodies and gas giants in the solar system (Lebonnois et al. 2012; Spiga et al. 2020, their appendix). Dynamical cores are not explicitly formulated to conserve axial angular momentum, and this can cause spurious variations of modeled angular momentum that can range from negligible to major, as was evidenced, e.g., in the case of exoplanet modeling by Polichtchouk et al. (2014). This question of angular momentum conservation is all the more critical in planets without a solid surface: a simple Rayleigh drag on horizontal winds is used as a bottom boundary condition in GCM studies of Jupiter and Saturn to emulate the closing of the angular momentum budget in simulated jets by a putative magnetic drag at depth (Liu & Schneider 2010; Young et al. 2019; Spiga et al. 2020). In the simulations of hot Jupiters, too, Rayleigh drag, used as a simple parameterization of surface drag on horizontal winds, has been used to capture the impact of magnetic drag in the deep atmosphere of hot Jupiters (Perna et al. 2010). However, in a similar fashion as the above-mentioned dissipation, its main use is for stability as by dragging the deep atmosphere to immobility where the radiative timescale is long, or indeed infinite (Iro et al. 2005), one can remove dependency of the simulated results on the initial conditions (see Mayne et al. 2014b; Amundsen et al. 2016; Tremblin et al. 2017; Sainsbury-Martinez et al. 2019 for various discussions on this issue). Early simulations with the MITGCM demonstrated a loss of global axial angular momentum without this inner boundary drag (Cho et al. 2015), and similar effects have been found for both UM and THOR. However, often the angular momentum conservation may vary with the spatial and temporal resolution adopted for model simulations; this opens the possibility of obtaining negligible changes in axial angular momentum with adjustment to the spatial and temporal resolution of the particular setup. The cause of the angular momentum conservation issue in dynamical cores is still not clearly understood.

Thermodynamics—In most of the existing GCMs for planetary atmospheres, molecular weight gradients, heat capacity, gravity, etc., are not taken into account. Those quantities are therefore assumed constant throughout the whole atmosphere and throughout time. This assumption can have several effects: when the thickness of the atmosphere becomes large with respect to the radius (such as for Titan or mini Neptunes), the mass of a given atmospheric cell should change through buoyancy, due to the change of gravity. A constant or variable value of the heat capacity at a constant pressure C_p would directly impact the stability profile of the atmosphere. Such effects play on atmospheric dynamics and therefore on the equilibrium between the different terms of the atmospheric equations, but are generally assumed to be second-order effects (except when the variations are really strong). Taking into account this variability into GCMs would require subsequent development and time, and this endeavor would benefit both solar system planets and exoplanets. Preliminary works on how to take into account variations of C_p with the temperature have already been presented in Lebonnois et al. (2010) and Mendonca & Read (2016).

Shocks—For some simulations the flow speed can approach, or exceed, a Mach number of 1, leading to some authors questioning whether shock-capturing solutions to the continuity equation are required (Li & Goodman 2010; Fromang et al. 2016). However, as shown by Fromang et al. (2016), shocks play a minimal role in atmospheres dominated by a large-scale superrotating jet. As is understood for the solar system's gas giants, flows of conductive material in the presence of a background magnetic field can lead to drag (as mentioned earlier) and heating through "ohmic dissipation" (see, for example, Ginzburg & Sari 2016). In hot Jupiters, the outer layers can become ionized, leading to magnetic drag, and the deeper layers potentially experience significant enough ohmic heating to alter the planetary radius. To date, the only simulations consistently capturing these impacts have been those of Rogers & Komacek (2014) and Rogers & Showman (2014), which revealed that ohmic heating is unlikely to be significant enough, for reasonable magnetic fields. Magnetic drag, aside from the work of Rogers & Komacek (2014) and Rogers & Showman (2014), has otherwise been captured through parameterized drag schemes (Perna et al. 2010).

Substellar objects—As more and more planets have been detected and other observable populations identified, either through discovery or improvements in instrumentation, the use of GCMs for gas-giant extrasolar objects has widened. Brown dwarfs, substellar objects that form similarly to stars, share parameter space with gas-giant planets in terms of bulk compositions (to some extent) and radii, leading to similarities in atmospheric circulations in the two kinds of objects (Showman et al. 2019). Additionally, due to high levels (relative to Jovian planets) of interior convection, these objects are self-luminous, and high-signal to noise data is available (Buenzli et al. 2015). GCMs have been applied to these objects exploring their flows and additionally the occurrence of clouds (e.g., Zhang & Showman 2014; Tan & Showman 2021). The same studies can be applied to young gas-giant planets, with high interior convection (e.g., directly imaged planets). The main challenges are in handling the strong interior fluxes from the convection and the extremely short rotation periods. Work is beginning on coupling models of the convective interior of gas-giant planets to atmospheric models to better capture the interaction between these two regimes. Irradiated brown dwarfs, with a partner star from either the main sequence or white dwarf (Casewell et al. 2018), have also recently been studied using an adapted GCM (Lee et al. 2020). Through studying this collection of gas-giant objects, from young self-luminous Jovian exoplanets to older short and long-period Jovian planets, and isolated or heavily irradiated brown dwarfs, a complete continuum of atmospheric regimes could be unraveled.

Adaptations related to exotic thermodynamics and chemistry—Observations of hot Jupiters have also begun to demarcate this subclass itself into further categories. In particular, ultra-hot Jupiters (with temperatures in excess of ∼2500 K) provide some real advantages while presenting new challenges. Although the relatively high temperatures result in the assumption of chemical equilibrium holding over most of the observable portion of the atmosphere, effects of the high temperature and photon fluxes, such as thermal and photodissociation, and the resulting H⁻ opacity must be included (Baxter et al. 2020). Most significantly, these objects span a temperature range in which hydrogen is present in both molecular and atomic forms (Bell & Cowan 2018; Tan & Komacek 2019). The variations this causes in the specific heat capacity are large enough to mean the standard assumption of a single value through the atmosphere, made within GCMs, may become problematic. The resolution of this issue requires a significant reworking of the dynamical cores developed under the assumption of constant heat capacity.

Specific challenges for Mini Neptunes—The drive in instrumentation is pushing toward the detection and characterization of smaller radii and longer orbital period planets, ultimately in the search for potentially habitable planets (the focus of this workshop). However, the next set of observational facilities and instruments will likely provide access to the subclass of planets discussed previously, termed mini Neptunes or super-Earths. The existing challenges listed in Section 4.2.1 clearly apply to those objects. Nevertheless, for these planets, initial work has shown that the standard primitive equations of motion often employed within a GCM may not be valid (Mayne et al. 2019), and/or elapsed simulation times must be significantly extended (Wang & Wordsworth 2020). Additionally, as temperatures cool, the role of photochemistry and condensation may become even more important, but for a range of species.

As highlighted by Leconte et al. (2017), a background gas lighter than the condensible gas (for example H₂O and H₂) can induce a mean molecular weight gradient that inhibits convection in the atmosphere. This process is stronger on giant planets such as mini Neptunes but could be observed on smaller rocky planets. Therefore, it may be interesting to explore this phenomenon with 3D simulations.

Convection, and moist convection especially, is an important driver of heat redistribution in planetary atmospheres. By forming clouds and depending on the boundary layer processes, moist convection is also a key part of complex feedback mechanisms in the climate system (e.g., Arakawa 2004). To accurately resolve convective plumes, a numerical climate model has to have a sufficiently high spatial resolution, making it extremely computationally expensive and thus infeasible for long simulations or multiplanet studies. Thus, all modern exoplanet GCMs rely on parameterizations to emulate the overall effect of subgrid-scale convective processes on large-scale atmospheric fields. These parameterizations always include a number of quasi-empirical parameters, usually inherited from Earth climate models or validated against observations and convection-resolving simulations on Earth, raising the question of their applicability to extraterrestrial atmospheres.

For Earth's atmosphere, an assumption of a dilute condensible is usually a good approximation (Pierrehumbert 2010). This is not the case for other planetary atmospheres in the solar system and beyond: the main condensible species can comprise a substantial portion of the atmosphere or have thermodynamic properties such that the convective mass-flux (Ooyama 1971) is sufficient to affect the large-scale dynamics, like on Pluto (Bertrand et al. 2018). In such a scheme, the depth of the convection is constrained by the distance the air, rising in convective updrafts, penetrates above its level of neutral buoyancy, and ceases to rise. To simulate these effects correctly, the LMDG is equipped with a convection parameterization accounting for nontrace condensible species (Pierrehumbert & Ding 2016). It is applicable to a wide range of atmospheric conditions and is described in Leconte et al. (2013b).

Whatever the convection parameterization is based on—the adjustment to reference profiles, subgrid-scale mass-flux, or other principles—it has to be validated against observations and convection-resolving simulations. In the absence of in situ observations for exoplanets, the next best option is to use high-resolution convection-resolving and cloud-resolving models (CRMs), which simulate convective processes explicitly. Targeted limited-area CRM simulations can be used to benchmark and improve convection parameterizations. For instance, Abbot (2014) has compared CRM simulations to GCM simulations in the case of a snowball Earth and found that they provide consistent results. This helped to confirm the hypothesis that clouds could provide a large warming on a snowball Earth and potentially lead to the deglaciation of the planet. Two talks at the THAI workshop presented work paving the way in this direction. Denis Sergeev demonstrated substantial differences between cloud profiles in a global coarse-resolution GCM experiment and a limited-area CRM experiment, conducted using the UK Met Office Unified Model (Appendix A.1.4) for the THAI Hab1 & Hab2 setups. Differences in convective cloud cover appear to be model and planet dependent because the opposite picture has been found by Maxence Lefèvre, who compared LMD-G GCM (Appendix A.1.2) to the (Weather Research and Forecasting) WRF model (see appendix in Skamarock & Klemp 2008) in a CRM mode for a case of convection on Proxima Centauri b. Further work will hopefully build on Sergeev's and Lefèvre's CRM simulations to explore convective processes in other atmospheric and planetary regimes, informing atmospheric modelers of parameterization biases and caveats.

It was a general consensus among the participants that a flexible convection parameterization based on the mass-flux approach is the best option for coarse-resolution 3D GCMs, while the adjustment scheme is usually too crude to represent convection. A shift toward fully resolved global convection simulations is not going to happen quickly, but limited-area CRMs should be used more actively in the exoplanet atmospheric modeling. As outlined in Section 4.3, one of the main applications of CRMs is to retune existing convection parameterizations for different extraterrestrial atmospheres. Such experiments are routinely performed for Earth weather and climate prediction models (Rio et al. 2010) and for Mars as well (Colaïtis et al. 2013), so the exoplanet GCM community should work closely with meteorologists and Earth model developers to benefit from their invaluable expertise. In practical terms, an important and already feasible project is building an archive of convection-resolving simulations of H₂, N₂, CO₂-dominated atmospheres, which then can serve as a standard benchmark suite for coarse-grid global models. This project can then evolve into an exoplanet CRM model intercomparison project and find its rightful niche under the CUISINES umbrella (Section 2.2.3). With a wide grid of CRM models at hand, certain aspects of convective processes can then become more tractable, such as whether the structure and dynamics of convective plumes and precipitation anvils in non-Earth planetary atmospheres might change. How might they be influenced by plume size, convective overshoot beyond the level of neutral buoyancy, and entrainment and detrainment as timescales of convection change when the composition of the atmosphere or stellar/planetary parameters (such as stellar spectrum and planetary gravity) are changed?

In addition, Earth science expertise is valuable for the development of generalized convection schemes from scratch. In this case, code developers should strive to make them flexible, modular, and portable so that it will be relatively easy to swap one convection parameterization in the GCM for another. The fact that most convection schemes usually operate column-wise without communicating with neighboring columns in a 3D GCM makes the issue of portability easier to tackle.

There are many ways in which a continental surface can affect a planet's climate, including (but not limited to) surface albedo, topography, thermal inertia, surface roughness, etc. Depending on the nature of the land surface, these parameters could change significantly and thus alter the planet's climate (Madden & Kaltenegger 2020).

The most extreme configuration in which continental surfaces play a major role is "land planets." Land planets, sometimes also known as dune planets, are desert rocky planets with limited surface water (Abe et al. 2011). It is thought that this type of planet is one of the most probable (along with water-rich ocean planets) around M stars, as a result of formation and escape processes (Tian & Ida 2015). These types of planets have already been studied with 3D GCM simulations (Abe et al. 2011; Leconte et al. 2013a; Menou 2013; Yang et al. 2014; Kodama et al. 2019; Way & Del Genio 2020) and specifically applied to the TRAPPIST-1 planets (Wolf et al. 2017; Turbet et al. 2018; Rushby et al. 2020). Rushby et al. (2020) recently provided a detailed analysis of how surface type/composition may affect the climates of TRAPPIST-1 assuming they are land planets. This work was also described in the presentation given by Aomawa Shields.

While the study of continental surfaces of Earth and other planets of the solar system (in particular Mars, characterized by its hypercontinental climate) is today our primary source of information on this matter, characterizing the climate of the planets of the TRAPPIST-1 system and other nearby planets (e.g., Proxima b) may provide us with crucial data on how continental surfaces are operating on alien worlds.

Ocean modeling is often overlooked by the exoplanet community, largely due to the large computational expense associate with spinning up dynamic ocean models, coupled with the challenges of observing an ocean on another planet (e.g., Robinson et al. 2010). Future exoclimate simulations of terrestrial worlds benefit from the use of a dynamic ocean component (Way et al. 2018; Yang et al. 2019b) coupled with a dynamic sea ice model. It has been shown that sea ice drift can alter the habitable zone limits in various cases (Hu & Yang 2014; Way et al. 2017, 2018; Del Genio et al. 2019). Near the inner edge of the habitable zone, Leconte (2018), Yang et al. (2019b), and Salazar et al. (2020) have shown that ocean heat transport is not always necessarily critical, especially when continents are present. Warmer planets (T_S > 300 K) tend to have more homogeneous surface temperatures and thus ocean currents may not cause a meaningful net change in ocean–atmosphere heat exchanges. Yang et al. (2019b, p. 29) show this to be the case and further state that "...ocean dynamics have almost no effect on the observational thermal phase curves of planets near the inner edge of the habitable zone. These results suggest that future studies of the inner edge may devote computational resources to atmosphere-only processes such as clouds and radiation." However, Yang et al. (2019b) also argue that ocean heat transport is critical for the climate and observables for "middle habitable zone" planets in M-dwarf systems. Still further, Way et al. (2018) demonstrated that planets with modern Earth-like land–sea masks show marked differences in mean surface temperature for fast rotators (spin less than eight Earth sidereal days) versus slow rotators (see Way et al. 2018; Figure 2) in their inner edge of the habitable zone studies.

In addition, ocean composition—specifically, ocean salinity—may affect the fraction of habitable surface area and phase curves of middle and outer habitable zone worlds (Cullum et al. 2016; Del Genio et al. 2018; Olson et al. 2020). Salinity (dissolved salt content) is a first-order control on the density of seawater, and it further modulates the relationship between temperature and density. Salinity thus influences ocean stratification, circulation, and heat transport. At the same time, salinity depresses the freezing point of seawater, potentially limiting the formation of sea ice in salty oceans with consequences for surface albedo. In sum, saltier oceans tend to result in warmer climates. Models that ignore ocean dynamics cannot currently simulate salinity impacts on OHT but may include freezing point depression. However, the relative contribution of heat transport versus freezing point depression to climate warming, and how the balance of these effects may differ under different climate states, is not well understood. It is thus unknown when/if simple hacks such as adjusting the freezing point of seawater in a model without a dynamic ocean is a reasonable strategy for simulating the climates of exoplanets with unknown ocean salinity or whether the likelihood that exo-oceans differ from Earth's ocean with respect to salinity requires the inclusion of a dynamic ocean.

To facilitate future studies, the exoplanet modeling community should consider working more closely with oceanographers. For example, there is a large parameter space of ocean tidal dissipation to be explored that affects planetary rotation rates over time (e.g., Green et al. 2019). As mentioned above, rotation rate has been demonstrated to affect climate. It must also be stressed that current GCMs used for exoplanetary studies have serious shortcomings in some cases. First, the putative thermodynamic oceans (also called q-flux (Miller et al. 1983; Russell et al. 1985) as a heat source "q," whose values are generally specified by a control run, is prescribed to represent seasonal deep water exchange and horizontal ocean heat transport) used in exoplanet GCMs have generally used zero horizontal heat transport or highly simplified parameterizations (e.g., Edson et al. 2012; Godolt et al. 2015; Kilic et al. 2017). In addition, current GCM dynamic ocean models are presently highly parameterized for modern-day Earth because they are the children of Earth parent GCMs. For this reason it is important to engage more closely with the oceanography community to better parameterize the current suite of dynamic oceans used in exoplanetary GCMs.

For ocean planets that have a low density, the ocean depth can reach tens or even hundreds of kilometers. At the bottom of the ocean, ice under high pressure may form. For the deep ocean, the equation of state for seawater is required to be changed. Moreover, the ocean-bottom ice can influence the friction and the exchange of heat and materials between the ocean and the solid planet. Key questions may be answered using ocean GCMs: How deep is ocean circulation (including both wind driven and thermal driven), and how does ocean circulation influence the concentrations of CO₂ and other greenhouse gases in the atmosphere (Checlair et al. 2019)?

Besides seawater oceans, another type that needs to be investigated are magma oceans. The lowest temperature of a magma ocean is about 1600 K. The density, viscosity, and diffusivity of the magma ocean are quite different from that of Earth's ocean. However, they may still of the same order. So, in this respect an Earth ocean GCM may be easily modified to simulate the circulation of a magma ocean. But, a key process for the magma ocean is silicate (or other materials) precipitation in the ocean, which acts like vertical convection and can significantly influence the heat and mass transports in the ocean.

Gravity waves—The challenge of current exoplanet GCMs is to overcome the expense of running simulations with the necessary horizontal and vertical resolution over a range of pressures that adequately resolves processes operating over a variety of spatial and temporal timescales. But if our focus is data driven, we ought to look at the processes important in Earth's (and other worlds') middle atmospheres. Earth, Venus, and Titan are all terrestrial worlds with significant superrotation in their middle atmospheres. Venus' slow (and counter) rotation makes it a hallmark case study for tidally locked exoplanets with substantial atmospheres. It has superrotation in the equator that is likely driven by the Gierasch–Rossow–Williams mechanism (Gierasch 1975; Rossow & Williams 1979), where planetary waves from high latitudes can transport angular momentum toward the equator and spin up superrotating jets. Titan too has been observed to have significant variability in winds at different altitudes spanning from the stratosphere to the lower thermosphere, and it too is a slow rotator at roughly 16 days. Earth has an alternating stratospheric jet oscillation, known as the quasi-biennial oscillation (QBO), which is the product of a complex interaction from a broad spectrum of waves. Stratospheric oscillations are also found to occur in giant planets of the solar system (Fouchet et al. 2008) and are a recent focus for climate modeling (Cosentino et al. 2017; Bardet et al. 2021). The point is that their upper atmospheric dynamics are driven by waves not easily resolved and their effects are mostly missing from current exoplanet models while potentially having strong observational effects.

Small-scale or high-frequency gravity waves have a large role in the upper atmospheric dynamics of terrestrial-size planets. Not only are they an important source of momentum to drive the QBO, they also affect midlatitude jet streams, the semiannual oscillation, and even the Brewer–Dobson circulation (Butchart 2014). They are an efficient means of transporting energy across latitudes and altitudes and effectively redistribute energy, either mechanical or thermal, away from its source to other areas until equilibrium or relaxed/balanced states are reached. Gravity waves try to keep the upper atmosphere of Earth in dynamical and radiative equilibrium by dispersing energy spatially over time. Gravity wave activity is also confirmed in Venus' and Mars's middle/upper atmosphere by several measurements and claimed to produce the observed variability in density, temperature, and cloud structure (e.g., Creasey et al. 2006; Garcia et al. 2009; Altieri et al. 2012). Linear wave theory has allowed parameterizations of the effects from gravity waves and their breaking action in Earth-based GCMs for some time. The typical approach is to assume some characteristics of the properties of the waves, amplitude or momentum flux, wavelength, phase speed, etc. that are inputs into the parameterization that is applied depending on the modeled local atmospheric environment, with a focus on horizontal and vertical wind shear and temperature gradients. Nevertheless, given the lack of systematic observations of gravity waves and the uncertainty of the source (for instance, on Mars and Venus), which are necessary to constrain model parameters, our experience with different GCM configurations let us conclude that the total zonal wind (i.e., averaged for all local times) value is very sensitive to many GCM quirks. Zonal wind in the middle/upper atmosphere can be either positive or negative, producing different circulation regimes.

Non-LTE effect in the upper atmosphere—The upper atmospheres of terrestrial-like planets in our solar system are similar in terms of physical processes, in spite of important differences in temperature, density, and composition (Gladstone et al. 2002). A basic property of these upper layers is the low gas density that is also responsible for situations of the breakdown of LTE, specific for each molecular species and each vibrational transition. These non-LTE effects result in populations of molecular energy states not dictated by Boltzmann statistics at the local kinetic temperature and occur when molecular collisions are so infrequent that other processes (e.g., radiative transfer) become important for the determination of those states' number population (Lopez-Puertas & Taylor 2001). In terrestrial planets, and for the main molecules and infrared emissions, those layers usually correspond to their mesosphere and thermosphere. At pressure layers above about 10⁻⁵ mbar, solar EUV heating and thermal conduction are the main processes controlling the energy balance, while in the mesosphere (on terrestrial planets between 1 and 10⁻⁵ mbar, approximately), absorption and emissions by atmospheric molecules with active rovibrational bands in the IR usually play a crucial role on the thermal structure (Gladstone et al. 2002).

These non-LTE processes have to be considered when interpreting strongly irradiating exoplanets. CO₂ and CO non-LTE fluorescence is common in telluric atmospheres, but CO emission has also been detected in Neptune (Fletcher et al. 2010). Furthermore, non-LTE radiative transfer modeling helped to explain unexpected observed features around 3.3 um on the hot-Jupiter HD 189733b from ground measurements, reported to be CH₄ non-LTE emission (Swain et al. 2010; Waldmann et al. 2012). Future detection by JWST, LUVOIR, together with ground-based measurements of the upper atmosphere of hot Jupiters by IR spectrographs (CRIRES/VLT, METIS/E-LTE) will make it possible to test composition and temperature models of warm and hot Jupiters. However, it is still challenging for terrestrial exoplanets.

Due to its expected significant effect on exoplanet atmospheric characterization, the improvement of the mid- and upper atmosphere processes in exoplanet atmospheric models should therefore be a priority in the era of JWST and ELTs.

The transport of photochemically produced gaseous species and hazes could have a significant impact on the characterization of exoplanet atmospheres (Carone et al. 2018). Different circulation regimes can result in different global distributions of important atmospheric species such as ozone (Yates et al. 2020; Chen et al. 2019). The community should first think about what kind of composition is the most interesting and imperative for near-term observations (e.g., exo-Venuses). Very large uncertainties remain for many deposition fluxes. For instance, CO and relatively minor tracers like Cl have an important catalytic activity for both Venus and Mars but are not accurately constrained. Along a similar vein, the surface emission fluxes of various biogenic compounds such as dymethylsulfide (DMS) are unconstrained, but different assumptions can dramatically affect their resultant global distributions (Chen et al. 2018).

Exoplanet observations of many regimes show that clouds and hazes significantly affect planetary spectra. Current GCMs provide the spatial mapping of water clouds, which can have strong effects on transmission spectra and thermal phase curves (Wolf et al. 2019). However, the majority of GCMs do not include a self-consistent treatment of photochemical haze. As an integral part of climate modeling, aerosols need to be included—either from the ground up (production rates from chemical networks) or decoupled (aerosols and photochemical hazes would be separated). For instance, in modeling Titan (and Titan-like exoplanets; Lora et al. 2018), the production rate is fixed to reproduce observations, and then the photochemistry is separated. For exoplanets, this would require the use of very detailed models to capture monomer formation self-consistently, which will be used in a simple photochemical model favoring the approach of a lower-resolution model but with a higher-complexity component. As a follow-up to Chen et al. (2019), one possibility is through adapting CARMA (Larson et al. 2015), a state-of-the-art microphysical model that can be used here to simulate the evolution of hazes. Note that CARMA is already coupled to ExoCAM with a nominal fractal aggregate haze model (Wolf & Toon 2010), with funded plans to use it for studying hazy iterations of the habitable zone planet TOI-700 d. Haze production rates will be sourced from offline 1D Atmos calculations (e.g., Arney et al. 2017). Lastly, the 3D temperature structure of a planet could also be important for the chemistry itself (exo-/endothermal reactions), even in the absence of photochemistry.

It has to be noted that for exoplanet photochemistry, the use of precalculated photolysis rate tables is much too specific. A rate table of reasonable size is unlikely to be generic. For an atmosphere in which we do not know the profile of the species, we cannot build a coherent table. It is necessary to go through an online calculation of the photolysis rates that leads to a reconsideration of the radiative transfer within the GCM. Because most GCMs use the correlated-k method, it cannot be used for photolysis calculations. In short, it is necessary to calculate the photolysis online. However, simulating the photochemistry and chemistry with 3D models is computationally very expensive. It may require prior work to reduce the number of gaseous families, which would depend on the bulk composition, temperature, and instellation. Such an approach has been started by Olivia Venot's group but it is uncertain how 3D transport will affect this optimization, which is performed in 1D.

Another caveat is that the majority of terrestrial 3D chemical models are restricted to present-day Earth compositions. This is due to the fact that these models largely inherited components from their model supersets originally developed for Earth-based research. For instance, Chen et al. (2019) adapted the National Center for Atmospheric Research model (Marsh et al. 2013) by deploying subroutines from ExoCAM with the Whole Atmosphere Community Climate Model. Thus, another goal is to extend 3D photochemistry models to nonoxygenated reducing and weakly oxidized (H₂-, N₂- and CO₂-rich) atmospheres. Such anoxic conditions dominate the atmospheric evolution histories of Earth, Mars, Venus, and Titan. This suggests that anoxic atmospheres are the "default" state of a planet's eventual fate in the absence of biological activity. As shown by previous 1D work (Lincowski et al. 2018), assessing non-Earth-similar atmospheres is important to gauge the detectability of photochemical byproducts on the myriad of potential rocky exoplanet compositions.

While out of the scope of this report, it is important to mention that star–planet interaction is fundamental to understand a range of key processes (stellar-sourced charged particles, EUV heating, photodissociation, and photoionization, etc) that could shape the upper layers of rocky exoplanet atmospheres. Collaboration with stellar physicists and observers will be crucial.

Characterization of nonwater condensibles is extremely important for understanding both the atmosphere and surface processes of other worlds.

In the low-temperature range, gases such as CO₂, H₂S, and SO₂ can condense (Fray & Schmitt 2009) for planets with hydrogen-dominated atmospheres and volcanic activity, while CH₄, NH₃, and N₂ are expected or found for worlds like Venus, Titan, Pluto, and other solar system moons. This condensation can have large consequences for the planet's atmosphere by removing greenhouse gases, forming clouds, and modifying surface albedo. For instance, Turbet et al. (2017b) showed that CO₂ condensation can strongly reduce the deglaciation of terrestrial planets as CO₂ condensation leads to the accumulation of surface CO₂ ice that can get permanently trapped under water ice. In the high-temperature range, a variety of condensibles would potentially be observable by JWST or future ELTs for highly irradiated exoplanets, but data on properties (e.g., microphysical properties) of these condensible species are sorely lacking.

Spectral information would be needed from laboratory experiments and missions regarding optical properties of exotic clouds and albedo properties for nonwater condensibles to improve model response to a broad variety of scenarios.

In addition to acquiring new data, model developments are needed within GCMs in order to include condensible species other than water and to precipitate condensibles with appropriate albedo and grain-size properties. In this area, the LMD-G GCM already has extensive capabilities for CO₂ condensation on Mars (e.g., Forget et al. 1998), early Mars (e.g., Forget et al. 2013), and exoplanets (e.g., Wordsworth et al. 2011; Turbet et al. 2017b); N₂ condensation on early Titan (e.g., Charnay et al. 2014); and N₂, CH₄ and CO condensation on Pluto (e.g., Forget et al. 2017).

Overall, much more data are needed to make significant progress on the question of condensible species in exoplanet's atmospheres and surfaces (in particular, for high-temperature condensibles), along with potentially complicated and time-consuming model development. Initial work coupling various schemes with multiple condensate species has been done for hot Jupiters (e.g., Lee et al. 2016; Lines et al. 2018a, 2018b, 2019).

Aerosols are present in every atmosphere of our solar system planets and moons. Clouds have also been observed in exoplanet's atmospheres such as the super-Earth GJ-1214b (Kreidberg et al. 2014), the gaseous giant WASP-12b (Wakeford et al. 2017), and WASP-31b (Sing et al. 2016). Hazes have been observed on WASP-6b (Nikolov et al. 2015) and HAT-P-12b (Sing et al. 2016). They have not been observed for terrestrial-size exoplanets yet but simulations using GCMs and PSG have shown that they dramatically flatten the transmission spectrum, preventing an exhaustive atmospheric characterization from space observatories (Fauchez et al. 2019; Komacek et al. 2020; Suissa et al. 2020).

However, the detailed aerosol microphysics is uncertain for exoplanet atmospheres. Yet, changes in the microphysical and optical properties can have a very large impact on the climate simulation and simulated spectra. To improve our understanding of aerosol properties in exoplanet atmospheres, we need data from the lab but also from JWST and ARIEL. If a statistical study is performed on such data it may allow us to discriminate between cloud particles and hazes that differ in terms of size and microphysical/optical properties. Also, linear polarization is a powerful tool with which to retrieve cloud microphysical/optical properties as is performed for Earth using, for instance, POLDER/PARASOL data (Goloub et al. 2000). Finally, we have to improve our connection with other communities such as Earth scientists, paleoclimatologists, and solar system planetary scientists to better share data (remote sensing + in situ) and methods of applying data to exoplanet atmospheres.

In addition data modeling studies can also be very helpful. For instance, a sensitivity study on cloud microphysics can be performed by varying cloud particle size and amount of cloud condensation nuclei to simulate their impact on simulated spectra. Also, modeling work could allow us to better understand how spatial heterogeneity (horizontal and vertical resolution, cloud gap fraction, and overlap) affects both transmission and reflection spectra. These improvements may necessitate the inclusion of biological coupling for both haze (microbial, plant) and fires and therefore to develop computational physics packages, chemistry packages, and increased capabilities for capturing spatial and temporal heterogeneity.

More direct collaboration needs to be done with the climate modelers and observation simulators to bring consistency between assumptions about radiative transfer. More observations of Earth transmission spectra are therefore needed, for example, the Atmospheric Chemistry Experiment—Fourier Transform Experiment (ACE-FTS) on board the SCISAT-1 satellite provides observation in the 2.2–13.3 μm window. This emphasizes the benefit of Earth and Planetary Science synergies.

The HEXTOR energy balance model (EBM) was used to conduct the THAI scenarios, using both a latitudinal and longitudinal mode (J. Haqq-Misra & B.P. Hayworth 2021, in preparation). HEXTOR is a one-dimensional EBM based on the model by Williams & Kasting (1997). The model is typically run in a latitudinal mode, which reproduces Earth's mean annual climate. The model can also be used to explore changes in Earth's climate due to past and future orbital variations and possible feedback from anthropogenic forcing (Haqq-Misra 2014). Prior versions of this model have represented radiative transfer with a basic linear relationship (Haqq-Misra 2014) or with a polynomial fit of 1D radiative-convective climate calculations (Williams & Kasting 1997; Batalha et al. 2016; Haqq-Misra et al. 2016; Hayworth et al. 2020). The current version of HEXTOR attempts to improve the accuracy of the radiative transfer in the model by using a lookup table, which conducts a nearest-neighbor interpolation for outgoing longwave radiation (OLR) and albedo using a database containing thousands of 1D radiative-convective climate calculations. This provides an advantage in accuracy at the cost of added computational expense. During the THAI workshop, Dr. Haqq-Misra showed that HEXTOR in a latitudinal configuration either underestimates or overestimates the global average temperature in the THAI simulations, because the hemispheric differences between the day and night sides cannot be represented with a single dimension in latitude; however, he also showed that HEXTOR can also be configured as a longitudinal EBM through a coordinate transformation, which places the substellar point at the north pole and allows the day-to-night side contrast to be represented more accurately (Fortney et al. 2010; Koll & Abbot 2015; Checlair et al. 2017, J. Haqq-Misra & B.P. Hayworth 2021, in preparation). Longitudinal EBMs, either along the equator like HEXTOR or with full lat-lon resolution, can provide constraints on climate across broad parameter spaces or for long time integrations, which can be useful in identifying specific problems to study further with GCMs.

VPLanet (Barnes et al. 2020) includes an EBM called POISE (Planetary Orbit-Influenced Simple EBM), a one-dimensional seasonal EBM that reproduces Earth's annual climate as well as its Milankovitch cycles (see North & Coakley 1979; Huybers & Tziperman 2008; Deitrick et al. 2018). Though the model lacks a true longitudinal dimension, each latitude is divided into a land portion and a water portion, with distinct heat capacities and albedos, and heat is allowed to flow between them. Ice can accumulate on land at a constant rate when temperatures are below 0 °C while melting/ablation occurs when ice is present and temperatures are above 0 °C. Sea ice forms when a latitude's temperature drops below −2 °C (accounting for salinity) and melts when higher. To account for ice-sheet flow, bedrock depression, lithospheric rebound, and ice-sheet height, they employ the formulations from Huybers & Tziperman (2008). The bedrock depresses and rebounds locally in response to the changing weight of ice above, always seeking isostatic equilibrium. POISE is thus a self-consistent model for ice-sheet growth and retreat due to instantaneous stellar radiative forcing, orbital elements, and rotational angular momentum.

One-dimensional radiative-convective climate and photochemical models are 1D models representing a vertical atmospheric column assuming plane-parallel waves in hydrostatic equilibrium. In the photochemical models, the vertical transport takes into account molecular and eddy diffusion and are able to represent a complex photochemistry. One-dimensional models have been widely used by the community to determine the edges of the habitable zone (Kopparapu et al. 2013) and to study the ancient Earth (Arney et al. 2016, 2017) and various exoplanets (Lincowski et al. 2018; Meadows et al. 2018). In this workshop, THAI simulations with the Atmos 1D model (Wunderlich et al. 2020) were presented by Andrew Lincowski following a two-column approach. Dr. Lincowski has shown that two 1D radiative-convective atmospheric columns are able to reproduce the day–night temperature contrast simulated by GCMs while keeping an advantage in terms of computational time (Lincowski et al. 2021, in preparation in this focus issue).

GCMs are very complex models that require significant time to converge. Lower-dimensional models such as EBMs or 1D radiative-convective climate and photochemical models, while ideal to explore large parameter sweeps, lack a representation of atmospheric dynamics, surface heterogeneity, and clouds. The computational efficiency of EBMs enables them to simulate climates on much longer timescales of thousands or millions of years to explore orbital and rotational effects on climate (e.g., Spiegel et al. 2009; Deitrick et al. 2018). EBMs are typically 1D in latitude and solve a single partial differential equation for surface temperature. Temperature then depends on incoming stellar flux (instellation), heat diffusion, albedo, and the outgoing longwave radiation (OLR). The OLR and albedo are parameterized with simple formulations (North & Coakley 1979; Spiegel et al. 2009; Rose et al. 2017; Palubski et al. 2020), though several studies have made advancements by fitting polynomials to radiative-convective models (Williams & Kasting 1997; Haqq-Misra et al. 2016). The chief challenge of these models comes from the parameterization of atmospheric dynamics in terms of a heat diffusion term and accuracy in parameterizing the radiative transfer. For this reason, synergy with GCMs is necessary to ensure some measure of accuracy and predictive power from EBMs.

One-dimensional EBMs and radiative-convective climate models coupled with photochemical models can explore a very large parameter space (i.e., star and planet properties, instellation, rotation and orbital periods, eccentricity, atmospheric properties, etc.) and can identify key points of interest in the parameter space that 3D models can then investigate. For instance, 1D models can be used to determine the likely chemical state as input to a GCM. GCMs can also be used to determine cloud coverage percentage and dynamical model as input to 1 and 1.5D models. GCMs can be run with simple tracer chemistry with haze precursors and 1D photochemical models can be used to figure out what happens next with haze formation/chemistry etc. Using both 1D and 3D models simultaneously would allow one to get a more complete picture of chemistry, clouds, and observables.

It has also been highlighted during the THAI workshop that interactions between the atmosphere and interior of terrestrial planets require more attention than currently given. This would require improving the collaboration with geologists/geophysicists. Such coupling should probably first be developed in 1D following a "planet evolution model" approach based on an asynchronous coupling employing a mixture of (short term) climate calculations and long-term simulations (for glaciers, but also longer processes as well).

Finally, it is important to predict in advance, with a hierarchy of models, what we might see and have the models ready to interpret the data. Ideally, upstream modeling work should not be constrained by anticipated observational sensitivity.

ExoCAM is an exoplanet branch of the Community Earth System Model (CESM) version 1.2.1. CESM is provided publicly by the National Center for Atmospheric Research in Boulder, CO (http://www.cesm.ucar.edu/models/cesm1.2/), and ExoCAM is freely available on GitHub (https://github.com/storyofthewolf/ExoCAM). To use ExoCAM, the user must first obtain CESM v1.2.1, and then ExoCAM is installed as a patch on top of the core CESM code. ExoCAM was developed by E.T. Wolf to facilitate accessible configurations for exoplanet and planetary modeling and is now used by several different research groups in the community. The ExoCAM code package includes model configurations, initial condition files, source code modifications, and an accompanying flexible correlated-k radiative transfer model, ExoRT https://github.com/storyofthewolf/ExoRT). ExoRT can be run coupled to the 3D model or in a standalone 1D mode and has several supported gas absorption schemes. Typically, ExoCAM is run utilizing the cloud and convection physics from the Community Atmosphere Model (CAM) version 4 (Neale et al. 2010), and a finite-volume dynamical core (Lin & Rood 1996); however, ExoCAM can be configured to leverage other CESM-supported dynamical cores (e.g., spectral element cubed sphere) and physics routines as desired (e.g., CAM5, CARMA). Likewise, ExoCAM is most often run using a 4° × 5° horizontal resolution and 40 vertical atmospheric layers up to 1 mbar pressures; however, ExoCAM can easily be run with other supported model resolutions (e.g., Wei et al. 2020) and model tops (e.g., Suissa et al. 2020) with relative ease. For the THAI simulations, ExoCAM was run with 4° × 5° horizontal resolution, 51 vertical layers extending to 0.01 mbar pressures, configured with CAM4 cloud and convection physics, and with the ExoRT radiation scheme originally developed for Archean Earth atmospheres described in Wolf & Toon (2013). ExoCAM, coupled to ExoRT, has been used to study a variety of problems including deep paleoclimates for Earth (Wolf & Toon 2013, 2014), stellar and CO₂-driven moist greenhouse climates (Wolf & Toon 2015; Wolf et al. 2018), the climate of Earth-like exoplanets around solar-type stars (Wolf et al. 2017; Adams et al. 2019; Kang 2019a, 2019b, 2019c), tidally locked exoplanets around M-dwarf stars (Kopparapu et al. 2017; Komacek & Abbot 2019; Yang et al. 2019a; Komacek et al. 2019, 2020, 2020; Rushby et al. 2020; Suissa et al. 2020; Wei et al. 2020; Zhang & Yang 2020), and Earth-like planets in circumbinary systems (Wolf et al. 2020).

The LMD-G GCM—or the LMD Generic model—is a 3D Global Climate Model historically developed at the Laboratoire de Meteorologie Dynamique (LMD) in Paris, France. The model originally derives from the LMDz 3D Earth (Hourdin et al. 2006) and Mars (Forget et al. 1999) Global Climate Models, but it benefits from the lessons learned by developing GCMs for most atmospheres in the solar system, where models can be tested against a wide range of observations (Forget & Lebonnois 2013). It solves the primitive equations of geophysical fluid dynamics using a finite-difference dynamical core on an Arakawa C grid. The LMD-G GCM is equipped with flexible radiative transfer (based on the correlated-k method; Wordsworth et al. 2011) and thermodynamics/cloud microphysics packages with the objective of being able to simulate any cocktail of atmospheric gases (as long as spectroscopic data sets are available) and aerosols. In particular it can account for the condensation of both minor and major constituents of an atmosphere. Most planetary (planet size, mass, rotation period, topography, etc.) and stellar parameters (insolation, input spectrum) can be easily adapted for a wide range of planets. LMD-G has been used in many climate studies for the past and future climates of solar system planets (Forget et al. 2013; Wordsworth et al. 2013; Charnay et al. 2013, 2014; Leconte et al. 2013b; Turbet et al. 2017a, 2017b, 2020) and exoplanets in a wide range of conditions (Wordsworth et al. 2011; Leconte et al. 2013a; Charnay et al. 2015a, 2021; Bolmont et al. 2016; Turbet et al. 2016). It was specifically adapted and used for the TRAPPIST-1 planets in Turbet et al. (2018) and Fauchez et al. (2019). More information on the model (code, user manual, tools, publications) can be found on http://www-planets.lmd.jussieu.fr/. Recently, the flexible physical parameterizations of LMD-G have also been interfaced with LMD's next-generation icosahedral dynamical core DYNAMICO (Dubos et al. 2015), which is particularly suitable for massively parallel architectures. DYNAMICO has been used to perform high-resolution simulations of the solar system's giant planets (Cabanes et al. 2020; Spiga et al. 2020; Bardet et al. 2021) and carries many promising perspectives for exoplanet studies (Sainsbury-Martinez et al. 2019).

The Resolving Orbital and Climate Keys of Earth and Extraterrestrial Environments with Dynamics (ROCKE-3D) is a GCM developed at NASA Goddard Institute of Space Studies (GISS) (Way et al. 2017). ROCKE-3D is based on its parent Earth climate GCM GISS ModelE2 (Schmidt et al. 2014), which is used for the Coupled Model Intercomparison Project Phase 6 (CMIP), currently in its sixth version. For the THAI intercomparison ROCKE-3D version Planet_1.0 was used. ROCKE-3D Planet_1.0 was run at an atmospheric horizontal resolution of 4° × 5° with 40 vertical atmospheric layers. The atmospheric model top was 0.1 mbarc (∼60 km altitude). Because ROCKE-3D is an extension of its parent Earth model, it brings along many features of the parent model including river and underground runoff, ground hydrology for different soil types, and a dynamic lakes mode where lakes can either accumulate or dissipate depending upon the competition between evaporation and precipitation. The Planet_1.0 version of ROCKE-3D used in the THAI intercomparison is extensively documented in Way et al. (2017), where one can find a more detailed description of its capabilities, features, and limitations. One important area where ROCKE-3D differs from its parent model comes from its use of a completely different radiative transfer scheme called SOCRATES (see Appendix A.1.4) which offers far more flexibility than the default GISS scheme. At the same time, SOCRATES is more computationally demanding than the default GISS scheme. Whereas the GISS scheme was designed to be extremely fast, it is only available for use with modern Earth atmospheric pressures and gas mixing ratios. ROCKE-3D coupled with SOCRATES has been used in a variety of climate studies for solar system planets through time (Way et al. 2016; Del Genio et al. 2018, 2020; Way & Del Genio 2020) and beyond (Way & Georgakarakos 2017; Kane et al. 2018; Way et al. 2018; Aleinov et al. 2019; Colose et al. 2019; Del Genio et al. 2019; Olson et al. 2020).

The UM has been developed by the UK Met Office and UM Partnership over the last 30 yr with the aim of being able to use the same model for both operational weather forecasting and climate simulation. The UM can be run with a range of planetary parameters, spatial and temporal resolutions, in global (Walters et al. 2019) or regional (Bush et al. 2020) configurations. The UM's dynamical core solves the equations of motion using a semi-implicit, semi-Lagrangian method (Wood et al. 2014), with variables discretized on an Arakawa C grid in the horizontal and a staggered height-based terrain-following Charney–Phillips grid in the vertical. The dynamical core is capable of solving a range of dynamical equations from the most simplified primitive equations to those close to the full nonhydrostatic equations for a compressible fluid (see White et al. 2005; Mayne et al. 2014b). The UM includes sophisticated physical parameterizations for subgrid-scale turbulence, convection, water cloud, and precipitation, as well as radiative transfer, which is solved by the open-source, two-stream, correlated-k code SOCRATES, accessible at https://code.metoffice.gov.uk/trac/socrates and surface and subsurface processes. ³⁶ Note that through the SOCRATES radiative transfer code, the UM is capable of generating synthetic spectra for any given 3D simulation (e.g., Lines et al. 2018a; Boutle et al. 2020).

Adaptation and application of the UM to exoplanets have been led by the Exeter Exoplanet Theory Group (EETG; exoclimatology.com). The UM was initially benchmarked against a range of standard Earth-like planet tests (Mayne et al. 2014a), opening an avenue for studying temperate extraterrestrial atmospheres. Using the UM, various climate processes in the atmospheres of tidally locked rocky exoplanets have been studied, focusing on the impacts of planet eccentricity and atmospheric composition (Boutle et al. 2017), size and location of a substellar continent (Lewis et al. 2018), treatment of convection (Sergeev et al. 2020), host star spectrum (Eager et al. 2020), presence of mineral dust (Boutle et al. 2020), and an interactive ozone cycle (Yates et al. 2020).

The versatility of the UM allowed for its application to a range of gas-giant atmospheres, primarily H-/He-dominated hot Jupiters. After a successful adaptation of the radiative transfer code (Amundsen et al. 2014, 2017), the UM was used to explore flow structures in hot-Jupiter atmospheres (Mayne et al. 2017; Debras et al. 2019, 2020) and to demonstrate that for smaller planets with extended atmospheres (i.e., mini Neptunes; see Section 4.2), the often used primitive equations may not accurately capture the atmospheric dynamics (Mayne et al. 2019). Additionally, a flexible gas-phase chemistry scheme (Drummond et al. 2016) was coupled to the UM, allowing for 3D simulations of H-/He-dominated atmospheres with both equilibrium (Drummond et al. 2018) and kinetic (Drummond et al. 2020) gas-phase chemistry. To parameterize clouds in the atmospheres of hot Jupiters, the UM was also coupled with both a detailed high-temperature microphysics scheme (Lines et al. 2018a, 2018b) and a steady-state simplified cloud scheme (Lines et al. 2019).