# Modified L-distributions method for atmospheric remote sensing problems

Romuald Tapimo, Laurence Croize, Frédéric André, Jean Pierro

► To cite this version:

Romuald Tapimo, Laurence Croize, Frédéric André, Jean Pierro. Modified L-distributions method for atmospheric remote sensing problems. Applied optics, 2025, 64 (29), pp.8603. 10.1364/AO.561291. hal-05468743

HAL Id: hal-05468743

https://hal.science/hal-05468743v1

Submitted on 21 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Modified L-distributions method for atmospheric remote sensing problems

ROMUALD TAPIMO,<fnref n="1" /> LAURENCE CROIZE,<fnref n="1" />,* FREDERIC ANDRE,<fnref n="2" /> AND JEAN PIERRO<fnref n="1" />

Received 7 March 2025; revised 21 August 2025; accepted 24 August 2025; posted 25 August 2025; published 29 September 2025

Many applications involving atmospheric radiation need radiative transfer codes with an appropriate combination of accuracy, speed, and spectral coverage. In this paper, a new formulation, to our knowledge, of the L-distributions (LD) method is developed for modeling infrared sensor transmittances and radiances in inhomogeneous thermodynamic atmospheres containing a mixture of absorbing gases with variable concentrations. The formulation, called the modified L-distributions (MLD) method, overcomes the limitations of the LD method and establishes a balance between speed and accuracy. The atmosphere is considered a plane-parallel medium, limited at the top by the vacuum, while the kernel bidirectional reflectance distribution function is used to describe the anisotropy of the ground surface. The vertical structure of the plane-parallel atmosphere is subdivided into homogeneous layers of constant or variable discretization steps. It is demonstrated that the transmittance in an inhomogeneous atmosphere using the MLD method requires solving an equivalent homogeneous problem with the absorption coefficient equal to the weighted sum of absorption coefficients of all layers (above or below the observation point) and with thickness equal to the highest discretization step. The weighting factor is the ratio of the layer height to the highest discretization step. Compared to the standard version of the LD method, the MLD formulation is shown to be (1) independent of the direction of propagation along a non-uniform path, (2) more accurate, and (3) more computationally economic and efficient. The line-by-line (LBL) results are used as a benchmark for the comparison of the atmospheric transmittance. Daytime and nighttime top-of-atmosphere (TOA) radiances for a Lambertian ground surface are plotted and compared with the atmospheric radiative transfer code MATISSE outputs. The MLD method results match well with the MATISSE ones. Four locations with anisotropic surface configurations are also studied, and the daytime and nighttime TOA irradiance and albedo are estimated. For the atmospheres studied in this work, the MLD method yields absolute errors less than 1.1% transmittance with a CPU time of 0.2 s on an Intel Xeon W1250P 4.10 GHz computer. This means that the MLD method is about 2.5 times more accurate and 30 times faster than the standard LD method.

Published by Optica Publishing Group under the terms of the [Creative Commons Attribution 4.0 License](https://creativecommons.org/licenses/by/4.0/). Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
https://doi.org/10.1364/AO.561291

## 1. INTRODUCTION

The electromagnetic radiation from ultraviolet to microwave wavelengths is widely used by Earth observation sensors (aircraft, satellite, etc.). In the infrared spectral band, sensors capture the electromagnetic radiation resulting from the ground surface emission/reflection, atmospheric emission, and the Sun. The measured radiation is sent to a ground station where it is processed, analyzed, and interpreted using radiative transfer models [1]. In this context, the radiative transfer models serve as a bridge between the physical processes that generate signals and data collected by the sensors: forward radiative transfer models map the state of the atmosphere to the measurements, while the inverse radiative transfer models map measurements

to the state of the atmosphere [2]. Applications that involve the use of sensor data (radiative transfer) include Earth's energy budget [3-5], weather prediction, global climate change studies [6], land cover mapping (city, forest, etc.), monitoring of vegetation, armed conflicts [7,8], estimation of soil parameters [9-11], surface temperature [12], atmospheric gas concentration [4,13], and properties of the Earth's atmosphere and oceans [14-16]. All these applications and the increase in spatial and spectral resolution of satellites lead to an exponential increase in the amount of remote sensing data produced in the world [17-21]. Data explosion has created data-management problems, which in turn have an impact on data transmission, storage, processing, analysis, and interpretation of information [18]. Accordingly, this explosion of remote sensing data

<footnote n="1">¹ONERA, DOTA, Université Paris-Saclay, FR-91123 Palaiseau, France</footnote>

<footnote n="2">²Laboratoire d'Optique Atmosphérique, Département de Physique, Faculté de Sciences et Technologies, Université de Lille, Lille, France</footnote>

<footnote>*laurence.croize@onera.fr</footnote>
has become a very strong challenge to solve. To contribute to the advancement of atmospheric remote sensing research and applications, radiative transfer codes are developed in specific spectral regions to simulate changes in radiation intensity in the atmosphere [22]. However, due to the vast volume of remote sensing data, the radiative transfer calculations for estimating the state of the atmosphere should be extremely fast to enable the processing and the analysis in a reasonable time [18,19]. Moreover, for the transmitted information to be reliable and really exploitable, the codes must also be stable and accurate in a large number of channels in the electromagnetic spectrum [2]. To meet these constraints, researchers in the remote sensing community still need to develop novel and accurate radiative transfer techniques that significantly reduce the computational cost.

The most accurate radiative transfer method to simulate changes in radiation intensity is the line-by-line (LBL) method [23]. The LBL method is a high-resolution analytical modeling of absorption lines based on spectroscopic parameters available in databases such as HITRAN or GEISA [24,25], and its solution is mainly used as a benchmark to verify other radiative transfer method solutions. Several forward radiative transfer codes to simulate changes in radiation intensity in the atmosphere are based on an LBL approach [26,27]. Among the available codes, the line-by-line radiative transfer model (LBLRTM) [28], automatized atmospheric absorption atlas (4 A/OP) [29], a general line-by-line atmospheric transmittance and radiance model (GENLN2) [13], a fast line-by-line code for high resolution infrared molecular spectra (RADIS) [30], blue sky transmission and radiance atmospheric model (BTRAM) [31], and reference forward model (RFM) [32] can be cited. All these radiative transfer codes using the LBL approach spend huge time to examine radiation transfer in the atmosphere to retrieve the geophysical parameters of atmospheric gases [23,26,28]. Other radiative transfer codes, such as advanced modeling of the Earth for environment and scenes simulation (MATISSE) [14,33], moderate resolution atmospheric transmission (MODTRAN) [34,35], algorithms for modeling spectral transmittance and radiance (LINEPAK) [36], Santa Barbara DISORT atmospheric radiative transfer (SBDART) [37], and rapid radiative transfer model (RRTM) and variants [6,26,27] use the band, *k*-distributions, or correlated *k*-distributions techniques to calculate the atmospheric transmittance and radiance. However, the band, *k*-distributions, or correlated *k*-distributions techniques use the line-by-line procedure to compute the absorption coefficient of atmospheric gases in specific spectral regions [6,38–40], and the computational time may be influenced greatly as the number of narrow bands increases [26,41,42]. Although the diversity of radiative transfer models with different features, their performances (accuracy and speed) remain debatable [26]. Generally, for these codes, a higher resolution dataset (spectral and spatial) yields higher accuracy with lower speed, while a lower resolution dataset yields less accurate outcomes with increased speed [9]. It remains, therefore, a challenge to establish an appropriate balance between the speed and the accuracy of the radiative transfer models used in remote sensing analysis.

André [43] developed the L-distributions (LD) method and computed the transmissivity in thin homogeneous/

inhomogeneous gaseous media, where high temperatures occur. In his work, the total pressure, minimum temperature, and maximum thickness of the medium were set to 1 atm, 300 K, and 100 cm, respectively. Moreover, the author considered only the fact that the light, in the medium, propagates perpendicularly to the boundaries. These considerations do not really describe a medium like the atmosphere. Tapimo et al. [44] used the LD method to compute the transmissivity of the inhomogeneous atmosphere in the medium wave infrared (MWIR) portion of the electromagnetic spectrum, and efficient results were obtained. In the application of the LD method for inhomogeneous atmospheres, the vertical structure of the atmosphere is subdivided into uniform layers parallel to the ground surface. Adjacent layers of constant thermophysical properties are successively transformed into homogeneous layers with equivalent path lengths. The look-up tables involved in the method to construct the set of equivalent lengths are created and stored for all thermophysical states of the gas, spectral bands, and a set of values of the mapping function over [0,1]. The mapping functions associate values of transmissivities or absorptivities of a fast model with the corresponding values of the LBL calculations [43]. Tapimo et al. [44] considered 5 km inhomogeneous atmospheres and computed the atmospheric transmittance for tropical, midlatitude summer, midlatitude winter, and the US standard atmospheres composed of 28 gaseous molecules. They also compared results with the traditional LBL results and demonstrated that the LD method can be more accurate than 2.5% in absolute error from 2000 to 2500 cm⁻¹ spectral region with a CPU time less than 6 s while considering 10³ values of the mapping function over [0,1] to tabulate the method’s look-up tables. It should be noted that the accuracy of the method could be improved by increasing the size of the look-up tables. Although the LD method is fast with an accuracy close to the LBL method, it is memory expensive, especially when high-resolution look-up tables are needed. Moreover, the successive transformations of adjacent layers do not yield the same result when the direction of transformation is from the top-of-atmosphere (TOA) toward the ground and inversely [44].

The purposes of this work are (1) to propose a modified L-distribution (MLD) method that overcomes several limitations of the LD method and (2) to establish the balance between the accuracy and the speed of the method for remote sensing applications. In the development of the MLD method, the vertical structure of the atmosphere is discretized into parallel layers with constant or variable height (discretization step along the vertical axis). By an appropriate change of variables, equations are transformed into forms similar to the LD method equations in homogeneous atmospheres with an effective absorption coefficient equal to the weighted sum of the absorption coefficients of different layers and with a thickness equal to the highest discretization step. The weight factor is the ratio of the layer height to the highest discretization step. With this transformation of equations, look-up tables of the MLD method can be generated. It should be noted that the LD and MLD methods differ only in inhomogeneous media and are identical in homogeneous ones. In this work, it is assumed that the non-scattering atmosphere emits or absorbs radiation only in the vertical direction and that
each level of the atmosphere is in local thermodynamic equilibrium. It is also assumed that atmospheres of unit refractive indices are limited at the top by the vacuum, while the spectral variation in ground surface emissivity is negligible in the considered spectral range.

In Section 2, details on the LD method equations in a homogeneous atmosphere under local thermodynamic equilibrium are given. This section is important for introducing the derivation of the MLD method equations in inhomogeneous atmospheres and their similarities to those in homogeneous ones. In Section 3, the inhomogeneous atmospheres are subdivided into many layers of constant or variable thicknesses, and the equations of the MLD method are derived and compared to those given in Section 2. Specific remote sensing problems of land cover mapping and the retrieval of gas concentrations, surface temperature, and vertical temperature profiles in a planetary atmosphere from measured spectral radiances are studied in Section 4. The bidirectional reflectance distribution function (BRDF) of the ground is also given. The BRDF is a mathematical tool that describes the way light scatters in different directions when it interacts with a surface. Transmissivities are validated in Section 5 for the tropical, midlatitude summer, midlatitude winter, and US standard atmosphere profiles, providing the volume mixing ratio of 28 gaseous molecules. A comparison between the LD and MLD formulations in inhomogeneous atmospheres is discussed. The sensor radiance and albedo are estimated for different reflected surface types. Section 6 is devoted to the conclusion of the study and additional remarks.

## 2. L-DISTRIBUTIONS METHOD IN HOMOGENEOUS ATMOSPHERES

In this subsection, the vertical structure of the plane-parallel atmosphere and the absorption cross section of the absorbers are independent on the altitude of a point within the atmosphere. Let's consider an atmosphere of thickness *H* and spectral absorption coefficient *κ*<sub>ν</sub> (in m<sup>-1</sup>). In a narrow band Δν centered on $\bar{\nu}$, the averaged transmittance of a radiation propagating in an angular direction *μ* is [45–47]

$$T_{r,\bar{v}}(L_{\mu}) = \frac{1}{\Phi} \int_{\Delta\nu} \phi(v) e^{-\kappa_{\nu} L_{\mu}} dv, \quad (1)$$

where $\mu \in [0, 1]$ is the cosine of the viewing zenith angle, $L_{\mu} = H/\mu$ is the path traveled by the radiation in the atmosphere, $\phi(v)$ the spectral response of the measuring instrument and

$$\Phi = \int_{\Delta\nu} \phi(v) dv. \quad (2)$$

For an atmosphere composed of a mixture of *J* gases, $κ_ν$ is the weighted sum over all individual species absorption coefficients, $κ_{ν,j}$. That is [45,48]

$$\kappa_{\nu} = \sum_{j=1}^{J} \chi_{j} \kappa_{\nu, j}, \quad (3a)$$

where $\chi_j$ is the volume fraction of the jth gas molecule. Equation (3a) can also be written as follows:

$$\kappa_{\nu} = \sum_{j=1}^{J} \sigma_{\nu, j} N_{j}, \quad (3b)$$

where $\sigma_{\nu,j}$ and $N_j$ are, respectively, the cross section and number density of the attenuating species *j*. The number density is related to the pressure *p* and temperature *T* of the atmosphere as follows:

$$N_j = \frac{\chi_j p}{k_B T}, \quad (3c)$$

where $k_B$ is the Boltzmann constant. In Eqs. (3a) and (3b), the mixing ratio, $\chi_j$, is expressed in ppm.

To derive equations of the LD method in a homogeneous slab atmosphere, let us consider the function

$$\Omega(L_{\mu}) = \ln[1 - \alpha_{\bar{v}}(L_{\mu})], \quad (4)$$

where $\alpha_{\bar{v}}(L_{\mu}) \in [0, 1]$ is the absorptivity of the atmosphere given by

$$\alpha_{\bar{v}}(L_{\mu}) = 1 - T_{r,\bar{v}}(L_{\mu}) = 1 - \frac{1}{\Phi} \int_{\Delta\nu} \phi(v) e^{-\kappa_{\nu} L_{\mu}} dv. \quad (5)$$

Since $\alpha_{\bar{v}}(L_{\mu})$ can be expanded in a Taylor series in the neighborhood of $L_{\mu}^0 = 0$ and $\alpha_{\bar{v}}(L_{\mu}^0) = 0$, it follows that $\Omega(L_{\mu})$ has a series expansion in the neighborhood of $\alpha_{\bar{v}} = 0$, so that

$$\Omega(L_{\mu}) = - \sum_{n=1}^{\infty} \frac{\alpha_{\bar{v}}^{n}(L_{\mu})}{n}. \quad (6)$$

In the neighborhood of $L_{\mu}^0 = 0$, the function $\Omega(L_{\mu})$ can be expanded in a Taylor series of order two as follows:

$$\begin{aligned} \Omega(L_{\mu}) &= \Omega(L_{\mu}^{0}) + L_{\mu} \left. \frac{\partial \Omega}{\partial L_{\mu}} \right|_{L_{\mu}=L_{\mu}^{0}} + \frac{L_{\mu}^{2}}{2} \left. \frac{\partial^{2} \Omega}{\partial^{2} L_{\mu}} \right|_{L_{\mu}=L_{\mu}^{0}} \\ &= \omega_{1} L_{\mu} + \omega_{2} L_{\mu}^{2}, \end{aligned} \quad (7)$$

where

$$\omega_1 = -k_p = -\frac{1}{\Phi} \int_{\Delta\nu} \phi(v) \kappa_\nu \, dv \quad (8)$$

and

$$\omega_2 = \frac{1}{2} \left[ \frac{1}{\Phi} \int_{\Delta\nu} \phi(v) \kappa_{\nu}^2 dv - k_p^2 \right]. \quad (9)$$

Following [43], the inverse series expansion of Eq. (7) can be obtained as follows:

$$k_p L_{\mu} = -\Omega(L_{\mu}) + \frac{\pi}{2\beta} \Omega^2(L_{\mu}), \quad (10)$$

where

$$\beta = \frac{\pi \omega_1^2}{2 \omega_2}. \quad (10a)$$

The solution of Eq. (10) for $\Omega(L_{\mu})$ gives

$$\Omega(L_{\mu}) = -\frac{\beta}{\pi} \left( \sqrt{1 + \frac{2\pi k_p L_{\mu}}{\beta}} - 1 \right). \quad (11)$$
Then, Eq. (4) provides the approximate absorptivity of a uniform atmosphere as follows:

$$ \alpha_{\bar{v},*}(L_{\mu}) = 1 - \exp \left\{ -\frac{\beta}{\pi} \left( \sqrt{1 + \frac{2\pi k_p L_{\mu}}{\beta}} - 1 \right) \right\}. \quad (12) $$

The look-up tables of the LD method can now be generated for a set of prescribed values of $X \in [0, 1]$ and all spectral bands considered using the equation [46]

$$ G_r(X) = \frac{1}{\Phi} \int_{\Delta v} \phi(v) e^{-\kappa_v \Lambda(X)} dv, \quad (13) $$

where

$$ \Lambda(X) = - \frac{\ln[1-X]}{k_p} + \frac{\pi k_p}{2\beta} \left[ \frac{\ln[1-X]}{k_p} \right]^2. \quad (14) $$

The transmittivity of a given uniform path is deduced by interpolation from the look-up tables so that

$$ T_{r,\bar{v}}(L_{\mu}) = G_r[\alpha_{\bar{v},*}(L_{\mu})]. \quad (15) $$

Functions $G_r(X)$ are called mapping functions and establish a relationship between the values taken by the approximate model [Eq. (12)] and the LBL data. It should be noted that more detailed information on the derivation of Eqs. (13)–(15) are found in Ref. [43] and are not reported here for brevity.

### 3. MODIFIED L-DISTRIBUTIONS METHOD IN INHOMOGENEOUS ATMOSPHERES

Here, the absorption cross section of the absorbers depends on the thermophysical state of the atmosphere at altitude $z$. The vertical structure of the plane-parallel atmosphere is subdivided into $M$ layers of constant thermophysical properties and each layer height (discretization step) is $z_m - z_{m-1}$ ($1 \le m \le M$). In inhomogeneous atmospheres, the LD method consists of successive transformation of adjacent layers to an equivalent layer with an equivalent path length. Its equations are not reported in this work and details can be found in Refs. [43,44,46]. Only the MLD method, which consists of transforming equations in inhomogeneous atmospheres into its equivalent in homogeneous atmospheres is developed. In a narrow band, $\Delta v$, and for an altitude $z$ above the ground surface, the band averaged transmissivity between $z$ and the top of atmosphere (TOA) is written as follows [44]:

$$ \begin{aligned} T_{r,\bar{v}}(L_{\mu,m}, L_{\mu,m+1} \dots L_{\mu,M}) \\ = \frac{1}{\Phi} \int_{\Delta v} \phi(v) \exp \left[ -\left( \kappa_v^m L_{\mu,m} + \dots + \kappa_v^M L_{\mu,M} \right) \right] dv, \end{aligned} \quad (16) $$

where $L_{\mu,m} = (z_m - z_{m-1})/\mu = L_m/\mu$ is the length traveled by the radiation in the layer, $\mu \in [0, 1]$ is the cosine of the viewing zenith angle, $\kappa_v^m$, in m$^{-1}$, is the absorption coefficient of the $m$th layer (just above the observation point), and $\kappa_v^M$ is the absorption coefficient of the layer just below the TOA ($z = z_\infty$). The subscript $m$ runs from $M$ to 1, and $z_0 = 0$ corresponds to the ground surface (see Fig. 1). In Eq. (16), $\kappa_v^m$ are constant with respect to the altitude and are given by Eq. (3b). To derive

<figure><img src="image_4.png" /><figcaption>Fig. 1. 1D geometry of the formulation of the absorption and emission processes in the atmosphere. The view angle $\theta$ is measured with respect to $-z$ and $\mu = \cos \theta$. O is the observation point (the point where the radiance is measured).</figcaption></figure>

equations of the MLD method in inhomogeneous atmospheres, Eq. (4) is rewritten as follows:

$$ \begin{aligned} & \Omega (L_{\mu,m}, L_{\mu,m+1} \dots L_{\mu,M}) \\ &= \ln [1 - \alpha_{\bar{v}} (L_{\mu,m}, L_{\mu,m+1} \dots L_{\mu,M})]. \end{aligned} \quad (17) $$

To follow the idea described in Section 2, Eq. (17) is expanded in Taylor series at order two in the neighborhood of $L_{\mu,m}^0, L_{\mu,m+1}^0, \dots, L_{\mu,M}^0$. To facilitate the expansion, Eq. (17) is first transformed into an equation of a single variable, $\xi$, by setting $L_{\mu,m} = L_{\mu,m}^0 + b_m \xi$. Second, the resulting equation is expanded at a neighborhood of $\xi = 0$. Third, $L_{\mu,m}^0$ and $\xi$ are set, respectively, to 0 and 1. If $\Delta z$ denotes the highest discretization step, then the coefficients $\omega_1$ and $\omega_2$ of the series expansion are obtained at the neighborhood of $L_{\mu,m}^0 = 0$ ($1 \le m \le M$) as follows:

$$ \omega_1 = -k_p = -\frac{1}{\Phi} \int_{\Delta v} \phi(v) \bar{\kappa}_v dv \quad (18) $$

and

$$ \omega_2 = \frac{1}{2} \left[ \frac{1}{\Phi} \int_{\Delta v} \phi(v) \bar{\kappa}_v^2 dv - k_p^2 \right], \quad (19) $$

where $\bar{\kappa}_v$ is the weighted sum of the absorption coefficients of all layers above the observation point and is written as follows:

$$ \bar{\kappa}_v = \frac{L_m}{\Delta z} \kappa_v^m + \frac{L_{m+1}}{\Delta z} \kappa_v^{m+1} + \dots + \frac{L_M}{\Delta z} \kappa_v^M. \quad (20a) $$

In Eq. (20a), $\Delta z = \max\{L_m, L_{m+1}, \dots, L_M\}$. The look-up tables of the inhomogeneous atmosphere above the observation point can now be generated for a set of variables $X \in [0, 1]$ and all spectral bands considered using Eqs. (13) to (14), where $\kappa_v$ and $L_\mu$ are, respectively, replaced by $\bar{\kappa}_v$ and $\bar{L}_\mu = \Delta z/\mu$. Then, the transmittivity of the inhomogeneous atmosphere is easily calculated. If the radiation measured is the terrestrial radiation then $\bar{\kappa}_v$ is the sum of absorption coefficients of all layers below the observation point. For AFGL reference atmosphere models, the maximum discretization step is $\Delta z = 5$ km and the weighting factors (ratio of the layer height to the highest discretization step) in Eq. (20a) are given by

$$ \frac{L_m}{\Delta z} = \begin{cases} \frac{1}{5}, & 0 \text{ km} \le z \le 25 \text{ km} \\ \frac{1}{2}, & 25 \text{ km} < z \le 50 \text{ km} \\ 1, & 50 \text{ km} < z \le 120 \text{ km} \end{cases}. \quad (20b) $$
It should also be noted that the main advantage of the MLD method over the LD method is that it does not require the inverse $G_r^{-1}(X)$ of $G_r(X)$ as needed in LD method for the successive transformations of adjacent layers [43,44,46]. Consequently, the MLD method is more computationally efficient and memory saving than the LD method. Moreover, the formulation is reversible, that is independent on the direction of propagation (from $z_\infty$ to $z$ or from $z$ to $z_\infty$) along a non-uniform path. At an observation point within the atmosphere, the sky/terrestrial signal traces back its path when traveling in the reverse direction with the same angular direction, and so, the same amount of radiation from the local point $z$ is sensed by the instruments at sky/ground. It is important to recall that the LD and MLD methods differ only for inhomogeneous atmospheres. The algorithm of computing the MLD method is described in Appendix A.

## 4. REMOTE SENSING PROBLEMS

Consider the remote sensing problem of land cover mapping, and the retrieval of the vertical temperature profile, surface temperature, and gas concentrations. Let us assume that the TOA is transparent and is illuminated in the angular direction $\mu_0$ by the solar radiation, while the Earth's surface can emit and reflect radiation at the temperature $T_0$. For a plane-parallel atmosphere under thermodynamic equilibrium, the spectral upward and downward radiances at the sensor's aperture at an altitude $z$ (see Fig. 1) can be expressed in terms of relevant atmospheric parameters, respectively, as follows:

$$ I_v(z, -\mu) = I_v(z_b, -\mu)e^{-\tau_v(z_b, z)/\mu} + \int_{z_b}^{z} B_v[T(z')] \frac{d}{dz'} \left[ e^{-\tau_v(z', z)/\mu} \right] dz' \quad (21) $$

and

$$ I_v(z, \mu) = I_v(z_\infty, \mu)e^{-\tau_v(z, z_\infty)/\mu} - \int_z^{z_\infty} B_v[T(z')] \frac{d}{dz'} \left[ e^{-\tau_v(z, z')/\mu} \right] dz', \quad (22) $$

where $\mu$ is the cosine of the observation angle measured with respect to negative $z$, $\tau_v$ is the spectral optical depth of the atmospheric column between $z$ and $z'$, $z_b$ the altitude of the lower boundary (at the ground surface, $z_b = z_0 = 0$) and $B_v(T)$ the spectral black-body radiance. The spectral optical depth is given by

$$ \tau_v(z, z') = \int_z^{z'} \kappa_v(z'') dz''. \quad (23) $$

The spectral black-body radiance is given by [13,42]

$$ B_v(T) = \frac{2hc^2v^3}{ehc_v/k_bT - 1}. \quad (24) $$

Since the average distance between the Earth and the Sun is approximately thousands of million kilometers, beams at the TOA are assumed to be parallel, and the spectral extraterrestrial solar irradiance, $I_{v,\infty}$, is approximated by the spectral black-body radiance at the temperature $T_{sky} = 5800$ K multiplied

by the angle from which the Earth sees the Sun [49]. Then, the boundary conditions subject to Eqs. (21) and (22) can be written as follows [13,50]:

$$ I_v(z_b, -\mu) = \varepsilon_v(\mu) B_v(T_b) + \frac{1}{\pi} \int_0^{2\pi} \int_0^1 \rho_v(\mu, \mu', \varphi) I_v(z_b, \mu') \mu' d\mu' d\varphi \quad (25) $$

and

$$ I_v(z_\infty, \mu) = I_{v,\infty}\delta(\mu - \mu_0)\delta(\varphi - \varphi_0), \quad (26) $$

where $\delta(\cdot)$ is the Dirac function and $T_b$ the temperature of the lower boundary. $\rho_v(\mu, \mu', \varphi)$ is the bidirectional reflectance distribution function (BRDF) that captures the variability of natural surfaces and land cover, $\varepsilon_v(\mu)$ the associated directional emissivity and $\varphi$ the viewing azimuthal angle.

Substituting Eq. (22), for $z = z_b$, into Eq. (25) leads to the following expression of the ground radiance [50–52]:

$$ I_v(z_b, -\mu) = F_{v,\infty}(\mu, \mu_0)e^{-\tau_v(z_b, z_\infty)/\mu_0} + I_{v,1}(z_b, -\mu), \quad (27) $$

where

$$ F_{v,\infty}(\mu, \mu_0) = \frac{\mu_0}{\pi} I_{v,\infty} \rho_v(\mu, \mu_0), \quad (28) $$

$$ I_{v,1}(z_b, -\mu) = \varepsilon_v(\mu) B_v(T_b) + 2 \int_0^1 \rho_v(\mu, \mu') I_{v,d}(z_b, \mu') \mu' d\mu', \quad (29) $$

and

$$ I_{v,d}(z_b, \mu') = - \int_{z_b}^{z_\infty} B_v[T(z')] \frac{d}{dz'} \left[ e^{-\frac{\tau_v(z_b, z')}{\mu}} \right] dz'. \quad (30) $$

Equations (21) and (22) assume that the radiance is invariant with respect to the azimuthal angle. In order to compare the calculation and observation results, the radiance at the sensor aperture must be modified by convolution with the spectral response of the instrument [21,47], that is

$$ I_{\bar{v}}(z, \pm\mu) = \frac{1}{\Phi} \int_{\Delta v} \phi(v) I_v(z, \pm\mu) dv, \quad (31) $$

where $\phi(v)$ is the spectral response of the sensor. Under the assumption that the variation of the blackbody spectral radiance is small over a narrow band $\Delta v$, and that the second term in Eq. (29) is small compared to the first one, the spectral correlation between $e^{-\tau_v}$ and $I_{v,1}(z_b, -\mu)$ can be neglected [53] and the mean spectral radiances measured at the upper and lower boundaries of the atmosphere can be approximated, respectively, as follows [47,52,54,55]:

$$ I_{\bar{v}}(z_\infty, -\mu) = F_{\bar{v},\infty}(\mu, \mu_0) T_{r,\bar{v}}(z_b, z_\infty, \mu^*) + I_{\bar{v},1}(z_b, -\mu) T_{r,\bar{v}}(z_b, z_\infty, \mu) + \int_{z_b}^{z_\infty} B_{\bar{v}}[T(z')] \frac{d}{dz'} [T_{r,\bar{v}}(z', z_\infty, \mu)] dz' \quad (32) $$
and

$$I_{\bar{v}}(z_b, \mu) = I_{\bar{v}}(z_{\infty}, \mu) T_{r, \bar{v}}(z_b, z_{\infty}, \mu) - \int_{z_b}^{z_{\infty}} B_{\bar{v}}[T(z')] \frac{d}{dz'} [T_{r, \bar{v}}(z_b, z', \mu)] dz'. \quad (33)$$

Referring to the results of Hulley and Hook (cf. Fig. 7) [56], this approximation on the spectral correlation can be reasonable for infrared regions below 1000 cm⁻¹ and beyond 1250 cm⁻¹. In Eq. (32), $\mu^* = \mu_0\mu/(\mu_0 + \mu)$. The first term on the right-hand side of Eq. (32) contributes to the outgoing direct solar radiation reflected on the lower boundary, attenuated along the two-way paths through the absorbing atmosphere. The second term is the sum of the radiation emitted from the surface and downward atmospheric radiation reflected on the surface, attenuated along the path through the absorbing atmosphere. The third term contributes to the emission from the atmosphere. Equation (32) assumes that the spectral variation of the directional emissivity is small [54,56], and thus, the subscript can be dropped for simplicity. Under this assumption, the BRDF can be directly determined from the sensor observations (see Eq. (27)) if the ground radiance is known [55]. The directional emissivity is related to the BRDF thanks to the first Kirchhoff's law as follows [57,58]:

$$\varepsilon(\mu) = 1 - 2\pi \int_{0}^{1} \rho(\mu, \mu')\mu' d\mu', \quad (34)$$

where

$$\rho(\mu, \mu') = \frac{1}{2\pi} \int_{0}^{2\pi} \rho(\mu, \mu', \varphi) \, d\varphi. \quad (35)$$

In a general situation of vegetation canopy region over a ground material [57–60]

$$\rho(\mu, \mu', \varphi) = k_{iso} + k_{vol} f_{vol}(\mu, \mu', \varphi) + k_{geo} f_{geo}(\mu, \mu', \varphi). \quad (36)$$

In Eq. (36), $k_{iso}$ is an isotropic Lambertian contribution and is equivalent to the reflectance when $\mu = \mu'$ (specular reflection), $k_{vol}$ is related to the protrusion height, length, and width, and $k_{geo}$ is related to the facet reflectance and shadowing effects. In a plane-parallel atmosphere, remote sensing provides semi-empirical expressions of the Roujean's volumetric kernel $f_{vol}$ [59], and LiSparse-R geometric kernel $f_{geo}$ depending on the geometric configuration [58]. It should be noted that the measured radiances are linked to the state of the atmosphere due to the absorption coefficient. In Eqs. (32) and (33), it is assumed that the source term is isotropic and its variation over the narrow band $\Delta v$ is negligible. The physical state of the atmosphere can be estimated from a set of radiance measurements as function of the band center $\bar{v}$, the angular direction of measurements, $\mu$, and the solar illumination, $\mu_0$. The observed radiances at the TOA can also be converted into radiative flux by averaging Eq. (32) over the upper hemisphere. The ratio of radiant flux density leaving a surface to the incident radiant flux density defines the spectral albedo [61]. Note that the spectral albedo is crucial to control the energy budget of the Earth, surface spectral albedo, and the amount of CO₂, water vapor and ozone in the

atmosphere [62]. At the TOA, the spectral albedo is defined as [50,62]

$$R_{\bar{v}, \text{TOA}} = \frac{2\pi \int_{0}^{1} I_{\bar{v}}(z_{\infty}, -\mu) \, \mu \, d\mu}{\mu_{0} I_{v, \infty}}. \quad (37)$$

For results presented in this paper, $F_0$ is the radiant flux density leaving the ground surface and is expressed as follows:

$$F_0 = 2\pi \int_{0}^{1} I_{\bar{v}}(z_b = 0, -\mu) \, \mu \, d\mu. \quad (38)$$

## 5. RESULTS AND DISCUSSION

In this work, the numerical solution is run on an Intel Xeon W1250P 4.10 GHz computer. The problem treated by Tapimo et al. [44] is considered to evaluate the performance of the MLD method. In the medium wave infrared (MWIR) region from 2000 to 2500 cm⁻¹, the absorption coefficients of absorbers are extracted from HITRAN 2008 using the LBL procedure. The vertical profiles of the tropical (model 1), midlatitude summer (model 2), midlatitude winter (model 3), and US standard (model 6) atmospheres are taken from the AFGL reference atmospheres. The atmosphere is subdivided into layers parallel to the ground surface with a highest discretization step, $\Delta z = 5$ km, while the spectral region is subdivided into narrow bands of width 1 cm⁻¹. The high-precision line-by-line calculations of the atmospheric transmittance are used as reference benchmark results for comparison. It should be noted that the spectral region considered is very important in remote sensing to retrieve the surface temperature, atmospheric temperature profile, concentrations of CO₂, H₂O, O₃, CO, and CH₄ [63–65], and to observe the coldest part of the atmosphere [64]. In practice, to consider the temporal variation of the concentrations of some greenhouse gases such as the CO₂ and CH₄, the mixing ratio $\chi_j$ in Eq. (3c) should be multiplied by a scale factor $f_j$ defining the ratio of the current mixing ratio (concentration) to that reported in the AFGL models. For all results presented in this work, the spectral response of the sensor is considered constant and it is assumed that $f_j = 1.0$.

The validation of the MLD method is first presented. For 5 km thick atmospheres, the maximum spectral difference between the line-by-line, and the L-distributions and modified L-distributions methods transmissivities are plotted in Fig. 2. With $10^3$ values of the variable $X$ in [0,1] (i.e., $10^3$ values of the mapping functions) and various observation angles, Fig. 2 provides maximum absolute errors less than 2.5% for the LD method while errors provided for the MLD formulation are below 1%. Thus, the MLD method accuracy is quite closer to the LBL method than the LD method. The computational times are 6 s for the LD method and 0.2 s for the MLD method. The results of the MLD method indicate that a high accuracy is obtained with high speed (low CPU cost). As the MLD method is also memory saving, one can conclude that a balance is established between the accuracy and the speed of the method. It was indicated that the size of the $G_r$ look-up tables affects the accuracy of the LD method [44]. Its effect on MLD results is examined in Fig. 3. Figure 3 shows that, $2 \times 10^3$ values of $G_r$ can yield absolute errors of about 0.5% for 5 km thick atmospheres. It should be noted that increasing the size of $G_r$ to
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Maximum absolute error of L-distributions (left) and modified L-distributions (right) methods on atmospheric transmittance of the terrestrial radiation at 5 km above the ground surface. Various observation angular directions $\mu = \cos \theta$ are considered. Also, $10^3$ values of the variable $X$ in [0,1] are considered. The angle $\theta$ is in degrees.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Maximum absolute error of modified L-distributions method on atmospheric transmittance of the terrestrial radiation at 5 km (left) and 105 km (right) above the ground surface. Various observation angular directions $\mu = \cos \theta$ are considered. Also, $2 \times 10^3$ values of the variable $X$ in [0,1] are considered. The angle $\theta$ is in degrees.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Spectral atmospheric transmissivity of radiation at nadir for US standard atmosphere and for four columns' atmosphere just above (a) the ground, (b) 20 km, (c) 50 km, and (d) 85 km. $2 \times 10^3$ values of the variable $X$ in [0,1] are considered, and the total thickness of the atmosphere is 105 km.</figcaption></figure>

$2 \times 10^3$ does not increase the computational time. For 105 km thick atmospheres, Fig. 3 also shows that maximum spectral differences between the LBL and MLD methods are below 1.1%. Solution verification indicates that the MLD method is well implemented for all atmospheric models considered and accounts appropriately for the changes in altitudes and angular directions.

For further results, only the US standard atmosphere with a thickness of 105 km will be studied. Also, $2 \times 10^3$ values of the variable $X$ will be considered. Figure 4 presents the spectral atmospheric transmissivity of radiation at nadir for the US standard atmosphere. Figures 4(b)–4(d) show the fraction of the transmitted radiation that a space-based sensor flying, respectively, at 20, 50, and 85 km can measure from above. Since the
<table><caption>**Table 1.** Parameters *k*<sub>iso</sub>, *k*<sub>vol</sub>, and *k*<sub>geo</sub> for the BRDF [57]</caption><thead><tr><th>Locations</th><th>k<sub>iso</sub></th><th>k<sub>vol</sub></th><th>k<sub>geo</sub></th><th>Surface Types</th></tr></thead><tbody><tr><td>A (32.72°E, 30.27°N)</td><td>0.0945</td><td>-0.1699</td><td>0.0274</td><td>Bare soil</td></tr><tr><td>B (34.31°E, 31.24°N)</td><td>0.0034</td><td>-0.0034</td><td>-0.0574</td><td>Open shrubland</td></tr><tr><td>C (33.86°E, 28.78°N)</td><td>0.0450</td><td>-0.1474</td><td>0.0312</td><td>Barren</td></tr><tr><td>B (31.23°E, 30.61°N)</td><td>0.0187</td><td>-0.1351</td><td>0.0157</td><td>Croplands</td></tr></tbody></table>

atmosphere is limited at the top by a vacuum, Fig. 4(a) shows the fraction of the transmitted terrestrial radiation that a space-based sensor flying above 105 km can measure from below. The MLD method being reversible (independent of the direction of propagation along a non-uniform path), Fig. 4(a) also shows the fraction of the transmitted extraterrestrial radiation that a ground-based sensor can measure from above. From Figs. 4(c) and 4(d), most radiation leaving the top of stratopause escapes into space in this portion of the spectrum. From Figs. 4(a) and 4(b), trace gases in the lower atmosphere are warmed up by the radiation coming from the ground surface.

To extend the MLD method to the radiance calculations, four locations namely A, B, C, and D are studied and the kernel BRDF model of the surface is considered. The surface configurations of these locations are the bare soil, open shrubland, barren or sparsely vegetated, and croplands, respectively. The locations were already studied by Tang *et al.* [57]. The parameters *k*<sub>iso</sub>, *k*<sub>vol</sub>, and *k*<sub>geo</sub> for the BRDF where estimated from MODIS MWIR data in the period from July 12 to July 30, 2005 [57]. These parameters as well as locations and surface types are reported in Table 1. Figure 5 gives the directional emissivity of these surfaces. The viewing zenith angle of each location varies significantly from 0° to 60° [57]. In plotting Fig. 5, it was assumed that the spherical crowns are separated from the ground by half their diameter [58], and that the BRDF shapes in the visible and near infrared spectral regions are the same as the ones in MWIR regions [57,66]. In Fig. 5, the numerical integration of the directional emissivity is also compared to the fitting point method presented by Tang *et al.* [57] and a good agreement is found. Figure 5 indicates that the light reflects at the Earth's surface depending on the surface emissivity. The amount of the emitted radiation increases as the viewing zenith angle increases.

The first scenario examined to compute the radiance captured by the sensor consists of Lambertian ground surface with emissivity of 0.87. The surface temperature is $T_0 = 288.2$ K.

<figure><img src="image_8.png" /><figcaption>**Fig. 5.** Directional emissivity [Eq. (32)] for bare soil, open shrubland, barren or sparsely vegetated and croplands' surface types. The surface types correspond to locations A, B, C, and D, respectively (see Table 1). The fitting point method is described in Tang *et al.* [57]. The angle $\theta$ is in degrees.</figcaption></figure>

In the nighttime, the TOA is not exposed to an external radiation and $I_{v,\infty} = 0$, while in the daytime, the Sun is at nadir. This case corresponds to the nadir observation of location B and is investigated to validate the TOA radiance of the MLD method against MATISSE output. Figures 6(a) and 6(b) show the at-sensor radiance calculated using the MLD method in the nighttime and daytime, respectively. In Figs. 6(a) and 6(b), the radiance emitted by the surface and the total radiance leaving the ground are also represented.

In the nighttime, the radiance emitted by the surface and the total radiance leaving the surface are merely the same, and the total spectral radiance escaped from the absorbing atmosphere is inferior to the emitted surface radiance that will be escaping from the atmosphere if this latter was transparent. In other words, the terrestrial infrared radiance measured by a sensor at an altitude higher than 105 km is influenced by the atmosphere and a limited amount of this radiation escapes out into space. During the day, the contribution of the solar radiation attenuated in the two-way paths is very significant and causes the total radiance that escapes from the atmosphere to be somewhat superior to the total radiance leaving the surface for some wavelengths. From Eq. (27), it is clear that the difference $I_v(0, -1) - \varepsilon B_v[T(0)]$ is proportional to the ground BRDF indicating that the BRDF can be directly determined from the sensor observations if the ground radiance is known. In both cases, a small band of atmospheric infrared window can be observed where terrestrial radiation reaches the space without significant losses. Figures 6(c) and 6(d) compare the spectral radiance of the MLD method with MATISSE outputs during the nighttime and daytime, respectively. The MLD method results agree well with MATISSE ones. Once more, the plots of the surface emitted radiance and the total radiance leaving the surface indicate that the approximation used in Eq. (32) on the spectral correlation is checked in this spectral interval of study.

The second scenario examined to compute the radiance captured by the sensor is the case of anisotropic surface reflectance. Note that the anisotropic surfaces are surfaces, which scatter light depending on the angle they are viewing. They are described in this work using the kernel BRDF. Four locations A, B, C, and D, where the surface configurations correspond to the bare soil, open shrubland, barren or sparsely vegetated, and croplands, respectively, are considered. The parameters *k*<sub>iso</sub>, *k*<sub>vol</sub>, and *k*<sub>geo</sub> used in Eq. (36) for these locations are shown in Table 1. The surfaces are at the temperature $T_0 = 288.2$ K with the directional emissivities represented in Fig. 5. The viewing zenith angle varies from 0° to 60°. In the nighttime, the TOA is not exposed to an external radiation and $I_{v,\infty} = 0$, while in the daytime, the Sun is at nadir. Figure 7 shows the TOA radiant flux density of each location measured by the sensor at the night. The total terrestrial infrared irradiance is also represented and results reveal that gases prevent the terrestrial infrared radiation in the nighttime to escape out into space and the equilibrium energy of the Earth atmosphere is maintained. Figure 8 shows the TOA radiant flux density of each location measured by the sensor during the day. A small unbalance of the energy observed for some wavenumbers diminishes increasingly for very high directional surface emissivity where the contribution of the downwelling radiance reflected at the surface is small. That is,
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** TOA radiance of the MLD method (a) nighttime and (b) daytime. Comparison of the TOA radiance of the MLD method with MATISSE (c) nighttime and (d) daytime. The observer is at nadir and the Sun is at noon for the daytime. The ground surface is isotropic with emissivity $ε = 0.87$. The temperature $T_0 = 288.2$ K and $εB_v(T_0)$ is the spectral radiance emitted by the ground surface. Moreover, $I_v(0, μ = -1)$ is the spectral radiance leaving the ground surface. The radiance is expressed in W/(m² cm⁻¹ sr).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Nighttime TOA radiant flux density for US standard atmosphere and for bare soil, open shrubland, barren or sparsely vegetated and crop-lands surface types. The surface types correspond to locations A, B, C, and D, respectively (see Table 1). The TOA incident radiance is $I_{v,∞} = 0$. Note that $F_0$ given by Eq. (39) is the radiant flux density leaving the surface considered (terrestrial irradiance). The surface temperature is $T_0 = 288.2$ K. The radiance flux is expressed in W/(m² cm⁻¹).</figcaption></figure>

the energy balance is also held for surfaces with directional emissivity close unity while the energy is unbalance for lowest surface directional emissivity. This unbalance of the energy is characteristics of the greenhouse effect. As the spectral albedo is crucial to control the energy budget of the Earth, surface temperature

and the amount of CO₂ in the atmosphere, it should also be accurately estimated. Figure 9 shows the TOA spectral albedo for location A. Figures 7 and 8 show a small band of atmospheric infrared window where terrestrial irradiance reaches the space without significant losses. In the spectral region where the
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Daytime TOA radiant flux density for US standard atmosphere and for bare soil, open shrubland, barren or sparsely vegetated, and crop-lands' surface types. The surface types correspond to locations A, B, C, and D, respectively (see Table 1). The TOA is illuminated in angular direction $\mu_0 = 1.0$. Note that $F_0$ given by Eq. (39) is the radiant flux density leaving the surface considered (terrestrial irradiance). The surface temperature is $T_0 = 288.2$ K. The radiant flux density is expressed in W/(m$^2$ cm$^{-1}$).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Spectral albedo at the TOA for US standard atmosphere and for bare soil surface type (location A) with anisotropic surface reflectance. The illumination solar direction is $\mu_0 = 1.0$.</figcaption></figure>

atmospheric absorption from CO$_2$ is high, the sensor irradiance presents large departure from the surface irradiance. Thus, the decrease in albedo increases the surface temperature and so the concentration of the atmospheric CO$_2$, which is essential for evaluating economic development, and understanding climate change and global warming. In overall, a poor knowledge of surface emissivity and atmospheric radiation causes enormous error in the determination of surface temperature, atmospheric temperature profile and gas concentration from remote sensing measurements.

The performance of the MLD been investigated, it is of great interest to check the validity of the MLD to a wide spectral range for use in most remote sensing problems. In the spectral range from 1.1 to 50 µm, the absorption coefficients of absorbers are extracted from HITRAN 2012 using the LBL procedure.

The US standard atmosphere vertical profile is taken from the AFGL reference atmospheres. As in the case treated to study the performance of the method, the atmosphere of 105 km thick is subdivided in layers parallel to the ground surface with a highest discretization step $\Delta z = 5$ km while the spectral interval is subdivided in narrow bands of width 1 cm$^{-1}$. The spectral radiance reaching the sensor aperture is represented for the nighttime and daytime, and a comparison with MATISSE results is shown. Moreover, the total radiance leaving the ground is also represented. The effect of the large view angle on the radiance reaching the sensor is analyzed. For a sake of clarity, the spectral radiance for the day time scenario is represented separately from 1.1 to 3 µm and from 3 to 50 µm. Figure 10 shows the spectral radiance during the nighttime, while Figs. 11 and 12 show the daytime radiance in 1.1–3 µm and 3–50 µm spectral ranges, respectively. Results clearly demonstrate that the MLD method results agree well with MATISSE ones. However, a slight difference between the two models, which appears in the range from 1.1 to 2.5 µm is due to the fact that MATISSE considers stratospheric aerosols that are not included in the MLD model. In Figs. 10, 11, and 12, the sensor view’s zenith angles of 0° and 60° are considered. It should be noted that MATISSE is a radiative transfer code developed at ONERA and does not provide results beyond 14 µm. It is important to note that $I_v(0, \mu = -1) = I_v(0, \mu = -0.5)$ in Figs. 10 and 12. This result is justified by the fact that the Lambertian ground surface is considered, and the radiation scattered on the ground surface is the same in all angular directions. For large view zenith angles, the sensor radiation is attenuated due to the long path traveled in the atmosphere.
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Nighttime TOA radiance of the MLD method with sensor view zenith angles 0° and 60°. $I_v(0, \mu = -1)$ and $I_v(0, \mu = -0.5)$ are the spectral radiances leaving the ground surface in the angular directions $-\mu = 1$ and $-\mu = 0.5$, respectively. The radiance is expressed in W/(m²·µm·sr).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 11.** Daytime TOA radiance of the MLD method with sensor view zenith angles 0° and 60°. The spectral interval is from 1.1 to 3 µm and the Sun is at noon. The radiance is expressed in W/(m²·µm·sr).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Daytime TOA radiance of the MLD method with sensor view zenith angles 0° and 60°. The spectral interval is from 3 to 50 µm and the Sun is at noon. $I_v(0, \mu = -1)$ and $I_v(0, \mu = -0.5)$ are the spectral radiances leaving the ground surface in the angular directions $-\mu = 1$ and $-\mu = 0.5$, respectively. The radiance is expressed in W/(m²·µm·sr).</figcaption></figure>

# 6. SUMMARY AND CONCLUSION

The objectives of this work were to overcome the limitations of the L-distributions (LD) method and establish a balance between the accuracy and speed. A modified L-distributions (MLD) method was formulated. It was shown that the MLD is computationally efficient, memory saving and reversible (independent on the direction of propagation along a non-uniform path). In the medium wave infrared (MWIR) region from 2000 to 2500 cm⁻¹, the accuracy and the speed of the method were tested on atmospheric transmittance for 5 km thick atmospheres. Keeping the size of the variable $X \in [0, 1]$ to 10³, results indicated that an absolute error less than 2.5% could be reached in 6 s for the LD method. Using the MLD method, the absolute error and CPU time drop below 1.0% and 0.2 s, respectively. How the size of the parameter $X$ affects the accuracy of the MLD method was also investigated. With $2 \times 10^3$ values of this parameter, the absolute errors are below 0.5% and 1.1% for 5 km and 105 km thick atmospheres, respectively.

Also, the computational time was 0.2 s. From these results, a balance is established between the accuracy and the speed of the MLD method as the results yield high accuracy with increase speed (small CPU cost). The US Standard atmosphere with a thickness of 105 km was considered to investigate the TOA spectral radiance and irradiance of four locations A, B, C, and D. These locations correspond to bare soil, open shrubland, barren or sparsely vegetated, and croplands surface cover types, respectively. For each location, the directional emissivity was numerically calculated from the kernel bidirectional reflectance distribution function. Results were compared to the approximated analytical method described by Tang *et al.* [57] and good agreement was shown. The nighttime and daytime scenarios with the atmosphere limited at the ground by a Lambertian surface were considered for a comparison of the MLD calculations with MATISSE outputs. At nadir observation, the TOA radiance for location B reveals that the MLD method results match well with the MATISSE outputs. The nighttime and daytime scenarios with the atmosphere limited at the ground by an anisotropic surface was also investigated for viewing angles ranging from 0° to 60°. The TOA irradiance calculation demonstrated that in the nighttime, the radiant flux density captured by the sensor is inferior to the radiant flux density leaving the surface and that would be sensed if the atmosphere was transparent. Thus, in the nighttime, gases maintain the energy equilibrium of the atmosphere by preventing infrared radiation to escape out into space. In the daytime, this energy balance is also held for surfaces with directional emissivity close unity while a small energy imbalance can be observed for lowest surface directional emissivity at some wavenumbers. This imbalance is characteristic of the greenhouse effect. For most remote sensing applications, a wide spectral range that extended from 1.1 to 50 µm was also considered and the effect of large view angle on the results was examined. Once more, the MLD results were compared with MATISSE in the spectral interval 1.1–14 µm and good agreements were obtained. It should also be noted that from albedo measurement data, the BRDF of the
surface can be retrieved and so the directional emissivity. The comparison of the retrieved directional emissivity to Eqs. (34) and (36) helps to determine parameters $k_{iso}$, $k_{vol}$, and $k_{geo}$ for specific land types. Since this work was focused on the solution of forward radiative transfer method for use in remote sensing, a comparison of the present results with data of high-quality spectral radiation measurements such as satellite, aircraft and ground-based instruments will be investigated in the future. The atmospheric temperature profiles as well as greenhouse gas concentrations will be retrieved for the climate change study.

## APPENDIX A: ALGORITHM FOR COMPUTING THE MLD METHOD

The modified L-distributions method developed in this paper computes the atmospheric transmittance in a few steps.

Step 1. Divide the inhomogeneous atmosphere into *M* layers of constant physical properties $T_m$, $p_m$, and $\kappa_v^m$, and with the discretization step $L_m$ ($1 \le m \le M$).

Step 2. Divide the spectral interval into *N* narrow bands with the spectral resolution $\Delta v_i$ ($1 \le i \le N$).

Step 3. Find the thickness of the equivalent homogeneous atmosphere $\Delta z = \max(L_1, L_2, \dots, L_M)$.

Step 4. Calculate the absorption coefficient of the equivalent homogeneous atmosphere $\bar{\kappa}_v$ using Eq. (20a).

Step 5. For each narrow band, compute the MLD parameters $k_P$ and $\beta$ using Eqs. (8), (9), and (10a), where $\kappa_v$ is replaced by $\bar{\kappa}_v$ and store the results in a table.

Step 6. In each narrow band and for the parameter X ranging in [0,1], compute the mapping function $G_r(X)$ using Eq. (13), where $\kappa_v$ is replaced by $\bar{\kappa}_v$ and store the results in a table.

Step 7. Compute the transmissivity $T_{r,\bar{v}}(L_\mu)$ with $L_\mu = \Delta z/\mu$ and $\mu$ the view directional angle.

**Funding.** Agence Nationale de la Recherche (ANR-20-CE04-0005 ASGARD).

**Disclosures.** The authors declare no conflicts of interest.

**Data availability.** Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

## REFERENCES

1. X. Calbet, R. Kivi, S. Tjemkes, et al., "Matching radiative transfer models and radiosonde data from the EPS/Metop Sodankyla campaign to IASI measurements," Atmos. Meas. Tech 4, 1177-1189 (2011).

2. K. Salesin, K. D. Knobelspiesse, J. Chowdhary, et al., "Unifying radiative transfer models in computer graphics and remote sensing, Part I: a survey," J. Quant. Spectrosc. Radiat. Transf. 314, 108847 (2024).

3. J. S. Delamere, S. A. Clough, V. H. Payne, et al., "A far-infrared radiative closure study in the Arctic: application to water vapor," J. Geophys. Res. 115, D17106 (2010).

4. R. Tapimo, H. T. T. Kamdem, and D. Yemele, "A discrete spherical harmonics method for radiative transfer analysis in inhomogeneous polarized planar atmosphere," Astrophys. Space Sci. 363, 52 (2018).

5. C. Zhan and S. Liang, "Improved estimation of the global top-of-atmosphere albedo from AVHRR data," Remote Sens. Environ. 269, 112836 (2022).

6. E. J. Mlawer, S. J. Taubman, P. D. Brown, et al., "Radiative transfer for inhomogeneous atmospheres: RRTM, a validated correlated-k model for the longwave," J. Geophys. Res. 102, 16663-16682 (1997).

7. F. A. M. Garzón and I. Valánszki, "Environmental armed conflict assessment using satellite imagery," J. Environ. Geogr. 13, 1-14 (2020).

8. Q. Huang, G. Jin, X. Xiong, et al., "Monitoring urban change in conflict from the perspective of optical and SAR satellites: the case of Mariupol, a City in the conflict between RUS and UKR," Remote Sens. 15, 3096 (2023).

9. J. R. B. Fisher, E. A. Acosta, P. J. Dennedy-Frank, et al., "Impact of satellite imagery spatial resolution on land use classification accuracy and modeled water quality," Remote Sens. Ecol. Conserv. 4, 137-149 (2018).

10. R. Tapimo, C. C. Atemkeng, H. T. T. Kamdem, et al., "Bidirectional transmittance and reflectance models for soil signature analysis," Appl. Opt. 58, 1924-1932 (2019).

11. C. C. Atemkeng, R. Tapimo, H. T. T. Kamdem, et al., "Radiative transfer technique for retrieving the radiative properties of agricultural soils," J. Quant. Spectrosc. Radiat. Transf. 269, 107698 (2021).

12. Z. Qin and A. Karnieli, "Progress in the remote sensing of land surface temperature and ground emissivity using NOAA-AVHRR data," Int. J. Remote Sens. 20, 2367-2393 (1999).

13. D. P. Edwards, GENLN2: A General Line-by-line Atmospheric Transmittance and Radiance Model, Version 3.0 Description and Users Guide, NCAR/TN-367-STR (National Center for Atmospheric Research, 1992).

14. K. Caillaut, S. Fauqueux, C. Bourlier, et al., "Multiresolution optical characteristics of rough sea surface in the infrared," Appl. Opt. 46, 5471-5481 (2007).

15. A. Fokou, R. Tchinda, G. L. Ymeli, et al., "Estimating the radiance in atmosphere-ocean systems with different atmosphere models by discrete spherical harmonics method," ASME J. Heat Mass Transf. 145, 052802 (2023).

16. A. Fokou, R. Tapimo, G. L. Ymeli, et al., "Radiation distribution in inhomogeneous atmosphere-ocean system by discrete spherical harmonics method," J. Quant. Spectrosc. Radiat. Transf. 270, 107707 (2021).

17. A. del Aguila, D. S. Efremenko, and T. Trautmann, "A review of dimensionality reduction techniques for processing hyper-spectral optical signal," Light Eng. 27, 85-98 (2019).

18. M. Chi, A. Plaza, J. A. Benediktsson, et al., "Big data for remote sensing: challenges and opportunities," Proc. IEEE 104, 2207-2219 (2016).

19. R. Kosarevych, O. Lutsyk, B. Rusyn, et al., "Spatial point patterns generation on remote sensing data using convolutional neural networks with further statistical analysis," Sci. Rep. 12, 14341-14349 (2022).

20. Y. Ma, H. Wu, L. Wang, et al., "Remote sensing big data computing: challenges and opportunities," Future Gener. Comput. Syst. 51, 47-60 (2015).

21. F. André, C. Delage, L. Guilmard, et al., "Bridging physics and statistical learning methodologies for the accurate modeling of the radiative properties of non-uniform atmospheric paths," J. Quant. Spectrosc. Radiat. Transf. 320, 108961 (2024).

22. J. Vicent, J. Verrelst, N. Sabater, et al., "Comparative analysis of atmospheric radiative transfer models using the atmospheric look-up table generator (ALG) toolbox (version 2.0)," Geosci. Model Dev. 13, 1945-1957 (2020).

23. X. Zeng, C. Liang, L. Duan, et al., "A GPU-based line-by-line method for thermal radiation transfer of H2O, CO2, and H2O/CO2 mixture," Appl. Therm. Eng. 167, 114799 (2020).

24. I. E. Gordon, L. S. Rothman, R. J. Hargreaves, et al., "The HITRAN2020 molecular spectroscopic database," J. Quant. Spectrosc. Radiat. Transf. 277, 107830 (2022).

25. T. Delahaye, R. Armante, N. A. Scott, et al., "The 2020 edition of the GEISA spectroscopic database," J. Mol. Spectrosc. 380, 111510 (2021).

26. Y. Wang, Y. Zhao, W. Li, et al., "Using a GPU to accelerate a longwave radiative transfer model with efficient CUDA-based methods," Appl. Sci. 9, 4039 (2019).
27. S. A. Clough, M. W. Shephard, E. J. Mlawer, et al., "Atmospheric radiative transfer modeling: a summary of the AER codes," *J. Quant. Spectrosc. Radiat. Transf.* **91**, 233-244 (2005).

28. M. J. Alvarado, V. H. Payne, E. J. Mlawer, et al., "Performance of the line-by-line radiative transfer model (LBLRTM) for temperature, water vapor, and trace gas retrievals: recent updates evaluated with IASI case studies," *Atmos. Chem. Phys.* **13**, 6687-6711 (2013).

29. N. A. Scott and A. Chedin, "A fast line-by-line method for atmospheric absorption computations: the automatized atmospheric absorption atlas," *J. Appl. Meteorol. Climatol.* **20**, 802-812 (1981).

30. E. Pannier and C. Laux, "RADIS: a nonequilibrium line-by-line radiative code for CO₂ and HITRAN-like database species," *J. Quant. Spectrosc. Radiat. Transf.* **222-223**, 12-25 (2019).

31. I. M. Chapman, D. A. Naylor, B. G. Gom, et al., "BTRAM: an interactive atmospheric radiative transfer model," in *Bridging Excellence - The 30th Canadian Symposium on Remote Sensing* (2009), pp. 22-25.

32. A. Dudhia, "The reference forward model (RFM)," *J. Quant. Spectrosc. Radiat. Transf.* **186**, 243-253 (2017).

33. L. Croize, J. Pierro, L. Labarre, et al., "A wideband and narrow band radiative transfer model based on a fast-correlated k-distributions generation for MATISSE V3.0," *20th EGU General Assembly, EGU2018* Vienna, Austria, 4-13 April, 2018, paper 13910.

34. A. Berk, L. S. Bernstein, G. P. Anderson, et al., "MODTRAN cloud and multiple scattering upgrades with application to AVIRIS," *Remote Sens. Environ.* **65**, 367-375 (1998).

35. A. Berk and F. Hawes, "Validation of MODTRAN®6 and its line-by-line algorithm," *J. Quant. Spectrosc. Radiat. Transf.* **203**, 542-556 (2017).

36. L. L. Gordley, B. T. Marshall, and D. A. Chu, "LINEPAK: algorithms for modeling spectral transmittance and radiance," *J. Quant. Spectrosc. Radiat. Transf.* **52**, 563-580 (1994).

37. P. Ricchiazzi, S. Yang, C. Gautier, et al., "SBDART: a research and teaching software tool for plane-parallel radiative transfer in the Earth's atmosphere," *Bull. Am. Meteorol. Soc.* **79**, 2101-2114 (1998).

38. C. Wang, B. He, and M. F. Modest, "Full-spectrum correlated-k-distribution look-up table for radiative transfer in nonhomogeneous participating media with gas-particle mixtures," *Int. J. Heat Mass Transf.* **137**, 1053-1063 (2019).

39. H. Li, J. Zhang, Y. Cheng, et al., "Calculations of radiative intensity in one-dimensional gaseous media with black boundaries using the statistical narrow band model," *J. Quant. Spectrosc. Radiat. Transf.* **240**, 106691 (2020).

40. M.-D. Chou, J. C.-C. Yu, W.-L. Lee, et al., "A new k-distribution scheme for clear-sky radiative transfer calculations in Earth's atmosphere. Part I: thermal infrared (longwave) radiation," *J. Atmos. Sci.* **77**, 2237-2256 (2020).

41. Y. M. Govaerts, Y. Nollet, and V. Leroy, "Radiative transfer model comparison with satellite observations over CEOS calibration site Libya-4," *Atmosphere* **13**, 1759 (2022).

42. P. F. Baumeister and L. Hoffmann, "Fast infrared radiative transfer calculations using graphics processing units: JURASSIC-GPU v2.0," *Geosci. Model Dev.* **15**, 1855-1874 (2022).

43. F. André, "The *l*-distribution method for modeling non-gray absorption in uniform and non-uniform gaseous media," *J. Quant. Spectrosc. Radiat. Transf.* **179**, 19-32 (2016).

44. R. Tapimo, L. Croize, and F. André, "Fast and accurate infrared radiative transfer calculation method in gaseous atmospheres," in *Proceedings of the IEEE IGARSS (2024)*, pp. 6993-6997.

45. W. Li, F. Zhang, Y.-N. Shi, et al., "Efficient radiative transfer model for thermal infrared brightness temperature simulation in cloudy atmospheres," *Opt. Express* **28**, 25730-25749 (2020).

46. F. André, C. Cornet, M. Galtier, et al., "Radiative transfer in the O₂ A-band - a fast and accurate forward model based on the *l*-distribution approach," *J. Quant. Spectrosc. Radiat. Transf.* **260**, 107470 (2021).

47. F. Becker and Z.-L. Li, "Temperature-independent spectral indices in thermal infrared bands," *Remote Sens. Environ.* **32**, 17-33 (1990).

48. A. Mossi, M. M. Galarça, R. Brittes, et al., "Comparison of spectral models in the computation of radiative heat transfer in participating media composed of gases and soot," *J. Braz. Soc. Mech. Sci. Eng.* **34**, 112-119 (2012).

49. K. N. Liou, *An Introduction to Atmospheric Radiation* (Elsevier, 2002).

50. R. Tapimo, H. T. T. Kamdem, and D. Yemele, "Discrete spherical harmonics method for radiative transfer in scalar planar inhomogeneous atmosphere," *J. Opt. Soc. Am. A* **35**, 1081-1090 (2018).

51. R. Tapimo, M. Lazard, G. L. Ymeli, et al., "Radiative transfer model for ground surface irradiance estimation: clear sky," *J. Opt. Soc. Am. A* **38**, 1640-1646 (2021).

52. M. J. Mumma, K. Fox, and J. Hornstein, *Vibrational-Rotational Spectroscopy for Planetary Atmospheres I* (NASA, 1982).

53. A. S. Fisher and S. L. Rani, "A narrow band model based on the absorption coefficient and its application to the calculation of radiative transfer in one-dimensional enclosures," *J. Quant. Spectrosc. Radiat. Transf.* **277**, 107989 (2022).

54. Z. Wan and Z.-L. Li, "Physics-based algorithm for retrieving land-surface emissivity and temperature from EOS/MODIS data," *IEEE Trans. Geosci. Remote Sens.* **35**, 980-996 (1997).

55. B. Tang and Z.-L. Li, "Retrieval of land surface bidirectional reflectivity in the mid-infrared from MODIS channels 22 and 23," *Int. J. Remote Sens.* **29**, 4907-4925 (2008).

56. G. Hulley and S. Hook, *HyspIR Level-2 TIR Surface Radiance Algorithm Theoretical Basis Document* (NASA, 2011).

57. B.-H. Tang, Z.-L. Li, and Y. Bi, "Estimation of land surface directional emissivity in mid-infrared channel around 4.0 µm from MODIS data," *Opt. Express* **17**, 3173-3182 (2009).

58. G.-M. Jiang and Z.-L. Li, "Intercomparison of two BRDF models in the estimation of the directional emissivity in MIR channel from MSG1-SEVIRI data," *Opt. Express* **16**, 19310-19321 (2008).

59. J. L. Roujean, M. Leroy, and P. Y. Deschamps, "A bidirectional reflectance model of the Earth's surface for the correction of remote sensing data," *J. Geophys. Res. Atmos.* **97**, 20455-20468 (1992).

60. S. Y. Kotchenova and E. F. Vermote, "Validation of a vector version of the 6S radiative transfer code for atmospheric correction of satellite data. Part II. Homogeneous Lambertian and anisotropic surfaces," *Appl. Opt.* **46**, 4455-4464 (2007).

61. G. Picard, M. Dumont, M. Lamare, et al., "Spectral albedo measurements over snow-covered slopes: theory and slope effect corrections," *Cryosphere* **14**, 1497-1517 (2020).

62. M. Bilal, M. Nazeer, J. E. Nichol, et al., "A simplified and robust surface reflectance estimation method (SREM) for use over diverse land surfaces using multi-sensor data," *Remote Sens.* **11**, 1344-1367 (2019).

63. X. Xiong and W. Barnes, "An overview of MODIS radiometric calibration and characterization," *Adv. Atmos. Sci.* **23**, 69-79 (2006).

64. M. W. Shephard, S. A. Clough, V. H. Payne, et al., "Performance of the line-by-line radiative transfer model (LBLRTM) for temperature and species retrievals: IASI case studies from JAIVEx," *Atmos. Chem. Phys.* **9**, 7397-7417 (2009).

65. A. Wu, X. Xiong, D. R. Doelling, et al., "Characterization of terra and aqua MODIS VIS, NIR, and SWIR spectral bands' calibration stability," *IEEE Trans. Geosci. Remote Sens.* **51**, 4330-4338 (2013).

66. D. S. Kimes, "Dynamics of directional reflectance factor distributions for vegetation canopies," *Appl. Opt.* **22**, 1364-1372 (1983).