Atmospheric circulation of exoplanets II

Atmospheric circulation of tidally locked exoplanets:
II. Dual-band radiative transfer and convective adjustment

Kevin Heng, Dargan M.W. Frierson and Peter J. Phillipps
Zwicky Fellow, ETH Zürich, Institute for Astronomy, Wolfgang-Pauli-Strasse 27, CH-8093, Zürich, Switzerland
Department of Atmospheric Sciences, University of Washington, ATG Building, Box 351640, Seattle, WA 98195, U.S.A.
Geophysical Fluid Dynamics Laboratory, 201 Forrestal Road, Princeton, NJ 08540, U.S.A.
E-mail:, (KH)Email: (DMWF)Email: (PJP)
Submitted 2011 May 20. Re-submitted 2011 July 15 and 2011 August 16. Accepted 2011 August 18.

Improving upon our purely dynamical work, we present three-dimensional simulations of the atmospheric circulation on Earth-like (exo)planets and hot Jupiters using the GFDL-Princeton Flexible Modeling System (FMS). As the first steps away from the dynamical benchmarks of Heng, Menou & Phillipps (2011), we add dual-band radiative transfer and dry convective adjustment schemes to our computational setup. Our treatment of radiative transfer assumes stellar irradiation to peak at a wavelength shorter than and distinct from that at which the exoplanet re-emits radiation (“shortwave” versus “longwave”), and also uses a two-stream approximation. Convection is mimicked by adjusting unstable lapse rates to the dry adiabat. The bottom of the atmosphere is bounded by a uniform slab with a finite thermal inertia. For our models of hot Jupiters, we include an analytical formalism for calculating temperature-pressure profiles, in radiative equilibrium, which accounts for the effect of collision-induced absorption via a single parameter. We discuss our results within the context of: the predicted temperature-pressure profiles and the absence/presence of a temperature inversion; the possible maintenance, via atmospheric circulation, of the putative high-altitude, shortwave absorber expected to produce these inversions; the angular/temporal offset of the hot spot from the substellar point, its robustness to our ignorance of hyperviscosity and hence its utility in distinguishing between different hot Jovian atmospheres; and various zonal-mean flow quantities. Our work bridges the gap between three-dimensional simulations which are purely dynamical and those which incorporate multi-band radiative transfer, thus contributing to the construction of a required hierarchy of three-dimensional theoretical models.

planets and satellites: atmospheres – methods: numerical
pagerange: Atmospheric circulation of tidally locked exoplanets: II. Dual-band radiative transfer and convective adjustmentApubyear: 2011

1 Introduction

Much of the previous research on three-dimensional simulations of exoplanetary atmospheric circulation has focused on the atmospheric dynamics (e.g., Cooper & Showman 2005, 2006; Koskinen et al. 2007; Showman et al. 2008; Dobbs-Dixon, Cumming & Lin 2010; Rauscher & Menou 2010; Thrastarson & Cho 2010; Heng, Menou & Phillipps 2011), while others (Showman et al., 2009; Lewis et al., 2010; Selsis, Wordsworth & Forget, 2011) have developed sophisticated schemes which combine dynamics with multi-band radiative transfer. Despite making progress, our understanding of atmospheric circulation on extrasolar planets (or “exoplanets”) remains incomplete. To move forward, we need to construct a hierarchy of three-dimensional simulations of varying sophistication, so as to elucidate the isolated effects of various pieces of physics as well as the complex interplay between them (Peixóto & Oort, 1984; Held, 2005; Pierrehumbert, 2010). The purpose of the present study is to bridge the gap between these two bodies of work by elucidating the details of a computational setup of intermediate sophistication, so as to serve as a foundation for future work. Our work may be partially regarded as the three-dimensional generalizations of Hubeny, Burrows & Sudarsky (2003), Hansen (2008) and Guillot (2010).

As natural extensions to purely dynamical simulation work, the effects of radiative transfer and convection need to be considered. As a first step in adding radiative transfer, we simplify the treatment by assuming that most of the stellar emission peaks at a wavelength which is shorter than and distinct from that at which the exoplanet re-radiates the absorbed heat (Hansen, 2008; Langton & Laughlin, 2008; Guillot, 2010; Pierrehumbert, 2010). For example, the Sun-like host star of an exoplanet has an emission spectrum which peaks in the optical, since


using Wien’s law, where is the effective temperature of the star. For M stars (red dwarfs), the stellar spectrum instead peaks at . The equilibrium or effective temperature of the exoplanet is111This expression is considered to be order-of-magnitude because there is a factor of order unity associated with the exoplanetary albedo and the efficiency of heat redistribution from the day to the night side.


where is the stellar radius and denotes the spatial separation between the star and the exoplanet. As the stellar photons travel downward into the atmosphere, they are eventually absorbed and re-emitted. The emission spectrum of an exoplanet is expected to peak in the infrared, since


Thus, to a first approximation one can regard the star as emitting “shortwave” radiation onto the exoplanet, which then re-emits in the infrared (“longwave”). The conversion of a shortwave photon into multiple longwave photons serves to increase entropy. Our treatment of the radiative transfer then requires the specification of the shortwave () and longwave () optical depths, as well as the flux associated with stellar irradiation (). Also as a first step, we mimic convection by adjusting unstable atmospheric lapse rates to the stable dry adiabat (Manabe, Smagorinsky & Strickler, 1965). The key advantages of our simpler setup, compared to simulations with multi-band radiative transfer, are that it allows for a more efficient sampling of parameter space and clean comparisons to analytical models.

In §2, we describe our computational setup, including our schemes for stellar irradiation, shortwave/longwave absorption, radiative transfer and convection. Prior to simulating the atmospheric circulation on hot Jupiters, we test our computational setup using Earth-like models (§3). In §4, we describe our method for simulating hot Jovian atmospheres — including an analytical formalism for calculating temperature-pressure profiles which includes the effect of collision-induced absorption — and examine results from several models based on the values of parameters appropriate to HD 209458b. The implications of our results, as well as a concise summary of the present study, are discussed in §5. Table 1 summarizes the list of adopted parameter values as well as the commonly-used symbols. Appendix A contains the technical details of the atmospheric boundary layer scheme used only for our Earth-like models.

Symbol Description Units Earth Hot Jupiter
stellar irradiation constant W m 938.4
meridional temperature gradient parameter 1.4
power law index for shortwave optical depth 2 1
power law index for longwave optical depth 4 2
surface optical depth of shortwave absorbers 0.2 1401
surface optical depth of longwave absorbers at equator 6
surface optical depth of longwave absorbers at poles 1.5
strength of well-mixed longwave absorbers 0.1 1/2000
specific heat capacity at constant pressure (of the atmosphere) J K kg 1004.64 14308.4
ideal gas constant (of the atmosphere) J K kg 287.04 4593
adiabatic coefficient J K kg 2/7 0.321
reference pressure at bottom of simulation domain bar 1 220
areal heat capacity of lower atmospheric boundary J K m
acceleration due to gravity m s 9.8 9.42
radius of (exo)planet km 6371
rate of rotation of (exo)planet s
critical bulk Richardson number 1
roughness length m
von Kármán constant 0.4
surface layer fraction 0.1
hyperviscous time day 0.1
time step s 1200 120
initial temperature K 264 1824
vertical resolution 20 33
latitude degrees - -
longitude degrees 0— 0—
dry adiabatic lapse rate K km
potential temperature K –500
Eulerian mean streamfunction kg s

: expressed in terms of a (exo)planetary day (i.e., rotational period).
: we term the “adiabatic coefficient”, following Pierrehumbert (2010), and avoid calling it
the “adiabatic index” in order to not confuse it with the ratio of specific heat capacities.
Note: one “Earth day” is exactly 86400 s (as used in the text and simulation/analysis codes).

Table 1: Table of parameters, symbols and their values

2 The GFDL-Princeton Flexible Modeling System

In this section, we describe our computational setup, which is based on the Lima release of the FMS, developed by the Geophysical Fluid Dynamics Laboratory (GFDL) at Princeton University (Gordon & Stern, 1982). As the FMS is set up to work in MKS (metres, kilogrammes and seconds) units, most of the discussion will follow suit. Unless specified otherwise, the term “day” refers to an Earth day (exactly 86400 s).

2.1 Dynamics

We implement the spectral dynamical core, which solves the three-dimensional primitive equations of meteorology — these equations assume hydrostatic equilibrium and other consistent approximations such as a small aspect ratio for the atmosphere (e.g., see Heng, Menou & Phillipps 2011 and references therein). Held & Suarez (1994) originally proposed the comparison of purely dynamical simulations with simplified thermal forcing and performed via different methods of solutions, known as the “Held-Suarez benchmark” for Earth. Heng, Menou & Phillipps (2011) extended the Held-Suarez benchmark to tidally-locked exoplanets by examining three additional cases: a hypothetical exo-Earth, a shallow model for hot Jupiters and a deep model of HD 209458b.222The adjectives “shallow” and “deep” refer to the fact that one and about 5 orders of magnitudes in vertical pressure are simulated, respectively. The spectral and finite difference (B-grid) dynamical cores of the FMS were both subjected to these tests. In the case of the deep model of HD 209458b, the predictions for the wind speeds are found to depend upon the chosen magnitude of the hyperviscosity, although the qualitative features of the simulated wind and temperature maps remain largely invariant.

Within the spectral dynamical core, the (linear) fluid dynamical quantities are expressed as a truncated sum of spherical harmonics — the larger the value of , the higher the horizontal resolution. The fiducial resolution we will adopt in the present study is T63L (), which corresponds to a horizontal grid of 192 by 96 points in longitude versus latitude. (See Table 2 of Heng, Menou & Phillipps 2011 for a list of and their corresponding horizontal resolutions.) A finite difference solver with grid points is used for the vertical coordinate, which is cast in terms of the “-coordinate” (Phillips, 1957),


where denotes the vertical pressure and is the time-dependent pressure at the surface. By contrast, the reference pressure at the bottom of the simulation domain, , is a constant. Analogous to , we define . For Earth-like simulations, the reference pressure is usually set to bar = Pa = dyn cm.

Numerical noise accumulates at the grid scale and has to be damped via a “hyperviscous” term with a damping order of , such that the operator acts on the relative vorticity (see §3.1 of Heng, Menou & Phillipps 2011). The pragmatic aim is to apply horizontal dissipation on a time scale that is a small fraction of an exoplanetary day. For example, in the Held-Suarez benchmark for Earth, the adopted time scale is Earth day. For our hot Jupiter simulations, we use HD 209458b day, consistent with the value used in Heng, Menou & Phillipps (2011); this is also, to within an order of magnitude, the highest value needed. It is important to keep in mind that the application of hyperviscosity and the specification of is unsupported by any fundamental theory — horizontal dissipation is solely a numerical tool meant to prevent a simulation from failing due to spectral blocking.

As in Heng, Menou & Phillipps (2011), the simulations are started from an initial state of windless isothermality with denoting the initial temperature, in contrast to using initial temperature-pressure profiles obtained from one-dimensional radiative-convective calculations (e.g., Showman et al. 2009; Dobbs-Dixon, Cumming & Lin 2010; Lewis et al. 2010). The possible dependence of our results on more exotic initial conditions (Thrastarson & Cho, 2010) is deferred to a future study. There is an initialization period for each simulation which we disregard, during which the temperatures and velocities, averaged over the entire (exo)planet, are decreasing/increasing from some initial values to their quasi-steady values. For Earth-like simulations, we find the initialization period to be 200 days. For hot Jupiters, a longer initialization period of 500 days is needed. After quasi-equilibrium is attained, we execute the simulations for a further 1000 days and temporally average the ouput. Minor asymmetries in the physical quantities between the northern and southern hemispheres are artifacts of averaging over a finite period of time (i.e., in this case, 1000 days).

2.2 Stellar irradiation and shortwave absorption

At the top of the atmosphere, stellar irradiation is specified via the function,


where is the stellar irradiation constant,333The Solar constant is typically taken to be about 1367 W m.


where denotes the Stefan-Boltzmann constant and the dimensionless stellar irradiation function is in general a function of both the latitude and the longitude . Since hot Jupiters orbit at distances from their host stars, we have W m for . By contrast, we have W m for Jupiter, which is a factor less than for hot Jupiters. In the case of HD 209458b, we use the values of , and measured or inferred by Mazeh et al. (2000) and Brown et al. (2001) to obtain the value shown in equation (6); we note that . In general, the influence of the stellar irradiation on a (exo)planet and its dependence on geometry (through ) is known as the “thermal forcing”.

The shortwave optical depth determines how the stellar irradiation is absorbed as it travels downward into the atmosphere (Frierson, 2007a; O’Gorman & Schneider, 2008; Merlis & Schneider, 2010),


Specifying the power law index as is equivalent to having a uniformly mixed shortwave absorber, since (Guillot, 2010; Heng et al., 2011b)


where is the shortwave opacity and is the surface gravity. When is constant, we have . When , the shortwave absorber is less well-mixed and resides lower down in the atmosphere. The received flux as a function of vertical pressure is then


Finally, we note that it is possible to specify a mean (exo)planetary albedo, but this only serves to diminish the value of . As such, one only needs to explore the effects of varying the stellar irradiation constant.

2.3 Radiative transfer and longwave absorption

The radiative transfer scheme is based on the two-stream approximation (Mihalas, 1978; Hansen, 2008), which assumes that the radiation can be separated into longwave, upward- () and downward-propagating () fluxes (see §4.2 of Pierrehumbert 2010),


where is the blackbody flux. At the bottom of the simulated atmosphere, the longwave boundary condition is where is the surface temperature; at the top, it is . The Newtonian relaxation term in the temperature equation is replaced by a radiative source term (Frierson, Held & Zurita-Gotor, 2006),


where is the specific heat capacity at constant pressure, is the mass density and is the vertical coordinate. The shortwave, upward-propagating flux is assumed to be zero, i.e., . Furthermore, substituting the Newtonian relaxation scheme for the dual-band radiative transfer scheme alleviates the need to rely on one-dimensional radiative-convective models (e.g., Iro, Bézard & Guillot 2005; see also Ramanathan & Coakley 1978) for calculations of the equilibrium temperature-pressure profile, thereby allowing self-consistency to be achieved. From the perspective of equation (10) alone, the radiative transfer scheme is strictly speaking grey (as already noted by Frierson, Held & Zurita-Gotor 2006) — the shortwave, downward-propagating flux is merely diluted on its way down into the atmosphere via equation (9) — but we term it “dual-band” to emphasize the dual-wavelength nature of our treatment. The effects of scattering are neglected.

The bottom of the atmosphere (or the “surface”) is an idealized slab with a specifiable heat capacity, described by a single temperature governed by the equation,


where is the (areal) heat capacity of the slab. Horizontal transport of fluid within the slab is not modelled. In the (terrestrial) atmospheric sciences literature, this idealized slab mimics the “mixed layer ocean”, which is the m-thick layer of ocean intermediate between the atmosphere and the deep ocean, characterized by rigorous mixing induced by wave motion and turbulence. For simplicity, we do not implement an internal heat flux, typically parametrized by the equivalent blackbody temperature . For example, Iro, Bézard & Guillot (2005) and Guillot (2010) use K in their models of HD 209458b, while Liu & Schneider (2010) adopt 5.7 W m as the internal heat flux in their simulation of Jupiter (which is equivalent to K).

As the radiation is absorbed and re-emitted, the (infrared) photons encounter a longwave optical depth near the surface, described by


where the longwave optical depths at the equator and the poles are given by and , respectively. As already noted by Merlis & Schneider (2010), equation (13) ignores the longitudinal dependence of introduced by water vapour feedback and is thus inadequate for treating tidally locked, Earth-like aquaplanets. Throughout the vertical extent of the atmosphere, the longwave optical depth is a linear combination of well-mixed and segregated atmospheric absorbers,


where and generally depend on geometry (i.e., and ) as well as atmospheric chemistry. The second term in equation (14) represents species of absorbers where the associated pressure scale height is that of the well-mixed species. For example, the linear term may represent well-mixed species like carbon dioxide, while setting approximates the structure of water vapour in the terrestrial atmosphere (Frierson, Held & Zurita-Gotor, 2006). As another example, setting approximates the process of collision-induced absorption (Herzberg, 1952) within the atmospheres of the giant planets in our Solar System (Liu & Schneider, 2010). Following Frierson, Held & Zurita-Gotor (2006, 2007) and Frierson (2007a), we use


with being a dimensionless parameter that controls the strength of the well-mixed longwave absorber.

In Frierson, Held & Zurita-Gotor (2006, 2007), Frierson (2007a) and O’Gorman & Schneider (2008), a key assumption is that the (water) moisture content does not affect radiative transfer, such that the dynamical effects of the latent heat release associated with the condensation of water vapour can be isolated from its radiative effects. Merlis & Schneider (2010) improved upon this treatment by allowing the optical depth associated with water vapour () to depend on the instantaneous column density, thereby providing a simple treatment of water vapour feedback. Since we are concerned only with dry atmospheres in the present study, we will ignore this issue.

2.4 Convection

In simulations of atmospheric circulation, convection is typically mimicked using simplified parametrizations known as “convective adjustment” schemes (e.g., Manabe, Smagorinsky & Strickler 1965; §3.6.9 of Washington & Parkinson 2005). For a dry atmosphere in hydrostatic equilibrium, the change in temperature of a parcel of atmosphere resulting from vertical motion is (e.g., §2.7.2 of Holton 2004)444Equation (16) may be generalized to describe an atmosphere containing moisture and with a corresponding moist adiabatic lapse rate.


where the lapse rate is defined as


with being the vertical spatial coordinate, and the dry adiabatic lapse rate is , which has a value of about 9.8 K km for the terrestrial atmosphere. The potential temperature,


where and is the ideal gas constant, is the temperature a parcel of atmosphere would have if it was compressed or expanded adiabatically from the pressure to the reference pressure at the bottom of the simulation domain. The adiabatic coefficient is related to the number of degrees of freedom of the atmospheric gas : (see §2.3.3 of Pierrehumbert 2010). For example, molecular hydrogen is experimentally found to have , close to the theoretical value for a diatomic gas (). The atmosphere is unstable if (i.e., super-adiabatic lapse rates).555The condition for stability is sometimes known as the Schwarzschild criterion. When the lapse rate is adjusted to a state of stability — subjected to the constraint of total energy conservation — convective adjustment is said to be performed. In other words, the tendency of convection is to minimize the magnitude of the potential temperature gradient and convective adjustment is an approximate way to mimic this process (Ramanathan & Coakley, 1978).

Within the FMS, the dry convective adjustment scheme follows the prescription of Manabe, Smagorinsky & Strickler (1965). At a given time, the temperature difference resulting from departures from the adiabatic lapse rate is computed,


where the integral is taken over the vertical extent of the unstable, model atmosphere layer. An implicit assumption made is that when the lapse rate of a layer is super-adiabatic, convection is rigorous enough to maintain a neutral lapse rate of potential temperature. The kinetic energy created by convection is then dissipated and instantly converted into heat, such that the total potential energy is invariant to convective adjustment. Within each model layer, the temperature is then adjusted by .

It is worth noting that there exists a body of work on moist convective adjustment schemes, largely motivated by the pioneering work of Manabe, Smagorinsky & Strickler (1965), Betts (1986) and Betts & Miller (1986), who proposed the convective adjustment to reference temperature and humidity profiles instead of . Frierson (2007a) adopted a simplified version of the Betts-Miller scheme, where the reference humidity profile is held at a constant threshold value.

2.5 Atmospheric boundary layer

Figure 1: Temporally-averaged, zonal-mean zonal wind (left column) and temperature (right column) profiles as functions of vertical pressure and latitude . Top row: Held-Suarez dynamical benchmark. Middle row: Earth-like model without convection. Bottom row: Earth-like model with convection. Contours are in units of m s (left column) and K (right column).
Figure 2: Temporally-averaged, zonal-mean potential temperature (left column) and the Eulerian mean streamfunction (right column) profiles as functions of vertical pressure and latitude . Top row: Held-Suarez dynamical benchmark. Middle row: Earth-like model without convection. Bottom row: Earth-like model with convection. Contours are in units of K (left column) and kg s (right column).
Figure 3: Globally-averaged temperature-pressure profiles for the trio of Earth-like simulations presented in §3.1.

The boundary layer is the part of the atmosphere where the flow is strongly influenced by interaction with the surface. It can be thought of as enforcing a “no slip” boundary condition, such that large velocity gradients — and hence shears — exist across a small vertical height, which induces turbulent mixing (e.g., §5 of Holton 2004). The boundary layer may become convectively unstable, which induces mixing. The end result is the same: eddies are produced which are smaller in size than those resulting from large-scale rotational flows, are typically unresolved in global simulations of atmospheric circulation and thus need to be accounted for using sub-grid models (e.g., Troen & Mahrt 1986). In other words, turbulence within the boundary layer may be driven mechanically or via buoyancy. The eddies tend to have comparable horizontal and vertical length scales — i.e., they are effectively three-dimensional — and are responsible for heat and momentum exchanges between the atmosphere and the surface. Consequently, the boundary layer affects the dynamics and thermodynamics of the terrestrial atmosphere and is responsible for a non-negligible fraction of kinetic energy loss by the flow (Garratt, 1994). On Earth, the height of the boundary layer is –3000 m and contains of the mass of the terrestrial atmosphere.

We find that the boundary layer scheme is needed if the bottom of the simulation domain is placed within the active part of the atmosphere, such as in an Earth-like simulation. However, this scheme is not required to successfully execute hot Jupiter simulations. Moreover, it is not entirely clear if the analogy with a terrestrial boundary layer can be carried over to mimic drag in the deep, inert layers of a hot Jupiter. As such, although we include a description of our boundary layer scheme for completeness (since it is used for our Earth-like simulations), we relegate the technical details to Appendix A. As described therein, our atmospheric boundary layer scheme requires the specification of four additional parameters: the roughness length , the critical bulk Richardson number , the von Kármán constant and the surface layer fraction .

2.6 Additional physical quantities

Figure 4: Held-Suarez mean flow quantities for the tidally-locked Earth simulation. Top row: zonal wind (m s). Middle row: temperature (K). Bottom row: potential temperature (K). Left column: without convection. Right column: with convection.
Figure 5: Wind maps from the tidally-locked Earth simulation with convective adjustment. Left column: zonal wind. Right column: meridional wind. The first, second and third rows show the maps at , 0.55 and 1 bar, respectively. Contours are in units of m s.

In the interest of clarity, we will define physical quantities, used in the present study, which may be unfamiliar within the astronomical/astrophysical literature.

The potential temperature (equation [18]) is related to the entropy by (e.g., §1.6.1 of Vallis 2006, equation [3.7] of Peixóto & Oort 1984)


The preceding expression is only valid if is constant (otherwise it involves an integration over ). Therefore, examining profiles of the potential temperature is equivalent to determining where the surfaces of constant entropy (“isentropes”) reside, which provides insight into the convective stability of the atmosphere.

Insight into the large-scale circulation patterns within the atmosphere may be obtained by examining the Euleran mean streamfunction, defined as (Peixóto & Oort, 1984; Pauluis, Czaja & Korty, 2008)


where is the temporally- and zonally-averaged meridional velocity. On Earth, we have kg s. We note that the “Sverdrup” (Sv), where 1 Sv kg s, is an alternative unit used in the oceanography community. In the case of hot Jupiters (§4), we have kg s. Our definition of the streamfunction follows that of Pauluis, Czaja & Korty (2008), such that for a circulation cell situated just above the equator in the northern hemisphere, there is southward flow at low altitudes (high values) and northward flow at high altitudes (low values). It is worth noting that the presence of vertical motions due to atmospheric circulation does not violate the assumption of hydrostatic balance provided the vertical accelerations are small compared to the acceleration due to gravity (see §2.5 of Pierrehumbert 2010).

In Figures 6, 8, 12 and 14, we opt to split up the Eulerian mean streamfunction into the day and night sides of the exoplanet. Strictly speaking, this implies that our use of the word “streamfunction” to describe the relevant panels in these figures becomes invalid, since streamfunctions are normally associated with circulations which are non-divergent in the plane being considered — an integration over all longitudes is strictly required. However, our intention is to capture the differences between the atmospheric circulation on the day and night sides of a tidally locked exoplanet. The reader should therefore be mindful of the way we use the term “streamfunction” in these figures.

3 Earth-like Models

3.1 Earth

The first step is to establish a Earth-like model, both as a consistency check and as an operational baseline from which to generalize to tidally-locked exoplanets (Heng, Menou & Phillipps, 2011; Heng & Vogt, 2011). The present model should be regarded as a dry, simplified version of the Frierson, Held & Zurita-Gotor (2006) model. We omit latent heating effects for simplicity and thus the (dry) Earth-like models presented here have less fidelity in simulating the terrestrial climate than the models of Frierson, Held & Zurita-Gotor (2006) in aspects which are strongly influenced by water condensation, such as lapse rates and meridional overturning circulation strengths. The Frierson, Held & Zurita-Gotor (2006) model with moisture has a significantly better fit to observations than the Held-Suarez model in aspects of the terrestrial circulation such as the strength of the Ferrel cell, meridional energy transports and radiative cooling rates. Following Frierson, Held & Zurita-Gotor (2006, 2007) and Frierson (2007a), we define the stellar irradiation function to be


The meridional temperature gradient can be adjusted by varying the value of the dimensionless constant , which we take to be 1.4. The Solar constant is taken to be W m, which corresponds to a global-mean albedo of about 30%. The seasonal cycle of insolation666In this study, we will use the terms “irradiation” and “insolation” interchangably. is not considered, meaning the exoplanet is implicitly assumed to reside on an orbit with zero eccentricity and obliquity. As an aside, we note that solar irradiation has been inferred to vary over much longer time scales (Berger & Loutre, 1991; Crowley, 2000). Other parameter values are listed in Table 1 and are largely adopted from Frierson, Held & Zurita-Gotor (2006) and Frierson (2007a).

A key difference between the Held-Suarez model and our dry Earth-like model is that the radiative equilibrium of the latter is strongly unstable to convection. Preliminary simulations without the boundary layer scheme (§2.5) implemented fail, suggesting that models in which the surface is located within the active — as opposed to inert — part of the atmosphere require a scheme to establish heat and momentum equilibrium between the atmosphere and the surface. This expectation is borne out in our hot Jupiter simulations (§4) — where the surface is located at bar, well into the inert part of the atmosphere — which do not require the boundary layer scheme to run to completion. We therefore include the boundary layer schemes only in our Earth-like models.

In Figure 1, we show the temporally- and zonally-averaged profiles of zonal wind and temperature as functions of the vertical pressure and the latitude , which we term the “Held-Suarez mean flow quantities” (Held & Suarez, 1994). For comparison, we include the Held-Suarez dynamical benchmark, as computed by Heng, Menou & Phillipps (2011). In the trio of simulations, the jet streams in the upper troposphere (–0.4 bar), at mid-latitudes, are clearly evident. There are quantitative differences between each of the temperature profiles, but the general trends are that the temperature decreases away from the equator and also with increasing altitude.

More insight on the convective stability of the model atmospheres may be gleaned by examining the temporally- and zonally-averaged potential temperature profiles, as shown in Figure 2. The potential temperature may be regarded as the “height-adjusted” temperature (§III of Peixóto & Oort 1984). If it monotonically increases with height and does not depend on latitude, its structure is termed “barotropic”. In this case, no convection occurs. If there is a gradient in the potential temperature across latitude, then its structure is said to be “baroclinic”. In such situations, horizontal or “slanted” convective motions occur. Tapping into the available potential energy of the atmosphere (Lorenz, 1955), the baroclinic eddies resulting from these motions are responsible for transporting sensible heat from the equator to the poles (Pierrehumbert & Swanson, 1995). The stratosphere and troposphere are thus regions of the atmosphere where the potential temperature structure is (largely) barotropic and baroclinic, respectively, at least within the temporally- and zonally-averaged context of our models. In reality, the terrestrial stratosphere is not strictly barotropic — along isobars, there are distinctly non-zero variations of temperature with latitude which are dynamically important. The height or vertical pressure at which the troposphere transitions into the stratosphere is the tropopause (e.g., see Figure 4 of Frierson 2007a), located at bar near the equator (the tropics) and –0.4 bar otherwise (the temperate and polar regions).

The middle row of Figure 2 shows the simulation with radiative transfer but without convective adjustment. Unlike in the Held-Suarez benchmark, the existence of a surface with a finite heat capacity (see §2.3) and its heating by solar irradiation (see §2.2) results in convectively unstable vertical columns, thus necessitating the treatment of convection. The bottom row of Figure 2 demonstrates how the dry convective adjustment scheme (see §2.4) brings the model atmosphere to convective stability in the vertical direction, while preserving the baroclinic instabilities (i.e., horizontal convection) produced in the simulation. Essentially, the role of convective adjustment is to straighten the isentropes such that they are vertically stable.

Atmospheric circulation of the terrestrial atmosphere can be generally regarded as having a three-cell structure (e.g., Peixóto & Oort 1984; §2.1.5 of Washington & Parkinson 2005). The Hadley and polar cells, located near the equator and poles, respectively, are direct circulation cells because they are driven by heating and cooling patterns of the Earth’s surface. Both cells are characterized by rising air in their warmer branches. By contrast, the mid-latitudinal Ferrel cells are indirect because they are driven by the presence of baroclinic eddies; cooler air is forced to rise. The presence of these cells may be revealed by examining the Eulerian mean streamfunction (e.g., Frierson 2007a; Merlis & Schneider 2010), as defined in equation (21). The Held-Suarez benchmark, which mimics heating and cooling patterns by forcing the temperature field to relax to an ad hoc “equilibrium temperature” profile on a specified “cooling” time scale, manages to produce the Hadley, Ferrel and polar circulation cells, although the circulation associated with the Ferrel cells is too strong by a factor of two (top right panel of Figure 2). As observed on Earth, the polar cells are relatively weak with having values about an order of magnitude lower than for the other two cells. The simulations with dual-band radiation (middle right and bottom right panels of Figure 2) manage to only produce the Hadley and Ferrel cells; there is a pair of cells near the poles, but they reside only at low altitudes (–1 bar). The absence of the polar cells in our simulations and also those of Frierson, Held & Zurita-Gotor (2006) is unsurprising since the details concerning the polar temperatures (e.g., sea ice, the presence of Antarctica) are not modelled. The Hadley cell is significantly too strong in the dual-band simulation due to the near-neutral convective stability in the tropics, while the Ferrel cell is confined to only high altitudes due to the presence of convective momentum fluxes in the lower troposphere (see Frierson, Held & Zurita-Gotor 2006 for more details) — these issues are ameliorated in simulations with moisture.

Figure 3 shows the globally-averaged temperature-pressure profiles for the trio of simulations. By construction, the Held-Suarez benchmark does not produce the temperature inversion empirically observed in the stratosphere, since the stratospheric temperature is set to be constant in this model. There are slight differences between the - profiles for the simulations with and without convective adjustment, where the latter model produces a small temperature inversion in the stratosphere.

Our dry model neglects the fact that the terrestrial tropospheric temperature profile is strongly influenced by latent heat release from the condensation of water vapour. In the tropics, the temperatures are approximately moist adiabatic (Xu & Emanuel, 1989), while in the extratropics they are determined by a coupling between latent heat release and baroclinic eddies (e.g., Frierson, Held & Zurita-Gotor 2006) with additional contributions from the land-sea contrast. On average, these factors assure that the tropospheric lapse rate on Earth is less than the dry adiabatic value.

In summary, we have generalized our Earth-like model away from the purely dynamical Held-Suarez benchmark and demonstrated the utility of examining profiles of the potential temperature and Eulerian mean streamfunction. We have also demonstrated the necessity of including the treatment of convection alongside radiative transfer, at least for our Earth-like models. We next extend the same simulation and analysis techniques to a hypothetical tidally-locked Earth.

3.2 Tidally-locked Earth

Figure 6: Held-Suarez mean flow quantities for the Eulerian mean streamfunction (in units of kg s) in the case of the hypothetical tidally-locked Earth. Top row: day side. Bottom row: night side. Left column: without convection. Right column: with convection. The irregular contour intervals are chosen to reveal the presence of weaker circulation cells.
Figure 7: Hemispherically-averaged temperature-pressure profiles for the day and night sides of a tidally-locked Earth (§3.2). Shown are the simulations with and without convective adjustment.

The case of a hypothetical tidally-locked Earth is a useful and computationally efficient case study for operationally transitioning between the simulation of Earth and hot Jupiters (Heng, Menou & Phillipps, 2011). The stellar irradiation function is (Merlis & Schneider, 2010),


where the substellar point is located at and . Our (arbitrary) choice is . The parameter values are the same as those adopted in the Earth-like simulation, but the rotational frequency is slowed down such that the exoplanet takes one year to rotate on its axis (Merlis & Schneider, 2010; Heng, Menou & Phillipps, 2011),


Our functional form for the normalization of the longwave optical depth, as described by equation (13), does not have a longitudinal dependence and is the same as that adopted by Frierson, Held & Zurita-Gotor (2006) and O’Gorman & Schneider (2008). Merlis & Schneider (2010) have noted that a more realistic treatment is to allow for to longitudinally vary due to longwave water vapour feedback. However, since our intention is to implement a simplistic, dry model for a tidally locked exo-Earth merely as a sanity check, we ignore this refinement while acknowledging its need if the intention is to model tidally locked, Earth-like aquaplanets.

The Held-Suarez mean flow quantities for zonal wind, temperature and potential temperature are shown in Figure 4. Our results for the zonal-mean zonal wind profile are similar to the top-left panel of Figure 6 of Merlis & Schneider (2010), who noted that the profile contains weak residuals of opposing contributions from various longitudes. It is therefore unsurprisingly that the simulations with and without convective adjustment have zonal-mean zonal wind profiles which are somewhat different. This phenomenon can be seen more clearly in the zonal and meridional wind maps as functions of longitude and latitude (Figure 5), which indicate the presence of large circulation cells. The potential temperature profile shows that the atmosphere is, on average, approximately barotropic down to bar, implying that the tropopause is located farther down. We show the zonal-mean potential temperatures without separating them into day and night side profiles as they are fairly similar; we will examine this issue in greater detail for the hot Jupiter simulations (§4). Rather, our main intention is to demonstrate the effect of convective adjustment.

The Eulerian mean streamfunction in Figure 6 reveals the presence of multiple circulation cells on the day and night sides. On the day side, a pair of cells exist in both the northern and southern hemispheres. The cells are able to extend from the equator to the poles, because the reduced rate of rotation of the exoplanet weakens the Coriolis deflection. On Earth, the relatively rapid rotation limits the Hadley cells to about . On the night side, a pair of large, direct cells reside at higher altitudes ( bar); a collection of four smaller circulation cells is evident at lower altitudes. These features are robust to the application of convective adjustment.

The hemispherically-averaged temperature-pressure profiles for the day and night sides of the tidally-locked Earth are shown in Figure 7. On average, the temperature differences between the two sides are barely 10 K and are only evident at bar. As in the Earth-like case (§3.1), convective adjustment has little effect on the - profiles.

3.3 Other studies

The computational setup described here has also been used for other studies of the terrestrial atmosphere, including the analysis of equatorial transients (Frierson, 2007b), the width of the Hadley circulation in different climates (Frierson, Lu & Chen, 2007), the vertical temperature profile in mid-latitudes (Frierson, 2008), how the cross-equatorial Hadley circulation is affected by remote heating (Kang et al., 2008) and the location of the jet streams (Lu, Chen & Frierson, 2010). It was also applied to the study of the atmosphere on Titan (Mitchell et al., 2009).

4 Hot Jupiter

Figure 8: Held-Suarez mean flow quantities for the T63L33 deep model of HD 209458b (see text). Top row: potential temperature (K). Bottom row: Eulerian mean streamfunction ( kg s). Left column: day side. Right column: night side. For the potential temperature profile, contour intervals are made uneven in the interest of clarity. For the Eulerian mean streamfunction profile, one contour interval is irregular so as to reveal the presence of weaker streamlines.

To facilitate comparison with our previous work, most of the parameter values we adopt for our hot Jupiter simulations are culled from the deep model of HD 209458b presented in Cooper & Showman (2005, 2006) and Heng, Menou & Phillipps (2011). As in Heng, Menou & Phillipps (2011), we use a T63L33 resolution with a range of pressures of bar. A noteworthy difference from the Earth-like models is the assumption that stellar irradiation impinges upon a well-mixed shortwave absorber of unspecified chemistry (). The adiabatic coefficient used is ().

4.1 Previous work: purely dynamical model

As a prelude to our results, we re-analyze the simulation output from the purely dynamical, deep model of HD 209458b as computed by Heng, Menou & Phillipps (2011). Figure 8 shows the Held-Suarez mean flow quantities for the potential temperature (using equation [18]) and Eulerian mean streamfunction (using equation [21]) separated into the day and night sides of the hot Jupiter. Near the top of the atmosphere (–10 mbar) on the day side, the atmosphere is slightly baroclinic but with flow converging towards the equator (unlike in the case of the terrestrial troposphere). On the night side, the atmosphere is noticeably more baroclinic. On both the day and night sides, the atmosphere becomes barotropic for bar by construction, since the Newtonian relaxation time is defined to be infinite at these pressures. The relative weakness of Coriolis deflection again allows for a pair of circulation cell extending from the equator to the poles, as indicated by the Eulerian mean streamfunction. Unlike in the Earth-like (§3.1) and tidally-locked Earth (§3.2) simulations, the atmosphere sinks at the equator and rises at the poles on the day side; the reverse happens on the night side. As demonstrated by Showman & Polvani (2010, 2011) using a hierarchy of analytical models and simulations, there is downward transport of eddy momentum at the equators of hot Jovian atmospheres, which causes counter-rotating, equatorial flow. However, the horizontal convergence of eddy momentum yields super-rotating, equatorial flow and tends to dominate the net flow (as shown in Figure 4 of Showman & Polvani 2011). Our simulations are consistent with this picture. It is somewhat counter-intuitive that the circulation is stronger on the night side (i.e., higher values of ). There is virtually no meridional circulation in the inert ( bar) part of the atmosphere.

4.2 Setup

A number of conceptual challenges exist when generalizing our setup to simulate hot Jovian atmospheres. Firstly, as discussed in Showman et al. (2009), it is presently unclear how to physically describe the frictional dissipation resulting from the differences in flow structures between the atmosphere and its deeper counterparts within a hot Jupiter. Therefore, while the boundary layer scheme described in §2.5 offers a tempting option to simulate frictional dissipation, we choose to switch it off for our hot Jupiter simulations. Essentially, our models assumes “free slip” lower boundary conditions. A different approach is adopted by Liu & Schneider (2010), who simulated the zonal wind structures of Jupiter, Saturn, Uranus and Neptune by varying the (Rayleigh) drag at the lower boundary in order to better fit observations, with the justification that such a formulation mimics the Ohmic dissipation caused by the induced magnetic fields opposing the zonal flows. Since there are virtually no empirical constraints on the zonal wind profiles of hot Jupiters, unlike for the giant planets in our Solar System, we do not adopt the approach of Liu & Schneider (2010).

Secondly, the bottom of the simulation domain (§2.3) does not lend itself to easy physical interpretation, unlike in the Earth-like cases where it mimics the mixed layer ocean. A plausible approach is to demand that the radiative time scale (Showman & Guillot, 2002),


is continuous at the bottom. At and beneath the bottom of the simulation domain, the hot Jupiter has a finite thermal inertia if the heat capacity is non-zero ( J K m) — both longwave and intrinsic heat are not instantaneously (re-)emitted. The thermal inertia time scale is


where the thermal inertia is defined as (e.g., Palluconi & Kieffer 1981)


with , and being the thermal conductivity, mass density and specific heat capacity at constant pressure associated with the bottom, respectively. The thermal conductivity ranges from W K m for gaseous hydrogen (and helium) to W K m for metallic hydrogen (Hubbard, 1968); the latter is not expected to exist below pressures bar (Stevenson, 1982; Hubbard, Burrows & Lunine, 2002).

To ensure continuity between the simulated atmosphere and the bottom, we demand that . Furthermore, may be eliminated in favour of and via the ideal gas law (e.g., §2.3.1 of Pierrehumbert 2010),


where the specific gas constant is , J K kg is the universal gas constant and is the mean molecular weight of the atmospheric molecules. The specific gas constant may be related to more familiar quantities via


where is the Boltzmann constant, is the mean mass of the atmospheric molecules and denotes the mass of a hydrogen atom. For example, if J K kg (as we are assuming for HD 209458b), then , close to the value for an atmosphere dominated by molecular hydrogen. By equating equations (25) and (26), we obtain a plausible estimate for the thermal inertia,


Besides the order-of-magnitude nature of the estimate, the main uncertainty lies in the value to assume for . We will see later that our results are insensitive to these uncertainties.

4.3 Analytical formalism for temperature-pressure profiles: generalization of Guillot (2010)

The most efficient way to initiate a given simulation is with a temperature-pressure profile in radiative equilibrium, as may be computed using the models of Hubeny, Burrows & Sudarsky (2003), Hansen (2008) or Guillot (2010). We simplify our initial condition even further: we assume a constant temperature which does not depend on latitude, longitude or pressure. Since the radiative time scale increases with pressure (Showman & Guillot, 2002; Iro, Bézard & Guillot, 2005), the temperatures at low pressures rapidly equilibrate to values consistent with dynamical and radiative equilibrium, whereas the temperatures at depth () remain close to radiative equilibrium. It is therefore plausible to select . In this sub-section, we generalize the analytical results of Guillot (2010) to obtain the temperature-pressure profile with both and 2, and subsequently use it to estimate the values for .

By retaining the assumption of a constant shortwave opacity , but allowing for an arbitrary functional form for the longwave opacity , equation (29) of Guillot (2010) may be generalized as


where is the blackbody temperature of the internal heat flux, , , is the mass density, and is the vertical spatial variable (with the centre of the exoplanet as a zero reference point). The distance from to the top of the atmosphere is ; we expect . The quantity is thus the column mass, per unit area, measured from the top of the atmosphere.

The quantity is a dilution factor meant to mimic stellar irradiation being the strongest at the substellar point, since equation (31) assumes isotropic irradiation. The equilibrium temperature of the hot Jupiter, assuming zero albedo and no heat redistribution, is


For the values of , and associated with HD 209458b (see equation [6]), we have K. The corresponding irradiation temperature is K.

If is constant, then equation (31) reduces to equation (29) of Guillot (2010). With K, it follows that the temperature at depth is




Note that the equality in the preceding definition only holds for the assumption of constant opacities. Equation (49) of Guillot (2010), which averages over latitude and longitude, provides a more accurate estimate for the (mean) temperature at depth,777If , equation (33) yields K.


Equating (33) and (35) allows for an estimate of the dilution factor required: . Following Guillot (2010), we have adopted .

Figure 9: Analytical temperature-pressure profiles calculated using our formalism in §4.3 (solid curves), compared against simulated ones with only one Earth day of integration (dotted curves). Note that the simulated profiles are hemispherically-averaged over the day side only. Left panel: the increases in the temperatures at depth, due to the effect of collision-induced absorption, are described by a single parameter . Right panel: varying the ratio of shortwave to longwave opacity normalizations for a fixed value of .

As noted by Guillot (2010), the assumption of becomes inadequate deep inside the atmosphere, because collision-induced absorption becomes a non-negligible effect. The longwave opacity then scales linearly as pressure, which implies that a correction term with needs to be included. We demand that




The quantity is the normalization for the longwave optical depth in the absence of collision-induced absorption, while the dimensionless, multiplicative factor accounts for the increase of the longwave optical depth at due to this effect. Collision-induced absorption becomes important when


We write the longwave opacity as


where , and the second equality follows from the assumption of hydrostatic equilibrium,


Evaluating equation (31) in conjunction with equation (39), we get


where we have defined . Note that


is defined such that it is constant and has no dependence on either or . When , we recover , , and equation (41) reduces to equation (29) of Guillot (2010).

It follows that the temperature at depth ( K, , ) becomes


The additional term in equation (43) represents the increase in temperature at depth due to the effect of collision-induced absorption. It is non-negligible only when . Specifically, when


Even with a wide range of values in , we get –1800 K. For example, setting and 2000 corresponds to and 1824 K, respectively. For these respective values of , collision-induced absorption dominates the long-wave optical depth at pressures greater than 24, 2.2, 0.22 and 0.11 bar. Figure 9 plots the temperature-pressure profiles corresponding to these values of , as well as three examples of profiles with different values of . We have chosen for the three cases with different values so as to emphasize the differences between them (which are enhanced with increasing ).

Note that in deriving the results in this sub-section, we have retained the values of the first and second (longwave) Eddington coefficients (1/3 and 1/2, respectively) as adopted by Guillot (2010), as well as the closure relation that the ratio of the second to the zeroth moments of the shortwave intensity is 1/3. Detailed modelling of the opacity sources present in hot Jovian atmospheres (e.g., Gustafsson & Frommhold 2003) will inform us about the actual value of , but for now we are content with being able to account for the effect of collision-induced absorption via a single parameter. Readers interested in broader generalizations of Guillot (2010) are referred to Heng et al. (2011b).

4.4 Baseline models

4.4.1 Constant opacities ()

Figure 10: Comparing analytical and simulated temperature-pressure profiles for the hot Jupiter models. The profiles labelled “global” and “isotropic” correspond to simulations with globally-averaged temperatures and uniform irradiation, respectively (see text). The dashed curve (“Simulated (full)”) is taken from a globally-averaged, 1500-day simulation where the first 500 days of initialization have been disregarded, while the dotted curves (“Simulated”) are taken from only 1 day of integration (with no output being disregarded).
Figure 11: Temperature (top panel), zonal wind velocity (middle panel) and meridional wind velocity (bottom panel), averaged over the entire , hot Jupiter simulation domain, as functions of the simulation time (in Earth days). For the velocities, we have computed the root mean square (RMS) values as functions of time. The dotted vertical line marks the end of the initialization (“spin up”) phase of the simulation.

As a first step, we test our computational setup against the analytical models of Guillot (2010), as stated in equations (29) and (49) of that study. (See also Hansen 2008). To be consistent with the constant shortwave and longwave opacities of cm g and cm g, respectively, as adopted by Guillot (2010), we use and such that . One can choose either or . In the absence of robust empirical constraints, we assume that the longwave optical depth does not depend on latitude,


an assumption we will retain for the rest of the study.

Figure 10 compares the analytical versus simulated temperature-pressure profiles. The isotropic profile is obtained from a simulation where the stellar irradiation is uniform everywhere on the surface. This is equivalent to assuming in equation (29) of Guillot (2010) — which corresponds to the situation where the exoplanet is uniformly irradiated with the substellar intensity everywhere — and accounts for the higher temperature ( K) at depth. For the profiles labelled “global”, the stellar irradiation is described by equations (5) and (23) and the resulting temperature-pressure profiles are averaged over the entire globe, thus resulting in a lower temperature ( K) at depth. The temperatures become isothermal at depths at which the stellar irradiation becomes completely absorbed (). We have extracted temperature-pressure profiles only from the first Earth day of integration, such that the zonal winds in the simulations have not had time to increase to their quasi-equilibrium values (labelled “Simulated”). Recall that our simulations are initiated with a constant temperature and not from a - profile in radiative equilibrium (see §4.3). Thus, the reasonable agreement between the dotted and solid curves in Figure 10 demonstrates that the code is implementing the initial radiative condition correctly — it is not a rigourous test of the radiative transfer scheme. The agreement is imperfect because the analytical expressions of Guillot (2010) use approximate closure conditions (i.e., the Eddington approximations) and also because the vertical grid points in the simulated profiles are taken to be the larger of the pressure half levels in a Simmons-Burridge scheme (see §3.1 of Heng, Menou & Phillipps 2011). The dashed curve in Figure 10 (labelled “Simulated (full)”) is the globally-averaged temperature-pressure profile where the simulation is executed for 1500 days with the first 500 days being disregarded. Its general agreement with the analytical profile constitutes a weak test of the radiative transfer scheme at lower pressures: at bar, the radiative time scale is days, which is more than an order of magnitude shorter than our nominal initialization period of 500 days, implying that radiative equilibrium has been attained at these pressures. At higher altitudes ( bar), the fact that the simulated profile (dashed curve) in Figure 10 is generally –100 K cooler than the analytical one (solid curve) suggests that thermal energy has been converted to mechanical energy in the form of winds, predominantly in the zonal and meridional directions.

Figure 11 shows the temperature, zonal wind velocity and meridional wind velocity, averaged over the entire simulation domain, as functions of the simulation time. The zonal wind is absent during the first few Earth days of simulation, but gradually builds up as the simulation attains quasi-equilibrium (sometimes termed “spin up” in the literature). When the zonal wind is accelerated from rest, its equilibrium value is reached when the accelerations due to vertical and horizontal eddy momentum convergence attain steady values, which is expected to happen within 100 Earth days (Showman & Polvani, 2011).888The spin up of the zonal wind is shown in Figure 11 of Showman & Polvani (2011), but this specific three-dimensional model was constructed for the case of HD 189733b. However, we do not expect dramatic changes in the spin up period for the case of HD 209458b. It is apparent that 500 days is a reasonable period for the initialization phase. The mean temperature plateaus at about 1430 K, which is approximately the equilibrium temperature of the exoplanet. It is possible that the very deepest atmospheric layers ( bar) have not completely “spun up” yet, but in simulations which we executed for longer than 1500 days we witness that the qualitative trends in our Held-Suarez mean flow quantities are invariant, although the associated quantitative details may change slightly.

An uncertainty of our technique is the value of to adopt: for W K m, we have J K m. We thus executed simulations with and J K m. The resulting mean temperature-pressure profiles are virtually identical with differences of at most K. We further verified that simulated profiles associated with J K m and J K m deviate substantially ( K) from the analytical ones at depth.

The zonal-mean zonal wind, temperature, potential temperature and Eulerian mean streamfunction are shown in Figure 12. It is re-assuring that the zonal wind profile resembles Figure 5 of Showman et al. (2009), where a super-rotating, equatorial jet reaches down to about 10 bar, flanked by counter-rotating jets at mid-latitudes. The super-rotating and counter-rotating jets have maximum speeds of about 5.9 km s and -0.6 km s, respectively. From examining the potential temperature profiles, one can infer that the atmosphere is baroclinic for bar (and barotropic at greater pressures), in agreement with the purely dynamical results presented in Figure 8. A difference from Figure 8 is that there are now two pairs of circulation cells in each hemisphere (as is evident from examining the Eulerian mean streamfunction). We note that convective adjustment is negligibly invoked, which is unsurprising given the deep radiative zone — which encompasses both the baroclinic and barotropic atmospheric components — of hot Jupiters.

4.4.2 Collision-induced longwave absorption (, )

Figure 12: Held-Suarez mean flow quantities for the baseline hot Jupiter simulation with (see §4.4.1). Top left panel: zonal wind (m s). Top right panel: temperature (K). Middle row: potential temperature (K). Bottom row: Eulerian mean streamfunction ( kg s).
Figure 13: Comparing analytical and simulated temperature-pressure profiles for the case of and . The solid curve is based on our generalization of the Guillot (2010) analytical profile (see §4.3). The dotted curve (“Simulated”) is taken from a hemispherically-averaged, 1-day simulation, while the dashed curve (“Simulated (full)”) is taken from a hemispherically-averaged, 1500-day simulation where the first 500 days of initialization have been disregarded. Note that both of the simulated profiles are extracted from the day side only.

Collision-induced absorption,999For example, the hydrogen molecule in isolation possesses no dipole moment, which implies it cannot absorb photons. However, it may form transient “super molecules” with other hydrogen molecules under high pressure and the resulting, weak dipole moment allows for absorption. a phenomenon first discovered by Herzberg (1952) to be at work within the atmospheres of Neptune and Uranus, becomes a non-negligible effect at large pressures. We account for this effect using the formalism derived in §4.3. As before, our fiducial model assumes , but we now have and where . We first check that our simulated and hemispherically-averaged (day side only) temperature-pressure profiles are consistent with the analytical ones as derived in §4.3 (Figure 9). The reasonable agreement between them for –2000 again demonstrates that our code is implementing the initial radiative condition correctly (but does not constitute a rigourous test of the radiative transfer scheme).

We next focus on a specific value of the correction factor to the longwave optical depth, at the bottom of the simulation domain, due to collision-induced absorption: . Figure 13 shows the temperature-pressure profiles from our analytical formalism (solid curve), a 1-day simulation meant to mimic the case of no atmospheric dynamics (dotted curve) and a 1500-day simulation with the first 500 days being disregarded (dashed curve). The simulated profiles are hemispherically-averaged over the day side only. The agreement between the solid and dotted curves is reasonable considering that equation (41) was formulated for isotropic stellar irradiation (with a dilution factor ). We have executed simulations with and J K m and witnessed differences of at most 1 K in the globally-averaged temperature-pressure profiles, again demonstrating that uncertainties in the value of do not significantly affect our results.

Figure 14: Same as Figure 12, but for a simulation with and .

The Held-Suarez mean flow quantities from the full simulation are shown in Figure 14, where it is apparent that they are qualitatively similar to those from Figure 12. The zonal-mean zonal wind profile resembles that obtained from the purely dynamical simulation of Heng, Menou & Phillipps (2011) (see the top left panel of their Figure 12). An equatorial, super-rotating jet with a maximum speed of about 6.1 km s is again flanked by slower, counter-rotating jets ( km s) at mid-latitudes. At bar, the zonal wind structure becomes uniform across latitude, which is consistent with the uniform temperature and barotropic structure present. Two pairs of circulation cells extend from the equator to the poles on both the day and night sides, in disagreement with the purely dynamical results, presented in Figure 8, where only a single pair of cells exists. Unlike in the terrestrial atmosphere, our baseline model shows a baroclinic upper atmosphere sitting atop a barotropic lower atmosphere. It is thus somewhat awkward to prescribe the terms “stratosphere” and “troposphere” to the upper and lower atmospheres, respectively. The analogue of the tropopause sits at bar.

Figure 15: Snapshots of the temperature field as functions of latitude () and longitude (), for the hot Jupiter simulation with and , at mbar (top left panel), 216 mbar (top right panel), 4.69 bar (bottom left panel) and 21.9 bar (bottom right panel). The snapshots are taken at 1500 days after the start of the simulation. Temperatures are given in K.

Figure 15 shows temperature maps as functions of latitude and longitude at mbar, 216 mbar, 4.69 bar and 21.9 bar. These pressure levels are chosen to match those shown in Figure 8 of Heng, Menou & Phillipps (2011). At mbar, the location where the temperature is at its maximum — which we term the “hot spot” — is shifted away from the substellar point, which is somewhat different from the top left panel of Figure 8 of Heng, Menou & Phillipps (2011) and the second panel (from the top) of Figure 16 of Showman et al. (2009). Otherwise, the characteristic chevron-shaped feature, which is seen in all of the three-dimensional simulations of HD 209458b (e.g., Showman et al. 2009; Rauscher & Menou 2010; Heng, Menou & Phillipps 2011), appears at mbar. Deeper in the atmosphere, the flow becomes dominated by advection and longitudinal differences in temperature become less apparent.

4.5 Comparison of different models

Figure 16: Hemispherically-averaged temperature-pressure profiles, separated into the day and night sides, for the hot Jupiter models presented in this study. Left panel: comparing and and 1000 cases. Note that the () and cases are almost indistinguishable in this plot. Right panel: comparing profiles with and different values of . The curves labelled by “HMP” are taken from the HD 209458b model of Heng, Menou & Phillipps (2011), which only considers atmospheric dynamics and uses Newtonian relaxation to mimic radiative cooling. The dashed, horizontal line indicates the approximate location of the longwave photosphere. In the right panel, the condensation curves for TiO from Spiegel, Silverio & Burrows (2009) are shown.

We next compare hemispherically-averaged temperature-pressure profiles from different models in Figure 16, separated into the day and night sides of the simulated hot Jovian atmospheres. We see that for a fixed value of , the models with –1000 essentially have the same temperature-pressure profiles, both on the day and night sides, at and above the longwave photosphere. The maximum zonal and meridional wind speeds are virtually unchanged: 5.9 km s and from to km s, respectively (at least for the fiducial magnitude of the hyperviscosity adopted in this study). The day-night temperature contrasts are almost indiscernible, even though the temperatures at depth are modified differently due to the effect of collision-induced absorption. The day-night temperature contrasts from the purely dynamical model of Heng, Menou & Phillipps (2011) are somewhat different, but we should keep in mind that this simulation employs Newtonian relaxation via an equilibrium temperature-pressure profile, which in turn contains a parameter that sets the (initial) day-night temperature contrast as a function of pressure (i.e., –1000 K from bar to mbar). In our improved simulations, the day-night temperature contrasts are computed self-consistently.

The model with () exhibits temperature inversions on both the day and night sides, consistent with the finding by Hubeny, Burrows & Sudarsky (2003), Hansen (2008) and Guillot (2010) that models should produce temperature inversions. (See Burrows & Orton 2010 for a review of the observational evidence for temperature inversions in hot Jovian atmospheres.) Curiously, the temperature at depth is K, which is lower than K. Thus, this model makes a prediction that infrared observations which are able to probe the deep, barotropic layers of a hot Jupiter should infer a equivalent-blackbody temperature which is lower than the equilibrium temperature of the exoplanet. For the model with , we note that the Held-Suarez mean flow quantities are qualitatively similar to the baseline hot Jupiter models previously presented (except for the temperature field), despite the differences present in the temperature-pressure profiles.

These differences in the temperature-pressure profiles are relevant to observations, because the longwave photosphere is, in the absence of clouds, located at a pressure of (Burrows & Liebert, 1993)


We conclude that the quantity controls the presence or absence of a temperature inversion, as well as the day-night temperature contrast. The key differences from the one- or two-dimensional models is that the redistribution of heat/energy and the depths at which shortwave/longwave absorption occurs are calculated self-consistently once the optical depths are specified. We also note that while we have shown hemispherically-averaged temperature-pressure profiles for the day and night sides, our calculations generally produce three-dimensional - profiles.

5 Discussion

5.1 Observational consequences

Figure 17: Temporally-averaged maps of the flux associated with the outgoing longwave radiation, as functions of latitude () and longitude (), for two different hot Jupiter models (, ): (left panel) and 2 (right panel). Fluxes are given in units of W m.
Figure 18: Latitudinally-averaged outgoing longwave radiation (OLR; top panel) as a function of the longitude and thermal phase curve (bottom panel) as a function of the phase angle , for the two models shown in Figure 17. In the bottom panel, the thermal phase curves peak before the secondary eclipses occur.

A direct consequence of our results concerns the vertical mixing of atmospheric fluid. As explained in equation §2.6, the strength of atmospheric circulation is quantified by the Eulerian mean streamfunction. On Jupiter, we have kg s (Liu & Schneider, 2010). As shown by our results in Figures 8, 12 and 14, the circulation in hot Jovian atmospheres has a strength of kg s, which is times stronger than on Jupiter. Since the circulation cells extend from mbar to bar in our models, the relevant length scale is , where is the pressure scale height. If we pretend that we may prescribe a “diffusion coefficient” to the atmospheric circulation, then


which is consistent with the required value of cm s, estimated by Spiegel, Silverio & Burrows (2009) for HD 209458b, in order to keep particles associated with TiO, with sizes of 0.1–10 m, aloft such that they may act as shortwave absorbers capable of producing a temperature inversion. Unfortunately, the Eulerian mean streamfunction is not a good representation of the vertical or horizontal mixing of particles embedded in the flow (termed “tracers” in the atmospheric science community; see §12 of Vallis 2006). So while the enhanced vertical mixing (relative to Jupiter), via large-scale circulation cells present in hot Jovian atmospheres, may be a tempting and plausible mechanism for TiO to be maintained at high altitudes, a more detailed study of the interaction of tracers with the atmospheric circulation is required. (See Youdin & Mitchell 2010 for toy models of vertical, turbulent mixing in hot Jovian atmospheres and their implications for the exoplanetary radii.)

An observational way of distinguishing between different hot Jovian atmospheres is to examine their thermal phase curves (Cowan & Agol, 2008, 2011). The flux associated with the outgoing longwave radiation (OLR), denoted by , is the emergent flux from the longwave photosphere. In Figure 17, we show maps of for two models (, ; versus 2), where it is apparent that the chevron-shaped feature seen at bar resides near the longwave photosphere. The OLR flux may be used to construct the thermal phase curves, as was done for exo-Earths by Selsis, Wordsworth & Forget (2011). One has to first average the OLR flux over latitude,


The flux associated with the thermal phase curve, , can then be constructed using equation (7) of Cowan & Agol (2008),


where , the phase angle is denoted by and the summation coefficients are related via the following expressions:




We perform the fits, using the function described in the first expression in equation (49), to our computed curves by considering terms up to . We then obtain by transforming the and coefficients, using equations (50) and (51), to and coefficients. As pointed out by Cowan & Agol (2008), the odd sinusoidal modes do not contribute to the thermal phase curve: we set when is odd except for . This implies that information is lost when the latitudinally-averaged OLR flux is converted to the thermal phase curve. For example, for the model, we have while