The radiative transfer model for the greenhouse effect

Radiative transfer is at the heart of the mechanism to explain the greenhouse effect based on the partial infrared opacity of carbon dioxide, methane and other greenhouse gases in the atmosphere. In absence of thermal diffusion, the mathematical model consists of a first order integro-differential equation coupled with an integral equation for the light intensity and the temperature, in the atmosphere. We revisit this much studied system from a mathematical and numerical point of view. Existence and uniqueness and implicit solutions of the Milne problem for grey atmospheres are stated. Numerical simulations are given for grey and non-grey atmospheres and applied to calculate the effect of greenhouse gases. In the context of a transparent atmosphere for sunlight, it is found that by doubling the absorption coefficient in the infrared absorption range of CO2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\texttt {CO}_2$$\end{document} the temperature decreases by 2%. On the other hand, the same changes but in the low infrared range of the sunlight leads to an increase of temperature in the atmosphere. Several computer codes were written to cross-validate the results. The authors conclude that the radiative transfer model without thermal diffusion for an atmosphere transparent to the incident sunlight is not capable of explaining the greenhouse effect due to the greenhouse gases. A decreasing temperature due to an increasing proportion of CO2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\texttt {CO}_2$$\end{document} has been observed in the high atmosphere (D.W.J. Thomson et al, nature11579). In the lower atmosphere thermal diffusion and convection cannot be neglected and since the absorption coefficient are highly dependent on the temperature, a full ocean–atmosphere–biosphere climate model is required. Hence, driving conclusions from this study on climate change should be cautiously avoided and a review of the hypothesis of the radiative transfer argument commonly found in textbooks should be revisited.


Introduction
The greenhouse effect is an important element of the current theory of climate change. Some gases in the Earth atmosphere like carbon dioxide C0 2 and methane CH 4 absorb infrared rays and thus contribute to a global warming of our planet. As explained in [9,14,23,26] the Sun radiates light with a heat flux Q = 1370 Watt/m 2 , in the frequency range (0.5, 20) × 10 14 Hz corresponding approximately to a black body at temperature of 5800 K; 74% of this light intensity reaches the ground because the atmosphere is almost transparent to this spectrum and about 36% is reflected back by the clouds or the ocean, snow, etc. (albedo).
The Earth behaves almost like a black body at temperature T e = 288 K and as such radiates rays of frequencies ν in the infrared spectrum (0.03, 0.6) × 10 14 Hz.
So both the Sun and the Earth are approximate black bodies. The Planck theory says that a black body at temperature T radiates electromagnetic waves in the entire frequency spectrum ν ∈ R + with intensities given by the Planck function: where is the Planck constant, c is the speed of light in the medium and k is the Boltzmann constant.
A major discrepancy between reality and the black body theory for Earth is shown on Fig. 1. It is due to the partial transparency of the atmosphere and the absorbing power of CO 2 , H 2 O, O 3 , CH 4 , etc., in the infrared range. Figure 2 gives the absorption coefficients κ ν for some gas (a transparent gas has κ ν = 0, and 1 if it is opaque.) Consequently, the infrared light emitted by Earth, seen from far, has a defect of radiance in its spectrum which is affected by the proportion of Green-House Gases (GHG): it is the greenhouse effect.
In this article we propose a mathematical and numerical investigation to quantify this phenomenon, using a popular assumption on the transparency of the atmosphere to the incoming sunlight ( [14], eq. (2.16), (2.18)).  Absorption coefficient κ ν of some gases of the atmosphere in the range of frequencies of interest, but versus wave length c/ν. (reprinted from wikipedia). Adding more CO 2 increases κ ν in the range (8,15)μm Photons travel at the speed of light; energy balance can thus be assumed instantaneous. The atmosphere is affected by wind, rain, chemistry, etc., but at a very different time-scale; it is believed -and to some extend asserted, (see [14]) -that even if all these other phenomena are ignored, still the greenhouse effect, is present in the equations, and sufficient to explain, partially, global warming. Consequently, in the article, we restrict the analyses to the energy conservation equations for the radiative intensity and the temperature, Eq. (3) below.
Radiative transfers have been studied by astronomers, nuclear physicists, combustionists and many other. Their work is summarized in [8,28]. Mathematical analyses are also numerous and we send the readers to [11,15,20].
More recently, for obvious reasons, there is a renewed interest in numerical simulations of radiative transfers. Among others the reader is sent to [10,18,19,22,24]. However, we are not aware of a simulation of the very specific greenhouse gases (GHG) effect, as presented here, namely: Compare an atmosphere in which the absorption coefficient is κ 0 at all frequencies with an atmosphere in which the absorption coefficient is κ 0 , except in an infrared frequency interval [ν 1 , ν 2 ], where it is 1 2 κ 0 . The results are unexpected: the temperature of the first case is 2% less than in the second case. This means that, with this model, greenhouse gases cool the Earth atmosphere! The authors are applied mathematicians, with a limited knowledge of atmosphere dynamics, yet competence in fluid mechanics and nuclear engineering, two fields concerned with radiative transfers. So climatologists were consulted. Apparently this cooling phenomenon is observed in the stratosphere [13,17]; however the changes in 0 3 may to play a role, something which is not taken into account here. In the lower atmosphere, neglecting thermal diffusion is not allowed because any change in temperature has a strong effect on the density of the gases and on the absorption coefficient, creating winds, clouds etc.
Consequently this study does not upset the body of knowledge for climate studies, but it invalidates the simple heuristic explanations of the green house effect found many books like [23] and on the internet like the Al Gore experiment www.climaterealityproject.org/video/ climate-101-bill-nye. In Sect. 2, the paper begins with a recall of the radiative transfer equations as derived from fundamental thermodynamic principles. Then, in Sect. 3, the various approximations are reviewed using isotropy and also the special case where the absorption and scattering coefficients do not depend on the light frequencies: Milne's problem for grey atmosphere. On the way the Chandrasekhar approximation is recalled, to account for the Earth radius. The Chandrasekhar correction is very small but the equations are less singular at the poles [20]. Being important in plasma physics Milne's problem has received many mathematical developments; these are reviewed in Sect. 4 and 5, applied to the grey atmosphere and compared with Fowler's [14]. The numerical study begins in Sect. 6. A finite element method with SUPG correction is proposed and compared with a semi-analytical method based on an integral representation of the light intensities. The later is computationally more expensive but it allows the validation of the first method and demonstrate that 1/ the effect of the Chandrasekhar correction can neglected numerically also, 2/ the (needed) SUPG correction does not degrade the precision of the result. Section 7 deals with the greenhouse effect; as explained above two cases are compared: one with a constant absorption coefficient κ = κ 0 , another with κ(ν) = κ 0 − 1 (ν 1 ,ν 2 ) δκ, δκ constant positive and less than κ 0 . Figure 8 shows that the two computer codes and a calculus of variation done on the first one, all give a decrease of temperature from case 2 to case 1: T (case 2)> T (case 1) at all altitudes when (ν 1 , ν 2 ) is in the infrared range and the opposite when (ν 1 , ν 2 ) is in the lower range of the sunlight. However the precision is poor.
To assert the results, in Sect. 7.2, an iterative method is proposed whereby the light intensity is eliminated and the resulting system is a set of coupled integral equations for the temperature, function of altitude. The numerical results confirm the previous ones and with this last formulation the results are independent of the mesh and time steps: the method is much more robust and does not require to compute singular integrals.
In Sect. 8, a brief asymptotic analysis shows that the thermal diffusion will not upset the results, so the defect of the current model lies in that κ is not a function of the temperature. The numerical methods can be extended to the temperature dependent case but in [27], the complexity of a model with temperature and altitude dependent coefficients is shown to require special numerical techniques in the family of level sets: the correlated-k methods.

The fundamental equations
Let I ν (x, ω) be the intensity of the radiation of frequency ν in the direction ω at point x of the physical domain . Let T (x) be the temperature. Energy balance (see [14,28]) yields, Here, S 2 is the unit sphere, ρ(x) is the density of the medium, κ ν is the mass-extinction coefficient, a ν is the single scattering albedo; 1 4π p(ω, ω ) is the probability that a ray in the direction ω scatters in the direction ω , (recall that 1 4π S 2 p(ω, ω )dω = 1); both κ ν and a ν usually depend on ν, and even x. The constant κ T is the thermal diffusion. Note that a ν and ρκ are dimensionless; the later is the absorption coefficient, i.e. the percentage of light absorbed per unit length.
As usual, boundary conditions have to be given. For the temperature we may prescribe its normal derivative to be zero for all x ∈ ∂ . The equation for the intensity being a first order equation, I ν should be given on − defined as where n is the outer unit normal of ∂ . However, we will deal also with cases which use on − some of the information arriving on + = {(x, ω) ∈ ∂ × S 2 : n(x) · ω > 0}.

Vanishing thermal diffusion
Proof It is shown by averaging (2) on S 2 .

Corollary 1
The temperature equation which is normally written with a flux of radiative energy (3), can be recast as (6): Corollary 2 If the thermal diffusion κ T is neglected in (6), then Remark 1 When κ ν and a ν are constant, (7) leads to the Stefan-Boltzmann law Note that the standard definition of the Stefan-Boltzmann constant has an extra π.
Some proofs concerning the existence, uniqueness and stability for solutions of simplified versions of this problem appear below. The most general case is discussed in the conclusion with relevant and updated references.

One dimensional approximations
Proposition 2 Consider (2), (7) in a vertical slab = (0, H )×R 2 . Assume that the boundary conditions at x = (x, y, z) are independent of y, z, and assume isotropic scattering ( p independent of ω and ω )). Let n be the outer unit normal at x = H . Then, the solution depends only on x and μ = cos φ = ω · n and I ν (x, ω) = I ν (x, μ) and T (x) are given by (1) and O ω r y x Sunlight P Atmosphere Earth θ Fig. 3 The Sun, in the far right, sends sunlight to point P on the Earth surface. A cross section of Earth and its atmosphere is shown in the plane defined by the axis Ox and the point P. The projection of the Earth rotation axis in that plane is shown (bold line) but plays no role. As the Earth radius R is large compared with the atmosphere thickness H we focus on a rectangle tangent to the Earth surface at P. In the end the radiative transfer equations are set on the line (P, r ), function of the angle ω where the observer observes the radiation intensity Proof Assume that B ν is a given function of x only. Both problems (2) and (9) have one and only one solution. Let us show that I ν (x, y, z, ω 1 , ω 2 , ω 3 ) = I ν (x, cos φ) is a solution of (2) when I ν is a solution of (9). Let t be the direction of the projection ω t of ω on the plane P of the slab boundary (Fig. 3).
Once I ν (x, μ) known, then B ν becomes a function of x only by (7).

Application to the Earth-Sun problem
Consider Fig. 3 and apply the invariance of Proposition 2 to the rectangle tangent to Earth at point P on its surface. The rectangle has width H , the thickness of the atmosphere and length L small compared to the Earth radius R.
Accordingly I ν depends only on the radial distance r to P, r ∈ (0, H ) and on μ, the cosine angle of the ray from Pr to the observer. So I ν (r , μ) is studied for r ∈ (0, H ) and μ ∈ (−1, 1).
The first condition applies only when μ > 0 because μ < 0 corresponds to the backside of the tangent plane.
The second condition says that at the top of the atmosphere no ray comes back into the atmosphere. Now I ν is uniquely defined by (9), (10), (11). As the equations are linear, I ν is proportional to cos θ . Hence to compute T at all points of planet Earth, one needs only to compute it at the point of intersection of the sphere and the Sun-Earth line and then multiply by cos θ . Reality is definitely more complex because this theory implies, in particular, that the Earth temperature at night is zero Kelvin! Note In this article we focus on methods rather than numbers; for clarity we neglect the scattering, i.e. a ν = 0, however, much of what is derived below applies also when isotropic scattering is present.

Spherical symmetry
Chandrasekhar showed in [8] that the one dimensionality argument of Proposition 2, using a tangent plane to Earth, can be extended to planets by using an osculatory spherical cap instead of a tangent plane, to take into account the radius of the planet.
For clarity, consider a spherical planet receiving parallel light rays from infinity. The planet's radius is R and its atmosphere thickness is H . Let the radial distance to the surface be r = |x|−R. Invariance with respect to the azimuthal and latitude angles, after suitable scalings by the appropriate cosines, and a similar argument as above, lead to the Chandrasekhar correction (see [8]): (12) and if thermal diffusion is neglected: Note that no additional boundary condition to (11) is needed because 1 − μ 2 is zero at μ = ±1.

Dimensionless variables
These "Chandrasekhar equations" can be de-dimentionalized by introducing a length scale λ, scaling factors for B, ν and ρ and set: r =r λ, R =Rλ and H =H λ, ρ = ρ 0ρ , ν = ν 0ν and B ν (T ) = B 0Bν (T ). Then we may rewrite the above and its boundary conditions as : withκν

The multi-group problem
In numerical computations one replaces the continuous map ν → I ν by a finite set of frequencies {ν k , I k } K 1 , and writes the system Recall that T , which is a function of τ only, couples all the {I k }. This formulation will be used in the numerical Sect. 7.

The Milne problem
When γ is neglected in (18), the problem is known as Milne's problem in = (0, Z ) × (−1, 1): When κ = κ ν ρ 0 , the integral with respect to ν of the solution of (16) is I times Q the integral ofQ ν with respect to ν.

More about the Milne problem
Emphasis on the Milne problem is motivated by the following facts: • it corresponds to a local in space description of the atmosphere say of height Z because for R large compared to Z the Chandrasekhar correction can be neglected; • one can introduce the point of view of functional analysis (cf. [11] chapter 21 Vol 9) without going into details but keeping things as explicit as possible; • one can also use very explicit computations, derived at a time when computers were not available, before the 1950.
We consider the abstract problem (22) .
the problem has a unique solution I ∈ L 2 ( ) which satisfies the estimate: Moreover when the data f , g 0 , g Z are non negative , the same is true for I , the solution of (22). With f = 0 one has: Proof The leading ideas are given below while details can be found in [11]. They are all based on the estimate of physical quantities which are translated into "mathematical" norms. When unambiguous, the symbol . will be used below to denote, the L 2 norm in (L 2 ( ), L 2 μ (−1, 1), L 2 μ (0, 1) and L 2 μ (−1, 0) ; the subscript μ indicates the variable of integration. Observe that the formula, gives the decomposition of I ∈ L 2 μ (−1, 1) in its orthogonal projection on the space of μ-independent functions and on its orthogonal (i.e. function of 0 μ-average).
One introduces for ≥ 0 the regularized equation: A priori estimate and uniqueness Let us multiply this equation by I and integrate with respect to τ and μ : Notice that Hence using the Cauchy Schwarz inequality: Since the problem is linear, denoting by I 1 − I 2 the difference of two solutions with the same boundary data g 0 and g Z and same external density f , one deduces from (29) the uniqueness because To extend this observation to the case = 0 one observes that in such case (30) gives: which gives the relation which implies I 1 = I 2 .
Existence of solutions for > 0 For clarity the proof of the estimate (29) and of the existence of a solution I are done in the absence of boundary source. Then, using the linearity of the problem it can be easily adapted to more general situations. First with > 0 for (29) one obtains a trivial stability estimate To prove the existence of the solution one considers with > 0 the Milne problem (25) in an iterative form, leading to the estimate which shows that the mapping I n → I n+1 is a strict contraction.
Then the same type of proof works also for the case f = 0 with non zero incoming data on − with the estimate: and the solution of the general problem with > 0 non zero, f and non zero (g 0 and g Z ), follows by linearity. The above construction will be used to prove convergence of the numerical method in the second part of the paper.
Existence of a solution for = 0 To let → 0, one proceeds with the following contradiction argument. If there would be no finite constant C for which holds the relation: that would imply the existence of a family of functions f of L 2 ( ) with norm equal to 1 while the corresponding solution of I would go to infinity in the same norm. Then it generates a solution to the problem: NowĨ converge weakly to a limit solution of the Milne problem with zero data, hence to 0 by the uniqueness of the solution. To complete (by contradiction) the proof one has to show the strong convergence ofĨ which is of norm 1. This follows from the so called averaging lemma (cf [11,16]) using the estimate.
Proof of the non negativity Positivity can be shown by the following standard intuitive arguments.
Denote by (τ + , μ + ) (resp. (τ − , μ − )) the point where I achieves its maximum (resp minimum). Then whenever the maximum (resp. minimum) is reached inside the open set (0, Z ) × (−1, 1) one has: And if it is reached on the boundary (τ + , μ + ) ∈ + and/or (τ − , μ − ) ∈ − one has: Consequently, the equation: implies that if the data f , g 0 and g Z are non negative and if the minimum is reached inside the domain it cannot be negative and if reached on + by (41) it cannot be negative either.
The only remaining case is the situation where this minimum is reached on − but then it coincides with g 0 or g Z which both are non negative. Hence in all cases one has In the same way for the solutions of the problem a positive maximum cannot be reached inside the domain because with (40) it should be negative which contradicts (45) , and it cannot be reached on + by the same argument since I coincides with g 0 (μ) or g Z (μ) . Since the above properties are independent of the proof of the positivity and of the estimate (23b) follow by letting → 0.

Remark 3
For the above problem L ∞ estimates and positivity for I or I are even more natural than L 2 estimates. However, to produce a complete mathematical proof one proceeds as follow: 1. Observe that the above estimates are fully valid for C 1 solutions. 2. Use the fact that smooth solutions with smooth data are dense in the set of solutions and that the above estimates remain valid under weak limit.
This approach is documented with details and used in [1].
To conclude this section it is convenient to recall the implicit formula for the solution of the problem: The solution of (22) with (21) satisfies This formula under different variants will be used below.

Milne Problem and "non explicit formula" for the temperature in terms of the albedo of the Earth
The fraction of the incoming solar energy scattered by Earth back to space is referred to as the planetary albedo and is an essential component of the Earth energy balance; cf. for instance [30]. In particular it can be combined with Milne problem to determine the temperature of the Earth as described below. Hence in (0, Z ) × (−1, 1) one considers an intensity of radiation which evolves according to the equation: Then on the upper atmosphere τ = Z an incoming boundary condition is given, for instance: with I ∞ representing the intensity of the radiation coming from the Sun. On the surface of the Earth, i.e. for τ = 0, the amount of radiation scattered back in the atmosphere I (0, μ)| μ>0 depends on the incoming radiation I (0, μ)| μ<0 . Therefore one assumes that it is given by the albedo operator A: Such operator may depend on many parameters in a very complex fashion (cf. the discussion in Sect. 2.2 of [30]). However, in the present setting of the Milne problem we assume that A -which is a data of the problem-is a linear contraction operator in the |μ| 1 2 weighted Sobolev spaces: Then one has the following: with the incoming data and the albedo data 1)) . 2. According to the Stefan Boltzmann law the temperature on the Earth is given by the formula: with C(Z , A) denoting a constant depending on the depth of the atmosphere Z and the albedo operator A .
Proof Existence follows by the same regularization as above as a consequence of the linearity of the boundary value problem and of the albedo operator. Uniqueness can be show by studying I = I 1 − I 2 where I 1 , I 2 are two distinct solutions. The operator I → I − 1 2 1 −1 I (τ, μ )dμ is the L 2 (−1, 1)-projection operator on the subspace of zero average functions and However the albedo operator is a contraction, so Consequently (56) implies that I = 0.
Remark 4 • Observe that the above analysis can be applied with almost no modification to the case where the hypothesis (52) for the incoming radiation is replaced by However, C(Z , A) is replaced by a coefficient C(Z , φ, A) which may depend on φ in a less explicit and more subtle way. • The case where no radiation is reemitted (in other world where all the radiation is absorbed by the Earth) fits simply in the above discussion with A = 0 . • The case where the Earth acts like a mirror reemitting all the incoming radiation fits also simply in the above discussion with • To describe a situation where a certain fraction α (Maxwell accommodation coefficient) of the radiation is reemitted while the rest is homogenized, maintaining the total intensity equal to 0 , one introduces the albedo operator.
which satisfies the hypothesis of the theorem 2 and leads to a constant C(Z , α) . In particular this accommodation coefficient may depend on the Earth temperature T and that would lead, for the atmosphere temperature to an even more implicit equation of the form (54) with a temperature dependent operator A(T ).

The half-space Milne problem
Let us study the case Z = +∞, the so-called half-space Milne problem. For the Milne equation (47) As a consequence to remain uniformly bounded with respect to τ for Z → ∞ any solution of (63) has to have I = 0. This justifies the following (cf.( [6] and [5]) Theorem 3 For any incoming data g 0 (μ) defined for τ = 0 and μ ∈ (0, 1) with g 0 (μ) ∈ L 2 (μ 1 2 , (0, 1)), there exists a unique uniformly bounded in τ solution of the half space Milne problem: This solution has zero flux and satisfies the estimates: This solution converges exponentially fast to a constant C(g 0 ); moreover the mapping g → C(g) is linear continuous from L 2 (|μ| Proof Once again the proof is only sketched below; for details see [5]. First one considers the solution on a double domain (0, 2Z )×(−1, 1) with incoming boundary data I (2Z , μ)| μ<0 = g 0 (−μ). This makes the solution of (47) unique, well defined and symmetric with respect to Z × (−1, 1) in the sense Then the decomposition of I into its average I a (τ ) = 1 2 1 −1 I (τ, μ )dμ and the orthogonal complement I or t = I − 1 2 1 −1 I (τ, μ )dμ gives, with the 0-flux property, the relation: Thus, one obtains the estimate: With (47) one can show that it gives the exponential convergence to a constant C(g) for Z → ∞.
The uniqueness of the solution is based on the same type of estimates.

Remark 5
The determination of g → C(g) and in the quest for an explicit or semi explicit formula has been in the last century the object of intensive activities involving in particular the Wiener-Hopf calculus (cf. [7]). However, the most explicit form is based on the introduction of the Chandrasekhar function H , defined by the implicit formula: which gives the constant C(g) by the relation In particular for g(μ) = μ one has ω 1 := C(μ) = 0.7014.

Approximate determination of the temperature on Earth
We present an approximation which yields a temperature on Earth based on the asymptotic behaviour of the half-space Milne problem.

Using Theorem 3
We return to climate dynamics where r ∈ (0, H ) is the altitude but not in an evanescent atmosphere, i.e. with (14) without the Chandrasekhar correction. We make the following change of variable and assume that r → ρ(r ) does not decrease too fast so that y(+∞) = +∞. Then we focus on the case H >> 1. If the atmosphere is grey, one observes that I (y, μ) = y − μ is solution of the Milne equation (63) for y ∈ R + , μ ∈ (−1, 1). Assume constant flux given bỹ Hence one introduces the solution e(y, μ) of the half space problem with 0 flux and equal to μ at y = 0 for μ > 0. As was proved in theorem 3 such solution exists, is unique, and converges exponentially fast to a constant ω 1 as y goes to ∞. As such, for some small function rem, provides a boundary layer approximation (i.e. for y > 0, y << 1) of the solution of the Milne problem with a flux given equal to 2 3 c and no incoming radiation for y = 0 and rem(y, μ) going exponentially fast to ω 1 when y → ∞. Hence for y large enough one has Now let us use the linearity of the solution with respect to the data and consider the same problem with incoming intensity without the term o(e −α Z ). As this is a small perturbation, we expect the solution at y = 0 to be Consequently, For a solar flux equal to the intensity I Earth is obtained after multiplication by 3 2 (see 70). This gives: and with the Stefan-Boltzmann law one obtains: Formula (76) with = Q(1 − a)/4 as in [14] eq. (2.2), p.66, gives T = 282 K.

Numerical analysis
Consider (16), dedimentionalized and without the Chandrasekhar correction. Given the physical constants of Table 1 and following (14), (15) (20) with g Z = 0 and T be given by Without dimensional scaling, with κ ν = κ constant, the total light intensity, is given by This comes from an integration in ν of withȊ (Z , μ)| μ<0 = 0. Indeed, an integration with respect to ν yields Hence we must have De-dimensionalization requires to sets I ν =Ȋ ν /B 0 , so the system becomes (79):

Numerical scheme
For (79), two numerical schemes are used. Both need the following fixed-point iterations: The first, referred below as "implicit", is based on a discretization of the exact solution of the first equation: Programming is straightforward; it needs only to be evaluated at all vertices of a triangulation of the rectangle . The integrals are approximated by a second order quadrature (trapezoidal rule) on a non uniform discretization of (0, Z ) to account for the fast variations in the tiny interval where infrared radiations occur. Updating T with the second Eq. of (80) is hard in general but simple in the grey case, κ constant, because the left side is κπ 4 T n+1 4 /15. In the grey case, Table 2 shows the error versus the mesh size for Ex. 6.2.1, below. Note however that the precision is weak: O(h), probably because the integrands are singular at μ = 0, τ = 0.
The second method is based on a finite element discretization of the PDE as in [19]. It uses a weak formulation of (80) discretized in V h , the space of Lagrangian-P 1 triangular elements. For stability a least square upwinding term (SUPG) is added, namely h SUPG (μ∂ τ I + κ I )(μ∂ τÎ +κÎ ) for a small h SUPG , whereÎ is the test function of the variational formulation. This means that at each iteration n of a fixed-point loop one must solve with I n+1 ∈ V h satisfying the boundary conditions and for allÎ ∈ V h withÎ (0)| μ>0 = I (Z )| μ<0 = 0 . The method has been implemented using FreeFem++ [21], which uses the library UMFPACK [12] to solve the linear systems. We found that the method works best when the triangulation is build first in the physical variables r , φ and then mapped to the rectangle of τ, μ. The automatic mesh refinement of FreeFem++, which is based on the Hessian of I n here, is also convenient to improve precision.

Example
The performance of both methods are reported in Table 2. The Finite element method (82) appears to be less precise than (81), but much faster. Adjustment of h SUPB is done once and for all proportionally to the number of points on ∂ .

Convergence of the iterative scheme
For clarity let κ = 1. Scheme (80) is modified slightly with a parameter As observed in [1,2] solutions of scalar kinetic equation with 0 incoming data are 0 viscosity limit of elliptic equations; therefore it is natural to introduce such regularization.

Results
In practice the convergence is much faster than predicted above, even with = 0, as shown by Table 3. A typical result is also shown, in the physical coordinates, on Fig. 4 for (78) with Q = 1. It shows the solution of the Milne problem computed on a grid 40 × 20. The temperature on Earth, given by (77), is T = 298 K.

The Chandrasekhar equation
For a grey atmosphere the model with the Chandrasekhar correction reduces to and (11). The two schemes above are easily modified to accommodate γ . For the FEM scheme one just adds to (82) γ (τ, μ)Î ∂ y I under the left integral of (82) plus an upwinding term like (85) below.
This scheme is consistent because I n+ 1 2 satisfies μ κ ∂ τ I n+ 1 2 + I n+ 1 2 = B n and adding this to the last equation above gives I n+1 + γ κ ∂ μ I n+1 + μ κ ∂ τ I n+ 1 2 = B n . In practice some additional artificial viscosity of amplitude δ should be added on the left in (84) When the above is discretized by a P 1 Finite Element Method, the convergence of the fixedpoint algorithm is equally fast; results are shown on Fig. 5 and illustrate the convergence of the solution of Chandrasekhar equations to the solution of the Milne equation when R increases.

Numerical simulation of the greenhouse effect
So much for the grey case. Our aim here is to compare the Earth surface temperature for two different functions ν → κ i ν (ν), i = 1, 2 and observe the relative change of temperature. The problem is defined in (79); the Chandrasekhar correction is not needed because R >> H . The values of κ ν are as follows.
The atmosphere is fairly, but not fully, opaque except in a region (ν 1 , ν 2 ) = (0.2, 0.3) where it is much less opaque. If a change in GHG proportion makes the atmosphere more opaque in this range then we may set We chose κ 0 = 1.225 because of the numerical value of the density of air (see (15)). We chose arbitrarily δκ = 0.5. The values for ν 1 and ν 2 are on the left side of the Boltzmann curve for Earth, shown in Fig. 6. In one computation the infrared clear window is narrow: The computer programs produce 4 temperatures τ → T j (τ ), j = 0, . . . , 3. According to Figs. 1 and 2, Green House gases increase κ and makes the window which is transparent to infrared radiations narrower. This will be measured by T 2 − T 1 and T 1 − T 3 .
The problem needs to be discretized in ν by choosing a grid in (ν m , ν M ), which means that we need to solve the multi-group problem introduced above with (16).
The following numerical scheme is used: .
Solve μ∂ τ I n+1 + κ ν I n+1 = κ ν B n ν (τ ),  (87) Finding T n by inverting the Planck function can be challenging. However, when κ ν = κ − δκ 1 (ν 1 ,ν 2 ) finding T n , solution of the first equation in (87) can be done by a fixed-point k-loop as follows: To assert the method we ask algorithm (87) to recover the solution of the Milne problem (κ ν constant). The results are shown on Fig. 7: a precision of 1% is obtained, but refining the mesh and the integration intervals did not improve the precision. This riddle will be solved in Sect. 7.2 Then we computed the relative change of temperature when κ ν is changed from κ 1 ν to κ 2 ν . The change is of the order of 10 −2 and negative (see Fig. 8). In view of the small magnitude of the change, for a verification we turned to a calculus of variations.

Solution by calculus of variations
Even though both the FEM-based code and the implicit one give the same results, evidently we are trying to observe a temperature variation which is at the limit of the precision of the computer codes (see remark 8 below). So let us try another method.
Let (86) be written as κ ν = κ + δκ ν with δk ν = −δk 1 ν∈(ν 1 ,ν 2 ) and let κ be constant. With the obvious notations of a calculus of variations: Let δĪ = ∞ 0 δ I ν dν and similarly for δB. Integrating the equations with respect to ν leads to Adding both gives and, knowing that δ B = δ( π 4 T 4 15 ) = 4 π 4 T 3 15 δT , the change in temperature is computed by (91) divided by κ: The numerical solution of (92) can be obtained both with FEM and the implicit method; results for δT agree roughly only as shown on Fig. 8. Yet here too we obtain a decrease of Earth temperature when κ ν increases in the infrared interval (ν 1 , ν 2 ) and the numerical values of the change obtained are of the same magnitude as those obtained by a direct simulation with κ 1 ν and κ 2 ν , as seen on Fig. 8.

A Formulation to compute only the temperature
These unexpected conclusions forced us to think again and track all potential precision losses. In doing so we turned to exponential integrals (see en.wikipedia.org/wiki/ Exponential_integral) to evaluate (81), an idea on which is based the computation of an analytic solution of the Milne problem in [14]-Appendix B. Function E 1 is not hard to program; it is also part of the Gnu Scientific Library gsl. One has: Observe that (81) can be integrated in μ so as to write everything in terms of We note also that to compute τ → T (τ ) when amounts to solve These give the following fixed-point iterative scheme to solve (97): .
. Then compute F n+1 by This turns out to be a very fast and easily programmable method, and so far, the least prone to precision difficulties because the only singular integrand is E 1 (|t|)| t∼0 ∼ − log(|t|); but as it appears under an integral, a dedicated quadrature rule can be used with t log t dt = t log t − t, which is not singular in R + .

Remark 6
Note that the scheme is a non-linear solver for the functional integral equation for τ → T (τ ): Discrete Fourier Transform could be used to convert the integral equations into a linear system for κ ν B ν , (still nonlinear in T ) but the fixed-point iteration process is so fast that it is not worth it.

Remark 7
Note that the method extends to constant isotropic scattering a . Consider It leads to a formulation involving a frequency dependent integral, G ν (τ ) := 1 −1 I ν (τ, μ)dμ: The method was tested with a = 0.3. It gives a temperature at ground level 10% higher and 10% less cooling effect due to the same changes in κ ν . It requires 15 iterations,instead of 10, to converge to 3 digits accuracy.

A new set of tests
All numerical tests where re-run using (98), (automatic differentiation included), and the same results were obtained (Figs. 9 and 10): at all altitudes, temperatures decrease by 1 or 2 percents with a narrowed infrared absorption interval, or when κ ν is multiplied by 2 in an infrared interval! An added set of frequencies were tested: ν ∈ (1.0, 1.2) and ν ∈ (1.0, 1.4).
These numbers correspond to absorption rays of the GHG in the lower frequency range of sunlight (see Figs. 1 and 6). With the new set (102), the temperatures increase when κ increases and/or when the partially transparent window decreases in size, yet the numbers are an order of magnitude smaller, about 0.17% near the ground level, i.e. 0.5C. We note also (Fig. 9) that we can obtain agreement to at least 3 digits between a direct solution of the Milne problem with constant κ and the same solved by the multi-group formulation (98) even though κ is constant. Figures 8 and 10 are comparable, but notice Calculus of Variations and Automatic Differentiation could differ by a factor of 2 from the direct simulations of the temperature differences; this is probably because Calculus of Variations linearizes the problem while the finite differences don't. It could be also an indication of a precision issue due to a very sharp variation of the derivative of T with respect to κ.

On the convergence and stability of scheme (99)
Schemme (98) cannot jam and cannot explode either . Indeed, although E 1 has a log singularity at 0, its integral is bounded; B n ν (T ) is bounded too, therefore F n+1 in (98) is always bounded and positive.

Altitude -km
Relative change  Hence the scheme will always generate finite numbers. Although we may not concluded that it converges, we can conclude that if |T n+1 − T n | ∞ → 0, then any accumulation point of T n is a solution of (99), i.e. a solution to the radiative transfer equations (79), for some appropriate function I ν (τ, μ), not given by the formulation but computable after convergence by (81).
A typical convergence behavior is reported in Table 4. The precision at one point near the Earth surface is also reported. Three digits are guaranteed after 6 to 10 iterations and 200 discretization points for the altitude and an integration step for the integrals equal to 0.005.

Discussion on the reliability of the numerical results
Greenhouse gases leading to a cooling of the atmosphere is counter intuitive. To check the codes we made a similar change of κ ν but in a bigger range: κ = 1 + 0.5 1 (1,1.5) where the precision should not be a problem. Results are shown on Fig. 11: increasing κ in the range of sunlight leads to 0.5% increase on Earth temperature. The direct simulation of the change agrees with the calculus of variations. So the computer program is probably correct. The change of κ ν in a small frequency interval (0.2,0.3) could be beyond the numerical precision of the method and indeed Fig. 7 shows that the precisions may not be sufficient.
We make another remark relevant to precision:

Remark 8
It is because the intersection of the Boltzmann curve for black body Earth with the Boltzmann curve for black body Sun is non zero (see Fig. 6) that there is an infrared re-emission due to sunlight on Earth.
As a final check of the results of Figs. 8 and 10 we wrote an entirely different computer program to implement (87), in C++, linked to an Automatic Differentiation library: a technique based on operator overloading which gives the exact values of derivatives of any variable in the program with respect to another variable, here the value of δκ in the frequency range of GHG absorption. The program uses a uniform grid in τ, μ, by opposition to the FreeFem++ program which uses a fairly uniform grid in the physical domain refined and adapted during computation to the Hessian ofĪ . It also produces a decrease of Earth temperature of the same magnitude! All the numerical methods, implicit, finite element, uniform mesh (C++) and direct computation of the temperature difference or calculus of variations or automatic differentiation, give correct values for the temperatures versus altitude, but imprecise values to their derivative with respect to κ ν ; nevertheless all of them predict a decrease of temperature due to GHG when the absorption κ ν increases in an infrared range (ν 1 , ν 2 ) and when the range (ν 1 , ν 2 ) is decreased.
While the precisions of the finite element method and of the implicit method can be doubted because I has singularities at the poles and because of singular exponentials, there is no reason to doubt the precision of the last method based on (99) for the temperature. The temperature known, the light intensity can be recovered from (81) but it has singular integrals.

Boundary layer near the Earth surface
Consider the Chandrasekhar equations with thermal diffusion: ∀r , μ, η ∈ (0, H ) × (−1, 1) 2 , where η = cos θ and assume that κ T = κ 0 , << 1. Then it is likely that with T 1 (r , μ) << 1 when r → ∞. This leads to the following cascade of equations with r = r . For clarity and without losing generality we assume R is large so as to reduce the above to The last line is also −∂ 2 The conclusion is that there is no strong variation of the temperature r → T (r ) near the surface (r=0) due to thermal diffusion, but there is a strong variation of the gradient.
To connect with the next section we notice that (107) can be rewritten as:

Boundary layer and Robin boundary condition
The temperature is a solution of an elliptic equation which requires a boundary condition on the entire boundary ∂ while the boundary condition for I needs to be given only on the incoming part of − . Observe that (107) involves two temperatures T 0 (r ) which could be expressed in terms of I by the Stefan-Boltzmann law and a temperature T (r ) which represents the "observed temperature" near the boundary (which is unknown ) and determined in terms of non explicit constants. Such fact was already observed in nuclear reactor technology, where it leads for the diffusion approximation to a Robin boundary condition and is explained in [31] p. 199, eq. (8.13).
Below, following [5,11] we propose a self contained derivation of this type of formula based on scaling analysis. Moreover for the sake of simplicity we consider the solutions I of a dependent half space 0-flux (cf. Sect. 5) Milne problem; one has the following I (0) independent of μ, converges to the μ independent solution of the diffusion equation with the Dirichlet boundary data I (0) = I (0), with a rate of convergence O( √ ) in L 2 (R + × (−1, 1)). However, the expression provides an approximation of order in L ∞ (R + × (−1, 1)).
One observes that I is uniformly bounded in L ∞ (R + ×(−1, 1)) hence by standard estimates related to the diffusion approximation , it converges to a μ-independent function I 0 (r ) solution of (110) with I 0 (0) = I (0). Then one observes also that is a solution with an error of the order of √ of the equation (108). This construction can be iterated giving a solution of any finite order of this equation. However, at r = 0 and μ > 0 one has: and this estimation concerns a boundary layer of size √ which can be only analyzed by the use of the zero flux solution e(τ, μ) of the half space problem: Therefore one introduces the functions: Constructed in such a way, I c (r , μ) enjoys the following properties.
• It is a solution of (108) with a remainder of order .
• I rem is the sum of two terms √ ω 1 ∂ r I 0 (r ) and the boundary layer term: According to the theorem 3 one has: sup As a consequence of these observations one has In an informal way the following can be derived:

Corollary 3 Assume that the intensity of radiation I of the above half-space Milne problem is coupled with the solution of a diffusion equation at the boundary of the domain by the Stefan-Boltzmann law; then the introduction of a Robin boundary condition of the type
in the diffusion approximation will improve it by an order of √ to .
Proof Starting from the relations one deduces from (116) that From which

Conclusion
To summarize, we recall that radiative transfer is an old topic, studied by astronomers and nuclear scientists and more recently for Earth science. Much of the ancient material can be discarded in view of the more powerful computer solutions. However, it turns out that for the simulation of the effect of sunlight on the atmosphere, the problem is numerically difficult, so that any mathematical and analytical properties gathered in the past are welcome. Over the last 50 years the mathematical approach enjoyed stimulus from a huge range of applications, and the introduction of functional analysis and computing. However, one observes that there is still room for progress on the full model, in particular to make the hypotheses needed for proofs much more in agreement with the case considered in any kind of physics. As underlined above, the Eqs. (2), (3) lead to the following comments • Concerning the time dependent equations, for sufficiently regular coefficients (κ, ρ and regular initial and boundary data), as it is expected, the problem has a unique well defined (for a finite time) solution, which can be extended on [0, ∞) when the volumic sources f = 0 (cf. [29] for instance, for proofs and recent references). One of the main observation used in this contribution is the fact that an estimate of the type 0 < m(0) ≤ T (τ, 0) < M(0) remains valid for later time with 0 < m(t) < T < M(t) . • In ( [29]) independent boundary conditions are assumed for I and T . It may be more realistic to include in the description some relation on the boundary. This would make use of the boundary layer analysis briefly described in the Sect. 8. • A more serious difficulty comes from the non-constant opacity κ ν (τ, T ), possibly not regular (cf. [27,28]) and very often only vaguely known. At least two directions have been proposed to deal with this issue. In one of the first contributions on the subject (cf [25]) it was assumed that with no other hypothesis on the dependance with respect to the frequency ν that the function T → κ ν (T ) was non increasing , while the function T → κ ν (T )B ν (T ) was non decreasing. Then some L 1 stability estimates lead to existence and uniqueness for the system. On the other hand a popular grey model, which also leads to a Milne problem, is based on the assumption that the opacity depends on the temperature T , yet independent of the frequency ν. With such hypothesis several stability results have been obtained with no constraint on the regularity of the mapping T → κ(T ) (see [1,3,4,15]). • Under some convenient scaling hypothesis, in particular large opacity with respect to the size of the media one may approximate the dynamics by a diffusion equation known as the Rosseland approximation. Once again mathematical results are well advanced for the grey model and some of its variants and more sparse in the general case. Such approximation is very well adapted to describe "interior problems" like fusion by laser confinement. It does not seem (to our knowledge) present in climatology. As a matter of fact, the height of the atmosphere being small with respect to the Earth radius it is by itself a boundary layer, except for the fact that the equations are considered only in the vertical direction. As sketched in the Sect. 8.1 this issue is closely related to the improvement of the accuracy of the Rosseland approximation and also well developed for grey model. It is worth mentioning that mathematically the problem is easier with the Chandrasekhar correction (12), [20,32].
Numerically, it is a mixed integro-differential problem for which a fixed-point approach works quite well, and for which a convergence proof is available for grey atmospheres. Four methods of discretization have been tried. A finite element method with upwinding, an implicit method based on the integral form of the solution of the equation for the light intensity at given temperature, a finite difference implementation of the implicit method on a uniform mesh, with automatic differentiation, and finally an integral formulation for the temperature only. The second one is more precise but slower; the third is just for checking that the programs have no bugs, and the fourth one is the most trustworthy. So we may use it for validation . Convergence with respect to grid refinement is fairly fast.
The present computations validate a decrease of Earth temperature due to CO 2 and other greenhouse gases responsible for a substantial change in the transmission coefficient κ ν in the lower part of the infrared frequency range emitted by Earth seen as a black body at temperature ∼ 300 K. Narrowing the infrared transparent window has also the same cooling effect. On the other hand, it gives a heating effect when the same changes on κ ν occur in the lower frequency range of sunlight! Consequently, the radiative transfer equations, used with sunlight supposedly unaffected by the atmosphere and without scattering, should not be presented as an explanation of the greenhouse effect using the infrared frequency range. In [13,17] the greenhouse effect is explained by radiative transfer principles but assuming that κ and B ν depend on T and air pressure. It may require non-isotropic scattering in the atmosphere to justify numerically their results.
Finally, recall that radiative transfer alone gives a very crude model for the Earth temperature: 76 C at noon, -273 C at midnight and 19 C on the average in Paris, computed by assuming that sunlight power is Q = 620/2). Climate change is indeed much more complex than just radiative transfer! In a forthcoming publication a full treatment of nu-dependent scattering and earth albedo will be presented. These terms dramatically alter the conclusion reached here.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.