A New Line-By-Line General Circulation Model for Simulations of Diverse Planetary Atmospheres: Initial Validation and Application to the Exoplanet GJ 1132b

TL;DR

This paper introduces a novel line-by-line 3D GCM for simulating diverse planetary atmospheres, validates it against existing models, and applies it to study the climate of the exoplanet GJ 1132b, revealing insights into atmospheric composition and circulation patterns.

Contribution

The paper presents the first 3D GCM using a line-by-line radiative transfer approach, enhancing accuracy and flexibility for diverse planetary atmospheres.

Findings

Minor CO2 can significantly alter circulation patterns.

Phase shift of hotspot can reveal superrotating jet profiles.

Oxygen-rich atmospheres' mass may be inferred from phase amplitude.

Abstract

Exploring diverse planetary atmospheres requires modeling tools that are both accurate and flexible. Here, we develop a three-dimensional general circulation model (3D GCM) that for the first time uses a line-by-line approach to describe the radiative transfer. We validate our GCM by comparing with published results done by different 1D and 3D models. To demonstrate the versatility of the model, we apply the GCM to the hot Earth-sized exoplanet GJ 1132b and study its climate and circulation assuming an atmosphere dominated by abiotic oxygen (O$_2$). Our simulations show that a minor CO$_2$ composition can change the circulation pattern substantially, intensifying the equatorial superrotation in particular. Computation of the phase-resolved spectroscopy indicates that the vertical profile of the superrotating jet could be inferred in future spectro-photometric observations by the phase shift of the hotspot in the CO$_2$ principle absorption band centered at 667 cm$^{-1}$. We also show that atmospheric mass could potentially be constrained by the phase amplitude in the O$_2$ vibrational fundamental band for planets with O$_2$-rich atmospheres, although further experimental and/or theoretical O$_2$-O$_2$ collision-induced absorption data at high temperatures is needed to confirm this. More physical schemes such as moist dynamics will be implemented in the GCM in the future so that it can be used to tackle a wide variety of planetary climate problems.

Tables

Table 1. Table 1: Comparison of GCM components between FMS PCM and LMD.

	FMS PCM	LMD
Dynamical core	finite volume	finite difference
Radiation	line-by-line	correlated-k distribution
PBL diffusion	K-profile	Mellor-Yamada 2.5 closure

Table 2. Table 2: Parameters used in the GCM simulations on GJ 1132b.

Parameter	Value
aa Data taken from Berta-Thompson et al. (2015). Planet mass $M\left(M_{\oplus }\right)$	1.62
aa Data taken from Berta-Thompson et al. (2015). Planet radius $r\left(r_{\oplus }\right)$	1.16
aa Data taken from Berta-Thompson et al. (2015). Surface gravity $g\left(g_{\oplus }\right)$	1.20
aa Data taken from Berta-Thompson et al. (2015). Orbital period $P$ (days)	1.63
aa Data taken from Berta-Thompson et al. (2015). Stellar flux $F(1367$ W m-2)	18.64
bb Data taken from Segura et al. (2005). Stellar spectrum	AD Leo
Orbital eccentricity	0.0
Obliquity	0.0
Surface roughness height (m)	$5\times {10}^{-2}$
Surface albedo	0.2
Experiment A	1 bar O2
Experiment B	1 bar O2 + 0.01 bar CO2

Equations

$$S_{ij}(T)=S_{ij}(T_0)\frac{Q(T_0)}{Q(T)}\frac{1-e^{-h{\nu }_{ij}/k_{B}T}}{1-e^{-h{\nu }_{ij}/k_B{T}_0}}\frac{e^{-h{\nu }_i/k_{B}T}}{e^{-h{\nu }_i/k_B{T}_0}}.$$

$$\frac{\partial }{\partial p}\left[\omega \left(c_{p}T+gz+\frac{u^2+v^2}{2}\right)\right]={\dot{Q}}_{dia}$$

$$ϵ\sigma T_a^4$$

$$OL{R}_n$$

Full text

More physical schemes such as moist dynamics will be implemented in the GCM in the future so that it can be used to tackle a wide variety of planetary climate problems.

methods: numerical — planets and satellites: atmospheres — planets and satellites: terrestrial planets

1 Introduction

On terrestrial-mass planets, the atmospheric mass and bulk composition is expected to be extremely diverse, because the external boundary conditions (particularly exchange with the planetary interior and escape to space) dominate the atmospheric evolution (Pierrehumbert, 2010; Catling & Kasting, 2017). The variety of planetary atmospheres in the present-day solar system is one example; the inferred evolution of the Earth’s climate over geological timescales is another (Holland & Turekian, 2014; Catling & Kasting, 2017). For rocky planets outside the solar system, atmospheric characterization will be possible in the foreseeable future and great diversity is also expected (Birkby, 2018; Roberge & Seager, 2018; Madhusudhan, 2018). Investigating these diverse planetary climates requires various tools, including flexible and accurate three-dimensional general circulation models (GCMs).

Probably the most significant obstacle to making flexible GCMs is the radiative transfer calculation in a realistic, vertically-inhomogeneous planetary atmosphere. Numerically computation of the radiative heating rates that drives the flow motion in GCMs is enormously computationally expensive due to the non-gray nature of the atmospheric absorption. This arises from the large number of spectral lines of the absorbing gases , as well as the dependence of the line absorption on the atmospheric state. As a result, statistical parameterizations of radiative transfer have previously been used to make long-term integration feasible (Stephens, 1984; Goody et al., 1989; Lacis & Oinas, 1991). The most popular of these tools, the correlated-k method, has been used for several years in studies of paleoclimate and exoplanet climates by various GCMs, including but not limited to the Laboratoire de Météorologie Dynamique Generic (LMDG) model (Wordsworth et al., 2011, 2013), the Substellar and Planetary Atmospheric Radiation and Circulation (SPARC)/MITgcm (Showman et al., 2009; Kataria et al., 2015), the Extraterrestrial Community Atmosphere Model (Wolf & Toon, 2014, 2015), and the Resolving Orbital and Climate Keys of Earth and Extraterrestrial Environments with Dynamics (ROCKE-3D) (Way et al., 2017). However, statistical tools have limitations , particularly when they are used outside the range for which planetary radiative transfer is well-constrained by observations. Errors and biases can arise if the model is applied to a situation for which it was not designed. Just as importantly when dealing with new temperature/pressure ranges or adding more species, a pre-processing step that relies on external line-by-line softwares to calculate the absorption coefficients is required. This is labor extensive and requires careful choice of approximations, which can limit the flexibility of GCMs based on this technique.

Here we take an alternative approach to this problem and present the development of a new GCM that for the first time uses a line-by-line approach to calculate the radiative transfer directly. The principal advantage of this type of GCM is that our direct algorithm allows high computational accuracy and versatility when varying the atmospheric state and composition. We describe the model framework, including the line-by-line approach to calculate the radiative transfer and validation, in Section 2. As an example application, we use this GCM to simulate the possible climates on an Earth-sized exoplanet GJ 1132b and discuss the associated spectroscopic observable features in Section 3. Future directions of the GCM are discussed in Section 4.

2 Model description and validation

2.1 1D radiative-convective model with the line-by-line radiative transfer calculation

Before implementing the line-by-line radiative transfer calculation in the 3D GCM, we first constructed a one-dimensional radiative-convective (RC) model. This model works in a similar way to the conventional radiative-convective models (e.g., Manabe & Strickler, 1964), except for the spectral radiative calculation. We use 2000-8000 points in wavenumber ; as we discuss later, we have tested that further increases in the spectral resolution have insignificant effects on the spectrally integrated flux and heating rates. Before performing the iteration towards the radiative-convective equilibrium state, the model first reads in HITRAN and/or HITEMP line data and computes the absorption cross-section at these discrete wavenumbers in a temperature-log pressure table for each absorbing species. The values of the temperature-log pressure grids are determined by the initial condition of the surface temperature and surface pressure so that the constructed table is made full use of for long-term integrations. As the model iterates and the atmospheric state changes, new absorption cross-sections are interpolated from this table to calculate the monochromatic optical thickness and the radiative fluxes.

Absorption cross-sections are calculated by first scaling line strengths from their reference values via the formula

[TABLE]

Here $i$ and $j$ denote the ground and excited states of the transition, respectively, $S_{ij}$ is the line strength, $\nu_{ij}$ is the line location in frequency unit, $\nu_{i}$ is the ground state frequency, $T$ is temperature, $T_{0}=296$ K is the reference temperature, $h$ is Planck constant and $k_{B}$ is Boltzmann constant (e.g., Rothman et al., 1998). Line broadening due to both collisional and Doppler effects are calculated via the usual formulae (Goody & Yung, 1995). The line shapes are then calculated using the Voigt profile, with the Humlíček algorithm used to compute the complex probability function efficiently (Humlíček, 1982; Schreier, 1992).

A version of this model using fixed temperature profiles was first described in Schaefer et al. (2016). The time-varying radiative convective model was presented in Wordsworth et al. (2017). To validate the 1D model, we first ensured that it reproduces analytical results when the atmosphere is gray or semi-gray. Next, for H2O runaway greenhouse calculations, we checked that the outgoing longwave radiation (OLR) computed by our model agrees closely with published results for the Earth’s atmosphere (Goldblatt et al., 2013). This inter-comparison is also discussed in Schaefer et al. (2016).

Here, for completeness we present a summary of the model setup and describe an additional validation test we performed – a comparison among various 1D radiative-convective models with different radiative transfer calculations. The simulated case is an Earth-like planet with 1 bar pure CO2 atmosphere orbiting the Sun at 1 AU. The surface has a gray albedo of 0.20. Here the line absorption coefficients for CO2 are calculated from the 2012 HITRAN line list (Rothman et al., 2013), and the lines are truncated at 500 cm*-1*. The spectral range of the longwave calculation is from 1 cm*-1* to 5 times the Wien peak wavenumber of the Planck function at the simulated surface temperature (3000 cm*-1* in this case). 8000 points in wavenumber are used so that the wavenumber increment is roughly 0.375 cm*-1*. The CO2 collision-induced absorption was included using the GBB parameterization (Borysow, 1998; Baranov et al., 2004; Wordsworth et al., 2010). Eight quadrature points (four upwelling and four downwelling) are used for both shortwave and longwave radiative flux calculations. Currently the eight-stream radiative calculation is only performed in the pure-absorption limit. For the scattering calculation, we simply impose the enhanced reflection by Rayleigh scattering at the surface and the analytical two-stream solution is used, because the molecular absorption and Rayleigh scattering occur in well-separated spectral regions. In situations where clouds or aerosols are present or when significant gas absorption extends into the visible, this approximation breaks down and a multiple scattering code is required. We will describe extensions of our model into the multiple scattering regime in subsequent papers.

We compare our 1D model to two others here for the same simulated case: the Clima model developed by Kasting and collaborators that uses the correlated-k distribution method to describe the radiative transfer (Kopparapu et al., 2013) and the linearized flux evolution (LiFE) model recently developed by Robinson & Crisp (2018) that uses a line-by-line approach but combines linear flux Jacobians to rapidly adapt radiative flux profiles to changes in atmospheric state.

A comparison of the radiative-convective equilibrium thermal structures is shown in Figure 1a. All three models yield similar surface temperatures, and hence the vertical temperature profiles in the troposphere where dry convection occurs are also similar. One notable difference is that in the upper atmosphere above 300 Pa, the two models with high spectral resolutions are consistent and both colder by $\sim 10$ K than the Clima model.

Unlike the LiFE and Clima models, which both compute the averaged radiative quantities over spectral intervals , our algorithm is more like a quasi-random sampling of the absorption lines, because we perform the radiative calculation at discrete regular wavenumbers, and the absorption line locations are extremely irregular (particularly for H2O). We show the dependence of our results on spectral resolution in Figure 3a. Interestingly, it is not necessary to use a very high spectral resolution in this radiative calculation. First, the surface temperature is very insensitive to the spectral resolution if the resolution is high enough. The surface is only warmed by 0.8 K when the number of spectral points increases from 1000 to 8000. For the upper atmosphere, which is dominated by narrow spectral bands with strong lines, a relatively high spectral resolution is necessary. For example, given 2000 spectral points, the temperature profile has relatively converged to the results with higher spectral resolutions from 300 Pa to the ground, but above 300 Pa the cooling rate may be underestimated. Slightly more efficient results could probably be obtained using an irregular spectral grid, which is something we plan to investigate in future. However, the convergence of the results in Figure 3a at 8000 points indicates the general robustness of our approach even using a regular spectral grid.

2.2 Three-dimensional planetary climate model (FMS PCM)

Our 3D GCM (hereafter referred to as FMS PCM) is developed based on the “vertically Lagrangian” finite-volume dynamical core (Lin, 2004) of the Geophysical Fluid Dynamics Laboratory (GFDL) Flexible Modeling System (FMS111https://www.gfdl.noaa.gov/fms/), which solves the atmospheric primitive equations in spherical coordinates. The model also includes a cubed-sphere gridding technique for atmospheric dynamics (Ronchi et al., 1996; Putman & Lin, 2007), which improves both the computational performance and accuracy compared to conventional latitudinal-longitudinal gridding. For example, the cubed-sphere dynamical core allows more processors to be used for parallel computation, given the same horizontal resolution, than spectral dynamical cores in which the prognostic variables are decomposed by spherical harmonics and domain decomposition along the circles of latitude for parallel computation is not allowed. This makes the implementation of the line-by-line radiative calculation in the dynamical core feasible given modern computational resources.

The FMS dynamical core only deals with the atmospheric transport, so to build a planetary climate model, additional physical schemes are required. The schemes we implement in FMS PCM are listed as follows:

• •

Our line-by-line radiative calculation, which was introduced in Section 2.1.

• •

A convection scheme. Our scheme is similar in spirit to that described in Manabe & Strickler (1964). The lapse rate in the unstable layer is adjusted towards a neutral state while conserving enthalpy. Currently the model only has a dry version, with no phase transitions of condensable substances or associated cloud formation.

• •

A planetary boundary layer (PBL) scheme. The surface momentum and buoyancy fluxes are computed by the Monin-Obukhov similarity theory from the lowest level winds, temperatures, and tracer mixing ratios. A fully backward time-step method is employed for the vertical diffusion calculation in the planetary boundary layer. The diffusivity is set by the K-profile scheme (Troen & Mahrt, 1986).

To validate this 3D model, we use it to simulate the atmospheric collapse on synchronously-rotating rocky planets. Our motivation for studying this setup is the simplicity of the climate dynamics – a single-component atmosphere dominated by an axisymmetric overturning circulation. We chose not to simulate the Earth’s present-day climate because it is a much more complicated system. The specific case we have chosen is CO2 atmosphere collapse on an Earth-like planet orbiting a M-dwarf (Joshi et al., 1997; Wordsworth, 2015). The detailed parameters of our simulation are given by Table 1 of Wordsworth (2015). In this simulation, the horizontal resolution is C48 ($48\times 48\times 6=13824$ grid boxes in total; similar to a $2^{\circ}\times 2^{\circ}$ Lat-Lon grid); the physical schemes are performed over 26 vertical layers in the hybrid sigma-pressure coordinate. The top pressure level is 220 Pa. As shown in Figure 3a, 2000 points in wavenumber is enough for the radiative calculation for our setup. The equilibrium temperature profile in 1D simulations has converged to the profiles derived from higher spectral resolution runs. Therefore we reduce the spectral resolution to 2000 points in wavenumber to save computation time. The radiative heating rate in the GCM is updated every four model hours. With these parameters, it takes about one hour for a 200-day integration using 384 CPU cores on our cluster. The results we present here are averages over the last 500 days of a 2000-day integration.

The total energy conservation of a GCM is essential for climate simulations. A simple way to check it is to calculate the global-mean top-of-atmosphere (TOA) net radiative flux. In this run, the global-mean outgoing longwave radiative flux (OLR) is 280.28 W m*-2* and the outgoing shortwave radiative flux (OSR) is 61.48 W m*-2*. Given that the global-mean insolation is the same as present Earth (341.77 W m*-2*), net TOA radiative flux is $\sim 0.01~{}$W m*-2*. This also confirms that all of the physical schemes we have implemented in the GCM conserve energy very well. Our simulation agrees closely with the published result simulated by the Laboratoire de Météorologie Dynamique (LMD) Generic Model (Wordsworth, 2015): the pure CO2 atmosphere collapses on the nightside surface when the surface pressure is less than 0.1 bar. Despite the completely different components in the FMS PCM and LMD model, from the dynamical core to the grid-box physical schemes (summarized in Table 1), the FMS PCM produces very similar surface temperature distributions and vertical air temperature profiles (compare Figure 4 with Figure 9 in Wordsworth (2015)).

As discussed in Section 1, broadband averaged parameterizations in radiative transfer calculations allow relatively efficient computation in conventional GCMs. But such models fail to produce flux data with sufficiently high resolution to be useful for interpretation of observations. Usually, offline radiative calculations are required to get the disk-integrated spectrum (e.g., Robinson et al., 2011). The FMS PCM provides another approach to calculate the disk-integrated spectrum because a relatively high-resolution radiative calculation is performed along with the model iteration. As an example, Figure 5 shows the spectrally resolved planetary albedo and outgoing longwave radiation (spectral irradiance) in this CO2 atmosphere collapse simulation. The left panel of Figure 5 is the spectral distribution of the global-mean planetary albedo. The Rayleigh scattering by the CO2 atmosphere becomes significant for wavenumbers higher than 15000 cm*-1* and is well separated from the CO2 absorption region, which validates our simple treatment of atmospheric scattering. The right panel of Figure 5 shows the spectral flux in the thermal infrared emitted from the dayside and nightside. The IR spectroscopic feature contains very useful information regarding the vertical temperature profiles and is consistent with the profiles shown in Figure 4. On the dayside, the air temperature decreases monotonically with height. As a result, absorption troughs emerge in the spectral region where CO2 absorbs the thermal radiation from the lower atmosphere and re-emits with a colder temperature. On the nightside the absorption troughs become emission peaks because of the strong temperature inversion of the lower atmosphere. In particular, the big absorption trough in the CO2 principal band centered at 667 cm*-1* becomes two emission peaks due to the strong CO2 absorption near 667 cm*-1*. The two peak emissions on both sides of 667 cm*-1* correspond to the temperature maximum at $\sim 4000$ Pa on the nightside shown in Figure 4. Note that there is also a clear double-peak feature near 1250 cm*-1* that is formed in a different way. It is related to the double peak of the CO2 collision-induced absorption feature (Baranov et al., 2004). The CO2 absorption in this spectral region is much weaker than that near 667 cm*-1*.

3 Example application of FMS PCM: possible climates on an oxidized GJ 1132b

In this section, we provide one demonstration of FMS PCM and use it to simulate the climates on a nearby exoplanet GJ 1132b, which was recently discovered by the MEarth ground-based transiting planet survey (Berta-Thompson et al., 2015). This rocky planet has many Earth-like parameters, including the radius, surface gravity, composition, and the spin rate assuming synchronous rotation. One major difference is that GJ 1132b receives $\sim$ 19 times more stellar insolation than the Earth and 10 times more than Venus. Schaefer et al. (2016) discussed the possibility of abiotic O2 buildup by the photolysis of water vapor and the subsequent atmospheric escape and mantle absorption during the magma ocean stage. The most common outcome from their simulations are tenuous atmospheres with at most a few bars of O2 and little to no steam remaining. The existence of carbon and nitrogen-bearing species is possible but not included by the chemistry model in Schaefer et al. (2016). The FMS PCM can be used to explore the possible climates on GJ 1132b with various surface pressures and atmospheric compositions. As a demonstration, we consider two simple scenarios: one with 1 bar pure O2 atmosphere; the other with an additional 1% CO2. Other parameters used in the GCM simulations are listed in Table 2.

3.1 Temperature and wind fields

Before running any FMS PCM simulation, we need to choose the appropriate parameters for radiative transfer calculations in the GCM. Figure 6 shows some sensitivity tests. Similar to the CO2 atmosphere collapse simulation discussed in Section 2.2, it is also not necessary to use a very high spectral resolution. For this specific case, even 1000 points in both the OLR and ASR calculation yields an accuracy better than 1%, which is good enough for long-term GCM integrations. Absorption by weak lines tends to be very important in this relatively hot atmosphere. We include all CO2 lines in the GCM runs (the minimal line strength in the HITRAN database is 10*-30* cm*-1*/(molecule cm*-2*)). If a line strength cutoff of 10*-27* cm*-1*/(molecule cm*-2*) is applied, the OLR flux will increase by $\sim$1.5%.

Our simulation setup is basically the same as the CO2 atmosphere collapse simulation in Section 2.2. The spectral resolution is also 2000 points in wavenumber, which is an appropriate choice to maintain both computation accuracy and efficiency.

The O2-O2 collision-induced absorption (CIA) data is calculated from the HITRAN database. The transition bands and spectral ranges are summarized in Table 1 of Richard et al. (2012). Note that for most of the bands for this dataset, the temperature range is around room temperature (296 K), while the equilibrium temperature ($T_{eq}$) of GJ 1132b is 548 K assuming a planetary albedo of 20%. The binary absorption coefficient does not vary linearly with temperature so extrapolation is not a good choice. Thus in order to simulate radiative transfer more accurately in a O2-rich atmosphere on such hot planet, further experimental and/or theoretical data at high temperatures is needed. Here, we simply use the binary absorption coefficients corresponding to their highest temperatures. For this simulation, we also run the model by 2000 days and present an average result over the last 500 days.

The equatorial surface temperature distributions and vertical temperature profiles in these two experiments are shown in Figure 7. While the peak surface temperatures are much higher, the general structure here is analogous to our CO2 collapse simulation shown in Figure 4, including the larger variations of surface temperatures on the dayside and the strong air temperature inversions in the lower atmosphere on the nightside. For the pure O2 simulation, the upper atmosphere is much colder than the skin temperature ($T_{skin}=T_{eq}/2^{1/4}=461$ K), while in 1D radiative-convective simulations (results not shown) the upper atmosphere is indeed nearly isothermal with the skin temperature due to the weak absorption by O2. Thus the atmospheric dynamics in the GCM plays a remarkable role in cooling the upper atmosphere efficiently. The energetics within a global-averaged atmospheric column can be simplified as (See Peixoto & Oort, 1992, Chap. 13, Eq. (13.21))

[TABLE]

where $\omega=Dp/Dt$ is the vertical pressure velocity, $c_{p}T,gz,(u^{2}+v^{2})/2$ are the enthalpy, potential energy and kinetic energy of air respectively. $\dot{Q}_{dia}$ is the diabatic heating rate, which in our GCM simulation is solely due to the radiative heating in the upper atmosphere. Because the temperature profile in this upper atmosphere deviates significantly from the radiative-equilibrium value (i.e., the skin temperature), the radiative heating rate is non-trivial and must be balanced by the convergence of vertical energy transport represented by the left-hand side of Eq. 2.

Figure 7 also shows that adding 1% CO2 to the atmosphere raises both the surface temperature and near-surface air temperature due to the greenhouse effect of CO2. However, the lower atmosphere on the nightside is cooled because of the strong temperature inversion. CO2 is only an IR absorber here whose radiative effect depends on the lapse rate of the atmosphere. In other words, the surface and lower atmosphere are more closely connected by the enhanced atmospheric emissivity. A simple one-layer atmosphere model for the nightside helps to understand this effect. We assume the surface is a blackbody with temperature of $T_{s}$, and the overlying atmosphere has a temperature of $T_{a}$ and a gray emissivity of $\epsilon$. Ignoring the surface sensible heat exchange between the atmosphere and the surface, the energy budgets at the surface and at the top of atmosphere are

[TABLE]

where $\sigma$ is the Stefan-Boltzmann constant, $OLR_{n}$ is the mean outgoing longwave radiation on the nightside that also equals to the heat transport from the dayside to the nightside. First, Eq. 3 shows that for a gray atmosphere $(\epsilon<1),T_{a}>T_{s}$ and the temperature inversion is necessary on the nightside. Then we can see that an increase of the atmospheric emissivity tends to warm the surface and cool the overlying atmosphere, assuming a small change in the $OLR_{n}$ from Eq. 4 (see also Section 3.3 of Wordsworth (2015) for a similar analysis assuming a uniform atmospheric temperature). However, the change of temperature in the GCM runs is much more complicated than this simple model with these assumptions. First, more energy is transported to the nightside in the O2+CO2 run. Second, the surface sensible heat is comparable to the radiative fluxes. Third, the atmosphere is not gray but only absorbs strongly in some spectral regions by CO2. Our qualitative explanation here does not take these effect into account. Koll & Abbot (2016) developed a two-column radiative-convective-subsiding model to study temperature structures of dry tidally-locked rocky exoplanets, which incorporates the atmospheric dynamics by constraining the large-scale wind speed using heat engine efficiency. Figure 10 in Koll & Abbot (2016) shows similar results to our GCM simulations: an increase in the longwave optical thickness warms the nightside surface but cools the nightside overlying atmosphere, when the atmosphere is not optically thick and has small wave-to-radiative timescale ratio (for our GCM runs on GJ 1132b it is $\sim 0.04$).

Note that GJ 1132b has an Earth-like spin rate assuming synchronous rotation. The effects of planetary rotation on the atmospheric circulation becomes significant compared to our CO2 collapse simulation shown in Figure 4. We can verify this by comparing the planetary radius with the equatorial Rossby deformation radius $L_{RO}^{2}\sim(\sqrt{RT_{eq}}/\Omega a)a^{2}$, where $R$ is the specific gas constant, $T_{eq}$ the equilibrium temperature, $\Omega$ the spin rate of the planet, $a$ the planetary radius (e.g., Koll & Abbot, 2015; Pierrehumbert & Ding, 2016). The non-dimensional parameter $(\sqrt{RT_{eq}}/\Omega a)$ for GJ 1132b is 1.14. As a comparison, this parameter is 0.62 for the Earth and 18 for our CO2 collapse simulation. As a consequence of planetary wave perturbations, the symmetry of the surface temperature around the substellar (270∘) and antistellar point (90∘) in Figure 4 breaks down in the two GJ 1132b simulations, as shown in Figure 7a.

Another consequence of the planetary wave perturbations is equatorial superrotation. The zonal-mean cross sections of the zonal wind are plotted in Figure 8. Equatorial superrotation is a common feature of hot-Jupiter simulations with speeds of 1–4 km s*-1* (Showman & Guillot, 2002; Dobbs-Dixon & Lin, 2008; Menou & Rauscher, 2009). It has also been observed indirectly by the eastward displacement of the hottest region from the substellar point in light curves of the hot Jupiter HD 189733b (Knutson et al., 2007) and by the the ellipsoidal/beaming amplitude discrepancy on Kepler-76b (Faigler et al., 2013). The equatorial super-rotating jets in the two experiments are of the order of 100 m s*-1*, much slower than the simulated jet speed on hot-Jupiters, probably because of the efficient surface momentum dissipation and the heavier molecular weight of the O2 atmosphere compared to the H2-dominated atmosphere on hot-Jupiters. The effects of molecular weight on the zonal wind speed was discussed in Zhang & Showman (2017). To illustrate the underlying mechanism of equatorial superrotation, we calculated the Eliassen-Palm (EP) flux in the latitude-pressure coordinate, which is a useful tool for diagnostics of wave activities in the Earth’s climate (e.g., Edmon et al., 1980). The horizontal and vertical components of the EP flux relate to the meridional eddy momentum and heat flux in the atmosphere, respectively. The divergence of the EP flux represents an internal forcing of the zonal jet by disturbances, so that the equatorial superrotation in Figure 8 can be explained by the divergence of the vectors at the equator. Comparing the two panels in Figure 8 shows that adding 1% CO2 to the atmosphere changes the circulation pattern substantially. It amplifies the contrast of radiative-equilibrium temperature between the dayside and nightside, and therefore increases the wave amplitude and the equatorial super-rotating jet speed. More angular momentum is transported equatorward from high latitudes, creating two easterly jets above 50∘. The position of the jet core also changes. In the pure O2 simulation, two jet cores form in the subtropics associated with the meridional mean angular momentum transport, although the air above the equator still super-rotates. With an additional 1% CO2, one well-formed jet core is found at 0.1 bar above the equator.

3.2 Thermal phase curves

Thermal phase curves, which are the variations of the apparent infrared emission of an planet with its orbital phase, have been observed for hot Jupiters and used to infer their temperature distribution and atmospheric circulation (e.g., Knutson et al., 2007). Furthermore, Selsis et al. (2011) addressed the ability to use the multiband phase curves to characterize the atmosphere on hot terrestrial exoplanets with next-generation observations. Here we compute the disk-integrated and observer-projected spectrum from GJ 1132b following Cowan & Agol (2008) and Koll & Abbot (2015, Appendix C), $F_{\nu}(\theta)$ where $\theta$ is the phase angle (we assume $\theta=0^{\circ}$ at primary transit and $\theta=180^{\circ}$ at secondary eclipse). Then we plot the spectral distribution of the amplitude of $F_{\nu}(\theta)$ and the phase angle corresponding to the highest emission $\theta_{max}$ in Figure 9a and 9b, respectively. Note that $F_{\nu}(\theta)$ here is the disk-integrated and observer-projected flux from the planet. Then the star-planet contrast seen by a distant observer should be $F_{\nu}r_{p}^{2}/F_{\nu,\ast}r_{\ast}^{2}$, where $F_{\nu,\ast}$ is the spectral flux from the M-dwarf GJ 1132 at the same wavenumber, $r_{p}$ and $r_{\ast}$ are the radius of the planet GJ 1132b and the star GJ 1132 respectively.

The spectral distribution of the amplitude of $F_{\nu}(\theta)$ is called the ‘variation spectrum’ in Selsis et al. (2011). In our simulations, it has a similar shape to the emission spectrum from the dayside atmosphere (e.g., Figure 5) because the temperature contrast between the dayside and nightside atmosphere decreases with height (shown in Figure 7b) and window regions exhibits larger variations than optically-thick spectral regions. Meadows (2017) performed a thorough review on the identification of an oxygenic photosynthetic biosphere and discrimination between biological and abiotic sources of O2, but only with a focus on transmission spectroscopy and direct imaging of exoplanets. As described in Selsis et al. (2011), the variation spectrum is an interesting possibility for molecular signature detections. In our simulations, O2 absorption in the variation spectrum is significant in its vibrational fundamental band near 1500 cm*-1*.

Other than the amplitude of the phase curve, our calculation indicates that the spectral distribution of the phase angle corresponding to the highest emission $\theta_{max}$ can also be used for molecular detections. Moreover, the $\theta_{max}$ distribution provides valuable information regarding the vertical profile of the hotspot shift in the atmosphere, which is associated with atmospheric dynamics. In the 1 bar pure O2 simulation, $\theta_{max}$ in the O2 fundamental band indicates that there is insignificant hotspot shift in the lower atmosphere relative to the surface. In the O2+CO2 simulation, $\theta_{max}$ in the same band indicates a 6*∘* eastward shift of the lower atmosphere hotspot, as a result of the circulation change discussed in Section 3.1. The absorption bands of CO2 have more complicated features. We can take the principle absorption band centered at 667 cm*-1* as an example. In both the window region with little absorption and the optically-thick region near the 667 cm*-1* center, $\theta_{max}\sim 180^{\circ}$, indicating small hotspot shifts at the surface and top of the atmosphere. But on the shoulder of this absorption band, $\theta_{max}$ can be as low as 75∘, indicating a significant eastward shift of the hotspot by 105∘ in the middle atmosphere.

Figure 10a shows the vertical cross section of the air temperature anomaly (obtained by subtracting the zonal mean) along the equator in the O2+CO2 simulation to confirm our inferences from the $\theta_{max}$ distribution. There is a remarkable eastward shift of the hotspot in the middle atmosphere near 104 Pa, which is also the altitude of the jet core shown in Figure 8b. Moreover, the longitudinal tilting of the hotspot correlates well with this equatorial super-rotating jet, because the hotspot shift is the consequence of atmospheric dynamics. Previous studies suggested that the eastward hotspot shift is caused by the equatorial super-rotating jet (e.g., Zhang & Showman, 2017; Hammond & Pierrehumbert, 2018), with Zhang & Showman (2017) focusing on the heat advection by the equatorial super-rotating jet and Hammond & Pierrehumbert (2018) focusing on the Doppler-shifting stationary wave response. Note that until now theoretical work has only been done in the shallow-water framework. Further investigation is needed to explore whether baroclinicity is important for wave-mean flow interactions and therefore the longitudinal tilting of the hotspot in this situation of large horizontal temperature variation.

4 Discussion

We have described a new 3D climate model, FMS PCM, which uses a line-by-line approach to describe the radiative transfer.

Before running the GCM, some preparation is required to determine the appropriate setup (most importantly the spectral resolution) and maintain the computational accuracy and efficiency for long-term integrations. Compared to conventional GCMs that use parametrized radiative calculations, however, our model is designed to be both more accurate and flexible for the study of diverse planetary atmospheres. Moreover, because all spectral quantities are calculated within FMS PCM, this model can help interpret next-generation exoplanet observations. For example, we show that the emission spectrum in the CO2 principle absorption band centered at 667 cm*-1* could be used to infer the near surface temperature inversion on the nightside of tidally-locked terrestrial planets.

We have used FMS PCM to study possible O2-dominated atmosphere on GJ 1132b. The simulation results show that a minor amount of CO2 in the atmosphere can change the climate substantially, notably by enhancing the near-surface temperature inversion on the nightside and intensifying the equatorial superrotation and high-latitude easterlies. Moreover, the vertical profile of the equatorial superrotation could be inferred from the phase shift in the CO2 principle absorption band centered at 667 cm*-1*. As an example of using GCM simulations to help interpret exoplanet atmospheric observations, the metallicity of a hot-Jupiter WASP-43b was constrained in a study of Kataria et al. (2015) by comparing a range of GCM results to the Hubble Space Telescope (HST)/WFC3 spectrophotometric phase curve measurements. While the phase-resolved spectroscopic feature on GJ 1132b we have discussed here lies in the mid-IR rather than the visible/near-IR, it should be possible to perform such measurements using next-generation instruments, such as JWST (Selsis et al., 2011; Beichman et al., 2014; Cowan et al., 2015; Koll & Abbot, 2016). Koll & Abbot (2016) estimated the amount of time that is required to measure the thermal phase curves of short-period terrestrial planets with JWST/MIRI, which is similar to the time required to detect molecules via stacked transit spectroscopy.

So far FMS PCM is only a dry model. Moist dynamics, including phase transitions of condensible substances and formation and transport of condensates, is not yet implemented. We also plan in future to extend the radiative transfer model to allow full treatment of multiple scattering by atmospheric particulates (clouds and hazes). Nonetheless, the dry runs discussed in the paper demonstrate the high potential of this approach. We believe that this new type of GCM model will allow the detailed study of a wide variety of important planetary climate problems in future, including the climate dynamics of condensible-rich atmospheres (Ding & Pierrehumbert, 2016; Pierrehumbert & Ding, 2016; Ding & Pierrehumbert, 2018), atmospheric characterization of exoplanets near the inner and outer edges of the habitable zone (Wordsworth et al., 2011; Yang et al., 2013; Kopparapu et al., 2016), and water loss and abiotic oxygen build-up in water-rich atmospheres (Wordsworth & Pierrehumbert, 2014).

Support for this work was provided by NASA grants NASA NNX16AR86G and 80NSSC18K0829. The GCM simulations in this paper were carried out on the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University. The authors would like to thank Tyler D. Robinson for insightful comments of this work and for facilitating model comparisons with the Clima and LiFE model, and Wanying Kang and Hannah Diamond-Lowe for helpful discussions.