Analysis of the Discrete Theory of Radiative Transfer in the Coupled “Ocean–Atmosphere” System: Current Status, Problems and Development Prospects

: In this paper, we analyze the current state of the discrete theory of radiative transfer. One-dimensional, three-dimensional and stochastic radiative transfer models are considered. It is shown that the discrete theory provides a unique solution to the one-dimensional radiative transfer equation. All approximate solution techniques based on the discrete ordinate formalism can be derived based on the synthetic iterations, the small-angle approximation, and the matrix operator method. The possible directions for the perspective development of radiative transfer are outlined.


Introduction
Radiative transfer theory is the principal method for modeling radiation propagation in the atmosphere and the ocean in the photometric ray approximation [1,2]. In this approximation, the radiation field is decomposed into a coherent part, which determines the optical characteristics of the medium, and an incoherent one, which is related to the processes of multiple light scattering and satisfies the radiative transfer equation (RTE).
The propagation of radiation in the medium is associated with numerous physical phenomena described by various physical theories. In this article, we will consider a physical model of a continuous medium in which light is absorbed and in which discrete scatterers are located. The spatial distribution of the medium and scatterers in space is generally arbitrary. This paper is devoted to numerical methods for solving the RTE in the discrete ordinate space. Statistical modeling techniques, such as the Monte Carlo method, follow their own logic which drops out from the presented scheme.
As will be shown in this paper, all numerical methods are based on one or another way of replacing the scattering integral with a finite sum, which replaces the desired continuous brightness distribution with discrete values or a set of coefficients for the expansion of this distribution over a system of functions. Thus, the transfer equation, its solutions, and all the corollaries from it acquire a discrete matrix form. In this case, the approximation is only the replacement of the integral by the sum, and all other conclusions can be made strictly analytically. In fact, we can talk about a discrete transfer theory. Kolmogorov A.N. [3] pointed out that with the development of modern computer technology in many cases, it is reasonable to study real phenomena, avoiding the intermediate stage of their stylization in the spirit of representations of mathematics of the infinite and continuous, moving directly to discrete models.
Today the radiative transfer theory is not only an integral branch of atmospheric physics but also a cornerstone of Earth observation technology. An analysis of modern literature over the past ten years on radiative transfer reveals a significant predominance of publications, in which different radiative transfer codes have been intercompared. In [4], the results of different radiative transfer codes agree up to the 5th digit. Such an agreement cannot be accidental and suggests that all codes presumably exploit the same method for solving the RTE.
In this regard, the state of the radiative transfer theory is very similar to the state of physics by the end of the 19th century, which was analyzed by W. Thomson (Lord Kelvin) [5]. He not only reviewed the state of physics at the turn of the century but most importantly, he highlighted two clouds (problems) that "obscure the beauty and clearness of the whole theory". These clouds were (a) the independence of the speed of light on the reference frame in the wave equation resulting from Maxwell's equations and (b) the inability to describe the dependence of blackbody radiation in the framework of classical thermodynamics. Later, these "clouds" led to the creation of quantum physics, which radically changed the whole physical paradigm of the world. Our goal in this paper is to analyze the current state of the discrete theory of radiative transfer and find its main problems, i.e., "Thomson clouds", which will determine its future development.

Boundary Value Problems of the Radiation Transfer Theory
Since the Earth radius is two orders of magnitude larger than the scales of the troposphere (which contains 75% of the atmosphere's mass) and the ocean, most radiative transfer models used in Earth remote sensing are one-dimensional. The atmosphere is treated as a set of homogeneous slabs. The horizontal inhomogeneity is considered by using the so-called independent pixel approximation (IPA) [6]. However, one-dimensional models cannot be applied at large viewing or solar zenith angles, as well as for limb observations. Moreover, in the ocean-atmosphere systems, the point sources are practically important, e.g., beacons, searchlights and lasers, which are widely used in navigation and communication systems. Cloud and surface inhomogeneities can lead to a bias in the radiances computed within one-dimensional radiative transfer models [7,8]. Such cases are the subjects for three-dimensional radiative transfer models.
Another important type of problem, which is not covered by traditional one-dimensional models is radiation transfer in stochastic media. Usually, because of mathematical difficulties, the excitements or stochastic processes in the medium are uncorrelated. Such simplifications do not consider many critical phenomena, such as stochastic lenses, which under certain circumstances, allow astronauts to see grains of sand on the ocean floor [9].

Architecture of the Boundary Value Problems in Radiative Transfer
We consider the slab of turbid medium illuminated by the plane-parallel solar beam. The direction of incidence is given by μ is the solar zenith angle cosine. The boundary value problem (BVP) reads as follows: where L(τ,μ,ϕ) is the radiance field in the viewing direction τ is the optical thickness of the slab, ˆ( , ) x ′ l l is the single scattering phase function, Λ is the single scattering albedo, while ρ(r) is the reflection coefficient of the surface. The inelastic scattering processes changing the wavelength of the light are not included. The BVP (1) is defined in the Cartesian coordinate system OXYZ, in which the OZ axis is directed downward perpendicular to the layer boundary, ẑ is a unit vector along OZ. The upper boundary of the layer is located at 0 z = . Throughout the paper, the unit vectors are denoted with a "^" sign, while column vectors, row vectors and matrices are marked with the right arrows, left arrows and double arrows respectively.
Essentially, the BVP (1) is a set of integral equations, which consists of the RTE itself, and boundary condition, in which the reflected radiance is defined via the transmitted one. These equations are linear with respect to L(τ,μ,ϕ). Hence, the solution can be expressed as the sum of the basic solutions according to the principle of superposition. Such a technique was considered in [10] and is referred to as the method of decomposition of BVPs. This method allows us to substitute the initial problem with boundary conditions with the hierarchy of BVPs that do not contain equations in boundary conditions. Also, the fundamental problem is formulated, through which the solutions to all problems from the hierarchy can be found. Our analysis then is focused on the fundamental problem.
Following [10], the solution to Equation (1) can be expressed as follows: where the "haze" component The surface reflection can be written as: is the reflection coefficient from the surface, while S is the surface area.
The BVP (4) can be analyzed in the framework of the perturbation theory, namely: where ξ is formally a small parameter, which in the final expressions is to be set ξ = 1. Substituting Equation (6) into Equation (4) gives: where index 2 shows that the radiance values ( ) 2 n L are taken at the surface level.
By equating the terms of Equation (7) at the same powers of ξ, we obtain: The signal component (0)ˆ( , , ) ( , ) L z L z ≡ r l l depends on the coordinate r and defines the radiance of the surface, which is due to the multiple reflection events from the surface ρ as well as multiple scattering events in the medium. All Equations (8) have the same structure, namely: which is a convolution type equation. It is solved using the Fourier transform: After the Fourier transform in Equation (10), Equation (9) takes the form: We define a new function ( , , ) z ψ l ν as such that the relation holds: The corresponding BVP for it has the form: which corresponds to the RTE solution for a point diffuse (PD) source: Substituting the obtained result in Equations (10) and (12), after transformations and the inverse Fourier transform, one reaches a solution of the BVP (Equation (9)): where is the mean hemispherical albedo of the turbid medium slab.
Accordingly, following [10] the solution of the BVP (Equation (1)) for an arbitrary layer with a reflecting Lambertian bottom can be reduced to a superposition of solutions: where ( , , ) L z r l  is the radiance related to the single reflection event by the object ( ) ρ r  and multiple scattering events from the surface, while ( , , ) L z ′ r l corresponds to the nonlinear part of the radiance, which is due to multiple reflections by the object and the surface. Thus, the solution of the BVP (Equation (1)) is reduced to the solution of the BVP for a slab with the plane parallel source, the plane diffuse source, and the point diffuse source. In [11], it was shown that any BVP of radiative transfer can be derived from the problem with the point unidirectional source.
An analysis of the solution of BVPs with concentrated sources was carried out [12]. It was shown that they always come down to solving the problem for a point unidirectional (PM) source or, in special cases, to a point isotropic source (PI). In [11,12], an analysis was performed regarding the properties of solving fundamental sources and it was shown that the angular distribution of radiance of only the source PD has no singularities, all the others have: the PM delta feature in the 0th-fold scattering term, TD, PI, and PM additionally have the inverse sine of the viewing angle in the first multiplicity, the logarithmic in the second term, and PM is the inverse sine in the third term.

The Boundary Value Problem for the Radiative Transfer Equation for a Slab
Let us analyze the numerical solution of the BVP in the case of the plane parallel source: Here we assume that the slab is vertically homogeneous.
To solve Equation (17) numerically, it must be discretized in the angular domain by replacing the scattering integral with a finite sum [13]. Such substitution is impossible if there are singularities in the angular distribution of the radiance field, as the singularity cannot be represented as a finite series in any basis. In [14], it was proposed to express the total radiation as a sum of the direct radiation of the source and the diffuse part. Such representation became an indispensable starting point in most solution techniques. However, all-natural objects (aerosols, clouds, etc.) contain suspended particles, which are significantly larger than the wavelength of the incident light. In accordance with the theory of G. Mie [15], the single scattering phase functions for such particles are highly forward-peaked. Consequently, in the case of highly anisotropic phase functions, replacing the integral with a finite sum can lead to significant errors [16].

Discretization of the Radiative Transfer Equation
A more advanced approach consists of the following representation of the radiance field [13]: where the source function ( , ) Δ τ l is due to the anisotropic part of the solution, namely: Since is a smooth function, it can be expressed on a finite basis. For instance, applying the discrete ordinate method (DOM) [17] gives: where C ( ) C ( , ), C( ) C ( ), C ( ) is the column-vector of radiance values in the Fourier space along the discrete ordinate directions j j ± μ = ζ ± , while ζj are the Gaussian quadrature points of the N/2 order. For the sake of simplicity, index m is further omitted. Such discretization allows substituting the integral with a finite sum. The BVP is transformed into the matrix nonlinear differential equation of the first order with constant coefficients: where wj are the Gaussian quadrature weights, P ( ) n l μ are the associated Legendre polynomials,

Propagators and Scatters
The general solution of the inhomogeneous Equation (22) is given as the sum of the particular solution of the inhomogeneous equation and general solution of the homogeneous equation: is the solution to the homogeneous equation. It is referred to as the propagator and relates the radiance fields at spatial coordinates t and τ. Note, the propagator comprises both negative and positive exponentials, which physically corresponds to the downward and upward radiances, respectively. The positive exponentials increase the condition number of the matrix P  making the numerical procedure unstable as the optical thickness increases. To get rid of this effect, the scaling transformation is applied [18]: where Note, that the BVP (17) and the corresponding Equation (22) are the two-point boundary problems, i.e., the boundary conditions are defined at two boundaries. Hence, solution (23) expressed through the propagators is not complete [19]. (25) correspond to the radiances defined by the boundary conditions. The column-vectors correspond to the transmitted and reflected radiances. They can be found from Equation (25): The solution in the form of scatters (see Equation (27)) relates the incoming and outgoing radiances and, thus, can be considered as a generalization of the radiance coefficients. The column-vectors F  correspond to the intrinsic radiation of the layer, while the matrices R  and T  are the discrete values of reflection and transmission coefficients, respectively [13]. Equation (27) can be transformed into the propagator-like form: In [20], such a kind of transformation was referred to as the stellar product. The entries of scatter matrices satisfy the Riccati equation, which corresponds to the one-point BVP [20,21].

Invariance Property of the Solution
Consider the medium consisting of two adjacent layers: where the bottom index is the layer index (i.e., 1 corresponds to the upper layer, while 2 corresponds to the bottom layer). Note that since the layers are adjacent and the solution is continuous, we have: By expressing transmitted 1 C −  and reflected 2 C +  radiances by incident radiances 1 C ↓  and 2 C ↑  , we obtain: We can see that solution (31) for two adjacent layers is equivalent to solution (27) for a single layer with effective parameters of the medium. Namely, the invariance property of the solution is expressed in the form of scatterers, which corresponds to the discrete setting of the invariant imbedding method of V.A. Ambartsumian [22]. Indeed, let us consider a system of two layers, where the upper layer 1 is infinitely thin with an optical thickness d τ , while the lower layer 2 is semi-infinite. As layer 1 is infinitely thin, its radiance field can be derived by using the single scattering approximation. For that, we consider the Peiperl's integral equation, which is the formal solution to the radiative transfer equation, given that the scattering integral is known. For a layer illuminated by radiance ( ) L ↑ l from the top, and by ( ) L ↑ l from the bottom, the integral radiative transfer equation reads as: x The single scattering approximation leads to the following solution: where the upper equation corresponds to the downwelling radiance, and the lower equation corresponds to the upwelling radiance. Performing the integration in Equation (33), while considering the Fourier expansion and applying the discrete ordinate method, we obtain a solution for the infinitely thin layer in the form of scatters: Since layer 2 is semi-infinite, we have: Moreover, adding the infinitely thin layer 1 to the semi-infinite layer 2 does not lead to the change in the reflection, namely: and Equation (31) is transformed into the following expression: The expressions for the transmission and reflection of layer 1 include an infinitely small layer thickness dτ. By considering the Taylor expansion we obtain: Then, substituting reflection matrices corresponding to Equation (34) into Equation (37) gives: which corresponds to the expression derived by V.A. Ambartsumian [22] but in a discrete form. The invariance principle can be applied to a layer of arbitrary vertical inhomogeneity, by breaking it into an arbitrary number of homogeneous layers. In this case, two adjacent layers can be replaced by one layer described by Equation (31). This approach in the radiative transfer theory is called the matrix operator method [23] (MOM). Note that in this paper Equation (31) is derived for the arbitrary anisotropic part of the radiance field (see Equation (18)).

The Anisotropic Part of the Solution
The accuracy of solution (11) as well as the computational time is governed by the number of discrete ordinates N, the number of azimuthal modes M and the number of phase function expansion coefficients, K. Generally, M ≈ N ≈ K [13]. However, by appropriate choice of the anisotropic part in Equation (18), the regular part is close to an isotropic function, and the number of discrete ordinates and the number of azimuthal modes can be substantially reduced. In this case, K >> N >> M. Since the computational time increases M and N, the reduction of M and N leads to the performance enhancement.
A question arises: how to choose the anisotropic part in the best way? To answer to this question, we expand the radiance field in the Legendre series: Note that the spectrum ( ) Next, we consider the Taylor expansion of ( , ) Z k τ : Substituting Equation (41) into the BVP for a L and applying the spherical harmonics method one can obtain the following expression [24]: which is referred to as the small angle modification of the spherical harmonics method. In [4] it was reported that such a technique accelerates computations of the radiance field by two orders of magnitude. The algorithm based on the allocation of the anisotropic solution on the basis of expression (42) received the name of the modified DOM (MDOM) [4]. We estimate the accuracy of the MDOM based on a comparison with the exact analytical solution of S. Chandrasekhar [25] for the radiance of reflected radiation from a semi-infinite medium with a Rayleigh scattering phase function. We take the functions necessary for calculating the exact solution from the tables [26]. The relative error of the MDOM relative to the exact analytical solution Δ for a medium with Λ = 0.999 and two incidence angles ϑ0 = 0° and 70° for different sight angles is shown in Figure 1. The observation scheme is adopted, in which the direction to the Nadir ϑ = 180°. As one can see, the error does not exceed 10 −5 in almost the entire range of sight angles. Here it is necessary to note the considerable role that accurate analytical methods play in the transfer theory, but this is a vast separate topic that goes beyond the scope of our article. Note that the development of the analytical approach did not stop at the classics of the transfer theory but is still actively continuing today, see, for example [27].

Reflection and Refraction at the Ocean-Atmosphere Interface
Let us consider a layer illuminated by a plane unidirectional source with a diffusely reflecting substrate at the bottom with the reflection coefficient ρ. In this case, the boundary condition at the bottom in Equation (17) reads as follows: For Fourier expansion coefficients in the discrete ordinate space, Equation (43) can be written in the matrix form as: However, this approach cannot be directly applied at the boundary with refraction, since the directions of the ordinates ("rays") change at the boundary according to Snell's law. The solution to this problem was proposed in [28]. We consider the coupled ocean-atmosphere system. For simplicity, we set the refractive index of the atmosphere to 1, while the reflective index of the ocean is no > 1. The discrete ordinates in the atmosphere are related to those in the ocean through Snell's law as follows: These oceanic ordinates belong to the refraction zone. For μ is the total reflection angle cosine) we have the total reflection zone, in which the rays cannot escape from the water body, but rather are reflected from the water-air interface. The boundary condition at the ocean-atmosphere interface reads as follows: The first and last integrals are related to the refraction zone, while the second integral corresponds to the total reflection zone. After performing the change of variables in the second integral: where t ′ μ = μ ν , the Gaussian quadrature formula can be applied with Nt points per hemisphere.
Consequently, we obtain the two column-vectors C , C where square brackets denote the concatenation of vectors, the matrix operator method gives the following equation: where R, T   are the Fresnel reflection and transmission matrices, respectively.

Synthetic Iteration Method
In [29,30] the performance assessment of MDOM was conducted. In particular, the computational time of MDOM was measured for computing the angular distributions of the reflected radiances with different values of N and M. At the same time, it was found that for all N and M the computed radiances coincide in the entire range of angles except for the small glory region. In other words, insufficient numbers of discrete ordinates and azimuthal harmonics N and M affects only a narrow range of angles. The calculation in the glory region required a 100 times increase of N and M, which led to the increase of the computational time by more than 300 times. Thus, the computations of the radiance field in the wide angular range is two orders of magnitude faster than computations of sharp peaks. Why is it so?
Let us consider the Legendre expansion of the radiance field: Legendre polynomials of an order not higher than N can be expressed through polynomials of order N + 1 by the discrete samples as follows: where μi are the roots of P N+1 (μ). Substituting Equation (52) into Equation (51) gives: This is the Lagrange interpolation formula for Lm (τ,μ). These relations can be regarded as analogous to the Nyquist-Shannon-Kotelnikov sampling theorem. The following comments are in order: 1) Since the separation of the anisotropic part of the solution provides the most accurate discrete representation of the scattering integral, the MDOM determines the mean convergence; 2) The convergence of the solution in the uniform metric is determined by the features of the angular distribution of the scattering phase function, then the convergence here of any method for isolating anisotropy will be equivalent; 3) To achieve good convergence in a uniform metric, the sampling interval must correspond to the angular size of the finest detail of the radiance distribution necessary to solve a practical problem. Thus, the Nyquist-Shannon-Kotelnikov sampling theorem put a lower constraint on the number of discrete ordinates N. Moreover, N is the main parameter that governs the accuracy and computational time. To accelerate the calculation, N should be kept as low as possible. The synthetic iteration method allows for the overcoming of the limit defined by the Nyquist-Shannon-Kotelnikov sampling theorem and achieves performance enhancement without compromising accuracy.
The method of synthetic iterations (SI) was initially proposed in nuclear physics [31]. It consists of two steps. In the first step, an approximate solution is sought that converges well in the average energy metric. In the second step, the iteration is performed, which refines the solution and significantly increases the convergence in a uniform angular metric. Since MDOM convergences fast in the average metric, one can expect a significant increase in convergence rate after a single iteration [32].
A numerical comparison of the reflected radiance after the first iteration of MDOM is given in [29,30]. To calculate the glory region in the MDOM program, N = 801 and M = 256 are required, which corresponds to a sampling step of less than 0.5°. To achieve the same accuracy in the synthetic iteration, only N = 11, M = 4 are required, which reduces the computation time by almost 60 times. Accordingly, the synthetic iteration from MDOM allows one to calculate the angular distribution of radiance in one second with the relative error in the uniform metric less than 1%.

Polarized Radiation Transfer
The considered scalar transfer theory allows a simple generalization to the vector case of polarization. In the ray approximation, the polarized radiation is described [13] l l l l z l l l z l l l l l l l l is the rotator with respect to the dihedral angle between the ˆˆ′ × l l and 0× l z planes, is the scattering matrix, 0 L  is the Stokes vector for the top boundary condition.
Discretization of the VRTE by DOM is based on reducing the double integral to the single one by applying the addition theorem. In the case of polarization, the scattering matrix is surrounded by the rotator matrices R( ) ⋅  which disturb the azimuthal symmetry and do not allow using the addition theorem for the surface harmonics. Kuščer-Ribarič [32] applied the circular basis to determine polarization. This form is called the CP-representation (circular polarization) and can be obtained from the real-number energetic SP-representation (Stokes polarization) using the following matrix transformation: where i is the imaginary number. Thus, CP-representation is the complex basis, but it makes possible to diagonalize the rotation matrix: It changes the rotator form and allows using generalized surface functions P (cos ) k mn θ for the scattering matrix representation, for which the special form of the addition theorem [33] can be applied: where ϕ and ϕ′ are the azimuth of l and ˆ′ l about an axis OZ, respectively.
As a result, all the coefficients in the VRTE become complex numbers. After applying the addition theorem for the generalized spherical function, one should return to the Stokes polarization (SP). For the column vector defined as: equation similar to Equation (22), can be derived (see [13,34]) but with all the matrices 4 times larger in each dimension.

Three-Dimensional Radiative Transfer Models
The analytical simplicity of solving a plane-parallel problem is determined by its symmetry. In the case of a three-dimensional medium (3D), such symmetry disappears and the BVP of the RTE in general cases must be solved numerically. In the numerical solution of the three-dimensional RTE, it is necessary to discretize the radiance field in the spatial and angular domains. To achieve the required accuracy of the RTE numerical solution, the fast convergence of the numerical calculation of the scattering integral over the total number of discrete angular directions is crucial. DOM is one of the gold standards for the angular discretization procedure [33]. Regardless of the specific scheme of the scattering integral representation for arbitrary 3D geometry of the medium, it is important to isolate the anisotropic part of the solution in order to reduce the number of discrete ordinates. As well as for the slab geometry, it is most appropriate to treat the anisotropic part of the solution by using the MSH, which has an analytical form for any arbitrary geometry of the medium [13], thereby reducing the number of ordinates required for the scattering integral representation.
However, like in the slab geometry case, synthetic iterations can be applied here as well: we find the smooth part of the solution in the diffusion approximation, and then make use of iterations. This approach was referred to as "quasi-diffusion" when the direct radiation was isolated [35]. For spatially localized sources, MSH is accurate in the entire front hemisphere of the viewing directions; therefore, only one iteration from MSH is required in many practical problems.
A special case of arbitrary geometry is stochastic media, when the optical characteristics of the medium are distributed in space according to a certain probabilistic law, which leads to a stochastic uncertainty of the spatial-angular distribution of the radiance field. The transition to the equation for averaged quantities is complicated by the covariance terms, i.e., mean values of the products of the spatial-angular distribution of radiance and a random function of the optical parameters [36]. Subsequent averaging procedures lead to the hierarchical system of equations. To solve this system, certain closure conditions are assumed (e.g., neglecting higher-order covariance terms). However, such simplified stochastic models cannot reproduce several phenomena observed in nature, e.g., statistical lenses formed by the oceanic surfaces.

Inverse Models
To retrieve a given parameter of interest by using non-linear least-square fitting, the corresponding forward model must be linearized. Now, we consider the computations of the derivatives of the radiance with respect to parameter of interest, ς . The analytical linearization method is one of the most efficient and generic techniques, as it can be applied to any implementation of the radiative transfer model. This approach is since the direct model is implemented as a sequence of derivatives of functions. The differential operator is applied, first, to the input parameters of the model. Then, by applying a chain rule, it propagates through the code until the derivatives of the radiances are computed. The linearization of the basic matrix operations is straightforward. The matrix multiplication is linearized by the Leibniz rule. The linearization of the eigenvalue problem (see Equation (26)) is not straightforward. As suggested in [37], we consider the eigenvalue problem with the normalization condition: (60) Applying the differential operator to Equation (60) and using the orthogonality condition for eigenvectors, one can reach the following equation: The derivatives for the initial eigenvalue problem for matrix B are expressed in terms of the solution of the system (61), which has to be solved for each layer and each azimuth harmonic. Later, by using the Leibniz rule, the derivatives for reflection and transmission matrices in the context of the matrix operator method (see Equations (29)-(31)) can be found.

Numerical Aspects
In practice, computing the anisotropic part of the radiance by using the small-angle approximation allows the number of discrete ordinates in the atmosphere to remain as low as two. For the ocean, at least four discrete ordinates are required. Note that the computations do not involve additional assumptions on the phase functions (i.e., no truncation procedures are required). The linearized model is just as accurate for the derivatives as the forward model for the radiances. Since all differentiations are performed "by hand", no additional assumptions are needed.
The coupled model MDOM has been compared with the coupled code LIDORT [37]. The results of the comparison are shown in Figure 2. The agreement up to the 6 th digit has been obtained. The computations have been performed for the coupled model with the flat oceanic surface. The ocean has been modeled as a two-layer system with the dissolved organic matter in the upper oceanic layer. The atmosphere was modeled as a 14-layer system containing the trace gases (O3, NO2) as well as the urban aerosol.

Conclusions
In this paper, our intension was to demonstrate modern advances in the discrete theory of radiative transfer. In a one-dimensional case, the theory provides an exact analytical solution as well as a numerical method for solving the radiative transfer equation. In this sense, a further comparison of different solvers reveals advantages of certain implementations (i.e., programming skills) or approximations rather than "new" methods behind the solvers [38]. Moreover, the comparison of various approximations looks misleading, because essentially, they are based on: 1) The synthetic iteration techniques equipped with the small-angle approximation, DOM or SHM or; 2) MOM with various degrees of assumptions. The accuracy and limits of applicability of all these approximations can be estimated based on the discrete theory of radiative transfer. From the user perspective, an open-source reference tool for the preparation of atmospheric scenarios and benchmark cases would be very useful.
In our opinion, unsolved problems of radiative transfer are exclusively outside one-dimensional models. Three-dimensional [39] and stochastic radiative transfer models [8,40] for retrieval of atmospheric and water constituents are those "Thomson clouds" in the radiative transfer theory. However, they are at the interdisciplinary boundaries between radiative transfer and related fields, such as cloud physics and data processing [41,42].