Free Access
Issue
A&A
Volume 564, April 2014
Article Number A59
Number of page(s) 16
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/201323169
Published online 07 April 2014

© ESO, 2014

1. Introduction

For Earth’s atmosphere, irradiation from the Sun is the primary source of energy. Any model of the Earth’s atmosphere therefore needs a robust and accurate treatment of radiation transport. Global circulation models (GCMs) of the Earth are used for both weather prediction and climate research, and they include a dynamical core that solves some variant of the Navier-Stokes equations and a radiation scheme that calculates the radiative heating rate. The dynamical cores are tested using benchmarks (see e.g. Held & Suarez 1994; Reed & Jablonowski 2011; Ullrich et al. 2013), and both dynamical cores1 and radiation schemes (Ellingson et al. 1991; Collins et al. 2006; Oreopoulos et al. 2012) are tested through intercomparison projects.

GCMs have also been successfully applied to other solar system planets such as Jupiter, Saturn, Mars, and Venus (see for example Yamazaki et al. 2004; Müller-Wodarg et al. 2006; Hollingsworth & Kahre 2010; Lebonnois et al. 2011, respectively). In the past decade, GCMs have been used to study large-scale circulation on hot Jupiters (Showman & Guillot 2002; Showman et al. 2009; Rauscher & Menou 2012; Cho et al. 2008; Thrastarson & Cho 2010; Dobbs-Dixon & Lin 2008), a class of extrasolar planets which are approximately the size of Jupiter but orbit less than 0.1ua from their parent star. These planets, thought to have tidally locked circular orbits due to strong tidal interactions between the planet and its parent star (Baraffe et al. 2010), experience intense irradiation yielding a significant temperature contrast between the (permanent) day-side and night-side. Winds in the atmosphere of these planets are therefore expected to transport heat from the day-side to the night-side.

Atmospheric properties of hot Jupiters are obtained by various observational techniques, mainly transmission spectroscopy and secondary eclipse measurements. Brightness maps (Knutson et al. 2007; Majeau et al. 2012) and wind velocities (Snellen et al. 2010) are now accessible, and constraints on the composition are becoming available, but large uncertainties remain. Observations indicate a hotspot shifted eastward of the substellar point (Knutson et al. 2007; Majeau et al. 2012) and temperature contrasts smaller than what is expected for these planets without winds (Knutson et al. 2007, 2009), indicating transport of heat from the day-side to the night-side (Watkins & Cho 2010; Perez-Becker & Showman 2013). HD 209458b appears to have a temperature inversion in its upper atmosphere (Knutson et al. 2008; Burrows et al. 2007) while HD 189733b does not (Charbonneau et al. 2008; Barman 2008; Knutson et al. 2009), indicating that, despite similar orbital properties, hot Jupiters may have very different circulation patterns that still need to be understood. GCMs are therefore very valuable when trying to understand the increasing amount of observations of these systems.

Benchmarking of the dynamical cores of GCMs applied to hot Jupiters has been used to investigate stability of the codes and discrepancies between them (Heng et al. 2011; Menou & Rauscher 2009; Bending et al. 2013; Mayne et al. 2013, 2014; Polichtchouk et al. 2014). These benchmarks, and early GCMs applied to hot Jupiters, used simple, parametrised radiation schemes termed “Newtonian cooling” or “temperature forcing”, where the temperature is relaxed towards assumed equilibrium pressuretemperature (PT) profiles on a given timescale (Showman & Guillot 2002). Pressuretemperature profiles and timescales can be estimated using one-dimensional time-dependent radiative transfer calculations (Iro et al. 2005), though such an approach has flaws: (i) the equilibrium PT profiles used in the forcing may have a limited accuracy; (ii) radiative timescales may also have a limited accuracy and will vary in a non-trivial way as a function of latitude, longitude and depth; (iii) the forcing parametrisation itself may not be physically realistic, though the use of time-averaged equilibrium states when analysing model results may make this less of an issue; (iv) the model flexibility is poor since for each new planet modelled, the forcing must be changed.

Later studies used more complicated schemes such as flux-limited diffusion (Dobbs-Dixon & Lin 2008) and the two-stream approximation (Showman et al. 2009; Rauscher & Menou 2012; Dobbs-Dixon & Agol 2013). For the opacity treatment, grey schemes (Rauscher & Menou 2012), binning and averaging of the absorption coefficients (Dobbs-Dixon & Agol 2013) and the correlated-k method (Showman et al. 2009) have been used. The correlated-k method has also been used for retrieval analysis and characterisation of hot Jupiter atmospheres (Irwin et al. 2008) and to model brown dwarf atmospheres (Burrows et al. 1997). Brown dwarfs atmospheres have many similarities with hot Jupiter atmospheres (e.g. temperature range and composition), but local conditions are very different due to the strong irradiation from the parent stars for hot Jupiters. There is a notable lack of analysis of the accuracy of these schemes when applied to hot Jupiter-like atmospheres and of details on how opacities have been calculated from line lists, preventing rigorous comparison with results previously published in the literature. These are serious shortcomings in a field of research which develops quickly and will deliver more and more accurate data requiring reliable tools for their interpretation.

Table 1

List of molecules included in our opacity database with associated line list and partition function sources.

Both the two-stream approximation and the correlated-k method (see e.g. Thomas & Stamnes 2002) are widely used in GCM simulations of the Earth, and the literature on the methods’ applicability to the Earth atmosphere and their accuracy is extensive (see e.g. Toon et al. 1989; Meador & Weaver 1980; Zdunkowski et al. 1980; Goody et al. 1989; Lacis & Oinas 1991; Mlawer et al. 1997). They have both been found to yield results with satisfactory accuracy when comparing to more accurate solutions obtained from e.g. discrete ordinate (DO), line-by-line (LbL) calculations and when different schemes are compared through intercomparison projects (Ellingson et al. 1991; Collins et al. 2006; Oreopoulos et al. 2012). They are, however, still under investigation (Goldblatt et al. 2009) and are still one of the limiting factors of the accuracy of both weather prediction and climate modelling.

We are currently adapting the UK Met Office GCM, the Unified Model (UM), to hot Jupiter-like conditions. The main advantage of this model is its dynamical core, which solves the full 3D Euler equations and is coupled to a radiation scheme based on the correlated-k method and two-stream approximation. This GCM is state-of-the-art for both Earth and hot Jupiter atmospheric dynamics modelling. In previous papers we have tested and confirmed the dynamical core’s suitability to model hot Jupiter-like atmospheres (Mayne et al. 2013, 2014). This paper presents the adaptation of the radiation scheme, the Edwards-Slingo (ES) radiation scheme (Edwards & Slingo 1996), to conditions prevailing in hot Jupiter atmospheres. Opacities from high-temperature line lists have been calculated for the dominant absorbers in these atmospheres and the radiation scheme, using k-coefficients calculated from these opacities, has been tested against more accurate DO LbL calculations.

Observations of absorbing and scattering species in hot Jupiter atmospheres have so far been limited to the detection of molecular absorbers (Huitson et al. 2013; Wakeford et al. 2013; Tinetti et al. 2007), with some observations suggesting Rayleigh/Mie scattering clouds (Pont et al. 2008; Sing et al. 2013). Due to the large uncertainties related to scatterers in hot Jupiter atmospheres and the complexity it adds to radiation transport, we limit the discussions in this paper to purely absorbing atmospheres and postpone the inclusion of scattering to a future work.

The motivation for the present work is the lack of accurate tests and analysis of these radiation schemes now widely used by the community. This work will help in the implementation of similar schemes in the future and provide some guidelines for further progress. In Sect. 2 we first discuss our opacity data including the high-temperature line lists we use and the calculation of cross-sections from these line lists by using estimates for the line widths. In Sect. 3 we briefly summarise the implementation of the correlated-k method and two-stream approximation in the UM and move on to testing this scheme for gradually more complicated hot Jupiter-like atmospheres in Sect. 4. These tests will be useful to other groups for comparison and benchmark purposes. Our conclusions follow in Sect. 5. Combination of the adapted dynamical core and radiation schemes is in progress and will be presented in a future work, which will result in a state-of-the-art GCM that can be applied to a variety of exoplanet atmospheres.

2. Opacities

In this section we present the calculation of mass absorption coefficients from high-temperature line lists. We first discuss our opacity data, including line lists and line broadening parameters, and then provide details of our numerical implementation.

2.1. Opacity data

The dominant absorbers in hot Jupiter atmospheres with solar metallicity are H2O, CO, CH4, NH3, TiO, VO and H2–H2 and H2–He collision induced absorption (CIA; Burrows & Sharp 1999; Baraffe et al. 2010). References for the line lists and partition functions we use are given in Table 1. Where available we used line lists from the ExoMol project (Tennyson & Yurchenko 2012; H2O, named BT2, and NH3, named BYTe), while HITEMP (Rothman et al. 2010; CO) and the CIA extension to HITRAN (Richard et al. 2012; H2–H2 and H2–He CIA) are also used. For line lists from the ExoMol project, partition functions are calculated explicitly by summing up the energy levels. For CO we use the HITRAN partition functions (Rothman et al. 2009), while for TiO and VO we use the polynomial fits in Sauval & Tatum (1984) as recommended by B. Plez (priv. comm.). We use isotopic abundances, Ia, from Asplund et al. (2009).

Our methane line list is calculated using the Spherical Top Data System (STDS; Wenger & Champion 1998) setting the cut-off at J = 60. This line list has several flaws, e.g. important bands in the absorption spectrum are missing. We are aware that a new ExoMol CH4 line list should soon be released, but this will not change the main conclusions of this paper. We use the 12CH4 partition function from Wenger et al. (2008), and 13CH4 partition function from HITRAN (Rothman et al. 2009).

Due to the very large number of lines in these line lists, local thermodynamic equilibrium (LTE) is assumed to calculate level populations. Non-LTE effects could be important in the upper part of irradiated planet atmospheres. Indeed some observational works invoke non-LTE effects to explain the detection of strong emission features in hot Jupiters (see e.g. Waldmann et al. 2012), and there has been some work on modelling non-LTE effects in these atmospheres (Schweitzer & Hauschildt 2004; Barman et al. 2002). The significance of these effects still needs to be proven, however, and it is beyond the scope of the present work to consider non-LTE effects.

The line intensity for transition i can then be calculated by using (Thomas & Stamnes 2002) where T is the local temperature, c is the velocity of light, is the wavenumber of the transition, gu is the degeneracy of the upper level, El is the energy of the lower level, kB is Boltzmann’s constant, Q(T) is the partition function evaluated at T, Ai is the Einstein A-coefficient of the transition, e is the electron charge in CGS-Gaussian units, me is the electron mass and glflu is the gf-value of the transition with gl and flu being the degeneracy of the lower level and the oscillator strength, respectively. The quantities , gu, El and Ai, or alternatively glflu, are given in the line lists. For line lists in the HITRAN format, line intensities Si(T0) evaluated at a reference temperature T0 = 296K are given, and can be converted to any other temperature by using (3)which follows from Eq. (1) by taking Si(T)/Si(T0). Both ionisation and molecular dissociation are ignored.

We calculate absorption coefficients including both Doppler broadening and pressure broadening from collisions with H2 and He, which are the dominant species in hot Jupiter atmospheres, and ignore other broadening processes such as natural, self- and turbulent broadening as they are relatively unimportant compared to the former in hot Jupiter atmospheres. The pressure broadened width αL depends on both pressure and temperature in a complex way, but the relationship is often approximated as (Thomas & Stamnes 2002; Sharp & Burrows 2007) (4)where T0 and P0 are the reference pressure and temperature, respectively, and z is the perturbing species with Pz the partial pressure of species z. The total pressure broadened width is the sum of the H2 and He broadened widths. Pressure broadening parameters are, however, not included in the line lists and must be estimated from other sources. Table 2 lists our pressure broadened width data sources, which are mostly gathered from experimental data, and partly overlap with those used by Bailey & Kedziora-Chudczer (2012). In these sources, is listed as a function of quantum numbers at a given temperature and pressure. The temperature dependence exponent nz also depends on the quantum numbers of the transition, but less so than . We use the same nz for all transitions for a given species and broadener, a mean of the values found in the relevant paper(s) if more than one value is given. This approach is similar to that used by Sharp & Burrows (2007).

Table 2

Overview of our line width sources for pressure-induced broadening by hydrogen and helium.

In Fig. 1, we have plotted H2 broadened line widths as a function of Jl, the total rotational quantum number of the lower level of the transition, for the sources listed in Table 2. Since we have only for a small fraction of the transitions in the line lists, we model it as a linear function of Jl using least squares fits. This is similar to previous line width studies (Voronin et al. 2010; Burrows et al. 1997)2. Unfortunately has to be extrapolated to high Jl values, which we do by keeping it constant at the value where data for the highest Jl is available to avoid introducing additional complexity (see Fig. 1). A different extrapolation scheme was tested, where widths were decreased further before flattening out. The effect on the local absorption coefficient was quite large, but averaging over many lines yields negligible difference between the two extrapolation schemes. Our final line widths as a function of Jl are shown in Fig. 1 as solid lines.

thumbnail Fig. 1

H2 pressure broadened line widths for H2O, CH4, CO and NH3 from the sources listed in Table 2 at PH2 = 105  Pa and T = 296  K. Blue dots are the raw data points, the solid line is our fit and extrapolation. Conversion to other temperatures and pressures are done by using Eq. (4).

Open with DEXTER

Note that and nz from the referenced papers have been obtained at approximately room temperature and pressure, while we are using Eq. (4) to extrapolate to both very high temperatures and pressures. Unfortunately there is currently no viable alternative. Any calculation of absorption coefficients including pressure broadening has to extrapolate the line width data, but we note that details on extrapolation schemes used in the literature is sparse.

Having obtained both line intensities and widths, the line profile cross section, , and mass absorption coefficient, , can be calculated using (Thomas & Stamnes 2002): where is the mean molecular weight of molecule z in kg/molecule and (7)Here, αD is the Doppler broadened width and is the Voigt profile. The high velocity winds expected in the atmospheres of these planets will induce small Doppler shifts of the line centres, . Non-zero pressure will also cause a shift of the line centres, but data on pressure-induced line-shifts is even more spare than that on line widths. The net effect of both Doppler and pressure-induced line shifts are uncertain, and we choose to ignore both.

We calculate mass absorption coefficients by summing the Voigt profiles, (8) on a fixed wavenumber grid using a grid spacing of 1 × 10-3  cm-1.

CIA is a continuous opacity source, and is consequently provided in a different format. In the HITRAN format, A–B CIA data is tabulated as the absorption per density of species A per density of species B as a function of temperature on a fixed wavenumber grid with 1  cm-1 spacing. Here, A is always H2, and we multiply the CIA data by the density of H2 to make it consistent with the molecular line data.

2.2. Numerical considerations

The use of high-temperature line lists raises computational issues due to their size. In the HITRAN 2012 database (Rothman et al. 2013), the NH3 line list has about 4.6 × 104 lines, while the ExoMol NH3 line list, BYTe (Yurchenko et al. 2011), has about 1.1 × 109 lines. Consequently, calculating absorption coefficients from high-temperature line lists is significantly more computationally expensive than using smaller line lists such as HITRAN. In the literature this problem is often overcome by ignoring all lines with line intensities smaller than some value, sometimes evaluated at a fixed temperature (Sharp & Burrows 2007). The line intensity is, however, a strong function of temperature, and knowing where to apply the cut-off may be difficult.

A second cut-off has to be made, both for physical and computational reasons. According to Eq. (6) lines are infinite in extent and follow a Voigt profile. Unfortunately, real line profiles are not perfectly Voigtian (Thomas & Stamnes 2002). The Voigt profile is fairly accurate provided interactions between molecules are weak, but for stronger interactions effects such as collisional narrowing may occur. Line wings are most affected, and to avoid overestimating the line wing absorption, it is common practice to apply a cut-off at some distance d from the line centre (Freedman et al. 2008; Sharp & Burrows 2007). This distance may be fixed or be a function of pressure and/or temperature. In addition, evaluation of the Voigt profile is computationally expensive, and computing the line profiles to distances where it can be neglected adds an unnecessary computational cost.

To cope with these problems we have developed a scheme to combine the line wing cut-off with an elimination of unimportant weak lines to decrease computation time. The cut-off distance d is calculated on-the-fly by estimating when the line mass absorption coefficient has reached some value, . This is done by approximating the line profile as Lorentzian with a width equal to the sum of the Doppler and pressure broadened widths to facilitate analytical treatment and ensure that the profile width is not underestimated. This yields the following formula for d: (9)For very weak lines, , i.e. the value of the line mass absorption coefficient at the line centre is smaller than , and consequently the line can be ignored completely.

Lines are added one-by-one to the total mass absorption coefficient spectrum. The value of is chosen to be some fraction fAK of the latest value of the total mass absorption coefficient at the line centre as line profiles are summed up: . We use the abbreviation AK (adaptive ) to denote this cut-off method. Some lines can become unrealistically broad, however, so we include an upper limit on d of 100  cm-1. Note that this cut-off scheme cannot be used if the water continuum is included since it requires a cut-off at a fixed distance from the line centre (Clough et al. 1989).

The main motivation for using this scheme is computational efficiency and not physical considerations, and the final absorption coefficient will depend slightly on the order in which lines have been added up. The advantages are, however, that weak lines can be discarded on-the-fly taking into account the line intensity at the current temperature, making a simple cut-off in line intensity unnecessary. It also ensures that strong lines are computed to larger distances from the line centres than weaker lines. The current lack of robust line broadening schemes for conditions characteristic of hot Jupiters forces us to use such artificial schemes.

We have, however, tried to limit the impact of our treatments by testing other schemes used in the literature. We compare our line profile cut-off scheme to two other schemes: (i) a cut-off at a fixed distance dFW from the line centre (fixed width, FW) and (ii) a cut-off at a distance from the line centre given by the sum of the Doppler and pressure broadened widths multiplied by some factor fFF (fixed factor, FF). The former scheme is similar to that used when including the water continuum from Clough et al. (1989), where all lines have to be cut-off 25  cm-1 from the line centre, while the latter is similar to that used by Sharp & Burrows (2007). Note that when using FF, we still apply the upper limit of 100  cm-1 on d.

In Fig. 2, we show both the average absorption coefficient between 1000  cm-1 and 1001  cm-1 at 105Pa, 1500K and the computation time required for the three schemes as a function of the cut-off parameters. The code has been parallelised and runs using an Intel Xeon X5660 processor with 12 cores at 2.8GHz. 3

thumbnail Fig. 2

Left panel: arithmetic mean of the H2O absorption coefficient between 1000  cm-1 and 1001  cm-1 calculated using the adaptive (AK), fixed width (FW), fixed factor (FF) cut-off schemes as a function of the cut-off parameters (fAK, dFW and fFF, respectively) at 105  Pa, 1500  K. The mean absorption coefficients have been normalised by the value obtained using AK with fAK = 10-8. Right panel: computation time required using 12 cores at 2.8  GHz as a function of the cut-off parameter. We see that our adaptive cut-off scheme is about two orders of magnitude faster than the two other methods for a given average absorption coefficient.

Open with DEXTER

The three cut-off methods reach approximately the same average absorption coefficient for the largest values of the cut-off distances (see Fig. 2). Comparing the levels in the left-hand panels to the computation times in the right-hand panel, however, clearly shows the advantage of the AK method. At fAK = 10-6, this method reaches approximately the same level as FW at dFW = 102  cm-1 and FF at fFF = 103. The computation time is, however, more than two orders of magnitude smaller. Due to the uncertainties in line widths and the significant decrease in computation time, we have decided to adopt this scheme for all our molecules using fAK = 10-6 except for CO.

CO lines are divided into several clearly separate, narrow bands. Consequently, the absorption coefficient will vary by many orders of magnitude on the scale of the bands. The AK scheme is unsuited to such situations since it will tend to produce large cut-off distances at the beginning of the cross-section calculation. These line wings will normally be hidden by stronger lines, but for CO they become non-negligible due to the lack of strong lines in certain wavelength regions. For CO we therefore use the FF method with fFF = 102.

Absorption coefficients are tabulated on a log P,log T grid with 30 pressure points between 10-1Pa and 108Pa and 20 temperature points between 70K and 3000K, uniformly distributed on a logarithmic scale.

2.3. Molecular abundances

For the present work and tests we use the closed-form expressions for the chemical equilibrium abundances of H2O, CO, CH4 and NH3 from Burrows & Sharp (1999). We use solar-like elemental abundances, listed in Table 3. We assume the gas to be ideal and H2 and He partial pressures are calculated by assuming the atmosphere to be pure hydrogen and helium with an atomic hydrogen number fraction of 0.91183. This yields a mean molecular weight of 2.3376g   mol-1. For a mixture of gases, k-coefficients are calculated from an effective mass absorption coefficient table obtained by summing mass absorption coefficients for each species weighted by the respective mass mixing ratios, similar to the approach in Showman et al. (2009). Alternatively, the random overlap method can be used to combine k-coefficients for individual gases, though this is much more computationally expensive (Lacis & Oinas 1991). In Sect. 4.1, however, we only include absorption by H2O and use k-coefficients calculated from the H2O mass absorption coefficient and multiply the k-coefficients by the mass mixing ratio in each atmospheric layer.

Table 3

Our adopted solar-like elemental abundances (Asplund et al. 2009).

For TiO and VO we use for the present tests a simple parametrisation scheme to prescribe their abundances. Ti and V are thought to sequester deeper than 106Pa to 107Pa, while at low temperatures Ti and V will be bound to condensates (Showman et al. 2009; Fortney et al. 2006). We therefore parametrise the abundance of TiO and VO by assuming no absorbing TiO and VO to be present in the atmosphere for temperatures below Tcrit or pressure above Pcrit. For temperatures above Tcrit and pressures below Pcrit, however, we assume TiO and VO to be present, with partial pressures estimated by assuming all Ti bound in TiO and similarly all V bound in VO. This is a reasonable assumption since the abundances of both Ti and V are much smaller than that of oxygen, see Table 3. TiO and VO will therefore have a negligible effect on the availability of oxygen in the atmosphere. We use Tcrit = 1500K, Pcrit = 107Pa.

Non-equilibrium chemistry, due to photochemical or mixing processes, might be important in hot Jupiter atmospheres (Cooper & Showman 2006). We are currently developing a non-equilibrium chemistry scheme, but for the sake of present tests we limit our analysis to equilibrium chemistry. Neither the exact form of the TiO and VO abundance nor the inclusion of non-equilibrium chemistry will, however, alter the main conclusions of this paper.

3. The Edwards-Slingo radiation scheme and Atmo

Before discussing the tests we have performed to investigate the accuracy of the UM radiation scheme, we briefly summarise the formulation of the two-stream approximation (see e.g. Thomas & Stamnes 2002) and correlated-k method (Lacis & Oinas 1991; Goody et al. 1989) implemented in the Edwards-Slingo (ES) radiation scheme (Edwards & Slingo 1996) currently used in the UM. This scheme has been widely used in the meteorological community, see e.g. Sun (2011). We also describe the line-by-line discrete ordinate code Atmo used for comparisons. Both models approximate the atmosphere as plane parallel, a standard approximation for radiation schemes in GCMs, but we note that sphericity my be important particularly near the terminator.

3.1. The two-stream approximation

The ES radiation scheme uses the two-stream approximated version of the radiative transfer equation as formulated by Zdunkowski et al. (1980, 1982); Zdunkowski & Korb (1985) to obtain fluxes and heating rates, with details available in Edwards & Slingo (1996). Without scattering, the equations reduce to (10)where and are the diffuse 4 upward and downward fluxes, respectively, is the optical depth, is the Planck intensity at temperature T and D is the diffusivity. Per definition, the downward diffuse flux is always zero at the upper boundary. The direct component 5 of the flux, , is given by (11)where is the solar flux at the top of the atmosphere and μ0 = cosθ0 where θ0 is the solar zenith angle. The total upward and downward flux is the sum of the diffuse and direct components: (12)and we define the total flux, , as the upward flux minus the downward flux, i.e. (13)In the following, the explicit wavenumber dependence will be dropped to simplify the notation. Compared to the two-stream equation in Toon et al. (1989) this formulation is slightly different in that the thermal source function is πB and not 2πB/D. The choice πB ensures the correct thermal source flux independent of the choice of D, and is consistent with the formulation previously used by Dobbs-Dixon & Agol (2013); Rauscher & Menou (2012). Note that, for the thermal component, the ES radiation scheme solves the two-stream equations for the differential fluxes (total component of the flux less the Planck flux at the local temperature), see Edwards (1996). It is also worth noting that Showman et al. (2009) used the two-stream source function method (Toon et al. 1989) for the thermal component, which is exact in the no scattering case. We show that the two-stream approximation yields fairly accurate results for the thermal component without scattering, however, and plan to include scattering particles in the future. In the absence of scattering, the two-stream approximation yields the exact stellar component due to the separation of the flux into a direct (stellar) component and a diffuse (stellar) component, the latter being 0 without scattering6.

Since we ignore scattering, eliminating the single scattering albedo and backscattering coefficient, the only free parameter is the diffusivity factor D. It is related to the mean angle of the radiation by , where and is the mean zenith angle. In the present work we explore three different values for D: 1.66, originally due to Elasser (1942), from the discrete ordinate method with two quadrature points, and 2, obtained assuming an isotropic radiation field (Thomas & Stamnes 2002). See discussions in e.g. Thomas & Stamnes (2002); Edwards (1996) for more details.

The optical depth is given by (14)or in integrated form, assuming hydrostatic equilibrium: where and ζi are the mass absorption coefficient and mass mixing ratio of species i, respectively. The total mass absorption coefficient is kρ(z), and the total absorption coefficient is k(z) = kρ(z)ρ(z).

The heating rate is given by (17) where g is the gravitational acceleration, the mean molecular weight in g mol-1, P the pressure, R the ideal gas constant and F the total flux integrated over the entire spectrum. Hydrostatic equilibrium is assumed and the ideal gas equation is used to derive the final expression.

3.2. The correlated-k method

Currently the ES radiation scheme uses a combination of the exponential sum fitting of transmissions (ESFT) technique (Wiscombe & Evans 1977) and the correlated-k method (Lacis & Oinas 1991; Goody et al. 1989) to obtain k-coefficients, some details of which can be found in Sun (2011). In each band the absorption coefficients from the line-by-line wavenumber grid are reordered according to strength in increasing order and divided into nk subintervals. The spacing of these subintervals must be the same for each PT and would ideally be spaced logarithmically in k. Rather than using logarithmic k-intervals defined at a particular PT point, an average absorption coefficient is calculated from the top of the atmosphere down to an optical depth of one: where ζi(z) is the mass mixing ratio of species i, is the column density of species i down to τ = 1 and hydrostatic equilibrium has been assumed. This is similar to the approach of Hogan (2010) and provides an optimal sorting for the PT where each part of the spectrum is most important. We use an isothermal PT profile at 1116K, one of the temperatures in our PT grid, for this calculation as a compromise between day-side and night-side PT profiles of hot Jupiters. These average absorption coefficients are then used for the initial calculation of k-coefficients. In the following, the species index i will be dropped for ease of notation.

The k-coefficient for subinterval l is found by fitting the transmission, , to the weighted transmissions of the line-by-line coefficients over a set of nu column densities, uj, i.e. (20)where w is the weighting function and is the optimal k-coefficient in subinterval l. The probability of finding an absorption coefficient between k and k + dk is dg = f(k)   dk, with g(k) being the cumulative probability distribution. In the above equation, gl is the g-coordinate corresponding to the beginning of the subinterval for k-term l and gnk    +    1 is the g-coordinate for the end of k-term nk. The error for k-term l is defined as the root mean square (RMS) of the difference between the fitted exponential and exact integral for all column densities: (21)The total error in a band, ϵ, is defined as (22)The weights are defined to be normalised over each band: (23)and can be either a black-body spectrum at the current temperature, the stellar spectrum or uniform (in wavenumber).

A tolerance is set on the total error in a band, ϵmax. The number of k-terms in a band is chosen to be the smallest satisfying the criterion ϵ < ϵmax. Once this division into subintervals has been made, the fitting of optimal k-coefficients is repeated for each PT. These subsequent fits use the same subinterval spacing but with a reordering appropriate for each PT.

In the tests presented here, we use two values for ϵmax: 5 × 10-3 and 10-4, where ϵmax = 10-4 is expected to reduce the error from the correlated-k method significantly. We note, however, that the correlated-k method will not completely converge to the LbL solution even with many k-coefficients (small ϵmax): g(k) is calculated at each PT independently. If the absorption coefficient decreases with height at one wavelength and increases with height at a different wavelength within the same band, the two wavelengths will no longer correspond to the same value of g. A pseudo-monochromatic calculation where g is kept constant is therefore not equivalent to a proper monochromatic calculation except in a very few special cases, see e.g. Goody et al. (1989) for more details.

The maximum column density over which k-coefficients are to be fitted must be determined. The column density of species i is given by (24)where we have assumed a constant mass mixing ratio. We set the maximum column density using Eq. (24) with the maximum pressure and mixing ratio in our PT table. In the ES radiation scheme, the mass mixing ratio used in Eq. (19) is calculated from the maximum column density using Eq. (24).

To obtain the total flux in band b, nk pseudo-monochromatic calculations are performed, each yielding a total flux Fl. The band-integrated flux is then given by (25)The weight for the thermal component should ideally be identical to wl evaluated at the local temperature. For simplicity, however, the weights wl at the temperature where τ = 1 are adopted. For the stellar component, is the stellar spectrum at the top of the atmosphere. We later compare these weighting schemes to a uniform weighting scheme, and show that, using the band structure adopted here, the exact weighting scheme does not affect fluxes and heating rates to a significant degree.

We use a band structure very similar to that used by Showman et al. (2009), and we list our bands in Table 4. Three small modifications have been made compared to the bands listed in Showman et al. (2009): (i) the upper limit of band 30 has been reduced from 38   314  cm-1 to 28   000  cm-1 to reduce the memory usage of our correlated-k code; (ii) Band 31 has been added to capture absorption up to the small wavelength limit of our line lists; and (iii) Band 32 has been added to capture most of the stellar flux at small wavelengths.

Table 4

Bands used.

3.3. Band-averaged absorption coefficients

Band-averaged absorption coefficients have been applied to exoplanet GCMs (Dobbs-Dixon & Agol 2013) and stellar/substellar atmosphere radiation hydrodynamical models (see e.g. Freytag et al. 2010). In Dobbs-Dixon & Agol (2013), an average absorption coefficient is calculated in each band as: (26)where is the lower bound of bin b. The upper bound of the last band is defined as , where nb is the number of bands. The fluxes Fb are obtained by performing nb pseudo-monochromatic calculations, the total flux being the sum of the individual fluxes in each band. In Dobbs-Dixon & Agol (2013), the weighting function is a black-body spectrum evaluated at the local temperature for the thermal component, i.e. is the Planck mean in each band, known to be applicable in the optically thin limit. The stellar spectrum at the top of the atmosphere is used as weights for the stellar component. The bands were selected as in Showman et al. (2009), see Table 4.

Improved schemes utilising mean absorption coefficients exist, but we limit our discussion to the approach used by Dobbs-Dixon & Agol (2013) and compare its accuracy to a full correlated-k treatment. For this purpose, we show in Sect. 4 results obtained with band-averaged absorption coefficients designed to replicate the treatment used in Dobbs-Dixon & Agol (2013).

3.4. Atmo

In order to investigate the accuracy of both the correlated-k method and various two-stream approximations, we compare the ES radiation scheme to our line-by-line discrete ordinate code Atmo. It follows the method of the MARCS code (for a description, see Gustafsson et al. 2008), with some modifications.

The 1D plane-parallel radiative transfer equation is solved using the discrete ordinates method (see e.g. Thomas & Stamnes 2002), i.e. solving the radiative transfer equation for discrete ray directions μi, which are selected according to Gauss-Legendre quadrature. In this paper we use 16 rays, and we have checked the convergence by using up to 32 rays. The discrete ordinate equations are solved iteratively using the integral method, and the code has been parallelised to facilitate a high wavenumber resolution. By successively increasing the resolution, we found that a resolution on the order of ~10-3  cm-1 was necessary for the solution to have converged, i.e. about 5 × 107 wavenumber points. In contrast, we use about 300 pseudo-monochromatic calculations in the ES scheme, illustrating the large gain in computational efficiency achieved by using the correlated-k method.

4. Testing the correlated-k and two-stream approximations

In this section we test the accuracy of the ES radiation scheme by comparing it to Atmo. The accuracy is investigated for several different scenarios designed to test a range of physical conditions representative for hot Jupiters.

To ease comparison between the different two-stream approximations, we list the L1 norm of the error, (27)given here for the heating rate, where ES and Atmo are the heating rates from the ES radiation scheme and Atmo, respectively, and Pmin (Pmax) is the minimum (maximum) pressure in our calculations. This is a convenient measure to use when comparing errors between different two-stream approximations and opacity treatments, and represents the relative error of some quantity, in this case the heating rate, weighted by the current value of that quantity integrated over all pressures.

In Appendix A, we describe a very simple test where analytical solutions to both the two-stream approximated and the full radiative transfer equation exist. A grey opacity is used to eliminate errors from the correlated-k method. This scenario is used to test both the accuracy of the numerical solvers and the two-stream approximation in isolation. The test confirms that the numerical solvers of both the ES radiation scheme and Atmo have satisfactory accuracies. The value of D giving the most accurate results (D = 1.66) yields an error of about 1% for the flux and 10% for the heating rate.

We have also performed a test designed to be more realistic, but still minimising the error caused by the correlated-k method by including only H2–H2 collision induced (continuum) absorption (CIA) in an isothermal atmosphere without irradiation. For brevity, we do not describe the test here, but only summarise the main results. We found that the value of D giving the most accurate results, D = 1.66, yielded an error in the flux of <1%, and about 7.5% in the heating rate using ϵ = 5 × 10-3 (1 to 2k-coefficients in each band). Using a mean absorption coefficient in each band yields similar errors.

In the next sections, we describe three other tests in detail. Test 1 includes only absorption by H2O, a very important absorber in the atmospheres of both the Earth and hot Jupiters. Tests 2 and 3 are designed to test the accuracy for a typical well-mixed hot Jupiter night-side and day-side, respectively, where the day-side includes irradiation and absorption by TiO and VO in the upper atmosphere. A gravitational acceleration of 9.42  m/s2 is used, suitable for HD 209458b. Unless otherwise stated, k-coefficients were calculated using ϵmax = 5 × 10-3, a Planckian weighting scheme for the thermal component and the stellar spectrum for the stellar component, unless stated otherwise.

thumbnail Fig. 3

Left panel: fluxes obtained with ES radiation scheme using the two-stream approximation and correlated-k method obtained with D = 1.66 and ϵmax = 5 × 10-3 (dash-dotted, green), ϵmax = 10-4 (dashed, cyan), and mean absorption coefficients (dotted, blue), for an isothermal atmosphere with pure H2O absorption. The Atmo LbL DO result is also shown in this panel (solid, black) and is used to calculate the relative errors shown in the right panel (except for the mean absorption coefficient case since errors are too large). Relative errors stay below 30% throughout the atmosphere.

Open with DEXTER

4.1. Test 1: Pure H2O absorption in a high-temperature isothermal atmosphere

This test includes only absorption by H2O. The temperature is fixed to 1500K and the atmospheric domain extends from 10-1Pa to 108Pa, using 100 pressure points on a logarithmic scale. Irradiation at the upper boundary is not included and the lower boundary emits as a black body at 1500K with zero albedo. A constant mass mixing ratio of 3.3477    × 10-3 is adopted, which corresponds to the smallest mass mixing ratio predicted by Eq. (A4) in Burrows & Sharp (1999). Adopting ϵmax = 5 × 10-3 and ϵmax = 10-4 yield ~10 and ~100 k-coefficients in each band, respectively. Note, however, that the number of k-coefficients can vary significantly between different bands.

Figures 3 and 4 show the fluxes and heating rates, respectively, with corresponding relative errors. Since the atmosphere is isothermal, the upward flux is the Planck flux throughout the atmosphere. At the upper boundary, the downward flux is 0, while at the lower boundary it is the Planck flux due to the high optical depth. The total flux is therefore 0 at high optical depths, while at low optical depths the planet radiates as a black body, as expected, with a heating rate peak at ~104Pa. Errors in the heating rate generally stay below about 10% at pressures where the heating is significant, while using ϵmax = 10-4 yields significantly more accurate results. Table 5 shows L1 errors, indicating that D = 1.66 and yield the most accurate results. There is no significant difference between a Planckian and a uniform weighting (UW) scheme.

Using mean absorption coefficients (see Figs. 3 and 4 and Table 5) yield very inaccurate fluxes and heating rates. The flux is underestimated at a given pressure, while the heating rate peak occurs at pressures about two orders of magnitude smaller than the peak of the LbL DO result. The fact that one k-coefficient per band is not sufficient to resolve the opacity is reflected by the need for ~10 k-coefficients in each band to achieve a tolerance of ϵmax = 5 × 10-3 for H2O. The failure of this method is discussed in more detail in Sect. 4.4.

thumbnail Fig. 4

Same as Fig. 3, but for heating rates. Relative errors stay below 30% throughout the atmosphere. Note that in the region where the heating rate (magnitude) is large, the error remains small.

Open with DEXTER

Table 5

Computed flux (F) and heating rate () L1 norms for test 1, showing that the most accurate fluxes and heating rates are obtained with D = 1.66 and .

4.2. Test 2: mixed night-side hot Jupiter atmosphere

We next consider conditions representative of a real hot Jupiter atmosphere. We use a polynomial fit (Heng et al. 2011) to the night-side PT profile from Iro et al. (2005) with the smoothing described in Mayne et al. (2014), shown in Fig. 5. The temperature varies from about 400K in the upper atmosphere to above 1600K at 108Pa, consistent with the literature (see e.g. Fig. 6 in Showman et al. 2009; and Fig. 7 in Baraffe et al. 2008). Irradiation at the upper boundary is not included and the lower boundary emits as a black body with Tlb = T(Plb), where Tlb and Plb are the temperature and pressure at the lower boundary, respectively. From the PT profile in Fig. 5, Tlb = 1662K. All molecules in Table 1 are included except TiO and VO, with abundances as described in Sect. 2.3.

thumbnail Fig. 5

PT profile used in test 2. From the polynomial fit (Heng et al. 2011) to the night-side profile of HD 209458b from Iro et al. (2005) with the smoothing described in Mayne et al. (2014). The temperature varies from about 400  K high in the atmosphere to above 1600  K in the deeper layers.

Open with DEXTER

Fluxes and heating rates with relative errors are plotted in Figs. 6 and 7, with L1 errors given in Table 6, from which it is clear that D = 1.66 yields the most accurate fluxes and heating rates overall.

thumbnail Fig. 6

Left panel: fluxes obtained in test 2 with the ES radiation scheme using the PT profile in Fig. 5. All opacity sources listed in Table 1 are included except TiO and VO. The Atmo LbL DO result is also shown in this panel and is used to calculate the relative errors shown in the right panel. For explanation of the different lines, see the caption of Fig. 3.

Open with DEXTER

thumbnail Fig. 7

Same as Fig. 6 for the heating rates. Note that when the heating rate is very small, relative errors may become large. The effect on the heating budget will be small, however, so we do not consider this as a problem.

Open with DEXTER

Table 6

Computed flux (F) and heating rate () L1 norms for test 2, showing that the most accurate fluxes and heating rates are obtained with D = 1.66.

When the flux or heating rate is close to zero, the relative error can become large. We do not consider this as a problem, however, as it does not have a large impact on the atmospheric heat budget. The overall results are similar to those obtained in Sect. 4.1. Note that calculations with mean absorption coefficients in each band significantly underestimates the flux and result in heating rate peaks with the wrong amplitude.

4.3. Test 3: mixed day-side hot Jupiter atmosphere

Our last test adopts conditions suitable to the day-side of hot Jupiters. We use the polynomial fit (Heng et al. 2011) to the day-side PT profile of HD 209458b from Iro et al. (2005) with the smoothing described in Mayne et al. (2014), plotted in Fig. 8, consistent with the literature (Showman et al. 2009; Baraffe et al. 2008). The thermal and stellar components of the flux are calculated separately and then summed to obtain the total flux and heating rate. For the thermal component, the lower boundary again emits as a black body at Tlb = T(Plb), i.e. 1998K using the PT profile in Fig. 8. Note that, due to the separation of the intensity into direct and diffuse components and the our neglect of scattering, the stellar component of the intensity will only be subject to errors caused by the correlated-k method.

thumbnail Fig. 8

PT profile used in test 3. From the polynomial fit (Heng et al. 2011) to the day-side profile of HD 209458b from Iro et al. (2005) with the smoothing described in Mayne et al. (2014). The temperature varies from about 1400  K in the upper atmosphere to about 2000  K in the deeper layers.

Open with DEXTER

We assume an orbital distance aorbit = 0.047ua and a parent star effective temperature = 5785K and radius Rstar = RSun. Using Stefan-Boltzmann’s law, the stellar irradiation at the top of the planet’s atmosphere is given by (28)We adopt a zero solar zenith angle and use a solar spectrum from Kurucz7. At smaller wavelengths than available, we set the stellar flux to zero, while at larger wavelengths we extrapolate using a black-body spectrum with the effective temperature of the Sun (T = 5785K).

Figures 9 and 10 show the thermal flux and heating rate with relative errors. The flux error is small in regions with non-negligible flux. The heating rate error also remains small in regions with significant cooling. This is confirmed by the computed L1 norms listed in Table 7. A diffusitivy of D = 1.66 yields the smallest error, and it is approximately halved by decreasing ϵmax to 10-4. Whether a black-body or uniform weighting scheme is used does not affect the accuracy significantly, but using a mean absorption coefficient does result in significant errors.

thumbnail Fig. 9

Left panel: thermal component of the flux as a function of total pressure for test 3. The Atmo LbL DO result is also shown in this panel and is used to calculate the relative errors shown in the right panel. For explanation of the different lines, see the caption of Fig. 3. At high pressures the relative error becomes large, similar to that seen in Fig. 7 for the heating rate. As the flux itself is small, however, this is not a problem.

Open with DEXTER

thumbnail Fig. 10

Same as Fig. 9 for the thermal component of the heating rate. Relative errors become unreasonably large only where the heating rate is very small.

Open with DEXTER

Table 7

Computed flux (F) and heating rate () L1 norms for the thermal component in test 3.

Table 8

Computed flux (F) and heating rate () L1 norms for the stellar component in test 3.

Figures 11 and 12 show the stellar flux and heating rate with relative errors. At the top of the atmosphere, the flux is –6.092 × 105  W/m2, as prescribed, and is subsequently absorbed. The heating rate is positive, i.e. the atmosphere is heated due to the absorption of stellar radiation, as expected. The accuracy is acceptable, the error in the flux stays below 10%, while the heating rate error also stays below 10% in the regions with strong heating. This is reflected in the L1 errors listed in Table 8. Using ϵmax = 10-4 significantly reduces the error from the correlated-k method, and changing the weighting scheme does not alter the results significantly. The use of an average absorption coefficient, however, is seen to still result in significant errors. See Sect. 4.4 for discussion and more details.

thumbnail Fig. 11

Left panel: stellar component of the flux as a function of total pressure for test 3 obtained with ϵmax = 5 × 10-3 (dash-dotted, green) and ϵmax = 10-4 (dashed, cyan). The Atmo LbL DO result is also shown in this panel (solid, black) and is used to calculate the relative errors shown in the right panel. Errors are small, and using ϵmax = 10-4 almost completely eliminates errors in the ES radiation scheme.

Open with DEXTER

thumbnail Fig. 12

Same as Fig. 11 for the stellar component of the heating rate.

Open with DEXTER

The total flux and heating rate, obtained by summing up the stellar and thermal components of the flux and heating rate, are shown in Figs. 13 and 14, respectively. The main region of heating and cooling, seen separately in Figs. 10 and 12, respectively, are still clearly distinguishable in Fig. 14. The atmosphere is heated at low pressures and cooled slightly at higher pressures. Note that errors remain satisfactory small for relevant (i.e. non zero) values of the heating rate, as also shown in Table 9. The error introduced by using mean absorption coefficients is significant.

thumbnail Fig. 13

Left panel: total flux as a function of total pressure for test 3. The Atmo LbL DO result is also shown in this panel and is used to calculate the relative errors shown in the right panel. For explanation of the different lines, see the caption of Fig. 3. Again, relative errors become unreasonably large only where the flux is very small.

Open with DEXTER

thumbnail Fig. 14

Same as Fig. 13 for the total heating rate. Relative errors become unreasonably large only where the heating rate is very small, with a negligible effect on the heating budget.

Open with DEXTER

Table 9

Computed flux (F) and heating rate () L1 norms for the total flux and heating rate in test 3.

4.4. Discussion of the failure of mean absorption coefficients

Inspection of Figs. 3, 4, 6, 7 and 11 to 14 suggests systematic deviations of results based on the band-averaged absorption coefficients. For a given pressure, the thermal and stellar fluxes are underestimated, often resulting in heating rate peaks occurring at lower pressures and having the wrong magnitude. In an attempt to explain this behaviour, we first consider the direct stellar component.

The band-integrated direct stellar component of the flux is given by (29)where the atmospheric slab has been assumed to be homogeneous where uρ is the mass column density down to some height z. Using a mean absorption coefficient instead, the corresponding flux is (30)For simplicity we assume the incoming stellar radiation at the top of the atmosphere is wavenumber independent within a given band. Using a mean absorption coefficient then implies (31)Within a band the absorption coefficient will vary by orders of magnitude, causing some regions in the band to have a small transmission and others to have a large transmission. The mean in Eq. (26) is an arithmetic mean, i.e. the largest values of will dominate . Regions with high transmission due to small will be overshadowed by this large mean, causing the overall transmission to be underestimated. This explains the deviation in the flux seen in Fig. 11.

A similar argument can be used in the thermal region, but upward and downward radiation need to be considered separately. The radiative transfer equation reads (32)where s is the path over which radiation travels. We first consider the isothermal case, where the upward radiation is constant and equal to the black-body flux throughout the atmosphere. At the top of the atmosphere, the downward flux is zero, i.e. the change in intensity will, according to Eq. (32), be dominated by thermal emission (). Using a band-mean absorption coefficient effectively increases , which in Eq. (32) yields a larger intensity at a given s or pressure. The downward radiation contributes negatively to the total flux, i.e. the total flux will be smaller for a given pressure, as seen in Fig. 3.

If the atmosphere has non-zero temperature gradients, the upward flux will also depend on pressure. In both PT profiles used here, the temperature decreases with height overall. At the lower boundary the upward intensity is simply the Planck intensity, i.e. the right-hand side of Eq. (32) is zero. As the temperature decreases, will generally decrease, causing . The upward flux is therefore dominated by absorption, and effectively increasing will cause the upward flux to become smaller for a given s or pressure. This explains why the total flux at the top of the atmosphere is underestimated when using a band-mean absorption coefficient, as seen in Figs. 3, 6 and 13.

These results show that large errors in both fluxes and heating rates may occur when the mean opacity scheme described in Dobbs-Dixon & Agol (2013) is applied in hot Jupiter GCMs. Improved mean opacity schemes have been developed by the stellar atmosphere community (see e.g. Nordlund 1982; Skartlien 2000), which may be applicable to hot Jupiter atmospheres. Further developments of these improved schemes may be needed, however, as they rely on correlations induced by strong vertical stratification, and longitude-latitude-dependent stellar heating has not been considered. Further discussion is, however, beyond the scope of the present work.

5. Conclusions

The accuracy of radiation schemes used in GCMs has been studied extensively for Earth-like conditions, but detailed analysis for hot Jupiter-like conditions are lacking. In this paper we have analysed the accuracy and uncertainties in state-of-the-art radiation schemes used in several GCMs applied to hot Jupiters. Opacity sources and calculation of absorption coefficients from high-temperature line lists have been discussed. We present a line profile cut-off scheme that decreases the computation time required to calculate absorption coefficients by a factor of ~100 compared to other methods used in the literature, while still giving accurate results. Both the two-stream approximation and correlated-k method’s applicability to hot Jupiter atmospheres have been analysed by comparing the Edwards-Slingo radiation scheme to discrete ordinate line-by-line calculations.

The ES radiation scheme’s performance in these tests shows that we have successfully adapted it to hot Jupiter-like atmospheres. Our main conclusions are:

  • Pressure broadening parameters for high-temperaturemolecular lines are very uncertain and usually extrapolated fromroom temperature and pressure and small quantum numbers.Improvements in this area will become important as higheraccuracy will be required to analyse results from future exoplanetcharacterisation projects (e.g. JWST, EChO, SPHERE andELT).

  • A diffusivity factor of D = 1.66, already widely used in both Earth and hot Jupiter GCMs, yields the smallest errors from the two-stream approximation, although is only slightly less accurate.

  • About 10k-coefficients in each band for molecular line absorption yield satisfactory accuracy. Using ~100 k-coefficients per band does improve the overall accuracy, but errors decrease by less than 50%, while the radiative transfer computation time increases by a factor of 10. We therefore choose to adopt the former as a balance between accuracy and computational cost.

  • Both the two-stream approximation and the correlated-k method contribute non-negligibly to the total error, with overall heating rate errors of 10% in regions with significant heating/cooling. Flux errors are similar or smaller.

  • Whether a black-body spectrum, solar spectrum or uniform (in wavenumber) weighting scheme is used has little effect on the overall accuracy given the band structure used here (Table 4). We therefore choose to adopt a uniform weighting scheme, enabling the use of the same k-coefficients in both the thermal and stellar spectral regions and for different irradiation spectra.

  • Using a mean absorption coefficient in each band, as in Dobbs-Dixon & Agol (2013), yields inaccurate fluxes and heating rates for molecular absorption. Heating rate errors can reach 100% or more, even in regions with significant heating. Band-averaged absorption coefficients should thus be used with caution.

Any radiation scheme applied to hot Jupiters should be checked against the tests we have presented here. These tests and the detailed descriptions of our methods and approximations will be useful for future adaptation of radiation schemes in other GCMs. Current observational constraints on exoplanets do not require the level of accuracy we have applied in this work. The field develops at an amazing pace, however, and modellers should now develop the best theoretical and numerical tools to tackle the challenges posed by the increasing accuracy expected from future large observational projects.


2

Note that for CO the line width is approximately constant as a function of Jl, and we therefore use a simple mean instead of a linear fit, as shown in Fig. 1.

3

For NH3, having the largest of our adopted line lists with >1.1 billion lines, computation of absorption coefficients take about 11 days using our adaptive cut-off scheme on a computer with an Intel Xeon X5660 processor with 12 cores at 2.8GHz.

4

Diffuse radiation is in general defined as thermal and scattered radiation. Since we ignore scattering, the diffuse radiation will simply be the thermal component.

5

The direct component is defined as the unscattered part of the stellar radiation.

6

Note that for very hot planets, or planets with a low mass parent star (i.e. with a flux shifted towards the infrared), the thermal and stellar components will overlap significantly in wavelength. We specifically use the terminology thermal and stellar components over infrared and visible wavelength regions to avoid confusion.

Acknowledgments

We would like to thank Derek Homeier, Jonathan Tennyson, Bernd Freytag, Bertrand Plez, France Allard and Travis Barman for insightful discussions. This work is supported by the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013 Grant Agreement No. 247060) and by the Consolidated STFC grant ST/J001627/1. This work is also partly supported by the Royal Society award WM090065. The calculations for this paper were performed on the DiRAC Facility jointly funded by STFC, the Large Facilities Capital Fund of BIS, and the University of Exeter.

References

Appendix A: Test 0

This test is based on a grey atmosphere without scattering and irradiation at the top of the atmosphere, and a lower boundary that emits as a perfect black body at a temperature Tc. These assumptions are consistent with the thermal component of the radiation. To facilitate analytical treatment, we make an additional assumption: the lower boundary is located at a constant optical depth τ = τ. This is done in both the ES radiation scheme and Atmo by explicitly keeping the total mass absorption coefficient, kρ, constant as a function of pressure and placing the lower boundary at a constant pressure.

Appendix A.1: Analytical solutions

Analytical solutions of the two-stream approximated and full radiative transfer equation are available under these specific conditions and are provided below. The analytical solutions are compared to the numerical solutions obtained by the ES radiation scheme and to the discrete ordinate solution from Atmo.

Appendix A.1.1: The two-stream approximation

The two-stream approximated radiative transfer equation in the thermal region, ignoring scattering, is given by Eq. (10). We now drop the diffuse flux subscript since stellar irradiation is ignored. Using Eq. (16) and assuming hydrostatic equilibrium, the optical depth can be related to the pressure by (A.1)since kρ(ν,P) = kρ is assumed to be independent of both frequency ν and pressure P. The optical depth is therefore proportional to pressure and substitution can be done using the equation above.

Integrating Eq. (10) with respect to frequency yields (A.2)where Stefan-Boltzmann’s law has been used. The above equation is a simple inhomogeneous linear first order differential equation in optical depth, τ, and can be solved using traditional techniques. The homogeneous solution, i.e. ignoring the Planck emission, is given by (A.3)where A± is determined by boundary conditions, while the particular solution in this case is given by (A.4)which yields the complete solution (A.5)At the upper boundary, i.e. τ = 0, we have (A.6)At the lower boundary, which we place at an optical depth of τ = τ, we have (A.7)The upwelling, downwelling and total fluxes are therefore The heating rate is given by Eq. (17), and using Eq. (14), we get

Appendix A.1.2: The angularly dependent radiative transfer equation

The full angular dependent (but still azimuthally averaged) radiative transfer equation without scattering is given by (Thomas & Stamnes 2002) (A.13)where u = cosθ. From above, the general solution is given by (A.14)Note that a discrete ordinate method gives the same equation and solution, except u is replaced by the quadrature points ui. No downward radiation at the upper boundary implies , which yields (A.15)Perfect black-body radiation in the upward direction at the lower boundary implies , which yields (A.16)Note that is anisotropic in the downward direction except in the limit of high optical depths, τ → ∞. The intensity at a given optical depth in the downward direction, u < 0, is dictated by the amount of atmosphere above it, in the direction of the radiation, emitting thermally. The “effective optical depth” or optical path is higher for smaller values of , and the intensity consequently cannot be isotropic. The exception is at high optical depths, where the atmosphere becomes optically thick in all directions.

The upward flux is (A.17)while the downward flux is This integral does not have a simple closed-form solution, but can be found numerically. Integrating over all wavenumbers and using Stefan-Boltzmann’s law, the total flux is given by (A.20)i.e. dictated by this integral and clearly not equivalent to Eq. (A.10). Note that the two-stream approximation effectively evaluates the integral in Eq. (A.20) using a single quaderature point which can be chosen using e.g. Gauss-Legendre quadrature or an empirical fit. The heating rate, Eq. (17), is similarly given by

thumbnail Fig. A.1

Left panel: fluxes obtained using the two-stream flux in Eq. (A.10) and exact flux in Eq. (A.20) obtained with (dotted, red), D = 1.66 (dash-dotted, green), D = 2 (dashed, cyan) and solving the fully angular dependent radiative transfer equation (solid, black). Right panel: calculated relative errors in the two-stream fluxes. Relative errors become unreasonably large only where the flux is very small.

Open with DEXTER

thumbnail Fig. A.2

Same as Fig. A.1 for heating rates. Relative errors become unreasonably large only where the heating rate is very small.

Open with DEXTER

Appendix A.2: Accuracy of the two-stream approximation

Fluxes and heating rates with errors are plotted in Figs. A.1 and A.2 using the solutions in Eqs. (A.10), (A.12), (A.20) and (A.22). At small optical depths, the flux is equal to the black-body flux while at large optical depths, the flux is zero, as expected. The heating rate is zero at both low and high optical depths, while at intermediate optical depths the atmosphere is cooled (the heating rate is negative). Interestingly, the relative error in both flux and heating rate approaches unity at large optical depths, a consequence of the two-stream solutions approaching zero faster than the full solution. We do not consider this a problem, however, since both the flux and heating rate are close to zero in this region.

In Table A.1 we show the calculated L1 norms using the analytical solutions derived above. The smallest errors are achieved with D = 1.66, but using only yields slightly larger errors. This is verified by looking at the relative error in Figs. A.1 and A.2.

thumbnail Fig. A.3

Relative error in the numerical solutions from the ES radiation scheme (dash-dotted, green) and Atmo (solid, black) calculated using the analytical solution in Eqs. (A.10), (A.12), (A.20) and (A.22). Relative numerical errors become large at large optical depths, caused by both fluxes and heating rates being very small.

Open with DEXTER

Table A.1

Computed flux (F) and heating rate () L1 norms for test 0 using the analytical solutions, thereby eliminating the errors from the numerical solution schemes.

It is worth noting that the different values for D yield different convergence towards zero heating rate at low optical depths, evident in the right-hand panel of Fig. A.2, caused by the factor D in Eq. (A.12). Comparing Eqs. (A.12) and (A.22), it is clear that only D = 2 will yield the correct behaviour of the heating rate at low optical depths. The effect on the heating rate itself is small, however, and D = 1.66 yields the most correct heating rate overall.

Table A.2

Computed flux (F) and heating rate () L1 norms for test 0 comparing the numerical and analytical solutions to check the accuracy of the numerical schemes.

Appendix A.3: Accuracy of the numerical scheme

We use the analytical expressions in Eqs. (A.10), (A.12), (A.20) and (A.22) to estimate the errors in the numerical solution schemes in both the ES radiation scheme and Atmo. We have plotted the numerical error in Fig. A.3 as a function of optical depth, and given the L1 errors in Table A.2. The errors are small for small optical depths, while at large optical depths the errors increase significantly. The error in the flux and heating rate reach 10% at about τ = 10 and τ = 4, respectively, for the ES radiation scheme, while Atmo is accurate to a somewhat larger optical depth. The L1 errors reflects this, keeping in mind that the error in Atmo is also caused by a finite number of rays in the Gaussian quadrature, which may become important at the accuracy level of the numerical solver. Both numerical schemes are seen to yield errors significantly smaller than errors caused by the two-stream approximation. This confirms that both numerical solvers yield satisfactory accuracy.

All Tables

Table 1

List of molecules included in our opacity database with associated line list and partition function sources.

Table 2

Overview of our line width sources for pressure-induced broadening by hydrogen and helium.

Table 3

Our adopted solar-like elemental abundances (Asplund et al. 2009).

Table 4

Bands used.

Table 5

Computed flux (F) and heating rate () L1 norms for test 1, showing that the most accurate fluxes and heating rates are obtained with D = 1.66 and .

Table 6

Computed flux (F) and heating rate () L1 norms for test 2, showing that the most accurate fluxes and heating rates are obtained with D = 1.66.

Table 7

Computed flux (F) and heating rate () L1 norms for the thermal component in test 3.

Table 8

Computed flux (F) and heating rate () L1 norms for the stellar component in test 3.

Table 9

Computed flux (F) and heating rate () L1 norms for the total flux and heating rate in test 3.

Table A.1

Computed flux (F) and heating rate () L1 norms for test 0 using the analytical solutions, thereby eliminating the errors from the numerical solution schemes.

Table A.2

Computed flux (F) and heating rate () L1 norms for test 0 comparing the numerical and analytical solutions to check the accuracy of the numerical schemes.

All Figures

thumbnail Fig. 1

H2 pressure broadened line widths for H2O, CH4, CO and NH3 from the sources listed in Table 2 at PH2 = 105  Pa and T = 296  K. Blue dots are the raw data points, the solid line is our fit and extrapolation. Conversion to other temperatures and pressures are done by using Eq. (4).

Open with DEXTER
In the text
thumbnail Fig. 2

Left panel: arithmetic mean of the H2O absorption coefficient between 1000  cm-1 and 1001  cm-1 calculated using the adaptive (AK), fixed width (FW), fixed factor (FF) cut-off schemes as a function of the cut-off parameters (fAK, dFW and fFF, respectively) at 105  Pa, 1500  K. The mean absorption coefficients have been normalised by the value obtained using AK with fAK = 10-8. Right panel: computation time required using 12 cores at 2.8  GHz as a function of the cut-off parameter. We see that our adaptive cut-off scheme is about two orders of magnitude faster than the two other methods for a given average absorption coefficient.

Open with DEXTER
In the text
thumbnail Fig. 3

Left panel: fluxes obtained with ES radiation scheme using the two-stream approximation and correlated-k method obtained with D = 1.66 and ϵmax = 5 × 10-3 (dash-dotted, green), ϵmax = 10-4 (dashed, cyan), and mean absorption coefficients (dotted, blue), for an isothermal atmosphere with pure H2O absorption. The Atmo LbL DO result is also shown in this panel (solid, black) and is used to calculate the relative errors shown in the right panel (except for the mean absorption coefficient case since errors are too large). Relative errors stay below 30% throughout the atmosphere.

Open with DEXTER
In the text
thumbnail Fig. 4

Same as Fig. 3, but for heating rates. Relative errors stay below 30% throughout the atmosphere. Note that in the region where the heating rate (magnitude) is large, the error remains small.

Open with DEXTER
In the text
thumbnail Fig. 5

PT profile used in test 2. From the polynomial fit (Heng et al. 2011) to the night-side profile of HD 209458b from Iro et al. (2005) with the smoothing described in Mayne et al. (2014). The temperature varies from about 400  K high in the atmosphere to above 1600  K in the deeper layers.

Open with DEXTER
In the text
thumbnail Fig. 6

Left panel: fluxes obtained in test 2 with the ES radiation scheme using the PT profile in Fig. 5. All opacity sources listed in Table 1 are included except TiO and VO. The Atmo LbL DO result is also shown in this panel and is used to calculate the relative errors shown in the right panel. For explanation of the different lines, see the caption of Fig. 3.

Open with DEXTER
In the text
thumbnail Fig. 7

Same as Fig. 6 for the heating rates. Note that when the heating rate is very small, relative errors may become large. The effect on the heating budget will be small, however, so we do not consider this as a problem.

Open with DEXTER
In the text
thumbnail Fig. 8

PT profile used in test 3. From the polynomial fit (Heng et al. 2011) to the day-side profile of HD 209458b from Iro et al. (2005) with the smoothing described in Mayne et al. (2014). The temperature varies from about 1400  K in the upper atmosphere to about 2000  K in the deeper layers.

Open with DEXTER
In the text
thumbnail Fig. 9

Left panel: thermal component of the flux as a function of total pressure for test 3. The Atmo LbL DO result is also shown in this panel and is used to calculate the relative errors shown in the right panel. For explanation of the different lines, see the caption of Fig. 3. At high pressures the relative error becomes large, similar to that seen in Fig. 7 for the heating rate. As the flux itself is small, however, this is not a problem.

Open with DEXTER
In the text
thumbnail Fig. 10

Same as Fig. 9 for the thermal component of the heating rate. Relative errors become unreasonably large only where the heating rate is very small.

Open with DEXTER
In the text
thumbnail Fig. 11

Left panel: stellar component of the flux as a function of total pressure for test 3 obtained with ϵmax = 5 × 10-3 (dash-dotted, green) and ϵmax = 10-4 (dashed, cyan). The Atmo LbL DO result is also shown in this panel (solid, black) and is used to calculate the relative errors shown in the right panel. Errors are small, and using ϵmax = 10-4 almost completely eliminates errors in the ES radiation scheme.

Open with DEXTER
In the text
thumbnail Fig. 12

Same as Fig. 11 for the stellar component of the heating rate.

Open with DEXTER
In the text
thumbnail Fig. 13

Left panel: total flux as a function of total pressure for test 3. The Atmo LbL DO result is also shown in this panel and is used to calculate the relative errors shown in the right panel. For explanation of the different lines, see the caption of Fig. 3. Again, relative errors become unreasonably large only where the flux is very small.

Open with DEXTER
In the text
thumbnail Fig. 14

Same as Fig. 13 for the total heating rate. Relative errors become unreasonably large only where the heating rate is very small, with a negligible effect on the heating budget.

Open with DEXTER
In the text
thumbnail Fig. A.1

Left panel: fluxes obtained using the two-stream flux in Eq. (A.10) and exact flux in Eq. (A.20) obtained with (dotted, red), D = 1.66 (dash-dotted, green), D = 2 (dashed, cyan) and solving the fully angular dependent radiative transfer equation (solid, black). Right panel: calculated relative errors in the two-stream fluxes. Relative errors become unreasonably large only where the flux is very small.

Open with DEXTER
In the text
thumbnail Fig. A.2

Same as Fig. A.1 for heating rates. Relative errors become unreasonably large only where the heating rate is very small.

Open with DEXTER
In the text
thumbnail Fig. A.3

Relative error in the numerical solutions from the ES radiation scheme (dash-dotted, green) and Atmo (solid, black) calculated using the analytical solution in Eqs. (A.10), (A.12), (A.20) and (A.22). Relative numerical errors become large at large optical depths, caused by both fluxes and heating rates being very small.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.