EDP Sciences
Free Access
Issue
A&A
Volume 598, February 2017
Article Number A97
Number of page(s) 10
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/201629322
Published online 09 February 2017

© ESO, 2017

1. Introduction

Rapid calculation of wavelength-integrated fluxes and heating rates are needed in most planetary and brown dwarf atmosphere models. As line-by-line approaches are too computationally expensive for practical use, the correlated-k method (Goody et al. 1989; Lacis & Oinas 1991; Thomas & Stamnes 2002) has been applied in retrieval models (e.g. Irwin et al. 2008), one-dimensional (1D) atmosphere models (e.g. Marley et al. 1996; Burrows et al. 1997; Fortney et al. 2005), and three-dimensional (3D) global circulation models (GCMs; e.g. Showman et al. 2009; Kataria et al. 2013; Amundsen et al. 2016). With the correlated-k 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 k-coefficients and corresponding weights. Pseudo-monochromatic calculations are performed using these k-coefficients, decreasing the required computation time by several orders of magnitude compared to line-by-line calculations.

The treatment of overlapping gaseous absorption, i.e. absorption by more than one gas in a single spectral interval, with the correlated-k 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 k-coefficients 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.

    Pre-mixedk-coefficients (PM, Goodyet al. 1989): k-coefficients for themixture are computed directly from the total line-by-line gasopacity.

  • 2.

    The random overlap method, both without (RO) and with (RORR) resorting and rebinning (Lacis & Oinas 1991): k-coefficients are computed for each gas and combined assuming their absorption cross-sections are uncorrelated.

  • 3.

    Equivalent extinction (EE, Edwards 1996): k-coefficients are computed for each gas and combined using an “equivalent grey absorption” for all minor absorbers and all k-coefficients for the major absorber in each band.

Pre-mixed k-coefficients 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 k-coefficients for different gases, but is inflexible as mixing must be assumed before k-coefficients are computed. Alternatively, gas mixing ratios can be added as dimensions to the look-up table of k-coefficients, 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 cross-sections of different gases are uncorrelated. The total number of k-coefficients in a band scales as the product of the number of k-coefficients for each overlapping gas, causing this method to become computationally expensive, but resorting and rebinning the resulting k-coefficients 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 pre-mixed k-coefficients, 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 k-coefficients individually for each gas and combine them on-the-fly in models using the current local mixing ratios. As all wavelength information is lost when the k-coefficients 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 cross-sections 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 line-by-line calculations.

In this paper we compare pre-mixing, random overlap and equivalent extinction in terms of computational efficiency and evaluate their accuracy by comparing to results from line-by-line calculations. In Sect. 2 we give a brief overview of the correlated-k 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 radiative-convective equilibrium atmosphere code ATMO (Tremblin et al. 2015, 2016) and our GCM radiation scheme SOCRATES1 (Edwards & Slingo 1996; Edwards 1996; Amundsen et al. 2014). We give our concluding remarks in Sect. 5.

2. The correlated-k method

As treating the wavelength-dependence of gaseous absorption explicitly is too computationally expensive to be performed in many atmosphere models, the correlated-k 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 one-dimensional (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 correlated-k method in detail here, but refer to e.g. Lacis & Oinas (1991), Goody et al. (1989) and Thomas & Stamnes (2002) for in-depth discussions. Note that we have previously verified the applicability of the correlated-k method in hot Jupiter and brown dwarf atmosphere models (Amundsen et al. 2014).

In the correlated-k method the opacity spectrum is divided into bands b. In each band k-coefficients 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.

Pseudo-monochromatic 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 nb is the number of bands.

The -coefficients are the k-coefficients 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 k-coefficients 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.

Table 1

Methods for treating overlapping gaseous absorption considered here and our adopted acronyms.

3.1. Pre-mixed

The total absorption coefficient can be calculated by summing line-by-line absorption coefficients for all absorbing species weighted by their relative abundances: (5)where the sum is over all Ns 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 ktotρ, where ρ is the total gas density. ktot can be used to compute and tabulate k-coefficients for the gas mixture as a function of temperature and pressure. This approach has several advantages: it is fast, requiring only one set of k-coefficients 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 k-coefficients. A potential solution would be to add gas mixing ratios as dimensions to the look-up table of k-coefficients, 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). k-coefficients are computed individually for each gas, and k-coefficients 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 (ux,uy) 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 line-by-line 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 k-coefficients for the individual gases x and y. The transmission through one layer is, using Eqs. (1), (2) and (6), Defining uxy = ux + uy, we can write the above transmission as (9)where and (12)For illustration we show schematics of the k-coefficients H2O and CO in Fig. 1a, with the combined k-coefficients in Fig. 1b calculated using Eq. (11) assuming ζH2O = 0.9 and ζCO = 0.1.

Running nk,xnk,y pseudo-monochromatic calculations using the kxy,lm-coefficients the total flux can be calculated as usual using Eqs. (3) and (4) with the weights wxy,lm. This procedure can be replicated for an arbitrary number of gases, however, the computation time increases by a factor of nk for each gas added. This method therefore quickly becomes too computationally expensive for practical use.

thumbnail Fig. 1

Schematic illustration of the random overlap method. The original k-distributions of H2O and CO are shown in panel a), with the combined k-coefficients shown in panel b) calculated using Eqs. (11) and (12) assuming ζH2O = 0.9 and ζCO = 0.1. Without resorting and rebinning (RO) the kxy,lm-coefficients with corresponding weights wxy,lm are used directly. With resorting and rebinning (RORR) the combined k-coefficients in panel b) are resorted, as shown in panel c), and then rebinned down to a smaller number of k-coefficients, 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 kxy,lm-coefficients and rebinning them to obtain a smaller number of k-coefficients , would circumvent the scaling issue. The procedure is illustrated in panels c and d of Fig. 1. First the k-coefficients of two gases are combined using Eqs. (11) and (12), as shown in Fig. 1b. These nk,xnk,yk-coefficients are sorted in increasing order, with the weights sorted using the same mapping, as shown in Fig. 1c2. The sorted kxy,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 kxy,lm. As kxy,lm will vary vertically, however, equal spacing in log kxy,lm in one layer will not correspond to equal spacing in log kxy,lm at a different level. Consequently a compromise must be made, with one possible solution being an equal spacing in log kxy,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 kxy,lm will change meaning the ideal spacing to change. In addition, particular care must be taken to treat cases where k-coefficients 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 Gauss-Legendre quadrature, while in ATMO we use uniform weights, both of which can have an arbitrary number of rebinned k-coefficients . The rebinned coefficients are found by computing a weighted average of all kxy,lm-terms belonging to each bin , where wxy,lm are used as weights. If a kxy,lm-term extends over more than a single bin, it is split over neighbouring bins such that the weights wxy,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 kxy,lm. This causes a significant increase in computation time, however, due to the need to compute the logarithm of many k-terms, and would also require particular care to treat cases where kxy,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 k-coefficients are used to compute the fluxes and heating rates for the atmosphere. This approach is consequently much more flexible than pre-mixing gases as gas abundances can be set at run-time.

3.3. Equivalent extinction

thumbnail 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 k-coefficients of the major gas (H2O 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 kx,l are the k-coefficients of the minor gas in the layer with corresponding weights wx,l, and Fv,l is the thermal flux in the layer including only absorption by k-term l of the gas. Pseudo-monochromatic calculations are performed for all nkk-coefficients of the major gas in each band, with all other absorbers included by using the equivalent grey absorption . This effectively reduces the number of pseudo-monochromatic calculations required to one per k-coefficient 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 non-zero if scattering is included, the equivalent extinction is defined by (14)Fs ∗ ,l is the direct flux at the lower boundary including only k-term l of the gas. The use of Fs ∗ ,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 k-coefficient 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 k-distributions are shown in Fig. 1a. From the k-coefficients 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 k-coefficients of the minor gas using the corresponding mixing ratios, here assumed to be ζH2O = 0.9 and ζCO = 0.1, with the combined k-coefficients 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 H2O assuming all O is in the form of H2O. 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 k-coefficients 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 k-coefficients 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 dwarf-like atmospheres by comparing them to line-by-line calculations. We have previously tested the applicability of the correlated-k 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 total3. Abundances are as in Amundsen et al. (2016), and we include TiO and VO opacity for the day side PT 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 correlated-k 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 correlated-k 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 PT profiles are shown in Fig. 3, which are derived from the equilibrium PT 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 line-by-line (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 k-terms 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 line-by-line 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 line-by-line fluxes and heating rates, obtained for the night side PT profile shown in Fig. 3. It is clear that both fluxes and heating rates obtained when using the correlated-k method with the random overlap method match the line-by-line results very well, with errors of a few percent. We note that these errors are both due to the use of the correlated-k method and the random overlap assumption, and in agreement with the errors found in Amundsen et al. (2014).

thumbnail Fig. 3

PT profiles used for the tests in the present work. Profiles are derived from the equilibrium PT 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

thumbnail Fig. 4

Fluxes (left) and absolute errors in fluxes (right) for the night side PT profile in Fig. 3. The line-by-line (LbL) calculation was run at a 0.001 cm-1 resolution corresponding to 5 × 107 monochromatic calculations, while we used 120 rebinned k-coefficients in each band corresponding to 3840 pseudo-monochromatic calculations. Both results were obtained using ATMO.

Open with DEXTER

thumbnail 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 PT profile shown in Fig. 3. As for the night side PT profile both fluxes and heating rates obtained with the correlated-k method using random overlap are in good agreement with the corresponding line-by-line results, with errors of a few percent. As for the night side these errors are both due to the use of the correlated-k method and the random overlap assumption, and in agreement with the errors found in Amundsen et al. (2014).

thumbnail Fig. 6

Fluxes (left) and absolute errors in fluxes (right) for the day side PT profile in Fig. 3.

Open with DEXTER

thumbnail 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 k-terms 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 two-stream 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 PT 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 k-terms 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 PT profile.

thumbnail Fig. 8

Fluxes (left) and absolute errors in fluxes (right) obtained with the night side PT 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 pre-mixed opacities (PM).

Open with DEXTER

thumbnail 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

Pre-mixed (PM) opacities are significantly less accurate than all the above overlap treatments, this stems from errors introduced by the interpolation in the pre-mixed opacity table. Changes in mixing ratios with temperature and pressure can cause large changes in the pre-mixed 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 PT profile, but using constant mixing ratios equal to the mixing ratios at P =104 Pa, T =1000 K both when computing the pre-mixed opacity table and when combining k-coefficients 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 k-coefficients, 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).

thumbnail Fig. 10

Fluxes (left) and heating rates (right) obtained with the night side PT 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 k-terms 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 k-coefficients, 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.

Table 2

Computation times of the thermal fluxes in SOCRATES for various overlap treatments using the night side PT profile in Fig. 3 not including TiO and VO opacity.

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 PT 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 k-terms 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.

thumbnail Fig. 11

Fluxes (left) and absolute errors in fluxes (right) obtained with the day side PT 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 pre-mixed opacities (PM).

Open with DEXTER

thumbnail 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 pre-mixed 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 PT-profile 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 PT points in the look-up k-coefficient 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 PT 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.

thumbnail Fig. 13

Fluxes (left) and heating rates (right) obtained with the day side PT 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 k-terms about 3 times slower than (A)EE.

Table 3

Computation times of the thermal fluxes in SOCRATES for various overlap treatments using the day side PT profile in Fig. 3 including TiO and VO opacity.

Table 4

Computation times of the stellar fluxes in SOCRATES for various overlap treatments using the day side PT profile in Fig. 3 including TiO and VO opacity.

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 k-terms.

Equivalent extinction (EE) is about three times faster than RORR with only 8 rebinned k-terms, 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 pre-mixed 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 pre-mixed 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 H2O, CH4 and CO may also lead to inaccuracies as seen in Sect. 4.2.1. A lower resolution can be tolerated in k-coefficient tables of individual gases than in pre-mixed tables as individual opacities vary more smoothly.

In our 1D atmosphere code ATMO we use RORR, usually with about 30 rebinned k-terms in each band. This gives us the flexibility to manipulate gas mixing ratios without recomputing or having to add additional dimensions to the k-table. It also uniquely allows us to treat non-equilibrium chemistry consistently, i.e. have non-equilibrium abundances feed back on the total opacity and consequently the calculated PT profiles (Tremblin et al. 2015, 2016; Drummond et al. 2016).

Unfortunately, RORR, even with only 8 rebinned k-coefficients, 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 short-wave scattering.


2

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.

3

We include the opacity of H2O (Barber et al. 2006), CO (Rothman et al. 2010), CH4 (Yurchenko & Tennyson 2014), NH3 (Yurchenko et al. 2011), H2–H2 and H2–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/2007-2013 Grant Agreement No. 247060-PEPS and grant No. 320478-TOFU). 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 Super-computer, a DiRAC Facility jointly funded by STFC, the Large Facilities Capital Fund of BIS and the University of Exeter.

References

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 fi(ki) 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 fi(ki). 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 fx(kx) and fy(ky) 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

Table 1

Methods for treating overlapping gaseous absorption considered here and our adopted acronyms.

Table 2

Computation times of the thermal fluxes in SOCRATES for various overlap treatments using the night side PT profile in Fig. 3 not including TiO and VO opacity.

Table 3

Computation times of the thermal fluxes in SOCRATES for various overlap treatments using the day side PT profile in Fig. 3 including TiO and VO opacity.

Table 4

Computation times of the stellar fluxes in SOCRATES for various overlap treatments using the day side PT profile in Fig. 3 including TiO and VO opacity.

All Figures

thumbnail Fig. 1

Schematic illustration of the random overlap method. The original k-distributions of H2O and CO are shown in panel a), with the combined k-coefficients shown in panel b) calculated using Eqs. (11) and (12) assuming ζH2O = 0.9 and ζCO = 0.1. Without resorting and rebinning (RO) the kxy,lm-coefficients with corresponding weights wxy,lm are used directly. With resorting and rebinning (RORR) the combined k-coefficients in panel b) are resorted, as shown in panel c), and then rebinned down to a smaller number of k-coefficients, as shown in panel d).

Open with DEXTER
In the text
thumbnail 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 k-coefficients of the major gas (H2O 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
thumbnail Fig. 3

PT profiles used for the tests in the present work. Profiles are derived from the equilibrium PT 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
thumbnail Fig. 4

Fluxes (left) and absolute errors in fluxes (right) for the night side PT profile in Fig. 3. The line-by-line (LbL) calculation was run at a 0.001 cm-1 resolution corresponding to 5 × 107 monochromatic calculations, while we used 120 rebinned k-coefficients in each band corresponding to 3840 pseudo-monochromatic calculations. Both results were obtained using ATMO.

Open with DEXTER
In the text
thumbnail Fig. 5

Same as Fig. 4 but for heating rates.

Open with DEXTER
In the text
thumbnail Fig. 6

Fluxes (left) and absolute errors in fluxes (right) for the day side PT profile in Fig. 3.

Open with DEXTER
In the text
thumbnail Fig. 7

Same as Fig. 6 but for heating rates.

Open with DEXTER
In the text
thumbnail Fig. 8

Fluxes (left) and absolute errors in fluxes (right) obtained with the night side PT 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 pre-mixed opacities (PM).

Open with DEXTER
In the text
thumbnail 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
thumbnail Fig. 10

Fluxes (left) and heating rates (right) obtained with the night side PT 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
thumbnail Fig. 11

Fluxes (left) and absolute errors in fluxes (right) obtained with the day side PT 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 pre-mixed opacities (PM).

Open with DEXTER
In the text
thumbnail 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
thumbnail Fig. 13

Fluxes (left) and heating rates (right) obtained with the day side PT 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

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

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

Initial download of the metrics may take a while.