M Herreras-Giralda

P Litvinov

O Dubovik

Y Derimian

T Lapyonok

D Fuertes

O Sourdeval

R Preusker

J Fischer

Thermal emission in the Successive Orders of Scattering (SOS) radiative transfer approach

. In comparison with the widely used code DISORT (DIScrete-ORdinatemethod for Radiative Transfer) [3], the developed SOS scheme achieves a mean accuracy of radiance calculation of -0.005 K (-0.003%) expressed in terms of brightness temperature. Under the same vertically inhomogeneous atmospheric conditions, GRASP SOS RT is approximately eight times faster than DISORT. The analysis of the sensitivity of GRASP TIR SOS scheme to the number of layers and the effect of polarization are also investigated in the paper.

Introduction

Recent progress in remote sensing has demonstrated the benefits of extended information content in advanced observations for the characterization of different atmospheric constituents. In particular, essential improvements in dust characterization are related to the extension of measurements from visible to Near Infrared (NIR) and short-wave infrared (SWIR) [4] spectral range. Recent studies showed new possibilities of dust and volcanic aerosol characterization using thermal infrared (TIR) spectral bands, where strong resonances can be observed for the complex refractive index of different volcanic and dust chemical components [5], [6], [7], [8]. In this regard a combination of VIS, NIR, SWIR and TIR measurements should provide higher sensitivity to both size distribution and composition of aerosol particles.

GRASP (Generalized Retrieval of Atmosphere and Surface Properties) algorithm [2] has been successfully used in several aerosol retrieval applications from observations covering ultraviolet, visible spectrum and near-infrared spectra for ground-based and satellite measurements as, for example, in [9], [10], [11], [12] and [13]. Further adaptation of GRASP forward and inversion modules for TIR spectral range opens new possibilities for a wide variety of synergistic applications. For example, this adaptation is expected to enable the combined retrieval of space borne missions as 3MI [START_REF] Marbach | Schlüssel in Polarization science and remote sensing VII[END_REF], [START_REF] Fougnie | [END_REF] and IASI (Infrared Atmospheric Sounding Interferometer) [START_REF] Blumstein | Infrared Spaceborne Remote Sensing XII[END_REF], [START_REF] Hilton | [END_REF]. Similarly, in the case of ground-based applications, the combination of AERONET (AErosol RObotic NETwork) sun-sky-photometer [18] and CLIMAT TIR radiometer [19], [20] observations is also expected to provide an enhanced aerosol mineral dust retrieval with an extended composition inversion.

The analytical solution of the RT equation does not exist for a particulate medium representing all complexity of Earth's atmosphere. There are several efficient and accurate numerical techniques for solving the RT problem [1], [21], [START_REF] Hovenier | Transfer of polarized light in planetary atmospheres: basic concepts and practical methods[END_REF]. In this paper we consider the implementation of atmosphere and surface thermal emission into the SOS method for the RT standard problem: plane-parallel vertically inhomogeneous particulate layer atmosphere, which is representative for a vast variety of remote sensing applications. Although there are no principle limitations to account for the emission processes, the existing numerical implementations of the SOS method (e.g. [1], [START_REF] Kotchenova | [END_REF], [24]) do not include thermal emission sources in the radiative transfer scheme. Hitherto, the SOS RT at TIR domain is not widely used yet. This is related to the fact that the radiative transfer in thermal infrared spectral range deals with problems which are not common at the Solar range. For example, the refractive index of atmospheric aerosol particles can produce scattering resonance effects at TIR range, while their spectral dependence in UV, VIS, and NIR/SWIR range is rather smooth [25], [26]. The resonances result in an essential growth of the required number of expansion terms for accurate calculations of optical characteristics of aerosol [START_REF] Van De Hulst | Rept Inst. Space Studies NASA[END_REF], [START_REF] Bohren | Absorption and scattering of light by small particles[END_REF] and strongly affect the computation speed of SOS RT. Moreover, accounting for thermal emission together with gas absorption [START_REF] Rothman | [END_REF] and multiple scattering requires a reliable procedure of discretisation for temperature, gases and aerosol vertical profiles ensuring an accurate RT solution.

This paper describes the application of RT SOS method for atmospheres with thermal emission and tests it through the inter-comparison with another existent RT scheme accounting for the thermal emitted component. The series of case studies are performed to analyze different approaches for vertical discretisation of temperature and aerosol properties in atmospheres for a wide range of optical thickness. Other minor effects, such as, the influence of polarization in the long-wave radiative transfer calculations are also discussed.

Presently, there is a number of RT codes based on different techniques accounting for thermal radiation emission, scattering and absorption in the atmosphere: the Matrix Operator Model (MOMO) ( [30], [31]), RTTOV ( [START_REF] Saunders | A fast radiative transfer model for assimilation of satellite radiance observations-RTTOV-5[END_REF], [START_REF] Saunders | [END_REF]) FASDOM ( [34]), or DISORT (DIScrete-ORdinate-method for Radiative Transfer) [3]. DISORT has already been widely used by the scientific community to validate radiative transfer methodologies, as for example in [START_REF] Kotchenova | [END_REF] and [35]. In this work the C# implementation (cdisort-2.1.3) of DISORT technique is used to test the developed SOS RT in TIR domain. In addition to the original developments of [3], cdisort-2.1.3 also includes the δ-M method [36], the correction to the intensity field [37] and the solution for a general source term [38].

The paper is organised as follows: first the basic general radiative transfer equations are presented (Section 2). In Section 3 the implementation of the SOS in the thermal infrared spectral range is discussed. A comparison of the performance of different methodologies to unify the optical depth and the temperature vertical discretization as well as the details about the implementation of the SOS approach in GRASP code can be found in the Section 4. The Section 5 is devoted to the validation and testing of the developed methodology and SOS code against DISORT. Analysis of the polarization effects in TIR is included in the Section 6. Finally, the main conclusions are summarised in the Section 7.

The radiative transfer equation in TIR

For a plane-parallel layer of particulate sparse medium with an arbitrary optical thickness, the directions of the scattered and incident light can be represented by the vectors k and k 0 . In the right handed coordinate system with axis z perpendicular to the boundaries of the medium, these directions are described by azimuth and zenith angles (φ, θ) and (φ 0 , θ 0 ) for the scattered and incident directions correspondingly. For the sake of clarity, the zenith angle in this paper is always expressed in terms of its cosine (µ = cos(θ)). Thus, µ > 0 stands for upward direction, µ < 0 for downward direction. Therefore, µ 0 < 0 in the chosen coordinate system. The vertical dimension is expressed in terms of the optical depth τ . Thus, a value of τ = 0 corresponds to the top of the atmosphere (TOA) and a value of τ = τ * to the bottom of the atmosphere (BOA), where the boundary surface is located. In general terms, the monochromatic vector radiative transfer equation can be expressed as follows [39]:

µ dL(τ, µ, φ) dτ = L(τ, µ, φ) -S(τ, µ, φ). ( 1
)
Here the subscripts corresponding to the wavelength are omitted because the interactions of radiation with different frequencies are not going to be considered. L(τ, µ, φ) in Eq.(1) corresponds to the vector radiance in the direction n defined by Stokes parameters vector [40]:

L(τ, µ, φ) = (I, Q, U, V ) T .
(

) 2
where I represents the total intensity, Q, U and V are the Stokes parameters describing linear and circular polarization [1] [41] [START_REF] Mishchenko | Radiative transfer and coherent backscattering[END_REF]. S(τ, µ, φ) in Eq.( 1) is the vector source function. If scattering and emission is taken into account it can be generally defined as follows:

S(τ, µ, φ) = ω(τ ) 4π P(τ, µ, φ, µ 0 , φ 0 )E 0 exp( τ µ 0 ) + ω(τ ) 4π 2π 0 +1 -1 P(τ, µ, φ, µ ′ , φ ′ )L(τ, µ ′ , φ ′ )dµ ′ , dφ ′ + ϵ λ (τ )B λ (T (τ )), (3)
where ω represents Single Scattering Albedo (SSA) of scatterers in the medium, P(τ, µ, φ, µ, φ) is the phase matrix of the scatterers, and B λ is the Planck function describing thermal emission of particles with temperature T (τ ) and emissivity ϵ λ (τ ). In other words, the first term on the right hand side of the equation stands for the single scattering, the second one accounts for the multiple scattering processes and the third one describes the thermal emission. Since the natural Sun radiation is unpolarized, incident irradiance to the atmosphere can be expressed as: E 0 = (E 0 , 0, 0, 0) T .

Eq.( 1) can be written in the integral form normally used for numerical solutions of the RT problem:

L(τ, µ, φ) = - τ 0 e -(τ ′ -τ )/µ S(τ ′ , µ, φ)dτ ′ /µ ∀ µ < 0 (4a) L(τ, µ, φ) = L(τ * , µ, φ)e -(τ * -τ )/µ + τ * τ e -(τ ′ -τ )/µ S(τ ′ , µ, φ)dτ ′ /µ ∀ µ > 0 (4b)
The integro-differential nature of Eq.( 1) resulting in the inter-dependencies of S(τ, µ, φ) and L(τ, µ, φ) in Eqs.(4a) and (4b) complicates general analytical solutions of the RT problem, especially when both scattering and emission processes are accounted. In order to avoid the complex numerical solution of the double integrals over azimuth and zenith angles in RT equations, expansions in Fourier series of the azimuthally dependent characteristics are used. In this way, the RT equations can be written as a system of linear independent integral equations for Fourier components independent of azimuth angle. The Fourier decomposition technique to solve RT equations is well known and described in detail in a number of papers, for example see [1]. For completeness this description is also provided in the Appendix A.

The next Section briefly revisits the numerical approach used in SOS method to solve RT equations for scattering media and extends it to the case when TIR emission is taken into account.

The SOS RT method accounting for thermal emission

The interaction of light with the atmosphere includes scattering, absorption and emission processes with its own methodological differences. The SOS RT approach is based on the method of successive approximation implemented via consequent solution of linear equations. Specifically, the solution is represented by an infinite series of "orders of scattering (approximation)" which represents different approximation terms in the series:

L(τ, µ, φ) = N n=1 L n (τ * , µ, φ), (5)
where the term with n = 1 corresponds to the First order approximation term. If thermal emission is accounted, it includes both emission and single scattering terms. The terms n > 2 represents multiple scattering contributions.

First order approximation

The usual way to account for thermal emission in the medium is based on the assumption that the particles can be modelled as a grey body in Local Thermodynamic Equilibrium (LTE) ( [3][31]). According to Kirchhoff's law [START_REF] Kirchhoff | Von Kirchhoff bis Planck[END_REF] particles emissivity at a certain wavelength corresponds directly to their absorptivity at the same wavelength. Therefore, the emission term of source function in Eq.(3) at each optical depth τ can be expressed via single scattering albedo (ω) as the following:

S 1,emis (τ ) = (1 -ω(τ ))B λ (T (τ )). (6)
For a horizontally isotropic atmosphere, the emission does not depend on azimuth angle φ.

The total source function in the first order approximation is a sum of the emission and single scattering parts expressed as follows:

S 1 (τ, µ, φ) = S 1,emis (τ ) + S 1,scat (τ, µ, φ), (7)
where

S 1,scat (τ, µ, φ) = 1 4π ω(τ )P(τ, µ, φ, µ 0 , φ 0 )E 0 e τ /µ0 . (8)
Thus, from the Eqs.(4a) and (4b) the scattering and emission parts of the radiation in the first order approximation can be calculated independently for the upward and downward directions:

L 1,i (τ, µ, φ) = - τ 0 e -(τ ′ -τ )/µ S 1,i (τ ′ , µ, φ)dτ ′ /µ ∀ µ < 0 (9a) L 1,i (τ, µ, φ) = L 1,i (τ * , µ > 0, φ)e -(τ * -τ ′ )/µ + τ * τ e -(τ ′ -τ )/µ S 1,i (τ ′ , µ, φ)dτ ′ /µ ∀ µ > 0, (9b)
where i stands for 'emis' or 'scat'. At the boundary level, τ = τ * , the upwelling (µ > 0) radiance in the first order approximation is defined by the surface scattering and emission:

L 1,scat (τ * , µ, φ) = - µ 0 π R(µ, φ, µ 0 , φ 0 )E 0 e τ * /µ0 ∀ µ > 0 (10a) L 1,emis (τ * , µ, φ) = ϵ sur B λ (T sur ). ∀ µ > 0 (10b)
Here R is the surface reflection matrix [1][44] [45] and ϵ sur is the surface emissivity at a given wavelength. Similarly to the consideration of the atmospheric layers, the surface emission is modelled as the emission of a grey body in LTE with a temperature T sur . If the surface is optically thick, the transmitted radiation through the surface is negligible in comparison to the scattering and absorption. In these conditions, the surface emissivity can be directly related to the surface albedo a sur :

ϵ sur = 1 -a sur . (11)
The total radiance in the first order approximation:

L 1 (τ ) = L 1,emis (τ ) + L 1,scat (τ ). ( 12
)

Multiple Scattering

In Eq.( 5) the terms with n > 1 describe the multiple scattering of the radiation in the medium. The total source function for the multiple scattering terms is expressed as the following:

S n (τ, µ, φ) = ω(τ ) 4π 2π 0 +1 -1 P(τ, µ, φ, µ ′ , φ ′ )L n-1 (τ, µ ′ , φ ′ )dµ ′ , dφ ′ ∀ n > 1 ( 13
)
The corresponding upward and downward radiation can be calculated following once again equations (4a) and (4b):

L n (τ, µ, φ) = - τ 0 e -(τ ′ -τ )/µ S n (τ ′ , µ, φ)dτ ′ /µ ∀ µ < 0 (14a) L n (τ, µ, φ) = L n (τ * , µ, φ)e -(τ * -τ )/µ + τ * τ e -(τ ′ -τ )/µ S n (τ ′ , µ, φ)dτ ′ /µ ∀ µ > 0 (14b)
At the bottom border of the medium, the radiance for the multiple scattering term can be written as follows:

L n (τ * , µ, φ) = 2π 0 0 -1 (-µ ′ )R(τ, µ, φ, µ ′ , φ ′ )L n-1 (τ * , µ ′ , φ ′ )dµ ′ dφ ′ /π ∀ n > 1 , µ > 0 (15)
A simple fact can be noted from Eqs.( 6)-( 15): accounting for the thermal emission in RT affects only the equation for the source function in the first order of approximation, whereas the equations for higher orders remain unchanged. Nevertheless, the accurate accounting for thermal emission, scattering and absorption in SOS solution faces some principle challenges in numerical implementation, such as the identification of a reliable procedure of discretisation for temperature, gases and aerosol vertical profiles.

4 Vertical discretization: scattering, absorption and temperature profiles

One of the most important factors for an accurate numerical solving of RT equations is the definition of an adequate vertical discretization of the atmosphere to account for the vertical variability of atmospheric scattering, absorption and emission properties. It is especially crucial for the SOS RT method, where the integral over source function should be calculated correctly for each discrete layer l (with optical thickness ∆τ l ) in both first order and multiple scattering approximations. Indeed, the atmosphere is divided into optically thin layers ranging from 1 to L max (which corresponds to levels ranging from 0 to L max ), where the optical properties are averaged. Thus, the radiance corresponding to an arbitrary level τ = τ L is expressed as:

L(τ L , µ, φ) = - L l=1 e -(τ l -τ L )/µ ∆τ l 0 e (τ ′ l -∆τ l )/µ S(τ l-1 + τ ′ l , µ, φ)dτ ′ l /µ ∀ µ < 0 (16a) L(τ L , µ, φ) = L(τ * , µ > 0, φ)e -(τ * -τ L )/µ + Lmax l=L+1 e -(τ l-1 -τ L )/µ ∆τ l 0 e -τ ′ l /µ S(τ l-1 + τ ′ l , µ, φ)dτ ′ l /µ, ∀ µ > 0 (16b)
where τ * = τ Lmax , ∆τ l = τ l -τ l-1 and S n (τ l-1 + τ ′ l , µ, φ) is the source function defined in the in the layer l:

S(τ l-1 + τ ′ l , µ, φ) = 1 4π 2π 0 +1 -1 ω(τ l-1 + τ ′ l )P(τ l-1 + τ ′ l , µ, φ, µ ′ , φ ′ )L(τ l-1 + τ ′ l , µ ′ , φ ′ )dµ ′ , dφ ′ (17)
In addition to the scattering and absorption, accounting for the thermal emission in the atmosphere also requires a correct representation of the temperature profile in each layer l (Eq.( 5)-( 17)). In general, a correct vertical discretisation of the atmosphere with scattering, absorption and thermal emission can be done by increasing the number of atmospheric vertical layers. In practice, this may essentially increase the computation time and make the SOS method less efficient for a number of applications. Therefore, this section is devoted to the analysis of optimal discretization of the temperature, as well as, the scattering and absorption profiles.

First of all, the current section presents a description of how the atmospheric layers are defined in the GRASP SOS RT code for optimum accounting for vertical variations of scattering properties in the presence of several optically distinct components. Then, the effect of vertical discretization over scattering and emission processes are independently analysed in detail.

Atmospheric Layer definition

Finding a numerical solution for the SOS RT equations at a certain location of the observations requires the subdivision of the integrals over optical thickness Eqs.(4a) and (4b) into optically thin layers Eqs.(16a) and (16b). The scattering and absorption properties within each of such layers can be represented by the layer average phase matrix and single scattering albedo. In case of presence of several scattering and absorbing atmospheric components these quantities can be calculated as follows:

ω(∆τ l ) = Nc i=1 ω i (∆τ l )∆τ (i) l Nc i=1 ∆τ (i) l , (18)
ω(∆τ l )P(∆τ l ) = Nc i=1 ω i (∆τ l )P i (∆τ l )∆τ

(i) l Nc i=1 ∆τ (i) l , ( 19
)
where N c is the total number of components in each layer (aerosol, gases, molecules or clouds).

The scattering source functions for each layer l can be written as follows:

S 1,scat (∆τ l , µ, φ) = 1 4π ω(∆τ l )P(∆τ l )E 0 e τ /µ0 (20)
As previously mentioned, the layers in the approach described here are considered to be in LTE conditions which requires all atmospheric components within the same layer to be at the same average temperature. Thus, each layer is represented just by one single Planck function. However, temperature variation through different layers is allowed. Consequently, the emission source term in the layer l and in the first order approximation represented in Eq.( 6) is expressed in terms of ω(∆τ l ) as the following:

S 1,emis (∆τ l ) = (1 -ω(∆τ l ))B λ (T (∆τ l )). ( 21
)
The representation of T (∆τ l ) in each layer l is described in the next sections.

GRASP SOS RT in high TOD conditions

Hitherto, the performance of GRASP SOS radiative transfer has only been tested for moderate TOD values normally existing for visible and SWIR spectral range applications. If strong gas absorption lines are accounted for, then the conditions can be far from those in visible and SWIR domains. For example, TOD can reach more than one order of magnitude higher. If a high enough number of layers is selected, the SOS RT approach can be successfully used under any TOD condition. However, as already mentioned the calculations with a high number of layers can be computationally expensive. The thermal part, in general, does not present special problems once that the discretization of temperature vertical variability has been properly resolved. In this section we therefore focus on the calculation of the scattering processes solely.

In order to explore the limits of GRASP SOS RT accuracy under high TOD conditions, a reference simulation of 600 layers is going to be compared with simulations performed with a smaller number of layers. In order to avoid the effects of temperature vertical discretization, four different scenarios with a high TOD based on an isothermal atmosphere at 250 K (including the surface) have been designed.

The single scattering albedo has been set to a value around 0.5 in all layers. This high single scattering albedo condition may not be very realistic in comparison with the zero value of this magnitude of pure gas absorption lines. However, in this analysis we want to assure that the scattering part plays a significant role in the calculations. The extreme characteristics of the atmospheric scenarios designed for this comparison mark a good reference for the accuracy limit of GRASP SOS RT to resolve scattering processes, because the real atmospheric situations present more favourable conditions to perform these calculations.

Figure 3 shows the average radiance difference between the reference simulation using 600 layers and the calculation with a smaller number of layers calculated by GRASP SOS RT. This average radiance difference includes a complete set of upward and downward geometries including the full azimuth range and a zenith range from 10 • to 70 • (it can be visualized in figure 5). Hereafter, the radiance results presented in this study will be expressed in terms of Brightness Temperature (BT). In the case of ground-based observations, there is no significant change in the final radiance even for a reduced number of layers. On the contrary, in the case of calculations for satellite geometries, the convergence is reached only with a higher number of layers (around 300). This dependency seems to be fully associated with multiple scattering effects, because if the same simulations are performed restricting the code to a first order approximation regime (not shown), the effect of the increase of the number of layers totally vanishes; only numerical noise around 1e-6 K is obtained without any significant tendency. The complete isothermal nature of the atmosphere and surface selected for this test shows that the origin of this observed dependency of multiple scattering on satellite measurements is purely geometrical, arising from the involved scattering angles. Therefore, even in the thermal infrared spectral range, where the non-directional nature of the Planck function dominates, the multiple scattering effects cannot be neglected, especially in modeling observations with the satellite-like geometries.

Emission source function and temperature vertical discretization: from levels to layers

Once the limits to adequately account scattering processes have been investigated in the previous section, the current section discusses the temperature vertical discretization and the calculation of the thermal source function within the layers. The re-gridding of the radiative magnitudes from levels to layers is not a trivial task. There are different approaches focused on the calculation of the most representative temperature of the layer to directly obtain the corresponding source function, while others take the source function of the levels as the starting point, in order to assume certain variation between them.

The GRASP SOS RT vertical discretization is based on the division of the full atmospheric column in thin layers of equal optical depth. No particular issues were observed for this methodology at the solar spectral range. However, at the thermal-infrared, the lack of a specific temperature discretization can raise some uncertainty when the vertical derivatives of the temperature and optical depth differ significantly. In order to illustrate the undesired interactions between the temperature profile and the TOD discretization, a TOD of 75 has been exponentially distributed along the atmospheric column in a US Standard atmosphere. This combination of atmospheric profile and optical depth distribution produces a fast temperature variation near the TOA which is not accompanied by an optical depth variation of the same magnitude, as it can be seen in fig. (4). 4) represents the temperature (T pl ) in the original grid discretized in pressure levels but expressed in terms of TOD (τ ′ ). The GRASP vertical discretisation of equal TOD levels (τ ) is marked with an orange line. No particular issues can be seen in most part of the profile, however, it is clear that in the last levels from a TOD of 0.5 to TOA the temperature variation is not adequately retained by the TOD levels.

Furthermore, figure 4 shows two different approaches to assign a temperature value to the center of the layers. In one case only a linear variation of temperature between layer limits (green) is assumed. Whereas in the other approach, the final layer temperature is calculated by a TOD weighted integration (T wt , marked in red) of all intermediate levels of the original pressure grid between each level of the GRASP TOD grid:

T wt (∆τ l ) = τ l τ l-1 T (τ ′ )dτ ′ ∆τ l . (22)
However, it is common to perform this change from level to layer grid on the emission source function instead of applying it directly on the temperature profile, like in the case of [31], [34], [38] and [46]. Analogously to the case of temperature, different assumptions of the source function variation between the levels can be used. Therefore, a linear variation from top of the layer and a TOD weighted integration of both temperature and emission source function will be also analysed here. The objective of such analysis is the comparison of the performance of the different methodologies to obtain thermal layer source function. A smaller change in the radiance values with the increase of the number of layers used in the radiative transfer calculations is considered as a sign of better correspondence between temperature and TOD discretizations. Thus, the methodology which presents a smaller sensitivity to the size of the step in the vertical discretization can be taken as more desirable. Figure 5 shows the change in average radiance, expressed as brightness temperature, of all geometries presented in figure 5 for a US-Standard atmosphere with an exponentially distributed TOD of 25. The reference radiance calculation corresponds to the calculation using 600 layers. This TOD value has been selected to minimize the effects of multiple scattering inconsistencies which have already been shown in figure 3 From the analysis of the fig.

(5) it can be concluded that ground-based measurements almost do not present any sensitivity to the number of layers. This can be explained by the fact that the possible disagreement between derivatives is located near the TOA. Therefore, in the case of the ground-based geometry the radiation reaching the ground level from layers near the TOA is very attenuated. On the other hand, in the case of the satellite geometry the problematic TOA layers are just in front of the sensor. Thus, much stronger sensitivity to discretization methodology is expected. Thus, in the ground-based case, the assumptions in the temperature discretization has a minor effect, i.e., the differences in radiance disappear with a very reasonable amount of layers (between 100 and 150).

From the curves corresponding to satellite geometry (figure 5), it is clear that the methodologies assuming linear variations present a higher dependency with the number of layers. Whereas in the case of the TOD weighted methodologies the convergence is achieved much earlier. No special difference, in this test, can be seen between the calculations based on temperature or on source function.

Thus, figure 5 serves also as a proof of the performance of the GRASP SOS RT approach to deal with elevated TOD values. At the same time, these results correspond to a specific scenario, while the peak performance of each methodology can vary among different applications.

GRASP framework enables the possibility to add extra information that is very valuable to choose between the different methodologies presented above. If only columnar information is available, then the sensitivity tests like the presented in figure 5 can be used to assess the retrieval performance. At the same time, LIDAR or another kind of vertical resolved observations (extinction, depolarization, attenuation, etc) can be used in GRASP as discussed by [47]. In such situations, the amount of necessary assumptions is reduced and the TOD weighted methodologies should be more consistent.

Tests and comparisons against DISORT

This section presents the efforts on testing and validating the radiance calculations by GRASP SOS RT code at the thermal-infrared electromagnetic spectrum. As mentioned earlier, the DISORT [3] community code will be used as the established reference. GRASP SOS and DISORT approaches have several similar features that makes direct inter-comparison of the codes performance rather transparent. For example, both codes provide RT simulation for a vertically inhomogeneous non-isothermal plane parallel atmosphere. Furthermore, the phase matrix is expressed in a Legendre polynomial expansion in both approaches. Nonetheless, despite these similarities, the radiative transfer equation solving principles are different. Some of these fundamental differences will be discussed and analyzed below.

The selected DISORT implementation is restricted to the scalar radiative transfer equation, while a full development of the vector equation has been realized in GRASP RT code. Therefore, henceforth, only total radiances will be compared without accounting for polarization effects. The calculations of atmospheric radiance were realized for a US-Standard atmosphere using one hundred atmospheric levels and a lambertian surface model at a temperature of 303.15 K. The Sun irradiance has been modeled as a perfect black body at 5250.0 K, solar direction is defined at 60 • zenith angle and at azimuth angle of φ = 0 • . These atmospheric conditions are used in all the tests presented in this section with no modifications.

The absolute difference in brightness temperature between GRASP RT and DISORT (∆BT = GRASP-DISORT) for all geometries at TOA and BOA can be found in figure 5. The upper part of each of the polar plots in this figure corresponds to observations with ground-based geometry, i.e. downward radiation. The lower part shows calculations for a satellite like geometry, i.e. upward radiation. The different panels in figure 5: A), B) and C) correspond to different amounts of total optical depth (TOD); 1.5, 0.15 and 0.015 correspondingly for each one. Each panel shows calculations at several wavelengths: 7, 8, 10, and 12 µm.

The mean difference for downward radiance averaging all selected scenarios shown in fig.( 5) is -0.008 K (-0.005%), and for upward radiance the value of the difference is -0.002 K (-0.001%). The total mean difference for all considered wavelengths, TOD values and geometries between GRASP and DISORT is -0.005 K (-0.003%). These values can be compared with the commonly expected levels of noise of the remote sensing instruments operating in the thermal infrared range, e.g.: the average accuracy for IASI observations [START_REF] Blumstein | Infrared Spaceborne Remote Sensing XII[END_REF], [START_REF] Hilton | [END_REF] of ∼0.1 K, for AIRS [48] is around 0.2 K, for CALIPSO IIR instrument measurements [48] is within 0.11 K, or for ground based observations by CLIMAT instrument [19], [20] within 0.5 K. Thus, considering these levels of observation accuracy, it can be concluded that the discrepancies between GRASP RT and DISORT are notably below of their limits of detection.

The conducted comparisons did not reveal any noticeable systematic deviation tendency or bias for any of the scenarios presented above for the diverse aerosol loads and for the different selected zenith or azimuth angles. However, in the case of downward radiation, where the influence of the phase matrix is more important, small differences in the azimuth angle can be found. This is an expected known behaviour, because in this geometry the directional nature of Sun irradiance becomes more important. However, for upward radiance, where the surface is the major contributor in this spectral range, a high degree of homogeneity in each scenario can be appreciated for all wavelengths.

Some very minor tendencies can be seen in respect to TOD variations if the radiance difference is analyzed in detail. Specifically, in the case of downward radiance, a slight shift to positive values of ∆BT is observed. The averaged difference of ∆BT at 10 µm changes from -0.013 K for the lowest TOD, to -0.004 K in the highest TOD case. The maximum change with respecto to TOD for this geometry can be found at 8 µm, where the change in ∆BT with TOD ranges from -0.001 K to 0.023 K. If an analogous comparison is performed for upward radiance, the tendency is opposite to the case of downward observations. Thus, the increase of TOD tends to slightly shift ∆BT to smaller values. At 7 µm, the change in ∆BT with TOD ranges from 0.041 K to -0.021 K. And at 8 µm, ∆BT varies from 0.043 K to 0.011 K. Despite the presence of these tendencies, its small value and the change of the sign of ∆BT between different wavelengths and geometries avoid the assessment of any significant bias in the radiance comparison.

The differences found between both RT models can be partially attributed to intrinsic characteristics of their numerical implementations. For example, these small differences can be originated by differences in the numerical rounding between FORTRAN90 and C#, the double precision of this DISORT implementation in comparison with the single precision of GRASP RT code, or the inconsistencies arising from the different analytical schemes used to solve the radiative transfer equation. Nonetheless, there are other factors that are interesting to be analyzed in detail in relation with the differences observed above. Specifically, our studies are focused on the thermal-infrared spectrum, thus the focus will be on the calculation of the Planck function of the source term for the different layers.

It has been shown that due to the non-isothermal nature of the atmosphere considered here, multiple possible implementations can be used to describe the temperature, the Planck function and the change of both magnitudes in the vertical dimension. In this respect, one of the main differences to calculate emission source function between GRASP RT and DISORT is that DISORT works with a spectrally integrated Planck function whereas GRASP calculates this magnitude in one single line. As a result, despite the smooth spectral behavior of Planck function, the differences between the two approaches are noticeable. Another main difference in the implementation of the two codes relates with how the vertical variation of temperature and Planck function is accounted for. Namely, DISORT operates with the exponential-linear interpolation method by [38] to calculate the corresponding B λ (τ ) of each layer, whereas the GRASP RT numerical approximation used in this comparison is based on assuming linear variation of the temperature profile between the levels for obtaining thermal layer source function.

In order to quantify these methodological differences between GRASP RT and DISORT, the S Emis (T (τ )) calculated by the two methodologies for the layers have been generated and extracted out of their respective complete radiative transfer schemes. Their absolute difference,∆S Emis = S Emis,GRASP (T (τ )) -S Emis,DISORT (T (τ )), is shown in fig. (6) at four wavelengths for different increases of temperature between the top and bottom limits of a single layer at a temperate around 270 K. Despite the comparison cannot be taken as an absolute reference, since both codes used other algorithmic adjustments and corrections in addition to the ones shown here, fig.( 6) is evidently useful for explaining some of the observed features of the radiance comparison shown in fig. (5).

The temperature difference between levels in the selected US-standard atmosphere used for all the comparisons shown in fig. (5) ranges between 0.5 and 1.5 K. Looking at the mean values of ∆S Emis in this temperature range, it can be noted that at 7 and 8 µm the GRASP RT methodology tends to overestimate the results of DISORT; in the case of 10 µm the differences are rather well centered around zero, and a slight underestimation can be seen at 12 µm. As fig. (5) shows very close tendencies as the differences in the calculations of S Emis (T (τ )) in the layers, it can be considered as one of the major sources of discrepancy between models. It is worth noting that even in fully isothermal conditions (∆T = 0.0 K), still small differences between GRASP and DISORT can be found due to the spectrally integrated nature of DISORT Planck function calculation, as it has been previously mentioned. Notwithstanding, there are some other minor features in the observed radiance difference that remain unexplained by the factors discussed above, for example, the differences related to the azimuth variation, or to changes of the surface reflectance.

Computation time efficiency

In order to access the computation efficiency of the developed code, the calculation speed of GRASP RT SOS radiative transfer was compared to that of DISORT. The radiances were calculated by both codes for the same scenarios and geometries as those used above illustrated in fig. (5). The execution time of all the scenarios and wavelengths have been averaged for each code for a different number of streams or guassian quadratures for Fourier expansion coefficients correspondingly for DISORT and GRASP. Figure 7 illustrates the comparison of the average execution times for GRASP RT SOS and DISORT. It is important to note that for this comparison the same configuration has been used in both RT codes including the same number of vertical layers. However, the amount of necessary atmospheric layers to achieve the same accuracy may be different between both schemes. Thus, the results shown here are an orientation. From the conducted comparison, it can be concluded that the GRASP code is significantly faster than DISORT in the whole considered range of streams. On average, the GRASP code computation time is around 8 times smaller than DISORT. It can be noted that the calculation time difference increases with the number of streams. Indeed, when less than 10 streams are used for the radiative transfer calculations, GRASP SOS runs in about 2 to 3 times faster than DISORT. These differences increase significantly when the number of streams increases. For example, for 40 streams, the GRASP SOS RT calculations are more than 30 times faster than calculations by DISORT.

6 Effect of polarization and multiple scattering in TIR At present, there is only a limited number of studies discussing the effects of polarization at TIR spectral range in the context of atmospheric remote sensing. For example, the analysis of polarimetric effects in TIR surface reflectance and in the observations of cirrus clouds can be found in [45], or in [49] respectively.

As mentioned earlier, the comparisons against DISORT have been limited to solving the scalar radiative transfer equation because of the chosen version of DISORT does not account for polarization effects. In general, the contribution of polarization into total radiances at TIR is not significant and the scalar approximation provides sufficiently accurate results. Indeed, the emitted light by the sun, atmosphere and surface in this spectral range is unpolarized. Therefore, in the first order approximation (Eq.( 6)) the radiance is unaffected by polarimetric effects and the degree of polarization (defined by the Stokes parameters ration Q/I and U/I) is close to zero if the emission from the Sun is negligible in comparison to the emission from the atmosphere and surface. However, due to multiple scattering in the atmosphere some polarization effect on the radiance and non-zero Q and U Stokes parameters of the scattered light still can arise. Although the conditions for the manifestation of the polarization seem to be quite specific, we have performed some analysis of the polarization effect in TIR for typical atmospheric aerosol conditions.

In this analysis, the TIR radiation was simulated for an scenario corresponding to a US-Standard atmosphere with a reference TOD value of 1.5 with different relative contributions of dust aerosol (AOD) and gas absorption (GOD). Figure 8 presents the absolute difference in radiance, expressed in terms of brightness temperature, calculated by GRASP SOS RT code using only the scalar equation and the complete one accounting for full polarization (∆BT Pol = Scalar -Vector). Upper parts of the polar plots correspond to a ground-based up-looking geometries, and the lower part is for satellite down-looking geometries. Panel A) corresponds to a gas dominated atmosphere (AOD=0.5 and GOD = 1.0) and panel B) to an aerosol dominated atmosphere (AOD=1.0 and GOD=0.5). The situations with less aerosol are not shown because the differences between the calculations of the two different radiative transfer schemes were hardly noticeable. Even in the situation where aerosol optical depth reaches 33% of total optical depth (panel A)), the average radiance difference for all geometries and the two selected wavelengths is only 3.09e-5 K, with a maximum difference of 2.06e-3 K. Whereas in the case where AOD is twice the GOD value (panel B)), the effects of polarization for some geometries can reach the same order of magnitude as the differences between GRASP and DISORT, with a maximum value of 1.02e-2 K. However, even in this case the average difference is 1.00e-4 K, which is far below the accuracy limits of common observations. Therefore, it can be concluded that accounting for the polarization effects in total radiances calculations at TIR is not significant for both satellite and ground-based observations.

The lack of azimuth dependence in the manifestation of polarimetric effects is related to the fact that the emission from the atmosphere and surface is dominant over the Sun incident radiation at TIR and there is no a predefined incident direction. At the same time, a dependence on the viewing zenith angle is noticeable. Namely, there is a slight increase of the radiance for smaller zenith angles if the scalar approximation is used in comparison with the values corresponding to higher angles. However, these observations can hardly be considered as a general tendency because this effect is strongly dependent on the structure of aerosol scattering matrix, aerosol vertical profile and the presence of the atmospheric gases. Thus, significant differences in this angular distribution could be found in another atmospheric scenario.

The results in fig. 8 suggest that only in the situations with a considerable load of aerosol particles and with low gaseous absorption the polarization may contribute noticeably to the total radiances at TIR spectral range.

Figure 9 shows the Q and U components of the radiance, both normalized by I, calculated for different aerosol loads. Panels A), B) and C) correspond to AOD of 1.0, 0.5 and 0.05 for a common GOD=0.5. In all the scenarios presented here the Q component is significantly more important than U, which is only noticeable at 8 µm and for medium with high aerosol loads. The azimuth structure that can be appreciated in this figure comes from the small influence of the solar radiance in TIR range. This also explain the more noticeable azimuth dependence at 8 µm. In contrary to VIS, NIR and SWIR spectral range, the degree of polarization in TIR appears predominately in the multiple scattering terms of RT solution which describes scattering of radiation emitted by atmosphere and surface. However, in all scenarios of fig. (9) both components of polarization at the spectral range considered here are very low.

Figure 10 shows the brightness temperature differences obtained in the full RT solution and in the first order approximation (∆BT Scat = Full RT -SS). Here the calculations were done for an atmosphere with equal aerosol and gas optical depths (AOD = 0.5 and GOD = 0.5). As it can be appreciated, despite the emission is originated in first order approximation of RT solution, the multiple scattering effects cannot be neglected for an accurate radiance computation. This multiple scattering effect is much more noticeable in ground-based geometries, where a maximum difference of 20 K can be found; whereas in the case of satellite, the strong influence of the surface diminishes these differences.

Conclusions

An extension of the Successive Orders of Scattering (SOS) radiative transfer technique to include thermal emission sources has been presented. Following the original radiative transfer solution of [1], the current work extends the methodology to the full electromagnetic spectrum. The numerical implementation has been performed in the framework of the GRASP code [2], which is based on the SOS radiative transfer model. Different methodologies for the discretization of the vertical variability of the optical depth and atmospheric temperature in the frame of SOS RT implementation have been analyzed. The results of this analysis suggest that the TOD weighted methodologies are preferable for realizing adequate temperature discretization in comparison to conventional linear approaches. On the other hand, any noticeable differences have not been found between the methodologies which work directly with the source function or with the temperature profile.

The validation of the GRASP SOS radiative transfer in thermal infrared has been done against DISORT reference code [3]. An average absolute difference of -0.005 K (-0.003%) was found between the radiances modeled by GRASP SOS and DISORT. This discrepancy can be considered as insignificant since these differences are far below the limit of detection of the actual atmospheric remote sensing instruments operating in the considered spectral range. Nevertheless, a comprehensive analysis of the methodological differences between the two approaches was conducted and showed that a significant part of the noted discrepancies can be explained by the use of different approaches to obtain the layer emission source function. Also the tests showed that within the outlined accuracy limit, the radiative transfer calculations by GRASP SOS are up to 30 times faster than those by DISORT if both codes are executed with the same number of streams and vertical layers. Furthermore, general insights on polarization and multiple scattering contribution in TIR have been discussed. In TIR, the main contribution to the polarized radiation (Stokes parameters Q and U) is originated in the interactions of multiple scattering regime with residual contribution of the single scattering interactions with the solar radiation in long wavelengths. However, in all the atmospheric scenarios considered here the polarimetric effects in the TIR are very low.

Appendix A

The Phase Matrix, P, is the optical characteristic which describes the directional distribution of Stokes parameters of the polarized light scattered by a particle.

For random media the phase matrix in Eq.( 3) is averaged over particle orientations and microphysical properties (size distribution, nonsphericity, etc). If Stokes parameters of the incident and scattered light are defined in the same scattering plane, the averaged over ensemble phase matrix for chaotically oriented and mirror symmetric particles depends only on the scattering angle Θ, and it can be represented by a 4x4 matrix with only 8 non-zero elements:

P(Θ) =    
For isotropic particles only six of the eight phase matrix elements are independent (p 21 = p 12 and p 43 = -p 34 ). This change between scattering and radiative transfer viewing angle can be expressed in matrix form as: P(τ, µ, φ, µ 0 , φ 0 ) = T(-χ)P(Θ)T(χ ′ ) (25) where

T(χ ′ ) =     1 0 0 0 0 cos(2χ ′ ) sin(2χ ′ ) 0 0 -sin(2χ ′ ) cos(2χ ′ ) 0 0 0 0 1     (26)
The angles χ and χ ′ refer to the angles between the scattering and the meridian planes for the scattered and incident radiance directions correspondingly. After some trigonometric transformations [START_REF] Hovenier | Transfer of polarized light in planetary atmospheres: basic concepts and practical methods[END_REF], these rotation angles can be written as: cos(χ ′ ) = cos(θ) + cos(θ 0 )cos(Θ) sin(θ 0 )sin(Θ ( 27) cos(χ) = cos(θ 0 ) + cos(θ)cos(Θ) sin(θ)sin(Θ) [START_REF] Bohren | Absorption and scattering of light by small particles[END_REF] The numerical solution of radiative transfer equation requires the expansion of the phase matrix elements in Legendre polynomials :

p ij (Θ) = L l=0 ξ i,j l P l (cos(Θ)). ( 29
)
Where p ij (Θ) are the different elements of P(Θ), and ξ i,j l are the corresponding expansion coefficients in the Legendre polynomials P l (cos(Θ)) series: ξ i,j l = 1 2(l + 1) π 0 p i,j (Θ)P l (cos(Θ))sin(Θ)dΘ (30) This expansion decouples complex angular dependencies and crucially simplifies the numerical solutions of RT equations over angular integrals.

In addition to this Legendre polynomial expansion, the necessary simplifications to solve the angular integrals appearing in the radiative transfer equation require a Fourier azimuth decomposition of the radiance and Phase matrix. The radiance can be written as:

On the other hand, the Phase Matrix is expressed as:

P(µ, φ, µ ′ , φ ′ ) = S s=0 (
A more detailed description of the expansion techniques as well as the numerical solutions can be found elsewhere [1][36] [50].

Figure 1 :

 1 Figure 1: Spherical coordinates which define radiative transfer viewing geometry.

Figure 2 :

 2 Figure 2: Schematic representation of the levels and layers discretisation.

Figure 3 :

 3 Figure 3: Radiance difference, expressed in brightness temperature, as a function of the number of layers used compare to the reference scenario with 600 layers for satellite (thick line) and ground-based (dashed line) geometries for four high TOD situations assuming a fully isothermal atmosphere at 250.0 K. Each point corresponds to the average radiance over upward and downward geometries including the full azimuth range, and a zenith range from 10 • to 70 • .

Figure 4 :

 4 Figure 4: Comparison of pressure (blue) and GRASP SOS RT equal TOD (orange) vertical discretizations for a US-Standard like atmosphere with a TOD value of 75 exponentially distributed over it. Linear temperature variation between levels (green) and TOD weighted calculations (red) to obtain layer temperatures are also shown.

Figure 5 :

 5 Figure 5: Radiance difference, expressed in brightness temperature, for different methodologies to obtain layer emission source function, compare to the reference the scenario with 600 layers shown as a function of the number of layers. Red lines show calculation for satellite-like geometry, blue lines correspond to ground-based like geometry. All the calculations are realized for a US-Standard atmosphere with a TOD of 25. Each point represents the average radiance of all geometries present in fig. (5).

Figure 6 :

 6 Figure 6: Absolute differences in layer emission source function between GRASP RT-like and DISORT-like methodologies shown as a function of the temperature difference between the top and the bottom limits of a layer at 270 K at four different wavelengths: 7, 8, 10 and 12 µm.

Figure 7 :

 7 Figure 7: Computation time in seconds as a function of the number of streams or gaussian quadratures of GRASP RT SOS and DISORT radiative transfer codes.

Figure 8 :

 8 Figure 8: Absolute difference of radiance expressed in terms of brightness temperature calculated with and without polarization effects (∆BT Pol = Scalar -Vector) at 8 and 10 µm, for a solar zenith angle of 60 • . Panel A) corresponds to a gas dominated atmosphere (AOD=0.5 and GOD = 1.0) and panel B) to an aerosol dominated atmosphere (AOD=1.0 and GOD=0.5). Upper part of the polar plots corresponds to downward radiation (ground-based like geometry) and the lower part to upward radiation (satellite-like geometry).

Figure 9 :

 9 Figure 9: Q and U radiance components normalized by I at 8 and 10 µm for different aerosol loads. Panels A), B) and C) corresponds to AOD equal to 1.5, 0.5 and 0.05, in this example the GOD is fixed to 0.5 for all of them. Upper part of the polar plots corresponds to downward radiation (ground-based like geometry) and lower part to upward radiation (satellite-like geometry).

Figure 10 :

 10 Figure 10: Absolute difference of radiance, expressed in terms of brightness temperature, for full GRASP SOS RT solution and the solution in first order approximation (∆BT Scat = Full RT -SS) at 8 and 10 µm, for an atmosphere with medium aerosol load (AOD=0.5 and GOD = 0.5) and a solar zenith angle of 60 • . Upper part of the polar plots corresponds to downward radiation (ground-based like geometry) and lower part to upward radiation (satellite-like geometry).

  cos(Θ) = -[cos(θ)cos(θ 0 ) + sin(θ)sin(θ 0 )cos(φ -φ 0 )]

L( 2 -( 2 -( 2 -

 222 n (τ, µ, φ) = S s=0 δ 0s )[cos[s(φ -φ 0 )]L s n,cos (τ, µ) + sin[s(φ -φ 0 )]L s n,sin (τ, µ)](31)withL s n,cos (τ, µ) = (I s n , Q s n , 0, 0) and L s n,sin (τ, µ) = (0, 0, U s n , V s n ) (32)whereI n (τ, µ, φ) = S s=0 δ 0s )cos(s(φ -φ 0 ))I s n (τ, µ) Q n (τ, µ, φ) = S s=0 δ 0s )cos(s(φ -φ 0 ))Q s n (τ, µ) U n (τ, µ, φ) = S s=1 2 sin(s(φ -φ 0 ))U s n (τ, µ) V n (τ, µ, φ) = S s=12 sin(s(φ -φ 0 ))V s n (τ, µ)

  2 -δ 0s )[cos[s(φ -φ ′ )]P s cos (µ, µ ′ ) + sin[s(φ -φ ′ )]P ssin (µ, µ ′ )] P 22 -sc ′ P 33 sP 43 sP 12 sc ′ P 22 -cs ′ P 33 c stands for cos(2χ), s for sin(2χ), and similarly for c' and s' with χ'. The source function in first order approximation and multiple scattering can be expressed as: -δ 0s )[cos[s(φ -φ 0 )]P s cos (µ, µ 0 )E 0 e τ /µ0 (36a)

				S 1 (τ, µ, φ) = (2 S n>1 (τ, µ, φ) = S ω(τ ) 4π s=0 S ω(τ ) s=0 2 (2 -δ 0s ) cos[s(φ -φ 0 )]	+1 -1	(P s cos L s n-1,cos -P s sin L s n-1,sin )dµ ′
					+1
				+sin[s(φ -φ 0 )]	(P s sin L s n-1,cos + P s cos L s n-1,sin )dµ ′
					-1
							(34)
		with
	   	P 11 cP 12 cc ′ P 22 + ss ′ P 33 c ′ P 12 0 0 0 0	0 0 ss ′ P 22 + cc ′ P 33 -cP 43 0 0 c ′ P 43 P 44	    =	S s=0	(2 -δ 0s )cos(s(φ -φ ′ ))P s cos (µ, µ ′ )
							(35a)
		0	0	s ′ P 12	0	
	  	0 0	0 -s ′ P 43	cs ′ 0 0	0 0	   =

S s=1 2 sin(s(φ -φ ′ ))P s sin (µ, µ ′ ) (35b)

where

Acknowledgements

This work was supported by Cifre/ANRT program 2019/0003 through a PhD thesis co-funding of M. Herreras-Giralda. The work was also supported by the CaPPA project (Chemical and Physical Properties of the Atmosphere) funded by the French National Research Agency (ANR) through the PIA (Programme d'Investissement d'Avenir) under contract ANR-11-LABX-0005-01 and by the Regional Council Hautsde-France and the European Funds for Regional Economic Development (FEDER).