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/00046361/201323169  
Published online  07 April 2014 
Accuracy tests of radiation schemes used in hot Jupiter global circulation models
^{1} Astrophysics Group, University of Exeter, Exeter EX4 4QL, UK
email: D.S.Amundsen@exeter.ac.uk
^{2} Met Office, Exeter EX1 3PB, UK
Received: 2 December 2013
Accepted: 3 February 2014
The treatment of radiation transport in global circulation models (GCMs) is crucial for correctly describing Earth and exoplanet atmospheric dynamics processes. The twostream approximation and correlatedk method are currently stateoftheart approximations applied in both Earth and hot Jupiter GCM radiation schemes to facilitate the rapid calculation of fluxes and heating rates. Their accuracy have been tested extensively for Earthlike conditions, but verification of the methods’ applicability to hot Jupiterlike conditions is lacking in the literature. We are adapting the UK Met Office GCM, the Unified Model (UM), for the study of hot Jupiters, and present in this work the adaptation of the EdwardsSlingo radiation scheme based on the twostream approximation and the correlatedk method. We discuss the calculation of absorption coefficients from hightemperature line lists and highlight the large uncertainty in the pressurebroadened line widths. We compare fluxes and heating rates obtained with our adapted scheme to more accurate discrete ordinate (DO) linebyline (LbL) calculations ignoring scattering effects. We find that, in most cases, errors stay below 10% for both heating rates and fluxes using ~10 kcoefficients in each band and a diffusivity factor D = 1.66. The twostream approximation and the correlatedk method both contribute nonnegligibly to the total error. We also find that using bandaveraged absorption coefficients, which have previously been used in radiativehydrodynamical simulations of a hot Jupiter, may yield errors of ~100%, and should thus be used with caution.
Key words: radiative transfer / opacity / planets and satellites: atmospheres / planets and satellites: gaseous planets
© 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 NavierStokes 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 cores^{1} 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üllerWodarg et al. 2006; Hollingsworth & Kahre 2010; Lebonnois et al. 2011, respectively). In the past decade, GCMs have been used to study largescale circulation on hot Jupiters (Showman & Guillot 2002; Showman et al. 2009; Rauscher & Menou 2012; Cho et al. 2008; Thrastarson & Cho 2010; DobbsDixon & Lin 2008), a class of extrasolar planets which are approximately the size of Jupiter but orbit less than 0.1 ua 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) dayside and nightside. Winds in the atmosphere of these planets are therefore expected to transport heat from the dayside to the nightside.
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 dayside to the nightside (Watkins & Cho 2010; PerezBecker & 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 pressure−temperature (P–T) profiles on a given timescale (Showman & Guillot 2002). Pressure−temperature profiles and timescales can be estimated using onedimensional timedependent radiative transfer calculations (Iro et al. 2005), though such an approach has flaws: (i) the equilibrium P–T profiles used in the forcing may have a limited accuracy; (ii) radiative timescales may also have a limited accuracy and will vary in a nontrivial way as a function of latitude, longitude and depth; (iii) the forcing parametrisation itself may not be physically realistic, though the use of timeaveraged 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 fluxlimited diffusion (DobbsDixon & Lin 2008) and the twostream approximation (Showman et al. 2009; Rauscher & Menou 2012; DobbsDixon & Agol 2013). For the opacity treatment, grey schemes (Rauscher & Menou 2012), binning and averaging of the absorption coefficients (DobbsDixon & Agol 2013) and the correlatedk method (Showman et al. 2009) have been used. The correlatedk 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 Jupiterlike 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.
List of molecules included in our opacity database with associated line list and partition function sources.
Both the twostream approximation and the correlatedk 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), linebyline (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 Jupiterlike 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 correlatedk method and twostream approximation. This GCM is stateoftheart 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 Jupiterlike atmospheres (Mayne et al. 2013, 2014). This paper presents the adaptation of the radiation scheme, the EdwardsSlingo (ES) radiation scheme (Edwards & Slingo 1996), to conditions prevailing in hot Jupiter atmospheres. Opacities from hightemperature line lists have been calculated for the dominant absorbers in these atmospheres and the radiation scheme, using kcoefficients 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 hightemperature line lists we use and the calculation of crosssections from these line lists by using estimates for the line widths. In Sect. 3 we briefly summarise the implementation of the correlatedk method and twostream approximation in the UM and move on to testing this scheme for gradually more complicated hot Jupiterlike 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 stateoftheart 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 hightemperature 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 H_{2}O, CO, CH_{4}, NH_{3}, TiO, VO and H_{2}–H_{2} and H_{2}–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; H_{2}O, named BT2, and NH_{3}, named BYTe), while HITEMP (Rothman et al. 2010; CO) and the CIA extension to HITRAN (Richard et al. 2012; H_{2}–H_{2} and H_{2}–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, I_{a}, from Asplund et al. (2009).
Our methane line list is calculated using the Spherical Top Data System (STDS; Wenger & Champion 1998) setting the cutoff 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 CH_{4} line list should soon be released, but this will not change the main conclusions of this paper. We use the ^{12}CH_{4} partition function from Wenger et al. (2008), and ^{13}CH_{4} 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. NonLTE effects could be important in the upper part of irradiated planet atmospheres. Indeed some observational works invoke nonLTE 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 nonLTE 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 nonLTE 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, g_{u} is the degeneracy of the upper level, E_{l} is the energy of the lower level, k_{B} is Boltzmann’s constant, Q(T) is the partition function evaluated at T, A_{i} is the Einstein Acoefficient of the transition, e is the electron charge in CGSGaussian units, m_{e} is the electron mass and g_{l}f_{lu} is the gfvalue of the transition with g_{l} and f_{lu} being the degeneracy of the lower level and the oscillator strength, respectively. The quantities , g_{u}, E_{l} and A_{i}, or alternatively g_{l}f_{lu}, are given in the line lists. For line lists in the HITRAN format, line intensities S_{i}(T_{0}) evaluated at a reference temperature T_{0} = 296 K are given, and can be converted to any other temperature by using (3)which follows from Eq. (1) by taking S_{i}(T)/S_{i}(T_{0}). Both ionisation and molecular dissociation are ignored.
We calculate absorption coefficients including both Doppler broadening and pressure broadening from collisions with H_{2} 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 T_{0} and P_{0} are the reference pressure and temperature, respectively, and z is the perturbing species with P_{z} the partial pressure of species z. The total pressure broadened width is the sum of the H_{2} 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 & KedzioraChudczer (2012). In these sources, is listed as a function of quantum numbers at a given temperature and pressure. The temperature dependence exponent n_{z} also depends on the quantum numbers of the transition, but less so than . We use the same n_{z} 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).
Overview of our line width sources for pressureinduced broadening by hydrogen and helium.
In Fig. 1, we have plotted H_{2} broadened line widths as a function of J_{l}, 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 J_{l} 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 J_{l} values, which we do by keeping it constant at the value where data for the highest J_{l} 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 J_{l} are shown in Fig. 1 as solid lines.
Fig. 1 H_{2} pressure broadened line widths for H_{2}O, CH_{4}, CO and NH_{3} from the sources listed in Table 2 at P_{H2} = 10^{5} 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). 
Note that and n_{z} 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, . Nonzero pressure will also cause a shift of the line centres, but data on pressureinduced lineshifts is even more spare than that on line widths. The net effect of both Doppler and pressureinduced 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 H_{2}, and we multiply the CIA data by the density of H_{2} to make it consistent with the molecular line data.
2.2. Numerical considerations
The use of hightemperature line lists raises computational issues due to their size. In the HITRAN 2012 database (Rothman et al. 2013), the NH_{3} line list has about 4.6 × 10^{4} lines, while the ExoMol NH_{3} line list, BYTe (Yurchenko et al. 2011), has about 1.1 × 10^{9} lines. Consequently, calculating absorption coefficients from hightemperature 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 cutoff may be difficult.
A second cutoff 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 cutoff 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 cutoff with an elimination of unimportant weak lines to decrease computation time. The cutoff distance d is calculated onthefly 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 onebyone to the total mass absorption coefficient spectrum. The value of is chosen to be some fraction f_{AK} 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 cutoff method. Some lines can become unrealistically broad, however, so we include an upper limit on d of 100 cm^{1}. Note that this cutoff scheme cannot be used if the water continuum is included since it requires a cutoff 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 onthefly taking into account the line intensity at the current temperature, making a simple cutoff 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 cutoff scheme to two other schemes: (i) a cutoff at a fixed distance d_{FW} from the line centre (fixed width, FW) and (ii) a cutoff at a distance from the line centre given by the sum of the Doppler and pressure broadened widths multiplied by some factor f_{FF} (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 cutoff 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 10^{5} Pa, 1500 K and the computation time required for the three schemes as a function of the cutoff parameters. The code has been parallelised and runs using an Intel Xeon X5660 processor with 12 cores at 2.8 GHz. ^{3}
Fig. 2 Left panel: arithmetic mean of the H_{2}O absorption coefficient between 1000 cm^{1} and 1001 cm^{1} calculated using the adaptive (AK), fixed width (FW), fixed factor (FF) cutoff schemes as a function of the cutoff parameters (f_{AK}, d_{FW} and f_{FF}, respectively) at 10^{5} Pa, 1500 K. The mean absorption coefficients have been normalised by the value obtained using AK with f_{AK} = 10^{8}. Right panel: computation time required using 12 cores at 2.8 GHz as a function of the cutoff parameter. We see that our adaptive cutoff scheme is about two orders of magnitude faster than the two other methods for a given average absorption coefficient. 
The three cutoff methods reach approximately the same average absorption coefficient for the largest values of the cutoff distances (see Fig. 2). Comparing the levels in the lefthand panels to the computation times in the righthand panel, however, clearly shows the advantage of the AK method. At f_{AK} = 10^{6}, this method reaches approximately the same level as FW at d_{FW} = 10^{2} cm^{1} and FF at f_{FF} = 10^{3}. 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 f_{AK} = 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 cutoff distances at the beginning of the crosssection calculation. These line wings will normally be hidden by stronger lines, but for CO they become nonnegligible due to the lack of strong lines in certain wavelength regions. For CO we therefore use the FF method with f_{FF} = 10^{2}.
Absorption coefficients are tabulated on a log P,log T grid with 30 pressure points between 10^{1} Pa and 10^{8} Pa and 20 temperature points between 70 K and 3000 K, uniformly distributed on a logarithmic scale.
2.3. Molecular abundances
For the present work and tests we use the closedform expressions for the chemical equilibrium abundances of H_{2}O, CO, CH_{4} and NH_{3} from Burrows & Sharp (1999). We use solarlike elemental abundances, listed in Table 3. We assume the gas to be ideal and H_{2} 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.3376 g mol^{1}. For a mixture of gases, kcoefficients 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 kcoefficients for individual gases, though this is much more computationally expensive (Lacis & Oinas 1991). In Sect. 4.1, however, we only include absorption by H_{2}O and use kcoefficients calculated from the H_{2}O mass absorption coefficient and multiply the kcoefficients by the mass mixing ratio in each atmospheric layer.
Our adopted solarlike 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 10^{6} Pa to 10^{7} Pa, 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 T_{crit} or pressure above P_{crit}. For temperatures above T_{crit} and pressures below P_{crit}, 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 T_{crit} = 1500 K, P_{crit} = 10^{7} Pa.
Nonequilibrium chemistry, due to photochemical or mixing processes, might be important in hot Jupiter atmospheres (Cooper & Showman 2006). We are currently developing a nonequilibrium 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 nonequilibrium chemistry will, however, alter the main conclusions of this paper.
3. The EdwardsSlingo 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 twostream approximation (see e.g. Thomas & Stamnes 2002) and correlatedk method (Lacis & Oinas 1991; Goody et al. 1989) implemented in the EdwardsSlingo (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 linebyline 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 twostream approximation
The ES radiation scheme uses the twostream 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 twostream 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 DobbsDixon & Agol (2013); Rauscher & Menou (2012). Note that, for the thermal component, the ES radiation scheme solves the twostream 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 twostream source function method (Toon et al. 1989) for the thermal component, which is exact in the no scattering case. We show that the twostream 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 twostream 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 scattering^{6}.
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 correlatedk method
Currently the ES radiation scheme uses a combination of the exponential sum fitting of transmissions (ESFT) technique (Wiscombe & Evans 1977) and the correlatedk method (Lacis & Oinas 1991; Goody et al. 1989) to obtain kcoefficients, some details of which can be found in Sun (2011). In each band the absorption coefficients from the linebyline wavenumber grid are reordered according to strength in increasing order and divided into n_{k} subintervals. The spacing of these subintervals must be the same for each P–T and would ideally be spaced logarithmically in k. Rather than using logarithmic kintervals defined at a particular P–T 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 P–T where each part of the spectrum is most important. We use an isothermal P–T profile at 1116 K, one of the temperatures in our P–T grid, for this calculation as a compromise between dayside and nightside P–T profiles of hot Jupiters. These average absorption coefficients are then used for the initial calculation of kcoefficients. In the following, the species index i will be dropped for ease of notation.
The kcoefficient for subinterval l is found by fitting the transmission, , to the weighted transmissions of the linebyline coefficients over a set of n_{u} column densities, u_{j}, i.e. (20)where w is the weighting function and is the optimal kcoefficient 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, g_{l} is the gcoordinate corresponding to the beginning of the subinterval for kterm l and g_{nk + 1} is the gcoordinate for the end of kterm n_{k}. The error for kterm 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 blackbody 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 kterms 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 kcoefficients is repeated for each P–T. These subsequent fits use the same subinterval spacing but with a reordering appropriate for each P–T.
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 correlatedk method significantly. We note, however, that the correlatedk method will not completely converge to the LbL solution even with many kcoefficients (small ϵ_{max}): g(k) is calculated at each P–T 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 pseudomonochromatic 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 kcoefficients 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 P–T 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, n_{k} pseudomonochromatic calculations are performed, each yielding a total flux F_{l}. The bandintegrated flux is then given by (25)The weight for the thermal component should ideally be identical to w_{l} evaluated at the local temperature. For simplicity, however, the weights w_{l} 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 correlatedk 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.
Bands used.
3.3. Bandaveraged absorption coefficients
Bandaveraged absorption coefficients have been applied to exoplanet GCMs (DobbsDixon & Agol 2013) and stellar/substellar atmosphere radiation hydrodynamical models (see e.g. Freytag et al. 2010). In DobbsDixon & 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 n_{b} is the number of bands. The fluxes F_{b} are obtained by performing n_{b} pseudomonochromatic calculations, the total flux being the sum of the individual fluxes in each band. In DobbsDixon & Agol (2013), the weighting function is a blackbody 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 DobbsDixon & Agol (2013) and compare its accuracy to a full correlatedk treatment. For this purpose, we show in Sect. 4 results obtained with bandaveraged absorption coefficients designed to replicate the treatment used in DobbsDixon & Agol (2013).
3.4. Atmo
In order to investigate the accuracy of both the correlatedk method and various twostream approximations, we compare the ES radiation scheme to our linebyline 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 planeparallel 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 GaussLegendre 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 × 10^{7} wavenumber points. In contrast, we use about 300 pseudomonochromatic calculations in the ES scheme, illustrating the large gain in computational efficiency achieved by using the correlatedk method.
4. Testing the correlatedk and twostream 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 twostream approximations, we list the L^{1} 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 P_{min} (P_{max}) is the minimum (maximum) pressure in our calculations. This is a convenient measure to use when comparing errors between different twostream 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 twostream approximated and the full radiative transfer equation exist. A grey opacity is used to eliminate errors from the correlatedk method. This scenario is used to test both the accuracy of the numerical solvers and the twostream 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 correlatedk method by including only H_{2}–H_{2} 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 2kcoefficients 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 H_{2}O, 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 wellmixed hot Jupiter nightside and dayside, respectively, where the dayside includes irradiation and absorption by TiO and VO in the upper atmosphere. A gravitational acceleration of 9.42 m/s^{2} is used, suitable for HD 209458b. Unless otherwise stated, kcoefficients 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.
Fig. 3 Left panel: fluxes obtained with ES radiation scheme using the twostream approximation and correlatedk method obtained with D = 1.66 and ϵ_{max} = 5 × 10^{3} (dashdotted, green), ϵ_{max} = 10^{4} (dashed, cyan), and mean absorption coefficients (dotted, blue), for an isothermal atmosphere with pure H_{2}O 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. 
4.1. Test 1: Pure H_{2}O absorption in a hightemperature isothermal atmosphere
This test includes only absorption by H_{2}O. The temperature is fixed to 1500 K and the atmospheric domain extends from 10^{1} Pa to 10^{8} Pa, 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 1500 K 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 kcoefficients in each band, respectively. Note, however, that the number of kcoefficients 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 ~10^{4} Pa. 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 L^{1} 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 kcoefficient per band is not sufficient to resolve the opacity is reflected by the need for ~10 kcoefficients in each band to achieve a tolerance of ϵ_{max} = 5 × 10^{3} for H_{2}O. The failure of this method is discussed in more detail in Sect. 4.4.
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. 
Computed flux (F) and heating rate (ℋ) L^{1} 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 nightside 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 nightside P–T profile from Iro et al. (2005) with the smoothing described in Mayne et al. (2014), shown in Fig. 5. The temperature varies from about 400 K in the upper atmosphere to above 1600 K at 10^{8} Pa, 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 T_{lb} = T(P_{lb}), where T_{lb} and P_{lb} are the temperature and pressure at the lower boundary, respectively. From the P–T profile in Fig. 5, T_{lb} = 1662 K. All molecules in Table 1 are included except TiO and VO, with abundances as described in Sect. 2.3.
Fig. 5 P–T profile used in test 2. From the polynomial fit (Heng et al. 2011) to the nightside 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. 
Fluxes and heating rates with relative errors are plotted in Figs. 6 and 7, with L^{1} errors given in Table 6, from which it is clear that D = 1.66 yields the most accurate fluxes and heating rates overall.
Fig. 6 Left panel: fluxes obtained in test 2 with the ES radiation scheme using the P–T 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. 
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. 
Computed flux (F) and heating rate (ℋ) L^{1} 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 dayside hot Jupiter atmosphere
Our last test adopts conditions suitable to the dayside of hot Jupiters. We use the polynomial fit (Heng et al. 2011) to the dayside P–T 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 T_{lb} = T(P_{lb}), i.e. 1998 K using the P–T 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 correlatedk method.
Fig. 8 P–T profile used in test 3. From the polynomial fit (Heng et al. 2011) to the dayside 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. 
We assume an orbital distance a_{orbit} = 0.047 ua and a parent star effective temperature = 5785 K and radius R_{star} = R_{Sun}. Using StefanBoltzmann’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 Kurucz^{7}. At smaller wavelengths than available, we set the stellar flux to zero, while at larger wavelengths we extrapolate using a blackbody spectrum with the effective temperature of the Sun (T = 5785 K).
Figures 9 and 10 show the thermal flux and heating rate with relative errors. The flux error is small in regions with nonnegligible flux. The heating rate error also remains small in regions with significant cooling. This is confirmed by the computed L^{1} 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 blackbody or uniform weighting scheme is used does not affect the accuracy significantly, but using a mean absorption coefficient does result in significant errors.
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. 
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. 
Computed flux (F) and heating rate (ℋ) L^{1} norms for the thermal component in test 3.
Computed flux (F) and heating rate (ℋ) L^{1} 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 × 10^{5} W/m^{2}, 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 L^{1} errors listed in Table 8. Using ϵ_{max} = 10^{4} significantly reduces the error from the correlatedk 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.
Fig. 11 Left panel: stellar component of the flux as a function of total pressure for test 3 obtained with ϵ_{max} = 5 × 10^{3} (dashdotted, 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. 
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.
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. 
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. 
Computed flux (F) and heating rate (ℋ) L^{1} 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 bandaveraged 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 bandintegrated 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 blackbody 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 bandmean 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 nonzero temperature gradients, the upward flux will also depend on pressure. In both P–T profiles used here, the temperature decreases with height overall. At the lower boundary the upward intensity is simply the Planck intensity, i.e. the righthand 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 bandmean 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 DobbsDixon & 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 longitudelatitudedependent 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 Earthlike conditions, but detailed analysis for hot Jupiterlike conditions are lacking. In this paper we have analysed the accuracy and uncertainties in stateoftheart radiation schemes used in several GCMs applied to hot Jupiters. Opacity sources and calculation of absorption coefficients from hightemperature line lists have been discussed. We present a line profile cutoff 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 twostream approximation and correlatedk method’s applicability to hot Jupiter atmospheres have been analysed by comparing the EdwardsSlingo radiation scheme to discrete ordinate linebyline calculations.
The ES radiation scheme’s performance in these tests shows that we have successfully adapted it to hot Jupiterlike atmospheres. Our main conclusions are:

Pressure broadening parameters for hightemperaturemolecular 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 twostream approximation, although is only slightly less accurate.

About 10kcoefficients in each band for molecular line absorption yield satisfactory accuracy. Using ~100 kcoefficients 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 twostream approximation and the correlatedk method contribute nonnegligibly 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 blackbody 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 kcoefficients in both the thermal and stellar spectral regions and for different irradiation spectra.

Using a mean absorption coefficient in each band, as in DobbsDixon & 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. Bandaveraged 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.
Note that for CO the line width is approximately constant as a function of J_{l}, and we therefore use a simple mean instead of a linear fit, as shown in Fig. 1.
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/20072013 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
 Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
 Bailey, J., & KedzioraChudczer, L. 2012, MNRAS, 419, 1913 [NASA ADS] [CrossRef] [Google Scholar]
 Baraffe, I., Chabrier, G., & Barman, T. 2008, A&A, 482, 315 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Baraffe, I., Chabrier, G., & Barman, T. 2010, Rep. Prog. Phys., 73, 016901 [NASA ADS] [CrossRef] [Google Scholar]
 Barber, R. J., Tennyson, J., Harris, G. J., & Tolchenov, R. N. 2006, MNRAS, 368, 1087 [NASA ADS] [CrossRef] [Google Scholar]
 Barman, T. S. 2008, ApJ, 676, L61 [NASA ADS] [CrossRef] [Google Scholar]
 Barman, T. S., Hauschildt, P. H., Schweitzer, A., et al. 2002, ApJ, 569, L51 [NASA ADS] [CrossRef] [Google Scholar]
 BelBruno, J. J., Gelfand, J., Radigan, W., & Verges, K. 1982, J. Mol. Spectr., 94, 336 [NASA ADS] [CrossRef] [Google Scholar]
 Bending, V. L., Lewis, S. R., & Kolb, U. 2013, MNRAS, 428, 2874 [NASA ADS] [CrossRef] [Google Scholar]
 Burrows, A., & Sharp, C. M. 1999, ApJ, 512, 843 [NASA ADS] [CrossRef] [Google Scholar]
 Burrows, A., Marley, M., Hubbard, W. B., et al. 1997, ApJ, 491, 856 [NASA ADS] [CrossRef] [Google Scholar]
 Burrows, A., Hubeny, I., Budaj, J., Knutson, H. A., & Charbonneau, D. 2007, ApJ, 668, L171 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Charbonneau, D., Knutson, H. A., Barman, T., et al. 2008, ApJ, 686, 1341 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Cho, J. Y.K., Menou, K., Hansen, B. M. S., & Seager, S. 2008, ApJ, 675, 817 [NASA ADS] [CrossRef] [Google Scholar]
 Clough, S. A., Kneizys, F. X., & Davies, R. W. 1989, Atmos. Res., 23, 229 [NASA ADS] [CrossRef] [Google Scholar]
 Collins, W. D., Ramaswamy, V., Schwarzkopf, M. D., et al. 2006, J. Geophys. Res. (Atmospheres), 111, 14317 [NASA ADS] [CrossRef] [Google Scholar]
 Cooper, C. S., & Showman, A. P. 2006, ApJ, 649, 1048 [NASA ADS] [CrossRef] [Google Scholar]
 DobbsDixon, I., & Agol, E. 2013, MNRAS, 435, 3159 [NASA ADS] [CrossRef] [Google Scholar]
 DobbsDixon, I., & Lin, D. N. C. 2008, ApJ, 673, 513 [NASA ADS] [CrossRef] [Google Scholar]
 Edwards, J. M. 1996, Journal of Atmospheric Sciences, 53, 1921 [Google Scholar]
 Edwards, J. M., & Slingo, A. 1996, Quart. J. Roy. Meteorol. Soc., 122, 689 [Google Scholar]
 Elasser, W. M. 1942, Heat transfer by infrared radiation in the atmosphere No. 6 (Harvard Meteorological Studies), 107 [Google Scholar]
 Ellingson, R. G., Ellis, J., & Fels, S. 1991, J. Geophys. Res., 96, 8929 [NASA ADS] [CrossRef] [Google Scholar]
 Fortney, J. J., Saumon, D., Marley, M. S., Lodders, K., & Freedman, R. S. 2006, ApJ, 642, 495 [NASA ADS] [CrossRef] [Google Scholar]
 Freedman, R. S., Marley, M. S., & Lodders, K. 2008, ApJS, 174, 504 [NASA ADS] [CrossRef] [Google Scholar]
 Freytag, B., Allard, F., Ludwig, H.G., Homeier, D., & Steffen, M. 2010, A&A, 513, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gamache, R. R., Lynch, R., & Brown, L. R. 1996, J. Quant. Spectr. Rad. Transf., 56, 471 [NASA ADS] [CrossRef] [Google Scholar]
 Goldblatt, C., Lenton, T. M., & Watson, A. J. 2009, Quarterly Journal of the Royal Meteorological Society, 135, 619 [NASA ADS] [CrossRef] [Google Scholar]
 Goody, R., West, R., Chen, L., & Crisp, D. 1989, J. Quant. Spectr. Rad. Transf., 42, 539 [NASA ADS] [CrossRef] [Google Scholar]
 Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hadded, S., Aroui, H., Orphal, J., Bouanich, J.P., & Hartmann, J.M. 2001, J. Mol. Spectr., 210, 275 [NASA ADS] [CrossRef] [Google Scholar]
 Held, I. M., & Suarez, M. J. 1994, Bulletin of the American Meteorological Society, 75, 1825 [NASA ADS] [CrossRef] [Google Scholar]
 Heng, K., Menou, K., & Phillipps, P. J. 2011, MNRAS, 413, 2380 [NASA ADS] [CrossRef] [Google Scholar]
 Hogan, R. J. 2010, J. Atmos. Sci., 67, 2086 [NASA ADS] [CrossRef] [Google Scholar]
 Hollingsworth, J. L., & Kahre, M. A. 2010, Geophys. Res. Lett., 37, 22202 [NASA ADS] [CrossRef] [Google Scholar]
 Huitson, C. M., Sing, D. K., Pont, F., et al. 2013, MNRAS, 434, 3252 [NASA ADS] [CrossRef] [Google Scholar]
 Iro, N., Bézard, B., & Guillot, T. 2005, A&A, 436, 719 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Irwin, P. G. J., Teanby, N. A., de Kok, R., et al. 2008, J. Quant. Spectr. Rad. Trans., 109, 1136 [NASA ADS] [CrossRef] [Google Scholar]
 Knutson, H. A., Charbonneau, D., Allen, L. E., et al. 2007, Nature, 447, 183 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Knutson, H. A., Charbonneau, D., Allen, L. E., Burrows, A., & Megeath, S. T. 2008, ApJ, 673, 526 [NASA ADS] [CrossRef] [Google Scholar]
 Knutson, H. A., Charbonneau, D., Cowan, N. B., et al. 2009, ApJ, 690, 822 [NASA ADS] [CrossRef] [Google Scholar]
 Lacis, A. A., & Oinas, V. 1991, J. Geophys. Res., 96, 9027 [NASA ADS] [CrossRef] [Google Scholar]
 LeMoal, M. F., & Severin, F. 1986, J. Quant. Spectr. Rad. Transf., 35, 145 [NASA ADS] [CrossRef] [Google Scholar]
 Lebonnois, S., Lee, C., Yamamoto, M., et al. 2011, in EPSCDPS Joint Meeting 2011, 144 [Google Scholar]
 Majeau, C., Agol, E., & Cowan, N. B. 2012, ApJ, 747, L20 [NASA ADS] [CrossRef] [Google Scholar]
 Mantz, A. W., Malathy Devi, V., Chris Benner, D., et al. 2005, J. Mol. Struct., 742, 99 [NASA ADS] [CrossRef] [Google Scholar]
 Margolis, J. S. 1993, J. Quant. Spectr. Rad. Transf., 50, 431 [NASA ADS] [CrossRef] [Google Scholar]
 Mayne, N. J., Baraffe, I., Acreman, D. M., et al. 2013, Geoscientific Model Development Discussions, 6, 3681 [NASA ADS] [CrossRef] [Google Scholar]
 Mayne, N. J., Baraffe, I., Acreman, D. M., et al. 2014, A&A, 561, A1 [Google Scholar]
 Meador, W. E., & Weaver, W. R. 1980, J. Atmos. Sci., 37, 630 [NASA ADS] [CrossRef] [Google Scholar]
 Menou, K., & Rauscher, E. 2009, ApJ, 700, 887 [NASA ADS] [CrossRef] [Google Scholar]
 Mlawer, E. J., Taubman, S. J., Brown, P. D., Iacono, M. J., & Clough, S. A. 1997, J. Geophys. Res., 102, 16663 [NASA ADS] [CrossRef] [Google Scholar]
 MüllerWodarg, I. C. F., Mendillo, M., Yelle, R. V., & Aylward, A. D. 2006, Icarus, 180, 147 [NASA ADS] [CrossRef] [Google Scholar]
 Nordlund, Å. 1982, A&A, 107, 1 [Google Scholar]
 Nouri, S., Orphal, J., Aroui, H., & Hartmann, J.M. 2004, J. Mol. Spectr., 227, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Oreopoulos, L., Mlawer, E., Delamere, J., et al. 2012, J. Geophys. Res. (Atmospheres), 117, 6118 [NASA ADS] [CrossRef] [Google Scholar]
 PerezBecker, D., & Showman, A. P. 2013, ApJ, 776, 134 [NASA ADS] [CrossRef] [Google Scholar]
 Pine, A. S. 1992, J. Chem. Phys., 97, 773 [NASA ADS] [CrossRef] [Google Scholar]
 Pine, A. S., Markov, V. N., Buffa, G., & Tarrini, O. 1993, J. Quant. Spectr. Rad. Transf., 50, 337 [NASA ADS] [CrossRef] [Google Scholar]
 Plez, B. 1998, A&A, 337, 495 [NASA ADS] [Google Scholar]
 Polichtchouk, I., Cho, J. Y.K., Watkins, C., et al. 2014, Icarus, 229, 355 [NASA ADS] [CrossRef] [Google Scholar]
 Pont, F., Knutson, H., Gilliland, R. L., Moutou, C., & Charbonneau, D. 2008, MNRAS, 385, 109 [NASA ADS] [CrossRef] [Google Scholar]
 Rauscher, E., & Menou, K. 2012, ApJ, 750, 96 [NASA ADS] [CrossRef] [Google Scholar]
 Reed, K. A., & Jablonowski, C. 2011, Mon. Weather Rev., 139, 689 [NASA ADS] [CrossRef] [Google Scholar]
 RégaliaJarlot, L., Thomas, X., von der Heyden, P., & Barbe, A. 2005, J. Quant. Spectr. Rad. Transf., 91, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Richard, C., Gordon, I. E., Rothman, L. S., et al. 2012, J. Quant. Spectr. Rad. Transf., 113, 1276 [NASA ADS] [CrossRef] [Google Scholar]
 Rothman, L. S., Gordon, I. E., Barbe, A., et al. 2009, J. Quant. Spectr. Rad. Transf., 110, 533 [NASA ADS] [CrossRef] [Google Scholar]
 Rothman, L. S., Gordon, I. E., Barber, R. J., et al. 2010, J. Quant. Spectr. Rad. Transf., 111, 2139 [NASA ADS] [CrossRef] [Google Scholar]
 Rothman, L. S., Gordon, I. E., Babikov, Y., et al. 2013, J. Quant. Spectr. Rad. Transf., 130, 4 [Google Scholar]
 Sauval, A. J., & Tatum, J. B. 1984, ApJS, 56, 193 [NASA ADS] [CrossRef] [Google Scholar]
 Schweitzer, A., & Hauschildt, P. H. 2004, in AIP Conf. Ser. 730, eds. J. S. Cohen, D. P. Kilcrease, & S. Mazavet, 111 [Google Scholar]
 Sharp, C. M., & Burrows, A. 2007, ApJS, 168, 140 [NASA ADS] [CrossRef] [Google Scholar]
 Showman, A. P., & Guillot, T. 2002, A&A, 385, 166 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Showman, A. P., Fortney, J. J., Lian, Y., et al. 2009, ApJ, 699, 564 [NASA ADS] [CrossRef] [Google Scholar]
 Sing, D. K., Lecavelier des Etangs, A., Fortney, J. J., et al. 2013, MNRAS, 436, 2956 [NASA ADS] [CrossRef] [Google Scholar]
 Skartlien, R. 2000, ApJ, 536, 465 [NASA ADS] [CrossRef] [Google Scholar]
 Snellen, I. A. G., de Kok, R. J., de Mooij, E. J. W., & Albrecht, S. 2010, Nature, 465, 1049 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Solodov, A. M., & Starikov, V. I. 2009, Mol. Phys., 107, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Steyert, D. W., Wang, W. F., Sirota, J. M., Donahue, N. M., & Reuter, D. C. 2004, J. Quant. Spectr. Rad. Transf., 83, 183 [NASA ADS] [CrossRef] [Google Scholar]
 Sun, Z. 2011, Quart. J. Roy. Meteorol. Soc., 137, 2138 [Google Scholar]
 Tennyson, J., & Yurchenko, S. N. 2012, MNRAS, 425, 21 [Google Scholar]
 Thomas, G. E., & Stamnes, K. 2002, Radiative Transfer in the Atmosphere and Ocean (Cambridge University Press) [Google Scholar]
 Thrastarson, H. T., & Cho, J. Y. 2010, ApJ, 716, 144 [NASA ADS] [CrossRef] [Google Scholar]
 Tinetti, G., VidalMadjar, A., Liang, M.C., et al. 2007, Nature, 448, 169 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Toon, O. B., McKay, C. P., Ackerman, T. P., & Santhanam, K. 1989, J. Geophys. Res., 94, 16287 [NASA ADS] [CrossRef] [Google Scholar]
 Ullrich, P. A., Melvin, T., Jablonowski, C., & Stanisforth, A. 2013, Quart. J. Roy. Meteorol. Soc., submitted [Google Scholar]
 Varanasi, P., & Chudamani, S. 1990, J. Quant. Spectr. Rad. Transf., 43, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Voronin, B. A., Lavrentieva, N. N., Mishina, T. P., et al. 2010, J. Quant. Spectr. Rad. Transf., 111, 2308 [NASA ADS] [CrossRef] [Google Scholar]
 Wakeford, H. R., Sing, D. K., Deming, D., et al. 2013, MNRAS, 435, 3481 [NASA ADS] [CrossRef] [Google Scholar]
 Waldmann, I. P., Tinetti, G., Drossart, P., et al. 2012, ApJ, 744, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Watkins, C., & Cho, J. Y.K. 2010, ApJ, 714, 904 [NASA ADS] [CrossRef] [Google Scholar]
 Wenger, C., & Champion, J. P. 1998, J. Quant. Spectr. Rad. Transf., 59, 471 [Google Scholar]
 Wenger, C., Champion, J. P., & Boudon, V. 2008, J. Quant. Spectr. Rad. Transf., 109, 2697 [NASA ADS] [CrossRef] [Google Scholar]
 Wiscombe, W. J., & Evans, J. W. 1977, J. Comput. Phys., 24, 416 [NASA ADS] [CrossRef] [Google Scholar]
 Yamazaki, Y. H., Skeet, D. R., & Read, P. L. 2004, Planet. Space Sci., 52, 423 [NASA ADS] [CrossRef] [Google Scholar]
 Yurchenko, S. N., Barber, R. J., & Tennyson, J. 2011, MNRAS, 413, 1828 [NASA ADS] [CrossRef] [Google Scholar]
 Zdunkowski, W. G., & Korb, G. J. 1985, Promet, 2/3, 26 [Google Scholar]
 Zdunkowski, W. G., Welch, R. M., & Korb, G. 1980, Beiträge zur Physik der Atmosphäre, 53, 147 [Google Scholar]
 Zdunkowski, W. G., Panhans, W.G., Welch, R. M., & Korb, G. J. 1982, Beiträge zur Physik der Atmosphäre, 55, 215 [Google Scholar]
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 T_{c}. 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 twostream 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 twostream approximation
The twostream 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 StefanBoltzmann’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 u_{i}. No downward radiation at the upper boundary implies , which yields (A.15)Perfect blackbody 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 closedform solution, but can be found numerically. Integrating over all wavenumbers and using StefanBoltzmann’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 twostream approximation effectively evaluates the integral in Eq. (A.20) using a single quaderature point which can be chosen using e.g. GaussLegendre quadrature or an empirical fit. The heating rate, Eq. (17), is similarly given by
Fig. A.1 Left panel: fluxes obtained using the twostream flux in Eq. (A.10) and exact flux in Eq. (A.20) obtained with (dotted, red), D = 1.66 (dashdotted, green), D = 2 (dashed, cyan) and solving the fully angular dependent radiative transfer equation (solid, black). Right panel: calculated relative errors in the twostream fluxes. Relative errors become unreasonably large only where the flux is very small. 
Fig. A.2 Same as Fig. A.1 for heating rates. Relative errors become unreasonably large only where the heating rate is very small. 
Appendix A.2: Accuracy of the twostream 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 blackbody 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 twostream 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 L^{1} 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.
Fig. A.3 Relative error in the numerical solutions from the ES radiation scheme (dashdotted, 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. 
Computed flux (F) and heating rate (ℋ) L^{1} 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 righthand 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.
Computed flux (F) and heating rate (ℋ) L^{1} 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 L^{1} 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 L^{1} 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 twostream approximation. This confirms that both numerical solvers yield satisfactory accuracy.
All Tables
List of molecules included in our opacity database with associated line list and partition function sources.
Overview of our line width sources for pressureinduced broadening by hydrogen and helium.
Computed flux (F) and heating rate (ℋ) L^{1} norms for test 1, showing that the most accurate fluxes and heating rates are obtained with D = 1.66 and .
Computed flux (F) and heating rate (ℋ) L^{1} norms for test 2, showing that the most accurate fluxes and heating rates are obtained with D = 1.66.
Computed flux (F) and heating rate (ℋ) L^{1} norms for the thermal component in test 3.
Computed flux (F) and heating rate (ℋ) L^{1} norms for the stellar component in test 3.
Computed flux (F) and heating rate (ℋ) L^{1} norms for the total flux and heating rate in test 3.
Computed flux (F) and heating rate (ℋ) L^{1} norms for test 0 using the analytical solutions, thereby eliminating the errors from the numerical solution schemes.
Computed flux (F) and heating rate (ℋ) L^{1} norms for test 0 comparing the numerical and analytical solutions to check the accuracy of the numerical schemes.
All Figures
Fig. 1 H_{2} pressure broadened line widths for H_{2}O, CH_{4}, CO and NH_{3} from the sources listed in Table 2 at P_{H2} = 10^{5} 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). 

In the text 
Fig. 2 Left panel: arithmetic mean of the H_{2}O absorption coefficient between 1000 cm^{1} and 1001 cm^{1} calculated using the adaptive (AK), fixed width (FW), fixed factor (FF) cutoff schemes as a function of the cutoff parameters (f_{AK}, d_{FW} and f_{FF}, respectively) at 10^{5} Pa, 1500 K. The mean absorption coefficients have been normalised by the value obtained using AK with f_{AK} = 10^{8}. Right panel: computation time required using 12 cores at 2.8 GHz as a function of the cutoff parameter. We see that our adaptive cutoff scheme is about two orders of magnitude faster than the two other methods for a given average absorption coefficient. 

In the text 
Fig. 3 Left panel: fluxes obtained with ES radiation scheme using the twostream approximation and correlatedk method obtained with D = 1.66 and ϵ_{max} = 5 × 10^{3} (dashdotted, green), ϵ_{max} = 10^{4} (dashed, cyan), and mean absorption coefficients (dotted, blue), for an isothermal atmosphere with pure H_{2}O 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. 

In the text 
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. 

In the text 
Fig. 5 P–T profile used in test 2. From the polynomial fit (Heng et al. 2011) to the nightside 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. 

In the text 
Fig. 6 Left panel: fluxes obtained in test 2 with the ES radiation scheme using the P–T 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. 

In the text 
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. 

In the text 
Fig. 8 P–T profile used in test 3. From the polynomial fit (Heng et al. 2011) to the dayside 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. 

In the text 
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. 

In the text 
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. 

In the text 
Fig. 11 Left panel: stellar component of the flux as a function of total pressure for test 3 obtained with ϵ_{max} = 5 × 10^{3} (dashdotted, 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. 

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

In the text 
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. 

In the text 
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. 

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

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

In the text 
Fig. A.3 Relative error in the numerical solutions from the ES radiation scheme (dashdotted, 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. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.