Analytical Models of Exoplanetary Atmospheres. IV. Improved Two-stream Radiative Transfer for the Treatment of Aerosols
Kevin Heng Daniel Kitzmann






1University of Bern, Center for Space and Habitability, Gesellschaftsstrasse 6, CH-3012, Bern, Switzerland.  Emails: kevin.heng@csh.unibe.ch, daniel.kitzmann@csh.unibe.ch
2Johns Hopkins University, Department of Earth and Planetary Sciences, 301 Olin Hall, Baltimore, MD 21218, U.S.A.
3Johns Hopkins University, Department of Physics and Astronomy, Bloomberg Center for Physics and Astronomy, Baltimore, MD 21218, U.S.A.


We present a novel generalization of the two-stream method of radiative transfer, which allows for the accurate treatment of radiative transfer in the presence of strong infrared scattering by aerosols.  We prove that this generalization involves only a simple modification of the coupling coefficients and transmission functions in the hemispheric two-stream method.  This modification originates from allowing the ratio of the first Eddington coefficients to depart from unity.  At the heart of the method is the fact that this ratio may be computed once and for all over the entire range of values of the single-scattering albedo and scattering asymmetry factor.  We benchmark our improved two-stream method by calculating the fraction of flux reflected by a single atmospheric layer (the reflectivity) and comparing these calculations to those performed using a 32-stream discrete-ordinates method.  We further compare our improved two-stream method to the two-stream source function (16 streams) and delta-Eddington methods, demonstrating that it is often more accurate at the order-of-magnitude level.  Finally, we illustrate its accuracy using a toy model of the early Martian atmosphere hosting a cloud layer composed of carbon-dioxide ice particles.  The simplicity of implementation and accuracy of our improved two-stream method renders it suitable for implementation in three-dimensional general circulation models.  In other words, our improved two-stream method has the ease of implementation of a standard two-stream method, but the accuracy of a 32-stream method.






§ INTRODUCTION



Two-stream solutions have been studied for decades in the context of atmospheres and come in various flavors <cit.>.  They originate from a neat mathematical trick: instead of solving the radiative transfer equation for the intensity, one solves for its moments.  Besides the loss of angular information, the two-stream solution performs poorly when aerosols reside in the model atmosphere.  A longstanding result, based on geomorphic evidence, that early Mars was able to harbor liquid water on its surface (see  for a review), due to the scattering greenhouse effect (e.g., ) mediated by carbon-dioxide ice clouds <cit.>, was called into question because the original two-stream calculation over-estimated the degree of warming <cit.>.  Mars teaches us the lesson that the choice of radiative transfer method may alter the qualitative conclusion of a study, and inspires us to improve the accuracy of the two-stream solution in order to apply it broadly to exoplanetary atmospheres.

The main source of error appears to be the over-estimation of the amount of infrared radiation reflected by aerosols, which leads to an over-estimation of the scattering greenhouse effect.  On Earth, this effect is subdued because water clouds are strong infrared absorbers but weak infrared scatterers <cit.>.  On Mars, it is pronounced because carbon-dioxide ice clouds scatter infrared radiation strongly <cit.>.  Figure <ref> illustrates these differences.  In general, we expect the two-stream method to perform poorly in the presence of medium-sized to large aerosols that have single-scattering albedos between 0.5 and 1 in the infrared range of wavelengths.  This shortcoming motivates us to design an improved two-stream method that calculates the amount of reflected radiation accurately.  Operationally, we accomplish this feat by revisiting the formalism surrounding the Eddington coefficients previously elucidated by <cit.>.  



Specifically, we relax the assumption that the first Eddington coefficients[As already noted in <cit.>, there is no consensus on how to number/order these Eddington coefficients, and we use the convention of <cit.>.] are equal and allow their ratio to depart from unity.  The top-left panel of Figure <ref> shows our calculations for this ratio, E.  We also show calculations for single atmospheric layers populated by aerosols with fixed values of the single-scattering albedo (ω_0) and scattering asymmetry factor (g).  We consider single atmospheric layers, because if one attains understanding (and accuracy) for a single layer, then it straightforwardly generalizes to an arbitrary number of layers in a model atmosphere.  For the sake of discussion, we refer to small, medium-sized and large aerosols as having ω_0=0.1 and g=0 (isotropic scattering), ω_0=g=0.5 and ω_0=g=0.9 (predominantly forward scattering), respectively.  We will explore other choices later.  To simplify terminology, we term the fraction of flux reflected and transmitted by an atmospheric layer the “reflectivity" and “transmissivity", respectively.  In the example of a layer populated by medium-sized aerosols (top-right panel of Figure <ref>), we see that the original, hemispheric two-stream solution (e.g., ) over-estimates the true solution, which is computed using a 32-stream discrete-ordinates method via the open-source  computer code <cit.>.  Our improved two-stream solution with E=1 matches the reflectivity computed by the hemispheric two-stream method well; deviations are due to modifications we have made to the transmission function, as we will discuss.  As E is varied from 1 to 1.4, we see that the reflectivity varies rather sensitively.  The true solution is matched by a value of E between 1.1 and 1.2.  This example illustrates that small variations of E from unity allow us to improve the accuracy of the two-stream solution drastically.

Following through on this property of E, the bottom-left panel of Figure <ref> shows calculations of the reflectivity for atmospheric layers with small, medium-sized and large aerosols.  For each calculation, the value of E has been chosen to match the reflectivity by construction; in <ref>, we will explain in detail how this is accomplished.  For all three calculations, the original, hemispheric two-stream method over-estimates the reflectivity by ∼ 10%.  For completeness, we also show the transmissivity associated with these three examples (bottom-right panel of Figure <ref>), where we see that the discrepancies between the hemispheric two-stream calculations and the true solutions are less pronounced.

The overaching goal of the present study is to elucidate the theory behind the improvement of the two-stream method and the calculation of E (presented in <ref>).  We further demonstrate that our improved two-stream method rivals or betters the two-stream source function method of <cit.>, which is widely implemented in the exo-atmospheres literature, in both accuracy and simplicity of implementation (in <ref>).  We discuss the implications of our findings in <ref>.  The present study is the fourth in a series of papers devoted to constructing analytical models for exoplanetary atmospheres to both aid in the development of intuition and provide algorithms for computation, following <cit.> (for shallow-water fluid dynamics), <cit.> (for two-stream radiative transfer) and <cit.> (for equilibrium chemistry). 






§ GENERALIZING THE TWO-STREAM FORMALISM



In the two-stream formalism, the reflectivity and transmissivity of a single layer are respectively <cit.>

    f_ T = ( ζ_-^2 - ζ_+^2 )  T/( ζ_-  T)^2 - ζ_+^2,  f_ R = ζ_- ζ_+ ( 1 -  T^2 ) /ζ_+^2 - ( ζ_-  T)^2,
 
where ζ_± are the coupling coefficients and T is the transmission function.  The coupling coefficients relate the relative strength of transmission versus reflection, and generally depend on ω_0 and g.  When the layer is transparent, we have T=1, f_ T=1 and f_ R=0.  When the layer is opaque, we have T=0, f_ T=0 and f_ R = ζ_-/ζ_+.  These asymptotic limits suggest that plausible improvements to the two-stream solution are accomplished by modifying the coupling coefficients.

In the present study, we focus on the reflectivity and transmissivity of a single layer, and not its emissivity (blackbody emission), because previous studies dealing with aerosols have shown that the largest sources of error originate from the reflectivity <cit.>.  It has been previously shown that the expressions for the coupling coefficients with the hemispheric closure conserve energy by construction <cit.>.

To modify the coupling coefficients, we need to understand their physical origin.  In sacrificing accuracy for simplicity, the two-stream solutions contain an ambiguity: ratios of the moments of the intensity (mean intensity, flux, radiation pressure) are assumed to be constants known as Eddington coefficients <cit.>.  Two of these Eddington coefficients are set to be equal by enforcing the condition of radiative equilibrium in the limit of pure scattering <cit.>.  This occurs when the single-scattering albedo is exactly unity.  However, the two-stream solutions formally break down at exactly ω_0=1 <cit.> and this limit is rarely reached in practice, which render this condition academic.  If we relax this condition, then the coupling coefficients have a more general form,

    ζ_±≡1/2[ 1 ±√(E - ω_0/E - ω_0 g)].

where E is the ratio of the Eddington coefficients.  In the original two-stream solutions with the hemispheric (or hemi-isotropic) closure, we have E=1.  The first improvement is to use

    E = ω_0 ( 1 - g r^2 )/1 - r^2

to compute E1 values for use in the coupling coefficients.  Here, we have r ≡ (1-R_∞)/(1+R_∞) and R_∞ is the asymptotic value of the reflectivity when a layer is opaque (τ≫ 1).  The preceding expression is derived by setting R_∞ = ζ_-/ζ_+ and using equation (<ref>).  It is worth emphasizing that this improvement ensures the asymptotic reflectivity matches the true solution by construction, as long as we have a way of computing R_∞.

In Figure <ref> (top-left panel), the grid of values for E obtains from computing R_∞ using the  code <cit.>, which uses the discrete-ordinates method of radiative transfer <cit.>.  We use 32-stream  calculations as the ground truth.  It should be emphasized that this is the entire parameter space of interest for aerosols embedded in atmospheres.  Part of the simplicity of the method is that the function E(ω_0, g) only needs to be computed once.  It may then be stored and used for all future calculations.

The second improvement we make is to modify the transmission function.  In the limit of pure absorption, the exact solution of the radiative transfer equation yields T = 2 E_3(τ), an expression that formally integrates over all angles <cit.>.  Here, E_3 is the exponential integral of the third order <cit.> and τ is the optical depth of the atmospheric layer.  For transmission, we use T = 2 E_3 ( τ^'), where τ^' = τ√((1-ω_0)(1-ω_0 g)) and the additional factor derives from the hemispheric two-stream solution <cit.>.  For reflection, we use T = 2 E_3( τ^'') and τ^'' = τ√((E-ω_0)(E-ω_0 g)).  



§ RESULTS





 §.§ Comparison to <cit.>






We now compare calculations of the transmissivity and reflectivity to those performed using other commonly used methods: hemispheric two-stream <cit.>[We specifically use the two-stream formalism written down by <cit.>; we are not claiming that <cit.> should be solely cited for the two-stream method.], delta-Eddington <cit.> and two-stream source function <cit.>.  The two-stream source function method is of particular interest, because it is widely implemented in the exo-atmospheres literature <cit.>.  It achieves a multi-stream solution using a clever mathematical trick: inserting the two-stream solution into the term of the radiative transfer equation involving the scattering phase function, which implies that the solution is, strictly speaking, not self-consistent.  We use the two-stream source function method with 8 streams in each hemisphere for a total of 16 streams, and sum up these streams numerically (weighted by the cosine of the polar angle) to construct the fluxes.

The delta-Eddington method uses the Eddington closure (see  or ), but includes an additional feature: it approximates the scattering phase function as consisting of a Dirac-delta function and a series expansion involving the cosine of the scattering angle <cit.>.  The motivation behind this approximation is to attain higher accuracy for radiative transfer with large aerosols, which tend to produce a strong forward peak in the scattered intensity.  However, the delta-Eddington method has been criticized as being ad hoc, as the relative weighting of the Dirac-delta function and series terms is chosen arbitrarily <cit.>.[<cit.> describes how the delta-Eddington method includes the procedures of truncation and renormalization, describes them as being “ad hoc" and remarks how “none of the proposed variants is demonstrably superior to any other."  <cit.> further adds that “the number of reasonable truncation and renormalization procedures is limited only by one's imagination."]  This criticism provided us with motivation to avoid the delta-Eddington method and its variants when constructing our improved two-stream method.  However, we include the delta-Eddington method in our comparisons because our emphasis is on reproducing Figures 2 and 3 of <cit.> as a benchmarking exercise.

We again use 32-stream  calculations of the transmissivity and reflectivity as the ground truth.  Figure <ref> shows the transmissivity and reflectivity for three different sets of values of ω_0 and g, chosen to faciliate comparison with the study of <cit.>.  Most of the absorption or scattering of radiation by an atmospheric layer occurs at τ∼ 1.  Dips in the error curves occur when the curve of the reflectivity or transmissivity intersects the  curve, such that the error between them formally drops to zero.  In practice, it drops to nearly zero, because the numerically-computed curves sampled at discrete points do not formally intersect.  The six sets of calculations in Figure <ref> suggest that our improved two-stream method achieves comparable or superior accuracy, compared to the other methods, often at the order-of-magnitude level, but with less computational effort.

Next, we focus on the error associated with the reflectivity, as it is known to be larger than for the transmissivity in atmospheric calculations involving aerosols <cit.>.  Figure <ref> quantifies the error incurred by using our improved two-stream method as a function of ω_0, g and τ.  When τ=10, the error is ∼ 0.01%.  This is unsurprising, as it is by construction.  When τ=1, the error is typically ∼ 1%, unless both ω_0 ≈ 1 and g ≈ 1, in which case the error approaches 10%.  In practice, the presence of gas, consisting of atoms and molecules, in the atmosphere reduces the value of ω_0 to below unity because they provide a source of absorption (which increases the total cross section), meaning that these large errors are rarely encountered.  Figure <ref> also shows the ratio of errors of the improved two-stream versus the two-stream source function methods.  For most combinations of ω_0 and g, the improved two-stream method is more accurate than the two-stream source function method by an order of magnitude, unless ω_0 ≳ 0.9.  This is remarkable, because it has the implementational simplicity of the two-stream method, but an accuracy that is superior to a 16-stream method.



 §.§ Toy model of early Martian atmosphere




So far, we have examined calculations with fixed values of ω_0 and g, because we were focused on benchmarking our improved two-stream method.  Real aerosols are associated with ω_0 and g that are functions of wavelength <cit.>.  To illustrate this behavior, we consider a toy model atmosphere of early Mars that hosts a cloud layer composed of medium-sized to large carbon-dioxide ice particles.  We do not compute a more realistic Mars model, because this has already been done in <cit.>.  Conceptually, understanding the interplay between aerosols and radiative transfer under Mars-like conditions is relevant to defining the outer edge of the classical habitable zone.

We imagine a scenario where starlight penetrates the atmosphere and heats up its surface (much like on Earth), which then re-emits the heat as infrared radiation.  The infrared radiation attempts to escape the atmosphere, but encounters the carbon-dioxide-ice cloud layer (assumed to be located at 0.1 bar), which reflects some of it back to the surface and heats up the atmosphere below the cloud layer.  The ice particles are assumed to follow a gamma distribution with an effective particle radius of 25 μm.  We assume that the atmosphere is dominated by gaseous carbon dioxide and that it has a temperature equal to the condensation temperature of carbon dioxide just below the cloud layer.  The absorptivity of the gas is parametrized by a grey opacity with a value chosen such that the optical depth contributed by the gas alone is 0.1.  We calculate the reflectivity of the cloud layer using realistic, wavelength-dependent single-scattering albedos and scattering asymmetry parameters ( and Figure <ref>).  

In Figure <ref>, our improved two-stream calculations yield errors at 1% or lower compared to 32-stream  calculations.  By contrast, the hemispheric two-stream and two-stream source function methods yield errors at the 5–10% level, depending on how opaque the cloud is.  As already eludicated by <cit.> using more realistic models of the early Martian atmosphere, these errors translate into an over-estimation of the surface temperature by about 40 K, enough to alter the qualitative conclusion.  Similarly, we anticipate that sophisticated simulations of the outer edge of the habitable zone using three-dimensional climate models would benefit from more accurate radiative transfer calculations using our improved two-stream method.



§ DISCUSSION



In the absence of scattering, the outgoing and incoming two-stream fluxes are decoupled.  The boundary condition at the bottom of the atmosphere (remnant heat from the formation of the exoplanet for gas giants and surface fluxes for rocky exoplanets) may be propagated upwards using the solution for the outgoing flux.  The boundary condition at the top of the atmosphere (stellar irradiation) may be propagated downwards using the solution for the incoming flux.  For each layer, one now has the values of the outgoing and incoming fluxes.  Taking the difference between these fluxes yields the net flux, which is then fed into the first law of thermodynamics to compute the change in temperature of each layer <cit.>.  Updating the temperature in each layer in turn alters the opacities and the fluxes.  This iteration is performed numerically for a model atmosphere with multiple layers, until convergence is attained and radiative equilibrium is reached, e.g., see the implementation of <cit.>.

When scattering is present, the pair of two-stream solutions feed into each other and an additional iteration is needed.  Physically, radiation may be scattered multiple times as it travels from a layer to its immediate neighbors and beyond.  It is worth emphasizing that this iteration is to enforce multiple scattering, and is distinct from the iteration for radiative equilibrium.  The alternative to such an iteration is to perform matrix inversion.  Instead of describing a pair of layers, the two-stream solutions may be used to describe the fluxes within a single layer bounded by two interfaces.  An atmosphere with a finite number of layers is represented by a set of equations, with each equation describing the fluxes at the layer center and interfaces.  Mathematically, this set of equations makes up a tridiagonal matrix, which may then be inverted using Thomas's algorithm <cit.>.  Since Thomas's algorithm is essentially a set of recursive algebraic relations, the procedure remains efficient.  This property renders the solutions feasible for implementation in three-dimensional general circulation models (e.g., ), which require radiative transfer to be highly efficient in order to simulate climates for ∼ 10^7 time steps or more.



We acknowledge financial support from the Swiss National Science Foundation, the PlanetS National Center of Competence in Research (NCCR), the Center for Space and Habitability (CSH) and the Swiss-based MERAC Foundation.  KH acknowledges a Visiting Professorship at Johns Hopkins University in the Zanvyl Krieger School of Arts and Sciences, held during the final revision and resubmission of this manuscript.



99

[Arfken & Weber(1995)]arfken Arfken, G.B., & Weber, H.J.  1995, Mathematical Methods for Physicists, fourth edition (San Diego: Academic Press)

[Cahoy et al.(2010)]cahoy10 Cahoy, K.L., Marley, M.S., & Fortney, J.J.  2010, Astrophysical Journal, 724, 189

[Chandrasekhar(1960)]chandra60 Chandrasekhar, S.  1960, Radiative Transfer (New York: Dover)

[Draine(2003)]draine03 Draine, B.T.  2003, Astrophysical Journal, 598, 1017

[Forget & Pierrehumbert(1997)]fp97 Forget, F., & Pierrehumbert, R.T.  1997, Science, 278, 1273

[Fortney et al.(2008)]fortney08 Fortney, J.J., Lodders, K., Marley, M.S., & Freedman, R.S.  2008, Astrophysical Journal, 678, 1419

[Goody & Yung(1989)]gy89 Goody, R.M., & Yung, Y.L.  1989, Atmospheric Radiation: Theoretical Basis, second edition (New York: Oxford)

[Hamre et al.(2013)]hamre13 Hamre, B., Stamnes, S., Stamnes, K., & Stamnes, J.J.  2013, AIP Conference Proceedings, 1531, 923

[Heng et al.(2012)]hhps12 Heng, K., Hayek, W., Pont, F., & Sing, D.K.  2012, Monthly Notices of the Royal Astronomical Society, 420, 20

[Heng & Workman(2014)]hw14 Heng, K., & Workman, J.  2014, Astrophysical Journal Supplements, 213, 27

[Heng et al.(2014)]hml14 Heng, K., Mendonça, J.M., & Lee, J.-M.  2014, Astrophysical Journal Supplements, 215, 4

[Heng & Tsai(2016)]ht16 Heng, K., & Tsai, S.-M.  2016, Astrophysical Journal, 829, 104

[Heng(2017)]heng17 Heng, K.  2017, Exoplanetary Atmospheres: Theoretical Concepts and Foundations (Oxford: Princeton University Press)

[Joseph et al.(1976)]joseph76 Joseph, J.H., Wiscombe, W.J., & Weinman, J.A.  1976, Journal of the Atmospheric Sciences, 33, 2452

[Kitzmann et al.(2013)]kitzmann13 Kitzmann, D., Patzer, A.B.C., & Rauer, H.  2013, Astronomy & Astrophysics, 557 A6

[Kitzmann(2016)]kitzmann16 Kitzmann, K.  2016 Astrophysical Journal Letters, 817, L18

[Malik et al.(2017)]malik17 Malik, M., Grosheintz, L., Mendonça, J.M., et al.  2017, AJ, 153, 56

[Marley & McKay(1999)]mm99 Marley, M.S., & McKay, C.P.  1999, Icarus, 138, 268

[Meador & Weaver(1980)]mw80 Meador, W.E., & Weaver, W.R.  1980, Journal of the Atmospheric Sciences, 37, 630

[Mihalas(1970)]mihalas70 Mihalas, D.  1970, Stellar Atmospheres, first edition (San Francisco: Freeman)

[Mihalas(1978)]mihalas78 Mihalas, D.  1978, Stellar Atmospheres, second edition (San Francisco: Freeman)

[Morley et al.(2013)]morley13 Morley, C.V., Fortney, J.J., Kempton, E. M.-R., et al.  2013, Astrophysical Journal, 775, 33

[Pierrehumbert(2010)]pierrehumbert10 Pierrehumbert, R.T.  2010, Principles of Planetary Climate (New York: Cambridge)

[Schuster(1905)]s05 Schuster, A.  1905, ApJ, 21, 1

[Showman et al.(2009)]showman09 Showman, A.P., Fortney, J.J., Lian, Y., et al.  2009, Astrophysical Journal, 699, 564

[Stamnes et al.(1988)]s88 Stamnes, K., Tsay, S.-C., Wiscombe, W., & Jayaweera, K.  1988, Applied Optics, 27, 2502

[Toon et al.(1989)]toon89 Toon, O.B., McKay, C.P., & Ackerman, T.P.  1989, Journal of Geophysical Research, 94, 16287

[Wiscombe(1977)]wiscombe77 Wiscombe, W.J.  1977, Journal of the Atmospheric Sciences, 34, 1408

[Wordsworth(2016)]wordsworth16 Wordsworth, R.D.  2017, Annual Review of Earth & Planetary Sciences, 44, 381





