Treatment of overlapping gaseous absorption with the correlatedk method in hot Jupiter and brown dwarf atmosphere models
^{1} Astrophysics Group, University of Exeter, Exeter, EX4 4QL, UK
^{2} Department of Applied Physics and Applied Mathematics, Columbia University, New York, NY 10025, USA
email: d.s.amundsen@columbia.edu
^{3} NASA Goddard Institute for Space Studies, New York, NY 10025, USA
^{4} Maison de la Simulation, CEACNRSINRIAUPSUVSQ, USR 3441, Centre d’Étude de Saclay, 91191 GifSurYvette, France
^{5} Met Office, Exeter, EX1 3PB, UK
^{6} Univ. Lyon, ENS de Lyon, Univ. Lyon 1, CNRS, CRAL, UMR 5574, 69007 Lyon, France
Received: 15 July 2016
Accepted: 4 October 2016
The correlatedk method is frequently used to speed up radiation calculations in both onedimensional and threedimensional atmosphere models. An inherent difficulty with this method is how to treat overlapping absorption, i.e. absorption by more than one gas in a given spectral region. We have evaluated the applicability of three different methods in hot Jupiter and brown dwarf atmosphere models, all of which have been previously applied within models in the literature: (i) random overlap, both with and without resorting and rebinning, (ii) equivalent extinction and (iii) premixing of opacities, where (i) and (ii) combine kcoefficients for different gases to obtain kcoefficients for a mixture of gases, while (iii) calculates kcoefficients for a given mixture from the corresponding mixed linebyline opacities. We find that the random overlap method is the most accurate and flexible of these treatments, and is fast enough to be used in onedimensional models with resorting and rebinning. In threedimensional models such as global circulation models (GCMs) it is too slow, however, and equivalent extinction can provide a speedup of at least a factor of three with only a minor loss of accuracy while at the same time retaining the flexibility gained by combining kcoefficients computed for each gas individually. Premixed opacities are significantly less flexible, and we also find that particular care must be taken when using this method in order to to adequately resolve steep variations in composition at important chemical equilibrium boundaries. We use the random overlap method with resorting and rebinning in our onedimensional atmosphere model and equivalent extinction in our GCM, which allows us to e.g. consistently treat the feedback of nonequilibrium chemistry on the total opacity and therefore the calculated P–T profiles in our models.
Key words: opacity / radiative transfer / methods: numerical / planets and satellites: atmospheres / brown dwarfs / planets and satellites: gaseous planets
© ESO, 2017
1. Introduction
Rapid calculation of wavelengthintegrated fluxes and heating rates are needed in most planetary and brown dwarf atmosphere models. As linebyline approaches are too computationally expensive for practical use, the correlatedk method (Goody et al. 1989; Lacis & Oinas 1991; Thomas & Stamnes 2002) has been applied in retrieval models (e.g. Irwin et al. 2008), onedimensional (1D) atmosphere models (e.g. Marley et al. 1996; Burrows et al. 1997; Fortney et al. 2005), and threedimensional (3D) global circulation models (GCMs; e.g. Showman et al. 2009; Kataria et al. 2013; Amundsen et al. 2016). With the correlatedk method the spectrum is divided into bands and, in each band, the opacity probability distribution is derived and described by a small number (usually 8 to 16) of kcoefficients and corresponding weights. Pseudomonochromatic calculations are performed using these kcoefficients, decreasing the required computation time by several orders of magnitude compared to linebyline calculations.
The treatment of overlapping gaseous absorption, i.e. absorption by more than one gas in a single spectral interval, with the correlatedk method is a difficult issue. Bands can be chosen such that absorption is dominated by a single gas in each band, however, this choice will be imperfect both because the relative strength of absorbers may change with temperature and pressure and due to spectral regions with significant overlap of the absorption of different gases. It is therefore necessary to take into account absorption by more than one gas in the same spectral interval. Several different schemes for deriving kcoefficients for gas mixtures have been developed for the Earth atmosphere (see e.g. Goody et al. 1989; Lacis & Oinas 1991; Fu & Liou 1992; Edwards 1996; Buchwitz et al. 2000; Yang et al. 2000; Li & Barker 2005; Shi et al. 2009; Hogan 2010; Sun 2011), each with advantages and drawbacks. The goal of this paper is not to review each of these, but to evaluate the accuracy and flexibility of three schemes that have previously been applied in hot Jupiter and brown dwarf atmosphere models:

1.
Premixedkcoefficients (PM, Goodyet al. 1989): kcoefficients for themixture are computed directly from the total linebyline gasopacity.

2.
The random overlap method, both without (RO) and with (RORR) resorting and rebinning (Lacis & Oinas 1991): kcoefficients are computed for each gas and combined assuming their absorption crosssections are uncorrelated.

3.
Equivalent extinction (EE, Edwards 1996): kcoefficients are computed for each gas and combined using an “equivalent grey absorption” for all minor absorbers and all kcoefficients for the major absorber in each band.
Premixed kcoefficients have been employed in solar system planet, exoplanet and brown dwarf atmosphere models (see e.g. Marley et al. 1996; Burrows et al. 1997; Fortney et al. 2005; Showman et al. 2009; Wordsworth et al. 2013). This method avoids problems related to combining kcoefficients for different gases, but is inflexible as mixing must be assumed before kcoefficients are computed. Alternatively, gas mixing ratios can be added as dimensions to the lookup table of kcoefficients, however, this leads to a very large number of dimensions in the table.
The random overlap method has been applied in retrieval models (Irwin et al. 2008), and 1D brown dwarf and hot Jupiter atmosphere models (Drummond et al. 2016; Mollière et al. 2015; Tremblin et al. 2015, 2016), and assumes that the absorption crosssections of different gases are uncorrelated. The total number of kcoefficients in a band scales as the product of the number of kcoefficients for each overlapping gas, causing this method to become computationally expensive, but resorting and rebinning the resulting kcoefficients can be used to circumvent this issue (Lacis & Oinas 1991). We have recently applied equivalent extinction in our GCM to study hot Jupiters (Amundsen et al. 2016). Like the random overlap method this method is more flexible than using premixed kcoefficients, but requires knowledge of which absorbers should be treated as the major and minor sources of opacity in each band.
It is clearly beneficial in terms of model flexibility to compute kcoefficients individually for each gas and combine them onthefly in models using the current local mixing ratios. As all wavelength information is lost when the kcoefficients are computed it is impossible to do this perfectly without loss of accuracy, and requires an assumption about the absorption of the different gases. The random overlap method assumes that the lines of different gases are randomly overlapping (or equivalently that the absorption crosssections are uncorrelated), while equivalent extinction assumes minor absorbers can be treated as grey. It is essential to verify the accuracy of these assumptions by comparing to linebyline calculations.
In this paper we compare premixing, random overlap and equivalent extinction in terms of computational efficiency and evaluate their accuracy by comparing to results from linebyline calculations. In Sect. 2 we give a brief overview of the correlatedk method and Sect. 3 describes the above overlap schemes in more detail. In Sect. 4 we apply them in hot Jupiter atmosphere models, compare them and evaluate their computational efficiency, by using our 1D radiativeconvective equilibrium atmosphere code ATMO (Tremblin et al. 2015, 2016) and our GCM radiation scheme SOCRATES^{1} (Edwards & Slingo 1996; Edwards 1996; Amundsen et al. 2014). We give our concluding remarks in Sect. 5.
2. The correlatedk method
As treating the wavelengthdependence of gaseous absorption explicitly is too computationally expensive to be performed in many atmosphere models, the correlatedk method is frequently used. It considers the probability distribution of the opacity in the spectral bands and assumes that the mapping between spectral regions and the probability distribution is vertically correlated. Originally developed for the Earth atmosphere (Lacis & Oinas 1991), it has since been adopted in both onedimensional (Marley et al. 1996; Burrows et al. 1997; Marley & Robinson 2015; Tremblin et al. 2015) and global circulation models (Showman et al. 2009; Kataria et al. 2013; Amundsen et al. 2016) of hot Jupiter and brown dwarf atmospheres. We do not discuss the correlatedk method in detail here, but refer to e.g. Lacis & Oinas (1991), Goody et al. (1989) and Thomas & Stamnes (2002) for indepth discussions. Note that we have previously verified the applicability of the correlatedk method in hot Jupiter and brown dwarf atmosphere models (Amundsen et al. 2014).
In the correlatedk method the opacity spectrum is divided into bands b. In each band kcoefficients and corresponding weights are computed from the probability distribution of the opacity, with where is the number of k coefficients within band b. The transmission through a homogeneous slab is given by where is the wavenumber, and are wavenumber limits of band b, is a weighting function, and and u are the opacity and column density of the gas, respectively. g(k) is the cumulative opacity probability distribution, where g(k) is the probability of having an opacity ≤k within the band.
Pseudomonochromatic fluxes are computed for each coefficient, with the integrated flux in band b given by (3)and the total spectral integrated flux given by (4)where n_{b} is the number of bands.
The coefficients are the kcoefficients for the gas mixture, i.e. taking into account all absorbers present. Spectral bands can be chosen such that absorption is dominated by only one gas, the major absorber, in each band. Other gases may still contribute significantly to absorption, however, which causes the need to treat overlapping absorption. In addition, in some spectral regions the major and minor absorbers may change depending on the gas mixing ratios. Consequently, there is a need to compute kcoefficients for a gas mixture.
3. Treatments of gaseous overlap
In this section we briefly discuss three different methods for treating overlapping gaseous absorption previously used in hot Jupiter and brown dwarf atmosphere models in the literature. For convenience we adopt a set of acronyms for these overlap methods, which we summarise in Table 1.
Methods for treating overlapping gaseous absorption considered here and our adopted acronyms.
3.1. Premixed
The total absorption coefficient can be calculated by summing linebyline absorption coefficients for all absorbing species weighted by their relative abundances: (5)where the sum is over all N_{s} species, and and ζ_{i}(P,T) are the absorption coefficient and mixing ratio of gas i at pressure P and temperature T, respectively. The total absorption coefficient at a given (P,T) is then given by k^{tot}ρ, where ρ is the total gas density. k^{tot} can be used to compute and tabulate kcoefficients for the gas mixture as a function of temperature and pressure. This approach has several advantages: it is fast, requiring only one set of kcoefficients for each temperature and pressure, and it is simple to implement. This technique has been used in 1D atmosphere models (e.g. Marley & Robinson 2015) and the SPARC/MITgcm (Showman et al. 2009). It is not particularly flexible, however, as the local mixing ratios ζ_{i}(P,T) must be determined before the time consuming calculation of kcoefficients. A potential solution would be to add gas mixing ratios as dimensions to the lookup table of kcoefficients, but the increased size of such a table is prohibitive for application in atmosphere models with many absorbing gases.
3.2. The random overlap method
The second method we discuss is the random overlap method (Lacis & Oinas 1991). kcoefficients are computed individually for each gas, and kcoefficients for different gases are combined assuming that the absorption coefficient of one gas x, is uncorrelated to that of a second gas y, i.e. that their lines are randomly overlapping. The total transmission of the gas mixture over some column density (u_{x},u_{y}) is then given by a simple scalar product, (6)We show in Appendix A that this is equivalent to convolving the opacity probability distributions of the different gases. The assumption that the absorption coefficients are uncorrelated will depend on the adopted bands and its applicability should be verified by comparing to linebyline calculations. We perform such a comparison in Sect. 4.1.
3.2.1. Without resorting and rebinning
Equation (6) can be rewritten in terms of the kcoefficients for the individual gases x and y. The transmission through one layer is, using Eqs. (1), (2) and (6), Defining u_{xy} = u_{x} + u_{y}, we can write the above transmission as (9)where and (12)For illustration we show schematics of the kcoefficients H_{2}O and CO in Fig. 1a, with the combined kcoefficients in Fig. 1b calculated using Eq. (11) assuming ζ_{H2O} = 0.9 and ζ_{CO} = 0.1.
Running n_{k,x}n_{k,y} pseudomonochromatic calculations using the k_{xy,lm}coefficients the total flux can be calculated as usual using Eqs. (3) and (4) with the weights w_{xy,lm}. This procedure can be replicated for an arbitrary number of gases, however, the computation time increases by a factor of n_{k} for each gas added. This method therefore quickly becomes too computationally expensive for practical use.
Fig. 1 Schematic illustration of the random overlap method. The original kdistributions of H_{2}O and CO are shown in panel a), with the combined kcoefficients shown in panel b) calculated using Eqs. (11) and (12) assuming ζ_{H2O} = 0.9 and ζ_{CO} = 0.1. Without resorting and rebinning (RO) the k_{xy,lm}coefficients with corresponding weights w_{xy,lm} are used directly. With resorting and rebinning (RORR) the combined kcoefficients in panel b) are resorted, as shown in panel c), and then rebinned down to a smaller number of kcoefficients, as shown in panel d). 

Open with DEXTER 
3.2.2. With resorting and rebinning
Lacis & Oinas (1991) suggested that ranking and reblocking, i.e. resorting the k_{xy,lm}coefficients and rebinning them to obtain a smaller number of kcoefficients , would circumvent the scaling issue. The procedure is illustrated in panels c and d of Fig. 1. First the kcoefficients of two gases are combined using Eqs. (11) and (12), as shown in Fig. 1b. These n_{k,x}n_{k,y}kcoefficients are sorted in increasing order, with the weights sorted using the same mapping, as shown in Fig. 1c^{2}. The sorted k_{xy,lm}coefficients are then binned down to coefficients, which we show in Fig. 1d.
The weights used in the rebinning, , must be the same for all layers, and should ideally be spaced equally in log k_{xy,lm}. As k_{xy,lm} will vary vertically, however, equal spacing in log k_{xy,lm} in one layer will not correspond to equal spacing in log k_{xy,lm} at a different level. Consequently a compromise must be made, with one possible solution being an equal spacing in log k_{xy,lm} defined where the optical depth in each band reaches unity. As the above procedure of resorting and rebinning is repeated to include more than two gases, however, log k_{xy,lm} will change meaning the ideal spacing to change. In addition, particular care must be taken to treat cases where kcoefficients of gases are zero.
For simplicity and ease of implementation we use a much simpler approach for determining the weights : in SOCRATES we use weights given by a GaussLegendre quadrature, while in ATMO we use uniform weights, both of which can have an arbitrary number of rebinned kcoefficients . The rebinned coefficients are found by computing a weighted average of all k_{xy,lm}terms belonging to each bin , where w_{xy,lm} are used as weights. If a k_{xy,lm}term extends over more than a single bin, it is split over neighbouring bins such that the weights w_{xy,lm} sum up to exactly in each bin. We use a linear average, but note that the accuracy can be improved somewhat by averaging in log k_{xy,lm}. This causes a significant increase in computation time, however, due to the need to compute the logarithm of many kterms, and would also require particular care to treat cases where k_{xy,lm}terms are zero.
After this resorting and rebinning, the process is repeated, adding one gas at a time, until all gases have been added. The final binned kcoefficients are used to compute the fluxes and heating rates for the atmosphere. This approach is consequently much more flexible than premixing gases as gas abundances can be set at runtime.
3.3. Equivalent extinction
Fig. 2 Schematic illustration of equivalent extinction (EE). An equivalent grey absorption is calculated for each minor gas (CO in our case), as shown in panel a), using Eq. (13) or 14. This grey absorption is added onto the kcoefficients of the major gas (H_{2}O in this example) using the corresponding mixing ratios, here assumed to be ζ_{H2O} = 0.9, ζ_{CO} = 0.1, as shown in panel b). 

Open with DEXTER 
The last method of treating gaseous overlap that we consider is equivalent extinction (Edwards 1996). It utilizes the fact that in most bands there is a primary (major) absorber, and includes additional absorbers through a grey “equivalent extinction”. In each layer and band an equivalent extinction is calculated for each minor gas, which for the thermal component is defined as (13)where k_{x,l} are the kcoefficients of the minor gas in the layer with corresponding weights w_{x,l}, and F_{v,l} is the thermal flux in the layer including only absorption by kterm l of the gas. Pseudomonochromatic calculations are performed for all n_{k}kcoefficients of the major gas in each band, with all other absorbers included by using the equivalent grey absorption . This effectively reduces the number of pseudomonochromatic calculations required to one per kcoefficient per gas.
The direct component of the stellar flux is readily included by calculating the transmission for each gas separately and then taking the product since, assuming random overlap, direct transmissions are multiplicative (see Eq. (6)). For the diffuse stellar beam, which will be nonzero if scattering is included, the equivalent extinction is defined by (14)F_{s ∗ ,l} is the direct flux at the lower boundary including only kterm l of the gas. The use of F_{s ∗ ,l} means that equivalent extinction in the current formulation is less suited for use in hot Jupiter atmosphere models as the direct stellar flux at the bottom boundary may be zero. In this case we use the smallest kcoefficient for the minor gas as . In this work, however, as we only consider Rayleigh scattering, the main stellar radiation is contained in the direct beam, making this a minor issue.
We show a schematic illustration of this overlap method in Fig. 2, where the original kdistributions are shown in Fig. 1a. From the kcoefficients of the minor absorber, in this case CO, an equivalent grey absorption is calculated according to Eq. (13) or 14, as illustrated in Fig. 2a. This grey absorption is then added onto the kcoefficients of the minor gas using the corresponding mixing ratios, here assumed to be ζ_{H2O} = 0.9 and ζ_{CO} = 0.1, with the combined kcoefficients used in the radiation calculations shown in Fig. 2b.
3.3.1. Determining the major absorber
We consider two approaches for determining the major absorber in each band:

1.
Calculating the transmission of each gas down to the lowerboundary using the maximum allowed mixing ratio for each gasgiven the elemental abundances, e.g. for H_{2}O assuming all O is in the form of H_{2}O. The gas with the smallest transmission in each band can be considered to be the major absorber. We label this approach EE, however, it is rather naive as local mixing ratios may significantly impact which gas dominates absorption in a band. Consequently, we have also considered a more sophisticated approach.

2.
In a given column the transmission of each gas is calculated from the top of the atmosphere down to the first layer where the total transmission is <e^{1}, i.e. where the total optical depth has reached 1 in the band. The major absorber is then defined as the gas with the smallest transmission at this level. We calculate the total transmission assuming random overlap and multiply individual gas transmissions according to Eq. (6). If the total optical depth does not reach 1 we use the gas transmissions down to the lower boundary of the model instead. We label this approach AEE for adaptive equivalent extinction.
To increase the efficiency and ease the implementation we use the average kcoefficients down to an optical depth of 1 for each gas as defined in Eq. (19) in Amundsen et al. (2014) in these calculations. Alternatively, the transmissions with AEE can be calculated using the local kcoefficients in each layer, however, we do not expect this to have a significant effect on the choice of major absorber.
4. Application to hot Jupiter and brown dwarf atmospheres
In this section we evaluate the accuracy of the overlapping gaseous absorption treatments discussed above when applied to hot Jupiter and brown dwarflike atmospheres by comparing them to linebyline calculations. We have previously tested the applicability of the correlatedk method to these atmospheres (Amundsen et al. 2014), so we limit the discussion here to the overlap treatments. The test atmospheres here are identical to those in Amundsen et al. (2014), however, in the present work we include all the updates to the opacities described in Amundsen et al. (2016) and include 13 opacity sources in total^{3}. Abundances are as in Amundsen et al. (2016), and we include TiO and VO opacity for the day side P–T profile, but not for the night side profile, in Fig. 3. We take the TiO/VO condensation curve from Fortney et al. (2008) and apply a small smoothing of the abundance as described in Amundsen et al. (2016), with a smoothing scale of 10 K.
An important aspect of the correlatedk method is the choice of spectral bands. Both accuracy and computational costs increase with increasing number of bands, and it is therefore necessary to make a compromise between accuracy and speed. We adopt the 32 spectral bands used in Amundsen et al. (2014) for all calculations using the correlatedk method presented here, as this is enough for the error to become acceptably small, and small enough to facilitate use in both 1D and 3D models. A study of how the error varies depending on the number and placement of the spectral bands is beyond the scope of the present work.
The adopted P–T profiles are shown in Fig. 3, which are derived from the equilibrium P–T profile of Iro et al. (2005) derived for HD 209458b as adopted by Cooper & Showman (2005, 2006), Rauscher & Menou (2010) and Heng et al. (2011) with the minor adjustment introduced by Mayne et al. (2014). The TiO/VO condensation curve is plotted as dashed black line.
4.1. Validity of the random overlap assumption
In this section we test the validity of the random overlap assumption, which is considered to be more accurate than equivalent extinction (Edwards 1996), by comparing it to linebyline (LbL) calculations using our 1D atmosphere code ATMO (Amundsen et al. 2014; Tremblin et al. 2015, 2016). Unfortunately the random overlap method without resorting and rebinning has not yet been implemented in ATMO. Instead we use 120 rebinned kterms in the tests presented in this section, which we have found to be more than sufficient for the solution to have converged (see Sect. 4.2 for convergence tests), in order to minimize errors caused by the rebinning. All linebyline calculations were run at a resolution of 0.001 cm^{1}.
4.1.1. Night side
In Figs. 4 and 5 we show the fluxes and heating rates, with corresponding errors calculated by comparing to linebyline fluxes and heating rates, obtained for the night side P–T profile shown in Fig. 3. It is clear that both fluxes and heating rates obtained when using the correlatedk method with the random overlap method match the linebyline results very well, with errors of a few percent. We note that these errors are both due to the use of the correlatedk method and the random overlap assumption, and in agreement with the errors found in Amundsen et al. (2014).
Fig. 3 P–T profiles used for the tests in the present work. Profiles are derived from the equilibrium P–T profile of Iro et al. (2005) as adopted by Cooper & Showman (2005, 2006), Rauscher & Menou (2010) and Heng et al. (2011) with the smoothing introduced by Mayne et al. (2014). We adopt the TiO/VO condensation curve from Fortney et al. (2008), plotted as a dotted line. 

Open with DEXTER 
Fig. 4 Fluxes (left) and absolute errors in fluxes (right) for the night side P–T profile in Fig. 3. The linebyline (LbL) calculation was run at a 0.001 cm^{1} resolution corresponding to 5 × 107 monochromatic calculations, while we used 120 rebinned kcoefficients in each band corresponding to 3840 pseudomonochromatic calculations. Both results were obtained using ATMO. 

Open with DEXTER 
Fig. 5 Same as Fig. 4 but for heating rates. 

Open with DEXTER 
4.1.2. Day side
In Figs. 6 and 7 we show the fluxes and heating rates, again with corresponding errors, for the day side P−T profile shown in Fig. 3. As for the night side P–T profile both fluxes and heating rates obtained with the correlatedk method using random overlap are in good agreement with the corresponding linebyline results, with errors of a few percent. As for the night side these errors are both due to the use of the correlatedk method and the random overlap assumption, and in agreement with the errors found in Amundsen et al. (2014).
Fig. 6 Fluxes (left) and absolute errors in fluxes (right) for the day side P–T profile in Fig. 3. 

Open with DEXTER 
Fig. 7 Same as Fig. 6 but for heating rates. 

Open with DEXTER 
Based on these results we conclude that the random overlap method is indeed sufficiently accurate to be applied to these atmospheres for the bands adopted here. We note that for a different choice of bands, particularly for wider bands, the random overlap assumption should be reevaluated to make sure it is still valid.
4.2. Comparison of gaseous overlap treatments
Unfortunately, using 120 rebinned kterms becomes too computationally expensive even in our 1D atmosphere model. In this section we evaluate the accuracy of more efficient overlap treatments, and compare them in terms of both accuracy and computational efficiency using our GCM radiation scheme SOCRATES (Edwards & Slingo 1996; Edwards 1996). We have previously presented the adaptation of this radiation scheme to hot Jupiters (Amundsen et al. 2014). The radiation scheme setup used here is identical to the one we use in our hot Jupiter GCM simulations presented in Amundsen et al. (2016), i.e. for the twostream approximation we use a diffusivity of D = 1.66 for thermal and D = 2 for stellar radiation. Rayleigh scattering is included for stellar radiation, and we ignore all other forms of scattering.
4.2.1. Night side
We show in Appendix A the thermal net upward fluxes and heating rates using the night side P–T profile in Fig. 3, with corresponding errors, for all overlap treatments considered here. Errors are calculated by comparing to results obtained using the random overlap method without resorting and rebinning (RO). First, it is clear that using the random overlap method with resorting and rebinning (RORR) with an increasing number of rebinned kterms significantly decreases the error in both fluxes and heating rates. Equivalent extinction (EE) is somewhat less accurate than RORR with , and with adaptive equivalent extinction (AEE) errors do not decrease significantly indicating the choice of major absorbers with EE was appropriate for this P–T profile.
Fig. 8 Fluxes (left) and absolute errors in fluxes (right) obtained with the night side P–T profile in Fig. 3 using SOCRATES. Fluxes obtained using the random overlap method without resorting and rebinning (RO) are used to calculate errors for the random overlap with resorting and rebinning (RORR) with , 16 and 32, equivalent extinction (EE), adaptive equivalent extinction (AEE) and premixed opacities (PM). 

Open with DEXTER 
Fig. 9 Same as Fig. 8 but for heating rates. L1 norms of the errors (see Amundsen et al. 2014, the average heating rate error weighted by the local heating rates) are 4.5 % for RORR with , 1.9 % for RORR with , 1.5 % for RORR with , 13 % for EE, 11 % for AEE and 38 % for PM. 

Open with DEXTER 
Premixed (PM) opacities are significantly less accurate than all the above overlap treatments, this stems from errors introduced by the interpolation in the premixed opacity table. Changes in mixing ratios with temperature and pressure can cause large changes in the premixed opacities which are not properly resolved by our opacity table. To illustrate this we have in Fig. 10 plotted fluxes and heating rates obtained again using the night side P–T profile, but using constant mixing ratios equal to the mixing ratios at P =104 Pa, T =1000 K both when computing the premixed opacity table and when combining kcoefficients using RO. This eliminates errors caused by the implicit interpolation in mixing ratio with PM. The very small differences remaining between RO and PM are mainly due to small differences in the precision of the kcoefficients, which for RO are derived for each gas separately while for PM for the mixture directly. As in Amundsen et al. (2014) we use an opacity table logarithmically spaced in temperature and pressure, with 20 temperature points between 70 K and 3000 K and 30 pressure points between 10^{1} Pa and 108 Pa, with the opacity interpolation performed linearly in temperature. This is similar to the resolution used in previous works (Showman et al. 2009).
Fig. 10 Fluxes (left) and heating rates (right) obtained with the night side P–T profile in Fig. 3 using constant mixing ratios equal to the mixing ratios at P =104 Pa, T =1000 K. This eliminates errors caused by the implicit interpolation of mixing ratios with PM which dominates the errors seen using this overlap method in Figs. 8 and 9. 

Open with DEXTER 
In Table 2 we give the relative computation times of the overlap treatments in Figs. 8 and 9. RO is, as expected, two to three orders of magnitude slower than the other overlap treatments. The quickest is PM, although (A)EE is only slightly slower. RORR, even with only 8 rebinned kterms is about a factor of 3 slower than (A)EE. We find that a significant fraction of the computation time with RORR is spent sorting the kcoefficients, and it is therefore important to use an efficient sorting algorithm. As mentioned in Sect. 3.2.2 we use a standard quicksort implementation, which we have found to consistently give good performance compared to shellsort and heapsort.
4.2.2. Day side
We show in Figs. 11 and 12 total (thermal plus stellar) net upward fluxes and heating rates obtained using the day side P–T profile in Fig. 3, with corresponding errors, for all overlap treatments considered here. Errors are, as for the night side, calculated by comparing to results obtained using RO. Results are overall similar to those obtained above for the night side, with errors being smallest for a large number of rebinned kterms with RORR. A significant improvement in the accuracy is seen when using AEE compared to EE, indicating that the appropriate major absorbers have changed compared to the night side profile.
Fig. 11 Fluxes (left) and absolute errors in fluxes (right) obtained with the day side P–T profile in Fig. 3 using SOCRATES. Fluxes obtained using the random overlap method without resorting and rebinning (RO) are used to calculate errors for the random overlap with resorting and rebinning (RORR) with , 16 and 32, equivalent extinction (EE), adaptive equivalent extinction (AEE) and premixed opacities (PM). 

Open with DEXTER 
Fig. 12 Same as Fig. 11 but for heating rates. L1 norms of the errors (see the caption of Fig. 9) are 7.6 % for RORR with , 3.0 % for RORR with , 1.8 % for RORR with , 7.0 % for EE, 2.2 % for AEE and 119 % for PM. 

Open with DEXTER 
Perhaps the most striking result is the large errors caused by using premixed opacities, which are significantly larger for the day side compared to the night side. The flux changes very rapidly between 103 Pa and 104 Pa, which causes a large increase in the heating rate. Looking at Fig. 3 this discontinuity occurs as the P−Tprofile crosses the condensation curve of TiO and VO. Both molecules are strong absorbers in the visible, and the presence of these molecules leads to a strong absorption of the incoming stellar radiation. The steep vertical gradient in the mixing ratios of TiO and VO when the temperature is near the condensation temperature causes a similarly steep gradient in the opacity. When using PM this transition is smoothed out as the resolving power is limited by the number of P–T points in the lookup kcoefficient table, thus reducing the accuracy of the interpolation. To illustrate this we have in Fig. 13 plotted fluxes and heating rates obtained again using the day side P–T profile, but with constant mixing ratios equal to the mixing ratios at P =104 Pa, T =1900 K, similar to Fig. 10 for the night side, thereby eliminating errors caused by the implicit interpolation in mixing ratio with PM.
Fig. 13 Fluxes (left) and heating rates (right) obtained with the day side P–T profile in Fig. 3 using constant mixing ratios equal to the mixing ratios at P =104 Pa, T =1900 K. This eliminates errors caused by the implicit interpolation of mixing ratios with PM which dominates the errors seen using this overlap method in Figs. 11 and 12. 

Open with DEXTER 
In Tables 3 and 4 we give the computation times for the different overlap treatments. Results are similar to those obtained for the night side in Table 2 with PM being the fastest, (A)EE slightly slower, and RORR with only 8 rebinned kterms about 3 times slower than (A)EE.
5. Conclusions
We have evaluated the applicability of several gaseous overlap treatments in hot Jupiter atmosphere models. We have shown that the random overlap method gives good accuracy and flexibility, but without resorting and rebinning (RO) it is too slow for practical use. With resorting and rebinning (RORR) it becomes much more computationally efficient, and the accuracy and speed can be adjusted by varying the number of rebinned kterms.
Equivalent extinction (EE) is about three times faster than RORR with only 8 rebinned kterms, and benefits from the same flexibility as RORR, however, it is clear that particular care must be taken when choosing the major absorber in each band. We present one way of determining the major absorber, which we term adaptive equivalent extinction (AEE), that benefits from better accuracy compared to a more naive choice without a major loss of computational efficiency.
The fastest overlap treatment considered here is premixed opacities, however, it lacks the flexibility of the random overlap method and equivalent extinction. In addition, particular care must be taken if there are regions of the atmosphere where the total opacity changes rapidly as a function of height to ensure steep variations in composition at important chemical equilibrium boundaries are adequately resolved. If these are not adequately resolved, as is the case in our premixed table, interpolation can cause such transitions to be significantly smoothed out. We have shown that TiO in particular can be a significant cause of inaccuracies, but other species such as H_{2}O, CH_{4} and CO may also lead to inaccuracies as seen in Sect. 4.2.1. A lower resolution can be tolerated in kcoefficient tables of individual gases than in premixed tables as individual opacities vary more smoothly.
In our 1D atmosphere code ATMO we use RORR, usually with about 30 rebinned kterms in each band. This gives us the flexibility to manipulate gas mixing ratios without recomputing or having to add additional dimensions to the ktable. It also uniquely allows us to treat nonequilibrium chemistry consistently, i.e. have nonequilibrium abundances feed back on the total opacity and consequently the calculated P–T profiles (Tremblin et al. 2015, 2016; Drummond et al. 2016).
Unfortunately, RORR, even with only 8 rebinned kcoefficients, is too slow for use in our GCM coupling dynamics and radiative transfer, and we consequently use equivalent extinction. The method for adaptively choosing the major absorber in each band (AEE) has not yet been implemented in the radiation scheme coupled to our GCM. In Amundsen et al. (2016) we therefore use the simple approach of determining the major absorber in each band instead (EE). Work is in progress to improve the treatment of overlapping gaseous absorption in our GCM, one possibility being using EE for bands with a clear major absorber and RORR for bands with more than one significant absorber. In addition, as the current definition of the equivalent extinction for the stellar component relies on direct stellar fluxes at the bottom boundary it will become important to improve this definition of the equivalent extinction when the diffuse stellar flux becomes more important, e.g. when including stronger shortwave scattering.
We have used quicksort, shellsort and heapsort, all available as standard library routines (e.g. Press et al. 2007), and found that quicksort is generally the fastest. We adopt quicksort in the current work.
We include the opacity of H_{2}O (Barber et al. 2006), CO (Rothman et al. 2010), CH_{4} (Yurchenko & Tennyson 2014), NH_{3} (Yurchenko et al. 2011), H_{2}–H_{2} and H_{2}–He collision induced absorption (Richard et al. 2012), TiO (Plez 1998), VO (B. Plez, priv. comm.), Li, Na, K, Rb and Cs (Heiter et al. 2008).
Acknowledgments
We thank the referee, Mark Marley, for comments that significantly improved the paper. This work is partly supported by the European Research Council under the European Community’s Seventh Framework Programme (FP7/20072013 Grant Agreement No. 247060PEPS and grant No. 320478TOFU). D.S.A. acknowledges support from the NASA Astrobiology Program through the Nexus for Exoplanet System Science. J.M. acknowledges the support of a Met Office Academic Partnership secondment. The calculations for this paper were performed on the DiRAC Complexity machine, jointly funded by STFC and the Large Facilities Capital Fund of BIS, and the University of Exeter Supercomputer, a DiRAC Facility jointly funded by STFC, the Large Facilities Capital Fund of BIS and the University of Exeter.
References
 Amundsen, D. S., Baraffe, I., Tremblin, P., et al. 2014, A&A, 564, A59 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Amundsen, D. S., Mayne, N. J., Baraffe, I., et al. 2016, A&A, 595, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Barber, R. J., Tennyson, J., Harris, G. J., & Tolchenov, R. N. 2006, MNRAS, 368, 1087 [NASA ADS] [CrossRef] [Google Scholar]
 Buchwitz, M., Rozanov, V. V., & Burrows, J. P. 2000, J. Geophys. Res., 105, 15247 [NASA ADS] [CrossRef] [Google Scholar]
 Burrows, A., Marley, M., Hubbard, W. B., et al. 1997, ApJ, 491, 856 [NASA ADS] [CrossRef] [Google Scholar]
 Cooper, C. S., & Showman, A. P. 2005, ApJ, 629, L45 [NASA ADS] [CrossRef] [Google Scholar]
 Cooper, C. S., & Showman, A. P. 2006, ApJ, 649, 1048 [NASA ADS] [CrossRef] [Google Scholar]
 Drummond, B., Tremblin, P., Baraffe, I., et al. 2016, A&A, 594, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Edwards, J. M. 1996, J. Atmos. Sci., 53, 1921 [NASA ADS] [CrossRef] [Google Scholar]
 Edwards, J. M., & Slingo, A. 1996, QJMRS, 122, 689 [NASA ADS] [CrossRef] [Google Scholar]
 Fortney, J. J., Marley, M. S., Lodders, K., Saumon, D., & Freedman, R. 2005, ApJ, 627, L69 [NASA ADS] [CrossRef] [Google Scholar]
 Fortney, J. J., Lodders, K., Marley, M. S., & Freedman, R. S. 2008, ApJ, 678, 1419 [NASA ADS] [CrossRef] [Google Scholar]
 Fu, Q., & Liou, K. N. 1992, J. Atmos. Sci., 49, 2139 [NASA ADS] [CrossRef] [Google Scholar]
 Goody, R., West, R., Chen, L., & Crisp, D. 1989, J. Quant. Spec. Radiat. Transf., 42, 539 [NASA ADS] [CrossRef] [Google Scholar]
 Heiter, U., Barklem, P., Fossati, L., et al. 2008, J. Phys. Conf. Ser., 130, 012011 [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]
 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. Transf., 109, 1136 [NASA ADS] [CrossRef] [Google Scholar]
 Kataria, T., Showman, A. P., Lewis, N. K., et al. 2013, ApJ, 767, 76 [NASA ADS] [CrossRef] [Google Scholar]
 Lacis, A. A., & Oinas, V. 1991, J. Geophys. Res., 96, 9027 [NASA ADS] [CrossRef] [Google Scholar]
 Li, J., & Barker, H. W. 2005, J. Atmos. Sci., 62, 286 [NASA ADS] [CrossRef] [Google Scholar]
 Marley, M. S., & Robinson, T. D. 2015, ARA&A, 53, 279 [NASA ADS] [CrossRef] [Google Scholar]
 Marley, M. S., Saumon, D., Guillot, T., et al. 1996, Science, 272, 1919 [NASA ADS] [CrossRef] [Google Scholar]
 Mayne, N. J., Baraffe, I., Acreman, D. M., et al. 2014, A&A, 561, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mollière, P., van Boekel, R., Dullemond, C., Henning, T., & Mordasini, C. 2015, ApJ, 813, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Plez, B. 1998, A&A, 337, 495 [NASA ADS] [Google Scholar]
 Press, W., Teukolsky, S., Vetterling, W., & Flannery, B. 2007, Numerical Recipes: The Art of Scientific Computing, 3rd edn. (Cambridge University Press) [Google Scholar]
 Rauscher, E., & Menou, K. 2010, ApJ, 714, 1334 [NASA ADS] [CrossRef] [Google Scholar]
 Richard, C., Gordon, I. E., Rothman, L. S., et al. 2012,J. Quant. Spectr. Rad. Transf., , 113, 1276 [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]
 Shi, G., Xu, N., Wang, B., Dai, T., & Zhao, J. 2009, J. Quant. Spectr. Rad. Transf., 110, 435 [NASA ADS] [CrossRef] [Google Scholar]
 Showman, A. P., Fortney, J. J., Lian, Y., et al. 2009, ApJ, 699, 564 [NASA ADS] [CrossRef] [Google Scholar]
 Sun, Z. 2011, QJMRS, 137, 2138 [NASA ADS] [CrossRef] [Google Scholar]
 Thomas, G. E., & Stamnes, K. 2002, Radiative Transfer in the Atmosphere and Ocean (Cambridge University Press) [Google Scholar]
 Tremblin, P., Amundsen, D. S., Mourier, P., et al. 2015, ApJ, 804, L17 [NASA ADS] [CrossRef] [Google Scholar]
 Tremblin, P., Amundsen, D. S., Chabrier, G., et al. 2016, ApJ, 817, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Wordsworth, R., Forget, F., Millour, E., et al. 2013, Icarus, 222, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Yang, S., Ricchiazzi, P., & Gautier, C. 2000, J. Quant. Spec. Radiat. Transf., 64, 585 [NASA ADS] [CrossRef] [Google Scholar]
 Yurchenko, S. N., & Tennyson, J. 2014, MNRAS, 440, 1649 [NASA ADS] [CrossRef] [Google Scholar]
 Yurchenko, S. N., Barber, R. J., & Tennyson, J. 2011, MNRAS, 413, 1828 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Relationship between convolution and random overlap
The transmission through a homogeneous slab is, from Eq. (1), where f(k) dk is the probability of the absorption coefficient being in the interval [k,k + dk] taking into account the weighting , and g(k) is the cumulative probability distribution (A.3)For each absorbing gas i, the opacity probability distribution f_{i}(k_{i}) can be derived, and need to be combined using the respective mixing ratios ζ_{i}. An exact procedure for doing this does not exist as the wavelength information is lost when deriving f_{i}(k_{i}). Assuming that random variables picked using these probability distributions are independent, however, they can be convolved to get the combined probability distribution f(k). We restrict the discussion here to combining the opacity distributions of two gases, as the procedure can easily be repeated to include an arbitrary number of gases.
In order to perform the convolution we need to take into account the mixing ratios of the gases, i.e. we need to find where . We know that (A.4)which yields (A.5)Having derived the opacity distributions f_{x}(k_{x}) and f_{y}(k_{y}) of two different gases x and y with mixing ratios ζ_{x} and ζ_{y}, the total opacity distribution will be given by the convolution of the two taking their respective mixing ratios into account: The total transmission becomes Changing the integration variable of the second integral to yields (A.10)Similarly changing the integration variable of the first integral to , and using(A.11)yields Equation (A.14) is identical to Eq. (6), i.e. the random overlap method is equivalent to convolving the opacity probability distributions.
All Tables
Methods for treating overlapping gaseous absorption considered here and our adopted acronyms.
All Figures
Fig. 1 Schematic illustration of the random overlap method. The original kdistributions of H_{2}O and CO are shown in panel a), with the combined kcoefficients shown in panel b) calculated using Eqs. (11) and (12) assuming ζ_{H2O} = 0.9 and ζ_{CO} = 0.1. Without resorting and rebinning (RO) the k_{xy,lm}coefficients with corresponding weights w_{xy,lm} are used directly. With resorting and rebinning (RORR) the combined kcoefficients in panel b) are resorted, as shown in panel c), and then rebinned down to a smaller number of kcoefficients, as shown in panel d). 

Open with DEXTER  
In the text 
Fig. 2 Schematic illustration of equivalent extinction (EE). An equivalent grey absorption is calculated for each minor gas (CO in our case), as shown in panel a), using Eq. (13) or 14. This grey absorption is added onto the kcoefficients of the major gas (H_{2}O in this example) using the corresponding mixing ratios, here assumed to be ζ_{H2O} = 0.9, ζ_{CO} = 0.1, as shown in panel b). 

Open with DEXTER  
In the text 
Fig. 3 P–T profiles used for the tests in the present work. Profiles are derived from the equilibrium P–T profile of Iro et al. (2005) as adopted by Cooper & Showman (2005, 2006), Rauscher & Menou (2010) and Heng et al. (2011) with the smoothing introduced by Mayne et al. (2014). We adopt the TiO/VO condensation curve from Fortney et al. (2008), plotted as a dotted line. 

Open with DEXTER  
In the text 
Fig. 4 Fluxes (left) and absolute errors in fluxes (right) for the night side P–T profile in Fig. 3. The linebyline (LbL) calculation was run at a 0.001 cm^{1} resolution corresponding to 5 × 107 monochromatic calculations, while we used 120 rebinned kcoefficients in each band corresponding to 3840 pseudomonochromatic calculations. Both results were obtained using ATMO. 

Open with DEXTER  
In the text 
Fig. 5 Same as Fig. 4 but for heating rates. 

Open with DEXTER  
In the text 
Fig. 6 Fluxes (left) and absolute errors in fluxes (right) for the day side P–T profile in Fig. 3. 

Open with DEXTER  
In the text 
Fig. 7 Same as Fig. 6 but for heating rates. 

Open with DEXTER  
In the text 
Fig. 8 Fluxes (left) and absolute errors in fluxes (right) obtained with the night side P–T profile in Fig. 3 using SOCRATES. Fluxes obtained using the random overlap method without resorting and rebinning (RO) are used to calculate errors for the random overlap with resorting and rebinning (RORR) with , 16 and 32, equivalent extinction (EE), adaptive equivalent extinction (AEE) and premixed opacities (PM). 

Open with DEXTER  
In the text 
Fig. 9 Same as Fig. 8 but for heating rates. L1 norms of the errors (see Amundsen et al. 2014, the average heating rate error weighted by the local heating rates) are 4.5 % for RORR with , 1.9 % for RORR with , 1.5 % for RORR with , 13 % for EE, 11 % for AEE and 38 % for PM. 

Open with DEXTER  
In the text 
Fig. 10 Fluxes (left) and heating rates (right) obtained with the night side P–T profile in Fig. 3 using constant mixing ratios equal to the mixing ratios at P =104 Pa, T =1000 K. This eliminates errors caused by the implicit interpolation of mixing ratios with PM which dominates the errors seen using this overlap method in Figs. 8 and 9. 

Open with DEXTER  
In the text 
Fig. 11 Fluxes (left) and absolute errors in fluxes (right) obtained with the day side P–T profile in Fig. 3 using SOCRATES. Fluxes obtained using the random overlap method without resorting and rebinning (RO) are used to calculate errors for the random overlap with resorting and rebinning (RORR) with , 16 and 32, equivalent extinction (EE), adaptive equivalent extinction (AEE) and premixed opacities (PM). 

Open with DEXTER  
In the text 
Fig. 12 Same as Fig. 11 but for heating rates. L1 norms of the errors (see the caption of Fig. 9) are 7.6 % for RORR with , 3.0 % for RORR with , 1.8 % for RORR with , 7.0 % for EE, 2.2 % for AEE and 119 % for PM. 

Open with DEXTER  
In the text 
Fig. 13 Fluxes (left) and heating rates (right) obtained with the day side P–T profile in Fig. 3 using constant mixing ratios equal to the mixing ratios at P =104 Pa, T =1900 K. This eliminates errors caused by the implicit interpolation of mixing ratios with PM which dominates the errors seen using this overlap method in Figs. 11 and 12. 

Open with DEXTER  
In the text 