Open AccessArticle

Linearizations of the Spherical Harmonic Discrete Ordinate Method (SHDOM)

Adrian Doicu

^* and

Dmitry S. Efremenko

Remote Sensing Technology Institute, German Aerospace Center (DLR), 82234 Oberpfaffenhofen, Germany

Author to whom correspondence should be addressed.

Atmosphere 2019, 10(6), 292; https://doi.org/10.3390/atmos10060292

Submission received: 23 April 2019 / Revised: 21 May 2019 / Accepted: 22 May 2019 / Published: 28 May 2019

(This article belongs to the Special Issue Radiative Transfer Models of Atmospheric and Cloud Properties)

Download

Browse Figures

Versions Notes

Abstract

Linearizations of the spherical harmonic discrete ordinate method (SHDOM) by means of a forward and a forward-adjoint approach are presented. Essentially, SHDOM is specialized for derivative calculations and radiative transfer problems involving the delta-M approximation, the TMS correction, and the adaptive grid splitting, while practical formulas for computing the derivatives in the spherical harmonics space are derived. The accuracies and efficiencies of the proposed methods are analyzed for several test problems.

Keywords:

3D radiative trasnfer; SHDOM; Jacobian; linearization; adjoint radiative transfer

1. Introduction

Accurate and efficient modeling of radiative transfer is important for the retrieval of earth atmospheric constituent information by means of remote sensing. Real clouds are an inhomogeneous three-dimensional scattering medium. For the new generation of satellite spectrometers with a relative high spatial resolution (such as Sentinel 5 Precursor, Sentinel 4 and Sentinel 5 [1]), it is important to account for the sub-pixel cloud inhomogeneities, or at least, to assess their effect on the radiances at the top of the atmosphere, and in particular, on the trace gas and cloud retrieval results. So far, the majority of operational retrieval algorithms are based on the one-dimensional (1D) radiative transfer model, which neglects all horizontal variability of the atmosphere. However, 1D approximation can lead to the significant errors in the calculated radiances and retrieved parameters [2,3,4,5,6]. In this regard, three-dimensional (3D) radiative trasnfer models are required.

The SHDOM method developed by Evans is a widely used deterministic method to solve the 3D radiative transfer equation [7]. The radiative transfer equation is solved iteratively by using the spherical harmonic and the discrete ordinate representations of the radiance field. Each iteration consists of four steps: (1) the transformation of the source function from spherical harmonics to discrete ordinates; (2) the integration of the source function along discrete ordinate directions to compute the radiance field; (3) the transformation of the discrete ordinate radiances to spherical harmonics; and (4) the calculation of the source function from the radiance in spherical harmonics space. An important feature of SHDOM is the adaptive grid technique. This technique improves the solution accuracy by increasing the spatial resolution in regions where the source function is changing more rapidly. The difference scheme employed in SHDOM is a characteristic scheme. In the so-called long characteristic method, the integration is performed along the whole characteristic of the differential operator, while in the short characteristic method, the integration is performed along a segment of the characteristic located inside a cell. SHDOM has been modfied for use on multiple processors [8], and to include polarization [9,10,11].

In atmospheric remote sensing and data assimilation, it is necessary to compute the partial derivatives of the outgoing radiances with respect to some atmospheric parameters of interest. For this purpose, a linearized forward approach, relying on an analytical computation of the derivatives, or a linearized forward-adjoint approach, relying on the application of the adjoint radiative transfer theory, can be used. For a plane-parallel geometry, scalar and vector linearized forward approaches were described in Refs. [12,13], while linearized forward-adjoint approaches were proposed in [14,15,16,17,18,19,20,21]. The forward-adjoint method is extremely efficient because only two radiative transfer calculations are required for derivatives calculations.

A theoretical framework for three-dimensional remote sensing of the atmosphere with adjoint methods was described by Martin et al. [22]. In a subsequent paper, Martin and Hasekamp [23] used the linearized-forward approach to retrieve cloud extinction fields from multi-angle measurements. Although, the analysis was restricted to two-dimensional problems (a two-dimensional domain and a set of directions confined to the unit circle and characterized by one angular variable), the efficiency of adjoint methods for multi-dimensional remote sensing of the atmosphere and surface was demonstrated.

In this paper, linearizations of SHDOM by using a forward and a forward-adjoint approach are presented. After formulating the problem in Section 2, we discuss in Section 3 the fundamentals of the linearized SHDOM with a forward-adjoint approach. The accuracy and efficiency of the linearized models is analyzed in Section 4, while the final section of our paper contains a few concluding remarks.

2. Problem Formulation

We consider the solar radiative transfer in a rectangular cuboid of lengths

L_{x}

L_{y}

and

L_{z}

as shown in Figure 1. The top and the bottom faces of the prism are denoted by

S_{t}

and

S_{b}

, respectively, while the lateral faces are denoted by

S_{1 x}

(

x = 0

S_{2 x}

(

x = L_{x}

S_{1 y}

(

y = 0

), and

S_{2 y}

(

y = L_{y}

). The boundary-value problem for the total radiance at point

r

in direction

Ω

consists of the inhomogeneous differential equation (with the thermal emission term neglected)

\begin{matrix} \frac{d I}{d s} (r, Ω) & = - σ_{ext} (r) I (r, Ω) + \frac{σ_{sct} (r)}{4 π} \int_{Ω} P (r, Ω, Ω^{'}) I (r, Ω^{'}) d Ω^{'}, \end{matrix}

(1)

the top-of-atmosphere boundary condition

I (r_{t}, Ω^{-}) = \frac{F_{0}}{| μ_{0} |} δ (Ω^{-} - Ω_{0}), r_{t} \in S_{t},

(2)

and the surface boundary condition

\begin{matrix} I (r_{b}, Ω^{+}) & = \frac{1}{π} \int_{Ω^{-}} ρ (r_{b}, Ω^{+}, Ω^{-}) |μ^{-}| I (r_{b}, Ω^{-}) d Ω^{-}, r_{b} \in S_{b} . \end{matrix}

(3)

At the horizontal boundaries we assume periodic boundary conditions, i.e.,

\begin{matrix} I (x = 0, y, z, Ω) & = I (x = L_{x}, y, z, Ω), \\ I (x, y = 0, z, Ω) & = I (x, y = L_{y}, z, Ω) . \end{matrix}

(4)

Here,

σ_{ext}

and

σ_{sct} = ω σ_{ext}

are the extinction and scattering coefficients, respectively,

ω

is the single-scattering albedo, P is the phase function,

ρ

is the bidirectional surface reflection function,

Ω_{0} = (μ_{0}, φ_{0})

with

μ_{0} < 0

, is the solar direction,

F_{0}

is the solar flux,

Ω^{+}

and

Ω^{-}

denote an upward and a downward direction, respectively,

Ω

is the unit sphere, and

Ω^{+}

and

Ω^{-}

stand for the upper and lower unit hemispheres, respectively.

Let us define the forward transport operator

L

\begin{matrix} (L I) (r, Ω) & = \frac{d I}{d s} (r, Ω) + σ_{ext} (r) I (r, Ω) \\ - \frac{σ_{sct} (r)}{4 π} \int_{Ω} P (r, Ω, Ω^{'}) I (r, Ω^{'}) d Ω^{'} \\ - \frac{1}{π} δ (z) H (μ) μ \int_{Ω} ρ (r, Ω, Ω^{'}) H (- μ^{'}) | μ^{'} | I (r, Ω^{'}) d Ω^{'}, \end{matrix}

(5)

and the forward source function

Q (r, Ω)

Q (r, Ω) = F_{0} δ (z - L_{z}) δ (Ω - Ω_{0}),

(6)

where

δ

is the Dirac delta function and H is the Heaviside step function. It can be shown that if the radiation field

I (r, Ω)

solves the inhomogeneous operator equation

(L I) (r, Ω) = Q (r, Ω),

(7)

with the boundary conditions (4),

I (r_{t}, Ω^{-}) = 0, r_{t} \in S_{t}, and I (r_{b}, Ω^{+}) = 0, r_{b} \in S_{b},

(8)

then

I (r, Ω)

solves the radiative transfer Equation (1) with the boundary conditions (2), (3), and (4). The converse result is also true.

The adjoint transport operator

L^{†}

is defined through the relation

〈L I, I^{†}〉 = 〈I, L^{†} I^{†}〉,

(9)

where the scalar product of the fields

I_{1}

and

I_{2}

is given by

〈I_{1}, I_{2}〉 = \int_{D} \int_{Ω} I_{1} (r, Ω) I_{2} (r, Ω) d Ω d V

(10)

and D is the domain of the rectangular cuboid. The expression of the adjoint operator, under the assumptions that (i) I satisfies the boundary conditions (4) and (8), and (ii)

I^{†}

satisfies the boundary conditions (4),

I^{†} (r_{t}, Ω^{+}) = 0, r_{t} \in S_{t}, and I^{†} (r_{b}, Ω^{-}) = 0, r_{b} \in S_{b},

(11)

is given by

\begin{matrix} (L^{†} I^{†}) (r, Ω) & = - \frac{d I^{†}}{d s} (r, Ω) + σ_{ext} (r) I^{†} (r, Ω) \\ - \frac{σ_{sct} (r)}{4 π} \int_{Ω} P (r, Ω^{'}, Ω) I^{†} (r, Ω^{'}) d Ω^{'} \\ - \frac{1}{π} δ (z) H (- μ) | μ | \int_{Ω} ρ (r, Ω^{'}, Ω) H (μ^{'}) μ^{'} I^{†} (r, Ω^{'}) d Ω^{'} . \end{matrix}

(12)

Note that the requirement of homogeneous and periodic boundary conditions for both I and

I^{†}

is essential when computing the adjoint of the differential operator

d / d s

by using the Gauss theorem over the domain D.

In our analysis we consider a situation which is typical for atmospheric remote sensing. Atmospheric radiation can be measured e.g., by a spectrometer, which in the following is called detector. The detector, which measures the radiance at the top of the atmosphere in the direction

Ω_{m}^{0} = (μ_{m}^{0}, φ_{m}^{0})

μ_{m}^{0} > 0

, is represented by a point with position vector

r_{D} = (x_{D}, y_{D}, z_{D})

and is characterized by a characteristic function of its field of view

χ_{FOV} = χ_{FOV} (Ω_{m} \cdot Ω_{m}^{0})

. The signal measured by the detector, can be written as

I = \int_{Ω} χ_{FOV} (Ω_{m} \cdot Ω_{m}^{0}) I (r_{t} (Ω_{m}), Ω_{m}) d Ω_{m},

(13)

where

I (r_{t} (Ω_{m}), Ω_{m})

is the radiance at the top of the atmosphere in direction

Ω_{m}

. For the characteristic function

χ_{FOV}

, we consider the representation

χ_{FOV} (Θ) = \{\begin{matrix} \frac{1}{Ω_{FOV}}, 0 \leq Θ \leq Θ_{FOV} \\ 0, otherwise \end{matrix},

(14)

where

Θ_{FOV}

and

Ω_{FOV} \approx π Θ_{FOV}^{2}

are the collecting polar and solid angles of the detector, respectively, and

cos Θ = Ω_{m} \cdot Ω_{m}^{0}

. Let

S_{tm} = S_{tm} (r_{D}, Ω_{m}^{0})

be the footprint of the detector on the top face

S_{t}

, i.e., an ellipse centered at

r_{t}^{0} = (x_{t}^{0}, y_{t}^{0}, L_{z})

with the semi-major and semi-minor axes

b = R_{D} Θ_{FOV} / μ_{m}^{0}

and

a = R_{D} Θ_{FOV}

, respectively, where

R_{D} = | r_{D} - r_{t}^{0} |

. Using the relation

d Ω_{m} = (μ_{m}^{0} / R_{D}^{2}) d S

, where

d S

is the surface element on the top face

S_{t}

, we express the measured signal as

I = \int_{S_{t}} h (r_{t}) I (r_{t}, Ω_{m} (r_{t})) d S_{t},

(15)

where

Ω_{m} (r_{t}) = \frac{r_{D} - r_{t}}{| r_{D} - r_{t} |} = (μ_{m} (r_{t}), φ_{m} (r_{t})),

(16)

is the unit vector along the line connecting the points

r_{t} \in S_{tm}

and

r_{D}

h (r_{t}) = \{\begin{matrix} \frac{1}{A_{tm}}, r_{t} \in S_{tm} \\ 0, otherwise \end{matrix}

(17)

is the characteristic function associated to the detector footprint, and

A_{tm} = π a b = R_{D}^{2} Ω_{FOV} / μ_{m}^{0}

is the area of the ellipse.

The radiative transfer in an inhomogeneous three-dimensional medium is usually solved by a numerical method. Considering a spatial discretization of the domain of analysis, the extinction coefficient

σ_{ext} (r)

, the single-scattering albedo

ω (r

), the Legendre expansion coefficients of the phase function

χ_{n} (r)

(Appendix A), are specified at every grid point of the so-called base grid and trilinearly interpolated within each cell. Furthermore, the bidirectional surface reflection function

ρ (r_{b}, Ω^{+}, Ω^{-})

is specified at every grid point on the bottom surface and in all upward and downward discrete ordinate directions.

3. Linearized SHDOM

In atmospheric remote sensing we are interested in the retrieval of say,

N_{p}

atmospheric parameters

ξ_{p}

p = 1, \dots, N_{p}

, which determine the fields

σ_{ext} (r)

ω (r)

χ_{n} (r)

, and

ρ (r_{b}, Ω^{+}, Ω^{-})

. The retrieval requires the knowledge of the weighting functions, i.e., the partial derivatives of the measured radiance with respect to the atmospheric parameters being retrieved, i.e.,

\partial I / \partial ξ_{p}

3.1. Linearized Forward Approach

The linearized forward approach relies on the analytical differentiation technique and is based on the fact that the radiance can be regarded as the composition of differentiable functions. By applying the chain rule repeatedly to these functions, derivatives of arbitrary order can be computed analytically. After choosing the independent variables with respect to which the derivatives of the radiance are required, the derivatives of each function are computed recursively, i.e., by repeatedly substituting the derivatives of the inner function in the chain rule.

The pertinent equations which are differentiated are given in Appendix A. Essentially, the analytical differentiation technique is applied to all steps of the algorithm including the transformation of the source function from spherical harmonics to discrete ordinates, the integration of the source function along discrete ordinate directions, the transformation of the discrete ordinate radiances to spherical harmonics, the calculation of the source function from the radiance in spherical harmonics space, the adaptive grid procedure, and the acceleration step. The benefit of this method is that no assumptions rather than those of the forward model have to be made. The disadvantage is that not only the radiance and the source function have to be stored at all grid points, but also their derivatives with respect to the atmospheric parameters of interest. As a result, the amount of required memory may exceed the size of RAM (several Gbytes).

3.2. Linearized Forward-Adjoint Approach

In view of Equation (15), the measured signal can be written as

I = 〈Q^{†}, I〉,

(18)

where

\begin{matrix} Q^{†} (r, Ω) & = F_{0}^{†} (r_{t}) δ (z - L_{z}) δ (Ω - Ω_{m} (r)), \end{matrix}

(19)

with

F_{0}^{†} (r_{t}) = h (r_{t}),

(20)

is the adjoint source function. Here,

Ω_{m} (r) = \frac{r_{D} - r}{| r_{D} - r |} = (μ_{m} (r), φ_{m} (r))

(21)

is the unit vector along the line joining the points

r

and

r_{D}

, and

r_{t} = r + s Ω_{m} (r) \in S_{tm}

with

s = (L_{z} - z) / μ_{m} (r)

and

r = (x, y, z)

, is the intersection point of this line with the top face

S_{t}

(Figure 1). Obviously, because the points

r

r_{t}

, and

r_{D}

are on the same line, we have

Ω_{m} (r_{t}) = Ω_{m} (r)

, and so,

μ_{m} (r_{t}) = μ_{m} (r)

The main result of the adjoint radiative transfer theory states that if (i) I solves the forward problem consisting in the operator equation

L I = Q

and the boundary conditions (4) and (8), and (ii)

I^{†}

solves the adjoint problem consisting in the operator equation

L^{†} I^{†} = Q^{†}

and the boundary conditions (4) and (11), then the measured signal

I

can be estimated either by the scalar product of the radiance field I (as determined by Q) and the adjoint source function

Q^{†}

, or by the scalar product of the adjoint radiance field

I^{†}

(as determined by

Q^{†}

) and the forward source function Q, i.e.,

I = 〈Q^{†}, I〉 = 〈L^{†} I^{†}, I〉 = 〈I^{†}, L I〉 = 〈I^{†}, Q〉 .

(22)

The partial derivative of the measured signal with respect to some atmospheric parameter

ξ_{p}

, can be computed as

\frac{\partial I}{\partial ξ_{p}} = 〈Q^{†}, \frac{\partial I}{\partial ξ_{p}}〉 = 〈L^{†} I^{†}, \frac{\partial I}{\partial ξ_{p}}〉 = 〈I^{†}, L \frac{\partial I}{\partial ξ_{p}}〉,

(23)

where in the first step we assumed that

\partial Q^{†} / \partial ξ_{p} = 0

, while in the last step we used the fact that

\partial I / \partial ξ_{p}

satisfies the boundary conditions (4) and (8). Taking the partial derivative of the operator equation

L I = Q

with respect to

ξ_{p}

and assuming

\partial Q / \partial ξ_{p} = 0

gives

\frac{\partial I}{\partial ξ_{p}} = - 〈I^{†}, \frac{\partial L}{\partial ξ_{p}} I〉 .

(24)

It should be pointed out that in general, the linearized forward-adjoint approach can be used to compute the weighted integral of the partial derivative of the radiance field

\partial I / \partial ξ_{p}

defined by

D_{p} = 〈Q^{†}, \frac{\partial I}{\partial ξ_{p}}〉 .

(25)

In this case, employing the same arguments as in the derivation of Equations (23) and (24), but without assuming that

\partial Q^{†} / \partial ξ_{p} = 0

, we get

D_{p} = - 〈I^{†}, \frac{\partial L}{\partial ξ_{p}} I〉 .

(26)

Thus, identifying the adjoint source function from Equation (25) and solving the corresponding adjoint radiative transfer equation

L^{†} I^{†} = Q^{†}

, we can compute the weighted integral

D_{p}

by means of Equation (26).

The forward and adjoint radiative transfer equations

L I = Q

and

L^{†} I^{†} = Q^{†}

, respectively, are related to each other. Replacing

Ω

- Ω

in the expression of the adjoint operator

L^{†}

, gives

L^{†} (- Ω) I^{†} (r, - Ω) = Q^{†} (r, - Ω)

. Defining the pseudo-forward field

{\hat{I}}^{†}

by the relation

{\hat{I}}^{†} (r, Ω) = I^{†} (r, - Ω)

and using the symmetry properties of the phase function

P (r, - Ω, - Ω^{'}) = P (r, Ω^{'}, Ω)

and of the normalized reflection function

ρ (r_{b}, - Ω, - Ω^{'}) = ρ (r_{b}, Ω^{'}, Ω)

, we find the identity

L^{†} (- Ω) I^{†} (r, - Ω) = L (Ω) {\hat{I}}^{†} (r, Ω)

. Thus, the pseudo-forward field

{\hat{I}}^{†}

solves the same type of radiative transfer equation as the radiance field I, i.e.,

L {\hat{I}}^{†} = {\hat{Q}}^{†}

, where the pseudo-forward source function, defined by

{\hat{Q}}^{†} (r, Ω) = Q^{†} (r, - Ω)

, is

\begin{matrix} {\hat{Q}}^{†} (r, Ω) & = {\hat{F}}_{0}^{†} (r_{t}) δ (z - L_{z}) δ (Ω - {\hat{Ω}}_{m} (r)), \end{matrix}

(27)

with

{\hat{Ω}}_{m} (r) = - Ω_{m} (r)

{\hat{μ}}_{m} (r) = - μ_{m} (r)

, and

{\hat{F}}_{0}^{†} (r_{t}) = F_{0}^{†} (r_{t}) = h (r_{t}) .

(28)

In other words,

{\hat{I}}^{†} (r, Ω)

solves the inhomogeneous differential equation

\begin{matrix} \frac{d {\hat{I}}^{†}}{d s} (r, Ω) & = - σ_{ext} (r) {\hat{I}}^{†} (r, Ω) + \frac{σ_{sct} (r)}{4 π} \int_{Ω} P (r, Ω, Ω^{'}) {\hat{I}}^{†} (r, Ω^{'}) d Ω^{'} \end{matrix}

(29)

with the top-of-atmosphere boundary condition

{\hat{I}}^{†} (r_{t}, Ω^{-}) = \frac{{\hat{F}}_{0}^{†} (r_{t})}{| {\hat{μ}}_{m} (r_{t}) |} δ (Ω^{-} - {\hat{Ω}}_{m} (r_{t})), r_{t} \in S_{t},

(30)

and the surface boundary condition

\begin{matrix} {\hat{I}}^{†} (r_{b}, Ω^{+}) & = \frac{1}{π} \int_{Ω^{-}} ρ (r_{b}, Ω^{+}, Ω^{-}) |μ^{-}| {\hat{I}}^{†} (r_{b}, Ω^{-}) d Ω^{-}, r_{b} \in S_{b} . \end{matrix}

(31)

We are now concern with the computation of the partial derivative

\partial I / \partial ξ_{p}

by means of Equation (24). Taking the variation of the forward transport operator

L

under the assumption that the phase function does not depend on the atmospheric parameter

ξ_{p}

, yields (cf. Equation (5))

\begin{matrix} (δ L I) (r, Ω) & = δ σ_{ext} (r) I (r, Ω) \\ - δ σ_{ext} (r) \frac{ω (r)}{4 π} \int_{Ω} P (r, Ω, Ω^{'}) I (r, Ω^{'}) d Ω^{'} \\ - σ_{ext} (r) \frac{δ ω (r)}{4 π} \int_{Ω} P (r, Ω, Ω^{'}) I (r, Ω^{'}) d Ω^{'} \\ - \frac{1}{π} δ (z) H (μ) μ \int_{Ω} δ ρ (r, Ω, Ω^{'}) H (- μ^{'}) | μ^{'} | I (r, Ω^{'}) d Ω^{'} . \end{matrix}

(32)

In this section we describe the computation of the partial derivative of the measured signal with respect to a parameter that determine the extinction field; the partial derivatives with respect to parameters that specify the single-scattering albedo and the bidirectional reflection function are given in Appendix B.

To compute the partial derivatives of the forward operator we need to find appropriate interpolation schemes for the extinction coefficient, source function, and radiance field within a grid cell. A major problem arises from the fact that SHDOM uses the integral form of the radiative transfer equation for radiance calculation, and in particular, the assumption that the extinction/source function product varies linearly along the characteristic. Unfortunately, when dealing with the differential form of the radiative transfer equation, this assumption can not be met; instead, we may assume that the extinction/source function product varies linearly within a cell. Thus, different interpolation schemes are used for radiance and derivative calculations.

In the following, we assume that the partial derivatives

\partial σ_{ext} (r_{i}) / \partial ξ_{p}

at all grid points

r_{i}

are known. Let c be the global index of a cell containing a grid point with global index j, and let

\bar{j}

be the local index of this point within cell c (Figure 2); we define the map g as

(\bar{j}, c) \overset{g}{⟶} j

, or equivalently,

j = g (\bar{j}, c)

. For

r \in D_{c}

, where

D_{c}

is the domain of cell c, we assume that the extinction and extinction/source function product vary linearly within the cell, i.e.

σ_{ext} (r) = \sum_{\bar{i}} L_{\bar{i}} (R) σ_{ext} (r_{g (\bar{i}, c)})

(33)

and

σ_{ext} (r) J (r, Ω) = \sum_{\bar{i}} L_{\bar{i}} (R) σ_{ext} (r_{g (\bar{i}, c)}) J (r_{g (\bar{i}, c)}, Ω),

(34)

respectively. Here,

J (r, Ω) = \frac{ω (r)}{4 π} \int_{Ω} P (r, Ω, Ω^{'}) I (r, Ω^{'}) d Ω^{'}

(35)

is the source function,

L_{\bar{i}}

are the first-order interpolation basis functions for a rectangular cuboid element (Appendix C),

R = r - ρ_{c}

is the position vector of a point in a local coordinate system attached to cell c,

ρ_{c} = (1 / 8) \sum_{\bar{i}} r_{g (\bar{i}, c)}

is the position vector of the center point of cell c,

r_{g (\bar{i}, c)} = r_{i}

is the position vector of the grid point with global index

i = g (\bar{i}, c)

, and the sum

\sum_{\bar{i}}

is taken over all grid points of a cell. Setting

σ_{ext g (\bar{i}, c)} = σ_{ext} (r_{g (\bar{i}, c)})

, we obtain

(\frac{\partial L}{\partial ξ_{p}} I) (r, Ω) = \sum_{\bar{i}} L_{\bar{i}} (R) [I (r, Ω) - J (r_{g (\bar{i}, c)}, Ω)] \frac{\partial σ_{ext g (\bar{i}, c)}}{\partial ξ_{p}},

(36)

implying

\frac{\partial I}{\partial ξ_{p}} = - \sum_{c} \sum_{\bar{i}} \int_{Ω} \int_{D_{c}} L_{\bar{i}} (R) {\hat{I}}^{†} (r, - Ω) [I (r, Ω) - J (r_{g (\bar{i}, c)}, Ω)] \frac{\partial σ_{ext g (\bar{i}, c)}}{\partial ξ_{p}} d V d Ω .

(37)

In the next step, we split the total radiance

I (r, Ω)

into the diffuse and direct radiances, i.e.,

\begin{matrix} I (r, Ω) & = I_{d} (r, Ω) + δ (Ω - Ω_{0}) T_{0} (r), \end{matrix}

(38)

where

T_{0} (r) = \frac{F_{0}}{| μ_{0} |} e^{- τ_{ext} (r, r_{0} | Ω_{0})}

(39)

is the transmission along the solar direction

Ω_{0}

from the starting point

r_{0} = r - s Ω_{0} \in S_{t}

to the end point

r

(the forward direct beam), and

τ_{ext} (r, r_{0} | Ω_{0}) = \int_{0}^{s} σ_{ext} (r_{0} + s^{'} Ω_{0}) d s^{'}

(40)

is the optical depth. Similarly, we split the pseudo-forward radiance

{\hat{I}}^{†} (r, Ω)

into its diffuse and direct components, i.e.,

\begin{matrix} {\hat{I}}^{†} (r, Ω) & = {\hat{I}}_{d}^{†} (r, Ω) + δ (Ω - {\hat{Ω}}_{m} (r)) {\hat{T}}_{m}^{†} (r), \end{matrix}

(41)

where

{\hat{T}}_{m}^{†} (r) = \frac{{\hat{F}}_{0}^{†} (r_{t})}{| {\hat{μ}}_{m} (r_{t}) |} e^{- τ_{ext} (r, r_{t} | {\hat{Ω}}_{m} (r))},

(42)

is the transmission along the detector direction

{\hat{Ω}}_{m} (r)

from the starting point

r_{t} = r - s {\hat{Ω}}_{m} (r) \in S_{tm}

to the end point

r

(the pseudo-forward direct beam). Approximating

Ω_{m} (r) = - {\hat{Ω}}_{m} (r) \approx Ω_{m} (ρ_{c}) : = Ω_{m c}, r \in D_{c},

(43)

we are led to the following representation for

\partial I / \partial ξ_{p}

in terms of diffuse radiances:

\frac{\partial I}{\partial ξ_{p}} = T_{A} + T_{B} + T_{C} + T_{D} + T_{E},

(44)

where

\begin{matrix} T_{A} & = - \sum_{c} \sum_{\bar{i}} \int_{Ω} \int_{D_{c}} L_{\bar{i}} (R) {\hat{I}}_{d}^{†} (r, - Ω) I_{d} (r, Ω) \frac{\partial σ_{ext g (\bar{i}, c)}}{\partial ξ_{p}} d V d Ω, \end{matrix}

(45)

\begin{matrix} T_{B} & = \sum_{c} \sum_{\bar{i}} \int_{Ω} \int_{D_{c}} L_{\bar{i}} (R) {\hat{I}}_{d}^{†} (r, - Ω) J (r_{g (\bar{i}, c)}, Ω) \frac{\partial σ_{ext g (\bar{i}, c)}}{\partial ξ_{p}} d V d Ω, \end{matrix}

(46)

and

\begin{matrix} T_{C} & = - \sum_{c} \sum_{\bar{i}} \int_{D_{c}} L_{\bar{i}} (R) {\hat{I}}_{d}^{†} (r, - Ω_{0}) T_{0} (r) \frac{\partial σ_{ext g (\bar{i}, c)}}{\partial ξ_{p}} d V, \end{matrix}

(47)

\begin{matrix} T_{D} & = - \sum_{c} \sum_{\bar{i}} \int_{D_{c}} L_{\bar{i}} (R) {\hat{T}}_{m}^{†} (r) I_{d} (r, Ω_{m c}) \frac{\partial σ_{ext g (\bar{i}, c)}}{\partial ξ_{p}} d V, \end{matrix}

(48)

\begin{matrix} T_{E} & = \sum_{c} \sum_{\bar{i}} \int_{D_{c}} L_{\bar{i}} (R) {\hat{T}}_{m}^{†} (r) J (r_{g (\bar{i}, c)}, Ω_{m c}) \frac{\partial σ_{ext g (\bar{i}, c)}}{\partial ξ_{p}} d V . \end{matrix}

(49)

The above integrals are computed in the spherical harmonics space by assuming the expansions

\begin{matrix} I_{d} (r, Ω) & = \sum_{m n} I_{m n} (r) Y_{m n} (Ω), \end{matrix}

(50)

\begin{matrix} {\hat{I}}_{d}^{†} (r, Ω) & = \sum_{m n} {\hat{I}}_{m n}^{†} (r) Y_{m n} (Ω), \end{matrix}

(51)

\begin{matrix} J (r, Ω) & = \sum_{m n} J_{m n} (r) Y_{m n} (Ω), \end{matrix}

(52)

where

Y_{m n} (Ω)

are the orthonormal real-valued spherical harmonic functions (Appendix A), and the sum

\sum_{m n}

should be understood as

\sum_{m n} = \sum_{m = - M}^{M} \sum_{n = | m |}^{N},

with N and M being the maximum expansion order and number of azimuthal modes, respectively. The linear approximation

f (r) = \sum_{\bar{j}} L_{\bar{j}} (R) f (r_{g (\bar{j}, c)}) = \sum_{\bar{j}} L_{\bar{j}} (R) f_{g (\bar{j}, c)},

(53)

where

f (r)

stands for

I_{m n} (r)

{\hat{I}}_{m n}^{†} (r)

T_{0} (r)

, and

{\hat{T}}_{m}^{†} (r)

, together with the identity

Y_{m n} (- Ω) = {(- 1)}^{n} Y_{m n} (Ω)

, then yield

\begin{matrix} T_{A} & = - \sum_{c} \sum_{m n} {(- 1)}^{n} \sum_{\bar{i}, \bar{j}, \bar{k}} {\hat{I}}_{m n, g (\bar{j}, c)}^{†} I_{m n, g (\bar{k}, c)} \frac{\partial σ_{ext g (\bar{i}, c)}}{\partial ξ_{p}} L_{\bar{i}, \bar{j}, \bar{k}}, \end{matrix}

(54)

\begin{matrix} T_{B} & = \sum_{c} \sum_{m n} {(- 1)}^{n} \sum_{\bar{i}, \bar{j}} {\hat{I}}_{m n, g (\bar{j}, c)}^{†} J_{m n, g (\bar{i}, c)} \frac{\partial σ_{ext g (\bar{i}, c)}}{\partial ξ_{p}} L_{\bar{i}, \bar{j}}, \end{matrix}

(55)

and

\begin{matrix} T_{C} & = - \sum_{c} \sum_{\bar{i}, \bar{j}, \bar{k}} {\hat{I}}_{d}^{†} (r_{g (\bar{j}, c)}, - Ω_{0}) T_{0 g (\bar{k}, c)} \frac{\partial σ_{ext g (\bar{i}, c)}}{\partial ξ_{p}} L_{\bar{i}, \bar{j}, \bar{k}}, \end{matrix}

(56)

\begin{matrix} T_{D} & = - \sum_{c} \sum_{\bar{i}, \bar{j}, \bar{k}} I_{d} (r_{g (\bar{j}, c)}, Ω_{m c}) {\hat{T}}_{m g (\bar{k}, c)}^{†} \frac{\partial σ_{ext g (\bar{i}, c)}}{\partial ξ_{p}} L_{\bar{i}, \bar{j}, \bar{k}}, \end{matrix}

(57)

\begin{matrix} T_{E} & = \sum_{c} \sum_{m n} Y_{m n} (Ω_{m c}) \sum_{\bar{i}, \bar{j}} {\hat{T}}_{m g (\bar{j}, c)}^{†} J_{m n, g (\bar{i}, c)} \frac{\partial σ_{ext g (\bar{i}, c)}}{\partial ξ_{p}} L_{\bar{i}, \bar{j}} . \end{matrix}

(58)

Here, the interpolations coefficients

L_{\bar{i}, \bar{j}, \bar{k}}

and

L_{\bar{i}, \bar{j}}

are given, respectively, by

\begin{matrix} L_{\bar{i}, \bar{j}, \bar{k}} & = \int_{D_{c}} L_{\bar{i}} (R) L_{\bar{j}} (R) L_{\bar{k}} (R) d V, \end{matrix}

(59)

\begin{matrix} L_{\bar{i}, \bar{j}} & = \int_{D_{c}} L_{\bar{i}} (R) L_{\bar{j}} (R) d V, \end{matrix}

(60)

while for

i = g (\bar{i}, c)

and

ω_{i} = ω (r_{i})

, the expansion coefficients of the source function in Equations (55) and (58) are computed as

J_{m n, i} = ω_{i} \frac{χ_{n, i}}{2 n + 1} I_{m n, i} + ω_{i} \frac{χ_{n, i}}{2 n + 1} Y_{m n} (Ω_{0}) T_{0 i},

(61)

where

χ_{n, i} = χ_{n} (r_{i})

are the Legendre expansion coefficients of the phase function (Appendix A), i.e.,

P (r, Ω, Ω^{'}) = 4 π \sum_{m n} \frac{χ_{n} (r)}{2 n + 1} Y_{m n} (Ω) Y_{m n} (Ω^{'}) .

(62)

Further comments and algorithm details are given below:

1.: The input parameters of the numerical algorithm are the partial derivatives $\partial σ_{ext} (r_{i}) / \partial ξ_{p}$ and $\partial ω (r_{i}) / \partial ξ_{p}$ at all grid points $r_{i}$ , and eventually, the partial derivatives $\partial ρ (r_{b i}, Ω_{i j}^{+}, Ω_{p q}^{-}) / \partial ξ_{p}$ at all grid points on the bottom surface $r_{b i}$ and in all upward and downward discrete ordinate directions $Ω_{i j}^{+}$ and $Ω_{p q}^{-}$ , respectively.
2.: The interpolation coefficients $L_{\bar{i}, \bar{j}, \bar{k}}$ and $L_{\bar{i}, \bar{j}}$ are computed analytically in a local coordinate system attached to the grid cell.
3.: The radiances ${\hat{I}}_{d}^{†} (r_{g (\bar{j}, c)}, - Ω_{0})$ and $I_{d} (r_{g (\bar{j}, c)}, Ω_{m c})$ in the expressions of $T_{C}$ and $T_{D}$ , respectively, are computed by the source integration method.
4.: For solar problems with the delta-M method [24], the radiative transfer equation is expressed in terms of the scaled quantities

${\bar{σ}}_{ext} = (1 - f ω) σ_{ext}, \bar{ω} = \frac{1 - f}{1 - f ω} ω,$

(63)

where f is the truncation fraction. At each grid point, we switch from the partial derivatives $\partial σ_{ext} / \partial ξ_{p}$ and $\partial ω / \partial ξ_{p}$ to $\partial {\bar{σ}}_{ext} / \partial ξ_{p}$ and $\partial \bar{ω} / \partial ξ_{p}$ , respectively, whereby the latter two are computed from Equation (63) by using the chain rule.
5.: The partial derivatives of the signal correction in the TMS method of Nakajima and Tanaka [25] are computed analytically (Appendix A).
6.: During the adaptive grid procedure, a base grid cell, also called a “parent cell”, is split into “child cells” to achieve higher spatial resolution. In this regard, the sums in Equations (54)–(58) are taken over all child cells.
7.: According to the main result of the adjoint radiative transfer theory, the measured signal can be computed by using Equation (15) or (22); in the latter case, the computational formula is

$\begin{matrix} I & = F_{0} \int_{S_{t}} {\hat{I}}_{d}^{†} (x_{t}, y_{t}, L_{z}, - Ω_{0}) d S_{t} . \end{matrix}$

(64)
8.: The code uses two alternative interpolation schemes for the source function: a linear variation of the extinction/source function product within a cell, i.e., Equation (34), and a linear variation of the source function within a cell, i.e.,

$\begin{matrix} J (r, Ω) & = \sum_{\bar{i}} L_{\bar{i}} (R) J (r_{g (\bar{i}, c)}, Ω), \end{matrix}$

(65)

for $r \in D_{c}$ . In the second case, the terms $T_{B}$ and $T_{E}$ are computed, respectively, as

$\begin{matrix} T_{B} & = \sum_{c} \sum_{m n} {(- 1)}^{n} \sum_{\bar{i}, \bar{j}, \bar{k}} {\hat{I}}_{m n, g (\bar{j}, c)}^{†} J_{m n, g (\bar{k}, c)} \frac{\partial σ_{ext g (\bar{i}, c)}}{\partial ξ_{p}} L_{\bar{i}, \bar{j}, \bar{k}}, \end{matrix}$

(66)

$\begin{matrix} T_{E} & = \sum_{c} \sum_{m n} Y_{m n} (Ω_{m c}) \sum_{\bar{i}, \bar{j}, \bar{k}} {\hat{T}}_{m g (\bar{j}, c)}^{†} J_{m n, g (\bar{k}, c)} \frac{\partial σ_{ext g (\bar{i}, c)}}{\partial ξ_{p}} L_{\bar{i}, \bar{j}, \bar{k}} . \end{matrix}$

(67)
9.: In the case of satellite remote sensing we may assume that (i) the distance from the top of the atmosphere to the detector $R_{D}$ is sufficiently large, so that we can approximate $I (r, Ω_{m} (r)) \approx I (r, Ω_{m}^{0})$ , and (ii) the footprint of the detector is a rectangle of lengths $2 a$ and $2 b$ centered at $r_{t}^{0}$ . The second assumption is made because the domain is discretized in rectangular cuboid elements. In this context, the measured signal is computed as

$I = \frac{1}{4 a b} \int_{S_{t}} \bar{h} (x_{t}, y_{t}) I (x_{t}, y_{t}, L_{z}, Ω_{m}^{0}) d S_{t},$

(68)

where the normalized characteristic function $\bar{h} (x_{t}, y_{t})$ takes the value 1 inside the footprint of the detector and 0 otherwise.
10.: The solution of the adjoint problem is a challenging task due to the spatial discontinuity of the pseudo-forward direct beam. In SHDOM, this type of problem is handled by the adaptive grid procedure. Essentially, the adaptive grid supplies extra spatial resolution along the boundaries of the pseudo-forward direct beam, so that there are not large discontinuities in the grid. In order to reduce the number of adaptive grid cells, the step characteristic function of the detector can be replaced by a smooth function, e.g., trapezoidal, Gaussian, cosine [23]. In the first case, the normalized characteristic function reads as

$\bar{h} (x_{t}, y_{t}) = \{\begin{matrix} {\bar{h}}_{0} (x_{t}, Δ h_{x}) {\bar{h}}_{0} (y_{t}, Δ h_{y}), (x_{t}, y_{t}) \in S_{tm} \\ 0, otherwise \end{matrix},$

(69)

where, for example, ${\bar{h}}_{0} (x_{t}, Δ h_{x})$ is given by

${\bar{h}}_{0} (x_{t}, Δ h_{x}) = c_{x} \times \{\begin{matrix} \frac{x_{t} - x_{t min}^{0}}{Δ h_{x}}, x_{t min}^{0} \leq x_{t} < x_{t min}^{0} + Δ h_{x} \\ 1, x_{t min}^{0} + Δ h_{x} \leq x_{t} \leq x_{t max}^{0} - Δ h_{x} \\ \frac{x_{t max}^{0} - x_{t}}{Δ h_{x}}, x_{t max}^{0} - Δ h_{x} < x_{t} \leq x_{t max}^{0} \end{matrix}$

(70)

with $c_{x} = 2 a / (2 a - Δ h_{x})$ , $x_{t min}^{0} = x_{t}^{0} - a$ , and $x_{t max}^{0} = x_{t}^{0} + a$ .
11.: In the linearized forward-adjoint approach, the forward and adjoint problems are solved successively on the same grid; in this way, the interpolation between different grids is avoided. Actually, in the first step, the adjoint problem is solved by using the adaptive grid procedure with a prescribed splitting accuracy $ε_{split}$ (Appendix A), and in the second step, the forward problem is solved on the resulting grid without splitting. In the linearized forward approach, the forward problem is solved first by using the adaptive grid procedure, and then, the partial derivatives are computed on the same grid; in this way, a less amount of memory for storing the derivatives of the radiance and the source function with respect to the atmospheric parameters of interest is required.

4. Numerical Simulations

In this section we analyze the accuracy and efficiency of the linearized versions of SHDOM.

4.1. Example 1

We consider a base grid with

N_{x} = N_{y} = 20

and

N_{z} = 11

points, and the discretization steps

Δ x = Δ y = Δ z = 0.1

km. For periodic boundary conditions, when the vertical boundaries in x and y have additional planes of grid points that duplicate the grid points in the planes

x = 0

and

y = 0

, respectively, the domain of analysis has the sizes

L_{x} = N_{x} Δ x = 2

km,

L_{y} = N_{y} Δ y = 2

km, and

L_{z} = (N_{z} - 1) Δ z = 1

km. The extinction field is Gaussian, i.e.,

σ_{ext} (x, y, z) = σ_{ext}^{\max} f (x, y, z),

(71)

f (x, y, z) = f_{G} (x, x_{0}, s_{x}) f_{G} (y, y_{0}, s_{y}) f_{G} (z, z_{0}, s_{z}),

(72)

f_{G} (u, u_{0}, s_{u}) = exp [- \frac{{(u - u_{0})}^{2}}{2 s_{u}^{2}}], u = x, y, z,

(73)

with

x_{0} = L_{x} / 2

y_{0} = L_{y} / 2

z_{0} = L_{z} / 2

s_{x} = L_{x} / 4

s_{y} = L_{y} / 4

, and

s_{z} = L_{z} / 4

. The parameters

σ_{ext}^{\max}

is 6

{km}^{- 1}

, in which case, the peak optical depth is

3.6

. The extinction field in the plane

y = L_{y} / 2

is illustrated in Figure 3. The number of discrete zenith and azimuthal angles are

N_{μ} = 16

and

N_{φ} = 2 N_{μ}

, respectively, the single-scattering albedo is

ω = 0.9

, and a Henyey–Greenstein phase function [26] with the asymmetry parameter

g = 0.8

is considered. In all our simulations, the solar zenith angle is

θ_{0} = 120^{\circ}

, and a Lambertian reflecting surface with the surface albedo

A = 0.2

is chosen. The coordinates of the center of the detector footprint are

x_{t}^{0} = L_{x} / 2

y_{t}^{0} = L_{y} / 2

, and

z_{t}^{0} = L_{z}

, the detector azimuthal angle is

φ_{m}^{0} = 0^{\circ}

, the lengths of the rectangular footprint are

a = 5 Δ x

, and

b = [a / (Δ y μ_{m}^{0})] Δ y

, where

[x]

is the integer part of x, and the lengths specifying the slope of the characteristic function in Equation (70) are

Δ h_{x} = 0.5 Δ x

and

Δ h_{y} = 0.5 Δ y

In the following, we compute the partial derivatives of the measured signal with respect to the extinction coefficients

σ_{ext p}

N_{p} = 11

base grid points (

x = x_{t}^{0}

y = y_{t}^{0}

z = \bar{p} Δ z

\bar{p} = 1, \dots, N_{p}

. Obviously, if p is the global index of the grid point with the local index

\bar{p}

in the derivative list, we have

\frac{\partial σ_{ext i}}{\partial σ_{ext p}} = δ_{i p}, \frac{\partial ω_{i}}{\partial σ_{ext p}} = 0 .

(74)

First, we compare the accuracies of the measured signal and of its derivatives. The accuracy of SHDOM radiances is determined by the angular and spatial resolution, i.e., the number of discrete ordinates in zenith angle

N_{μ}

and, for a specified base grid, the splitting accuracy

ε_{split}

. In Figure 4, we plot the partial derivatives computed with the linearized forward-adjoint approach and the linearized approach for different values of the splitting accuracy

ε_{split}

. Note that the influence of

N_{μ}

on the results accuracy is rather low. Taking the values corresponding to

ε_{split}^{0} = 5 \times 10^{- 4}

as a reference, we illustrate in Table 1, the relative error in the partial derivatives, expressed as an rms difference normalized by the mean of the quantity over the

N_{p}

comparison points, i.e.,

ε_{\partial I} = \frac{\sqrt{\sum_{\bar{p}} {[{(\frac{\partial I}{\partial σ_{ext p}})}_{ε_{split}} - {(\frac{\partial I}{\partial σ_{ext p}})}_{ε_{split}^{0}}]}^{2}}}{\sqrt{\sum_{\bar{p}} {(\frac{\partial I}{\partial σ_{ext p}})}_{ε_{split}^{0}}^{2}}},

(75)

and the relative error in the measured signal, i.e.,

ε_{I} = | I_{ε_{split}} - I_{ε_{split}^{0}} | / I_{ε_{split}^{0}} .

(76)

The results show that for the same

ε_{split}

1.: the accuracy of the measured signal is higher than that of its derivatives, and
2.: the accuracy of the linearized forward-adjoint approach is higher than that of the linearized forward approach.

However, in the next calculations, to avoid a large CPU time and computer memory, we take

ε_{split} = 10^{- 3}

. For this splitting accuracy, we show in Figure 5, the adaptive grids in the plane

y = L_{y} / 2

. From these plots we infer that

1.: for the linearized forward-adjoint approach, the number of adaptive grid cells is higher in the domain “seen” by the detector along the direction $θ_{m}^{0}$ and in particular, along the boundaries of the pseudo-forward direct beam,
2.: for the linearized forward approach, the number of adaptive grid cells is higher in the entire domain along the solar direction $θ_{0} = 120^{\circ}$ and when large discontinuities in the source function are present, and
3.: in general, the number of adaptive grid cells for the linearized forward approach is grater than that for the linearized forward-adjoint approach.

In Figure 6 we plot the partial derivatives

\partial I / \partial σ_{ext p}

for different values of the detector zenith angle

θ_{m}^{0}

. The results are computed with the linearized forward-adjoint approach (FAA), the linearized forward approach (FA), and the finite-difference approach (FDA). In the latter case, a centered finite difference scheme is used for derivative calculations. The partial derivatives

\partial I / \partial σ_{ext p}

N_{p} = 19

base grid points (

x = \bar{p} Δ x

y = y_{t}^{0}

z = L_{z} / 2

\bar{p} = 1 \dots ., N_{p}

, for different values of the solar azimuthal angle

φ_{0}

are plotted in Figure 7. Taking the values computed by the finite-difference approach as a reference, we illustrate in Table 2 the corresponding relative errors in the partial derivatives, defined, for example, in the case of the linearized forward-adjoint approach as

ε_{\partial I}^{FAA} = \frac{\sqrt{\sum_{\bar{p}} {[{(\frac{\partial I}{\partial σ_{ext p}})}_{FAA} - {(\frac{\partial I}{\partial σ_{ext p}})}_{FDA}]}^{2}}}{\sqrt{\sum_{\bar{p}} {(\frac{\partial I}{\partial σ_{ext p}})}_{FDA}^{2}}} .

(77)

The relative errors corresponding to the linearized forward-adjoint approach are smaller than

2.5 \times 10^{- 2}

, while the relative errors corresponding to the linearized forward approach are smaller than

4.3 \times 10^{- 3}

. The relative large errors of the linearized forward-adjoint approach can be explained as follows.

.1: As it can be inferred from Figure 5, the linearized forward-adjoint approach uses a rougher discretization grid (depending on the detector zenith angle $θ_{m}^{0}$ ) as compared to the forward and finite-difference approaches (which use the same grid).
2.: In the linearized forward-adjoint approach, the extinction/source function product is assumed to vary linearly within the cell, while in the forward and finite-difference approaches, the extinction/source function product is assumed to vary linearly along the characteristic.

4.2. Example 2

In the second example we consider in addition to the scattering and absorption by a cloud, molecular Rayleigh scattering, and the absorption by

{NO}_{2}

at the wavelength

λ = 443

nm. The difference to the previous geometry is that now we consider a base grid with

N_{z} = 16

points along the z axis, implying that the height of the domain of analysis is

L_{z} = (N_{z} - 1) Δ z = 1.5

km. The cloud extinction field is given by Equations (71)–(73), while the cloud single-scattering albedo and the phase function are computed by Mie theory [27] for a water-cloud model with a Gamma size distribution

P (a) \propto a^{α} exp [- α (\frac{a}{a_{\mod}})]

(78)

of parameters

a_{eff} = 10

μ m

a_{\mod} = 2 a_{eff} / 3

, and

α = 6

. Here, a is the particle radius, and the droplet size ranges between 0.02 and 50.0

μ

m. The computed value for the cloud single-scattering albedo is

ω_{cld} = 0.865

, the Rayleigh cross-section and depolarization ratios are calculated as in [28], while the pressure and temperature profiles correspond to the US standard model atmosphere [29]. The number density and absorption cross section profiles for

{NO}_{2}

n (z)

and

C_{abs} (z)

, respectively, are shown in Figure 8, and the absorption coefficient is calculated as

σ_{abs} (z) = n (z) C_{abs} (z)

. The peak optical depth is

5.4

In the following, we compute the partial derivative of the measured signal with respect to maximum value of the extinction field

σ_{ext}^{\max}

in Equation (71); in this case, we approximate

P \approx P_{cld}

and use

\frac{\partial σ_{ext i}}{\partial σ_{ext}^{\max}} = f (x_{i}, y_{i}, z_{i}) and \frac{\partial ω_{i}}{\partial σ_{ext}^{\max}} = \frac{ω_{cld} - ω_{i}}{σ_{ext i}} \frac{\partial σ_{ext i}}{\partial σ_{ext}^{\max}} .

(79)

We also compute the partial derivative of the measured signal with respect to the total column of

{NO}_{2}

X = (k_{B} T_{0} / p_{0}) \int_{0}^{L_{z}} n (z) d z

, where

k_{B}

is the Boltzmann constant, and

p_{0}

and

T_{0}

are the standard pressure and temperature. Assuming that the number density profile for

{NO}_{2}

is the a priori profile

n_{a} (z)

scaled by a constant, i.e.,

n (z) = α n_{a} (z)

, where

α

is the scaling factor, we get

\frac{\partial σ_{ext i}}{\partial X} = \frac{1}{X_{a}} σ_{abs} (z_{i}) and \frac{\partial ω_{i}}{\partial X} = - \frac{ω_{i}}{σ_{ext i}} \frac{\partial σ_{ext i}}{\partial X},

(80)

implying

\partial σ_{sct i} / \partial X = 0

In Table 3 and Table 4, we illustrate the partial derivatives

\partial I / \partial σ_{ext}^{\max}

and

\partial I / \partial X,

and the corresponding relative errors using the values computed by the finite-difference approach as a reference. Also shown are the relative errors of the measured signal computed by Equation (22) (or, equivalently, by Equation (64)) with respect to the measured signal computed by Equation (15). In these simulations, the number of discrete zenith angles is

N_{μ} = 20

, and the splitting accuracy is

ε_{split} = 3 \times 10^{- 3}

. The relative errors corresponding to the linearized forward-adjoint approach are now smaller than

2.3 \times 10^{- 2}

, while the relative errors corresponding to the linearized forward approach are smaller than

3.8 \times 10^{- 3}

. It should be pointed out that the small errors in the measured signal certify that the adjoint problem is correctly solved.

In Table 5 we compare the computational efficiency of all linearization methods. Note that the code is parallelized with the OpenMP interface. The simulations are performed on a server Intel(R) Xeon(R) CPU E5-2695 v3 @ 2.30 GHz with 56 cores, while the CPU times are given for the parallelized implementation. Clearly, the most efficient method is the linearized forward-adjoint approach, while the most time-consuming method is the finite-difference approach.

Finally, in Table 6, we show the numbers of adaptive cells and grid points for the test problems considered in this study. Because in the linearized forward-adjoint approach, the numbers of cells and points are smaller than in the linearized forward approach, we infer that the first method is not only faster but also less memory demanding than the second one.

5. Conclusions

Two linearization methods for SHDOM have been discussed. The first method is an analytical linearization approach, while the second method is based on the adjoint radiative transfer theory. In the latter case, practical formulas for computing the derivatives of the measured signal in the spherical harmonics space have been derived. Essentially, SHDOM has been specialized for derivative calculations and radiative transfer problems involving the delta-M approximation, the TMS correction, and the adaptive grid splitting.

Our numerical analysis leads to the following conclusions.

1.: The linearized SHDOM with analytical derivatives is an accurate approach. However, the method is time-consuming and demands a large computer memory. The reason is that not only the source function has to be stored as a spherical harmonic series at each grid point, but also its derivatives with respect to the atmospheric parameters of interest.
2.: The linearized SHDOM with a forward-adjoint approach requires less storage for derivatives calculation, are much faster, but relatively less accurate. The main reason for this lower accuracy is that different interpolation schemes are used for radiance and derivative calculations.

According to our numerical analysis, the relative errors of the linearized forward-adjoint approach are smaller than 2.5%. In this regard, it should be pointed out that in the framework of the Gauss-Newton method, errors of about 2–3% in the Jacobian do not play a significant role. The effect is somehow similar to the computation of the root of an equation by the secant method instead of the Newton method. However, if the errors in the Jacobian are important for the retrieval, the regularized total least squares, which accounts on the errors in both the Jacobian and the data, can be used.

In our analysis, we considered the case of a fixed detector position

r_{D}

and a single-angle measurement

Ω_{m}^{0}

. For multi-angle measurements and eventually, several detector positions, the solution of several adjoint problems (corresponding to each pair

(r_{D}, Ω_{m}^{0})

) can be avoided by employing the approach proposed in Ref. [23]. The main ideas of this approach, which can be applied for multi-dimensional remote sensing of cloud extinction fields from airborne radiometer measurements, are discussed in Appendix D. However, it should be pointed out that the computational cost for the retrieval of 3D extinction fields is extremely high. A reduction of the computation time can be achieved if a two-dimensional geometry (a constant extinction field along the y-axis, and a set of directions on the unit sphere) is considered. In this case, the application of the forward-adjoint approach requires that both the solar and detector azimuthal angles are zero (the solar and detector directions are in the

x z

-plane).

Potential application of the linearized SHDOM include the following.

1.: Retrieval of cloud model parameters. For broken clouds [30], the indicator function $f (x, y, z)$ in Equation (71) takes the values 1 inside the cloud and 0 inside the clear sky region. As in [31], the extinction field can be parametrized in terms of $σ_{ext}^{\max}$ (cf. Equation (71)), the cloud top height $h_{t}$ , and the cloud bottom height $h_{b}$ . The cloud optical thickness and the cloud geometrical parameters can be retrieved in the oxygen A-band by using the correlated k-distribution method for the broadband integration of the gaseous line absorption. Note that SHDOM is able to compute monochromatic and broadband radiative transfer (with a k distribution).
2.: Trace gas retrievals under cloudy conditions. The retrieval of total column of trace gases, e.g., $O_{3}$ , ${NO}_{2}$ , etc., can be performed by means of the differential optical absorption spectroscopy (DOAS) technique [32]. This approach requires the knowledge of the air mass factor (AMF), i.e., the partial derivative of the measured signal with respect to the total column of the trace gas at a specific wavelength. The main goal of the analysis is the computation of the air mass factor when clouds are outside the footprint of the instrument (effect of horizontal cloud edges on AMF calculation).

Author Contributions

Methodology, A.D.; Software, A.D. and D.S.E.; Writing—Original Draft Preparation, A.D. and D.S.E.; Writing—Review & Editing, A.D. and D.S.E.

Funding

This research received no external funding.

Acknowledgments

The authors are pleased to acknowledge very fruitful and stimulating discussions with Frank Evans during the completion of this work.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

In this appendix we summarize the SHDOM algorithm in order to emphasize the main equations which are differentiated in the framework of the linearized forward approach.

SHDOM uses both spherical harmonics and discrete ordinates to represent the radiance field during different parts of the solution algorithm. The spherical harmonics are employed for computing the source function including the scattering integral, while the discrete ordinates are used to integrate the radiative transfer equation spatially.

The set of discrete ordinates are a reduced Gaussian grid. There are

N_{μ}

Gaussian quadrature cosine zenith angles,

μ_{j}

, and

N_{φ}

evenly spaced azimuth angles,

φ_{k}

. The discrete ordinate set is reduced by having fewer azimuth angles at larger

{| μ |}_{j}

; thus,

N_{φ}

depends on j. The corresponding Gauss-Legendre quadrature weights and azimuthal integration weights (normalized appropriately) are

w_{μ j}

and

w_{φ j k}

, respectively.

The orthonormal real-value spherical harmonics are defined by

Y_{m n} (μ, φ) = P_{n}^{| m |} (μ) u_{m} (φ),

(A1)

where

P_{n}^{| m |} (μ)

are the normalized associated Legendre functions, and

u_{m} (φ) = \{\begin{matrix} (1 / \sqrt{π}) cos (| m | φ) \\ (1 / \sqrt{2 π}) \\ (1 / \sqrt{π}) sin (| m | φ) \end{matrix} \begin{matrix} m > 0 \\ m = 0 \\ m < 0 \end{matrix}

(A2)

are the Fourier harmonics. The source function

\begin{matrix} J (r, Ω) & = \frac{ω (r)}{4 π} \int_{Ω} P (r, Ω, Ω^{'}) [I_{d} (r, Ω^{'}) + δ (Ω^{'} - Ω_{0}) T_{0} (r)] d Ω^{'}, \end{matrix}

(A3)

\begin{matrix} T_{0} (r) & = \frac{F_{0}}{| μ_{0} |} e^{- τ_{ext} (r, r_{0} | Ω_{0})}, \end{matrix}

(A4)

where

τ_{ext} (r, r_{0} | Ω_{0}) = \int_{0}^{s} σ_{ext} (r_{0} + s^{'} Ω_{0}) d s^{'}

(A5)

with

s = | r - r_{0}

|, is the optical depth from the point

r

to the upper domain boundary point

r_{0}

along the solar direction

Ω_{0}

, is expressed in the spherical harmonic space as

J (r, Ω) = \sum_{m = - M}^{M} \sum_{n = | m |}^{N} J_{m n} (r) Y_{m n} (Ω),

(A6)

where

N = N_{μ} - 1

is the maximum expansion order and

M = N_{φ} / 2 - 1

is the maximum number of azimuthal modes.

By means of the addition theorem for the Legendre functions, the phase function, expanded in terms of unnormalized Legendre polynomials

\begin{matrix} {\tilde{P}}_{n} (cos Θ) & = \sqrt{\frac{2}{2 n + 1}} P_{n} (cos Θ), \end{matrix}

(A7)

P (r, cos Θ) = \sum_{n = 1}^{N_{rank}} χ_{n} (r) {\tilde{P}}_{n} (cos Θ),

(A8)

where

χ_{n} (r)

are the Legendre phase function coefficients and

N_{rank}

is the maximum expansion order of the phase function, is truncated in the spherical harmonic space as

\begin{matrix} P (r, cos Θ) & = 4 π \sum_{m = - M}^{M} \sum_{n = | m |}^{N} \frac{χ_{n} (r)}{2 n + 1} Y_{m n} (Ω) Y_{m n} (Ω^{'}), \end{matrix}

(A9)

where

cos Θ = Ω \cdot Ω^{'}

The solution method is based on Picard iterations. At the beginning of each iteration step, the expansion coefficients

J_{m n} (r_{i})

m = - M, \dots M

n = | m |, \dots N

, are assumed to be known at all grid points

r_{i}

. They are updated according to the following four computational steps.

Step 1.: The source function is transformed to discrete ordinates by means of the relation

$\begin{matrix} J (r_{i}, μ_{j}, φ_{k}) & = \sum_{m = - M}^{M} u_{m} (φ_{k}) \sum_{n = | m |}^{N} J_{m n} (r_{i}) P_{n}^{| m |} (μ_{j}) . \end{matrix}$

(A10)
Step 2.: The discrete ordinate radiance $I_{d} (r_{i}, μ_{j}, φ_{k})$ is computed from the source function $J (r_{i}, μ_{j}, φ_{k})$ by integrating the radiative transfer equation. Essentially, (i) the radiances are computed at all grid points in all downward discrete ordinate directions by integrating the radiative transfer equation with the top boundary condition $I_{d} (r_{t i}, μ_{j}^{-}, φ_{k}^{-}) = 0$ , where $(μ_{j}^{-}, φ_{μ j}^{-})$ and $(w_{μ j}^{-}, w_{φ j k}^{-})$ are the quadrature nodes and weights in the lower hemisphere, (ii) the radiances at the bottom boundary grid points in all upward discrete ordinate directions $(μ_{j}^{+}, φ_{k}^{+})$ are computed from the boundary condition

$\begin{matrix} I_{d} (r_{b i}, μ_{j}^{+}, φ_{k}^{+}) & = \frac{1}{π} ρ (r_{b i}, μ_{j}^{+}, φ_{k}^{+}, μ_{0}, φ_{0}) |μ_{0}| T_{0} (r_{b i}) \\ + \frac{1}{π} \sum_{p = 1}^{N_{μ} / 2} \sum_{q = 1}^{N_{φ p}} w_{μ p}^{-} w_{φ p q}^{-} ρ (r_{b i}, μ_{j}^{+}, φ_{k}^{+}, μ_{p}^{-}, φ_{q}^{-}) \\ \times | μ_{p}^{-} | I_{d} (r_{b i}, μ_{p}^{-}, φ_{q}^{-}), \end{matrix}$

(A11)

and (iii) the radiances are computed at all grid points in all upward discrete ordinate directions by integrating the radiative transfer equation with the boundary condition (A11).
Step 3.: The discrete ordinate radiance at each grid point is transformed to the spherical harmonic space according to

$\begin{matrix} I_{m n} (r_{i}) & = \sum_{j = 1}^{N_{μ}} w_{μ j} P_{n}^{| m |} (μ_{j}) \sum_{k = 1}^{N_{φ j}} w_{φ j k} I_{d} (r_{i}, μ_{j}, φ_{k}) u_{m} (φ_{k}) . \end{matrix}$

(A12)
Step 4.: At each grid point, the source function is computed from the radiance in the spherical harmonic space as

$\begin{matrix} J_{m n} (r_{i}) & = ω (r_{i}) \frac{χ_{n} (r_{i})}{2 n + 1} I_{m n} (r_{i}) + ω (r_{i}) \frac{χ_{n} (r_{i})}{2 n + 1} Y_{m n} (Ω_{0}) T_{0} (r_{i}) . \end{matrix}$

(A13)

Some comments are in order.

1.: The radiance and source function are initialized before the solution iterations with an Eddington radiative transfer solution on independent columns of the base grid.
2.: The explicit form of the transforms in Equations (A10) and (A12) illustrates how the azimuthal and zenith angle parts partially separate. For more than about 12 azimuthal angles, an FFT is used for the azimuthal Fourier transform.
3.: In Step 2, the radiative transfer equation is integrated backward from each grid point to a grid cell face that has known radiances at its bounding grid points. In the short-characteristic method, the integration is across just one cell, while in the long characteristic method, the characteristic $Ω_{j k} = (μ_{j}, φ_{k})$ is traced backward until the transmission falls below some minimum specified value. In the latter case, the error from interpolating the radiance at the grid cell face is reduced. According to the SHDOM difference scheme, the radiance at the exiting point A along a characteristic $Ω_{j k}$ , $I_{d} (r_{A}, Ω_{j k})$ , is computed from the integral form of the radiative transfer equation

$I_{d} (r_{A}, Ω_{j k}) = I_{d} (r_{B}, Ω_{j k}) e^{- \int_{0}^{s_{0}} σ_{ext} (r_{B} + s Ω_{j k}) d s} + S (s_{0}, J),$

(A14)

where $I_{d} (r_{B}, Ω_{j k})$ is the radiance at the entering point B, $s_{0} = | r_{A} - r_{B} |$ is the distance between the points A and B, and

$\begin{matrix} S (s_{0}, J) & = \int_{0}^{s_{0}} σ_{ext} (r_{B} + s Ω_{j k}) J (r_{B} + s Ω_{j k}, Ω_{j k}) \\ \times [e^{- \int_{s}^{s_{0}} σ_{ext} (r_{B} + s^{'} Ω_{j k}) d s^{'}}] d s \end{matrix}$

(A15)

is the integral of the source function along the characteristic. Setting $σ_{ext A, B} = σ_{ext} (r_{A, B})$ and $J_{A, B} = J (r_{A, B}, Ω_{j k})$ , and assuming that the extinction $σ_{ext}$ and the extinction/source function product $σ_{ext} J$ vary linearly along the characteristic, yields

$\begin{matrix} S (s_{0}, J) & = (1 - e^{- τ_{ext}}) \frac{1}{σ_{ext A} + σ_{ext B}} [σ_{ext A} J_{A} + σ_{ext B} J_{B} \\ + \frac{σ_{ext A} σ_{ext B}}{6} (J_{A} - J_{B}) s_{0}], \end{matrix}$

(A16)

for $τ_{ext} \leq 2$ , and

$\begin{matrix} S (s_{0}, J) & = (1 - e^{- τ_{ext}}) \frac{1}{σ_{ext A} + σ_{ext B}} [σ_{ext A} J_{A} + σ_{ext B} J_{B} \\ + \frac{2 σ_{ext A} σ_{ext B}}{σ_{ext A} + σ_{ext B}} (J_{A} - J_{B}) (1 - \frac{2}{τ_{ext}} + \frac{2 e^{- τ_{ext}}}{1 - e^{- τ_{ext}}})], \end{matrix}$

(A17)

for $τ_{ext} > 2$ , where

$τ_{ext} = τ_{ext} (r_{A}, r_{B} | Ω_{j k}) = \int_{0}^{s_{0}} σ_{ext} (r_{B} + s Ω_{j k}) d s$

(A18)

is the optical depth along the characteristic. Note that under the above assumptions, the representation (A16) follows from an expansion of the solution to the first order in the path distance $s_{0}$ , while the representation (A17) is a particular case of the general solution derived in [33]. The exiting and the entering values of the extinction and extinction/source function product $σ_{ext A}$ , $σ_{ext B}$ , $σ_{ext A} J_{A}$ , and $σ_{ext B} J_{B}$ are computed by bilinear interpolation from the grid point values of the faces pierced by the characteristic, while the radiance is also bilinearly interpolated between the surrounding grid points to give the initial radiance $I_{d} (r_{B}, Ω_{j k})$ for integration.
4.: To increase the convergence rate, an acceleration method based on geometrical convergence is applied. At the iteration step K, the “accelerated” source column vector is computed as

$J_{acc}^{(K)} = J^{(K)} + a Δ J^{(K)}$

(A19)

where $J^{(K)}$ is the column vector encapsulating the expansion coefficients $J_{m n}^{(K)} (r_{i})$ , $m = - M, \dots, M$ , $n = | m |, \dots, N$ , at all grid points $r_{i}$ , i.e., $J^{(K)} = [J_{m n}^{(K)} (r_{i})]$ ,

$\begin{matrix} Δ J^{(K)} & = J^{(K)} - J^{(K - 1)}, \end{matrix}$

(A20)

is the residual source vector at the iteration step K, and

$\begin{matrix} a & = \frac{1 - r cos ψ + r^{1 + π / (2 ψ)}}{1 + r^{2} - 2 r cos ψ} - 1 \end{matrix}$

(A21)

with

$\begin{matrix} r & = \frac{| | Δ J^{(N)} | |}{| | Δ J^{(N - 1)} | |}, \end{matrix}$

(A22)

$\begin{matrix} cos ψ & = \frac{{(Δ J^{(N - 1)})}^{T} (Δ J^{(N)})}{| | Δ J^{(N - 1)} | | | | Δ J^{(N)} | |}, \end{matrix}$

(A23)

is the acceleration parameter. The iterations are stopped when the solution criterion

$\frac{| | Δ J^{(N)} | |}{| | J^{(N - 1)} | |} \leq ε_{sol}$

(A24)

is satisfied, where $ε_{sol}$ is the solution tolerance.
5.: If the delta-M scaling method is applied, then the optical parameters $σ_{ext}$ , $ω$ , and $χ_{n}$ are scaled before their use. The scaled quantities ${\bar{σ}}_{ext}$ , $\bar{ω}$ , and ${\bar{χ}}_{n}$ are given, respectively, by (the dependency on $r$ is omitted)

$\begin{matrix} {\bar{σ}}_{ext} & = (1 - f ω) σ_{ext}, \end{matrix}$

(A25)

$\begin{matrix} \bar{ω} & = \frac{1 - f}{1 - f ω} ω, \end{matrix}$

(A26)

$\begin{matrix} \frac{{\bar{χ}}_{n}}{2 n + 1} & = \frac{1}{1 - f} (\frac{χ_{n}}{2 n + 1} - f), n = 0, \dots, N, \end{matrix}$

(A27)

where

$f = \frac{1}{2 N + 3} χ_{N + 1}$

(A28)

is the truncation fraction.
6.: The radiance at a specified direction $Ω_{m}$ and location $r$ is computed by integrating the source function through the medium, while the spherical harmonic representation of the source function is transformed to the desired viewing direction $Ω_{m}$ .
7.: For solar problems with the delta-M method, the TMS method is used to compute the source function. This method replaces the scaled, truncated Legendre phase function expansion for the singly scattered solar radiation by the full, unscaled phase function expansion, while the multiply scattered contribution still comes from the truncated phase function. The source function at grid point $r_{i}$ in direction $Ω_{m}$ is then computed as

$\begin{matrix} J (r_{i}, Ω_{m}) & = \sum_{m = - M}^{M} \sum_{n = | m |}^{N} J_{m n} (r_{i}) Y_{m n} (Ω_{m}) + Δ J (r_{i}, Ω_{m}), \end{matrix}$

(A29)

$\begin{matrix} Δ J (r_{i}, Ω_{m}) & = \frac{F_{0}}{| μ_{0} |} e^{- {\bar{τ}}_{ext} (r_{i}, r_{i 0} | Ω_{0})} [\frac{\bar{ω} (r_{i})}{1 - f} \frac{1}{4 π} \sum_{n = 1}^{N_{rank}} χ_{n} (r_{i}) {\tilde{P}}_{n} (cos Θ) \\ - \bar{ω} (r_{i}) \sum_{m = - M}^{M} \sum_{n = | m |}^{N} \frac{{\bar{χ}}_{n} (r_{i})}{2 n + 1} Y_{m n} (Ω_{0}) Y_{m n} (Ω_{m})], \end{matrix}$

(A30)

where $cos Θ = Ω_{0} \cdot Ω_{m}$ .
8.: The adaptive grid evolves from the base grid by splitting cells where more resolution is judged to be needed. The criterion for splitting cells is based on how much the source function times extinction changes across a cell. A cell may be split in half in either of the three Cartesian directions, depending on whether any of them exceed the splitting criterion. For cell c, the cell-splitting criterion is

$C_{c u j} = | Δ J_{u j} | (1 - e^{- τ_{u j}}),$

(A31)

where u is one of the three Cartesian directions x, y, and z,

$\begin{matrix} | Δ J_{u j} | & = \frac{1}{{\bar{σ}}_{ext u j}} \sqrt{\frac{1}{4 π} \int_{Ω} {[σ_{ext u j}^{+} J (r_{u j}^{+}, Ω) - σ_{ext u j}^{-} J (r_{u j}^{- .}, Ω)]}^{2} d Ω} \\ = \frac{1}{{\bar{σ}}_{ext u j}} \sqrt{\frac{1}{4 π} \sum_{m = - M}^{M} \sum_{n = | m |}^{N} {[σ_{ext u j}^{+} J_{m n} (r_{u j}^{+}) - σ_{ext u j}^{-} J_{m n} (r_{u j}^{-})]}^{2}}, \end{matrix}$

(A32)

$r_{u j}^{+} \in F_{u}^{+}$ and $r_{u j}^{-} \in F_{u}^{-}$ , $j = 1 \dots, 4$ , are two cell corners along direction u, $F_{u}^{+}$ and $F_{u}^{-}$ are the two cell faces crossed by direction u, $τ_{u j} = {\bar{σ}}_{ext u j} s_{u j}$ , ${\bar{σ}}_{ext u j} = (σ_{ext u j}^{+} + σ_{ext u j}^{-}) / 2$ , $σ_{ext u j}^{+} = σ_{ext} (r_{u j}^{+})$ , $σ_{ext u j}^{-} = σ_{ext} (r_{u j}^{-})$ , and $s_{u j} = | r_{u j}^{+} - r_{u j}^{-} |$ . The criterion is averaged over the four corners, i.e., ${\bar{C}}_{c u} = (1 / 4) \sum_{j = 1}^{4} C_{c u j}$ , and the cells are sorted by the maximum over the splitting directions $u = x, y, z$ of the averaged splitting criterion ${\bar{C}}_{c u}$ . At the iteration step K, those cells c with the highest criterion ${\bar{C}}_{c u}$ above a certain value $ε_{split}^{(K)}$ are split first. The new cells may themselves be divided during one solution iteration. As the solution iterations proceed, the cell splitting accuracy $ε_{split}^{(K)}$ is gradually lowered to the desired final cell-spliting accuracy $ε_{split}$ , and so more grid points are added at each iteration during this process.

Appendix B

The partial derivative of the measured signal with respect to a parameter

ξ_{p}

that determine the single-scattering albedo at every grid point

ω_{i} = ω (r_{i})

is computed in the spherical harmonic space as follows. For

r \in D_{c}

, we assume that the extinction/source function product varies linearly within cell c, i.e.,

\begin{matrix} σ_{ext} (r) J (r, Ω) = \sum_{\bar{i} = 1} L_{\bar{i}} (R) σ_{ext g (\bar{i}, c)} J (r_{g (\bar{i}, c)}, Ω) \\ = \sum_{\bar{i} = 1} L_{\bar{i}} (R) σ_{ext g (\bar{i}, c)} [\frac{ω_{g (\bar{i}, c)}}{4 π} \int_{Ω} P (r_{g (\bar{i}, c)}, Ω, Ω^{'}) I (r_{g (\bar{i}, c)}, Ω^{'}) d Ω^{'}] . \end{matrix}

(A33)

We obtain

\begin{matrix} (\frac{\partial L}{\partial ξ_{p}} I) (r, Ω) & = - \sum_{\bar{i}} L_{\bar{i}} (R) [\frac{σ_{ext g (\bar{i}, c)}}{4 π} \\ \times \int_{Ω} P (r_{g (\bar{i}, c)}, Ω, Ω^{'}) I (r_{g (\bar{i}, c)}, Ω^{'}) d Ω^{'}] \frac{\partial ω_{g (\bar{i}, c)}}{\partial ξ_{p}}, \end{matrix}

(A34)

yielding

\begin{matrix} \frac{\partial E}{\partial ξ_{p}} & = \sum_{c} \sum_{\bar{i}} \int_{Ω} \int_{D_{c}} L_{\bar{i}} (R) {\hat{I}}^{†} (r, - Ω) [\frac{σ_{ext g (\bar{i}, c)}}{4 π} \\ \times \int_{Ω} P (r_{g (\bar{i}, c)}, Ω, Ω^{'}) I (r_{g (\bar{i}, c)}, Ω^{'}) d Ω^{'}] \frac{\partial ω_{g (\bar{i}, c)}}{\partial ξ_{p}} d V d Ω . \end{matrix}

(A35)

Using

\int_{Ω} P (r, Ω, Ω^{'}) I_{d} (r, Ω^{'}) d Ω^{'} = 4 π \sum_{m n} \frac{χ_{n} (r)}{2 n + 1} I_{m n} (r) Y_{m n} (Ω)

(A36)

and

P (r, Ω, Ω^{'}) = 4 π \sum_{m n} \frac{χ_{n} (r)}{2 n + 1} Y_{m n} (Ω) Y_{m n} (Ω^{'}),

(A37)

we end up with

\frac{\partial I}{\partial ξ_{p}} = T_{A} + T_{B} + T_{C} + T_{D},

(A38)

where

\begin{matrix} T_{A} & = \sum_{c} \sum_{m n} {(- 1)}^{n} \sum_{\bar{i}, \bar{j}} \frac{χ_{n, g (\bar{i}, c)}}{2 n + 1} σ_{ext g (\bar{i}, c)} I_{m n, g (\bar{i}, c)} {\hat{I}}_{m n, g (\bar{j}, c)}^{†} \frac{\partial ω_{g (\bar{i}, c)}}{\partial ξ_{p}} L_{\bar{i}, \bar{j}}, \end{matrix}

(A39)

\begin{matrix} T_{B} & = \sum_{c} \sum_{m n} {(- 1)}^{n} Y_{m n} (Ω_{0}) \sum_{\bar{i}, \bar{j}} \frac{χ_{n, g (\bar{i}, c)}}{2 n + 1} σ_{ext g (\bar{i}, c)} T_{0 g (\bar{i}, c)} {\hat{I}}_{m n, g (\bar{j}, c)}^{†} \frac{\partial ω_{g (\bar{i}, c)}}{\partial ξ_{p}} L_{\bar{i}, \bar{j}}, \end{matrix}

(A40)

and

\begin{matrix} T_{C} & = \sum_{c} \sum_{m n} Y_{m n} (Ω_{m c}) \sum_{\bar{i}, \bar{j}} \frac{χ_{n, g (\bar{i}, c)}}{2 n + 1} σ_{ext g (\bar{i}, c)} I_{m n, g (\bar{i}, c)} {\hat{T}}_{m g (\bar{j}, c)}^{†} \frac{\partial ω_{g (\bar{i}, c)}}{\partial ξ_{p}} L_{\bar{i}, \bar{j}}, \end{matrix}

(A41)

\begin{matrix} T_{D} & = \sum_{c} \sum_{m n} Y_{m n} (Ω_{m c}) Y_{m n} (Ω_{0}) \sum_{\bar{i}, \bar{j}} \frac{χ_{n, g (\bar{i}, c)}}{2 n + 1} σ_{ext g (\bar{i}, c)} T_{0 g (\bar{i}, c)} {\hat{T}}_{m g (\bar{j}, c)}^{†} \frac{\partial ω_{g (\bar{i}, c)}}{\partial ξ_{p}} L_{\bar{i}, \bar{j}} . \end{matrix}

(A42)

Considering the partial derivative of the measured signal with respect to a parameter

ξ_{p}

that determine the bidirectional surface reflection function

ρ (r_{b}, Ω^{+}, Ω^{-})

, we use Equation (32), to obtain

\begin{matrix} \frac{\partial I}{\partial ξ_{p}} & = \frac{1}{π} | μ_{0} | μ_{m} (r_{b}) T_{0} (r_{b}) {\hat{T}}_{m}^{†} (r_{b}) \frac{\partial ρ}{\partial ξ_{p}} (r_{b}, Ω_{m} (r_{b}), Ω_{0}) \\ + \frac{1}{π} | μ_{0} | T_{0} (r_{b}) \int_{Ω} H (μ) μ {\hat{I}}_{d}^{†} (r_{b}, - Ω) \frac{\partial ρ}{\partial ξ_{p}} (r_{b}, Ω, Ω_{0}) d Ω \\ + \frac{1}{π} μ_{m} (r_{b}) {\hat{T}}_{m}^{†} (r_{b}) \int_{Ω} \frac{\partial ρ}{\partial ξ_{p}} (r_{b}, Ω_{m} (r_{b}), Ω) \\ \times H (- μ^{'}) | μ^{'} | I_{d} (r_{b}, Ω) d Ω \\ + \frac{1}{π} \int_{Ω} H (μ) μ {\hat{I}}_{d}^{†} (r_{b}, - Ω) \\ \times [\int_{Ω} \frac{\partial ρ}{\partial ξ_{p}} (r_{b}, Ω, Ω^{'}) H (- μ^{'}) | μ^{'} | I_{d} (r_{b}, Ω^{'}) d Ω^{'}] d Ω . \end{matrix}

(A43)

The integrals over the unit sphere are computed in the discrete ordinate space under the assumption that the partial derivatives

\partial ρ / \partial ξ_{p}

are known at all grid points on the bottom surface, and in all downward and upward discrete ordinate directions. For a Lambertian surface,

ρ (r_{b}, Ω^{+}, Ω^{-})

is simply the albedo A, and for

ξ_{p} = A

, the computational formula is

\begin{matrix} \frac{\partial I}{\partial A} & = \frac{1}{π} | μ_{0} | μ_{m} (r_{b i}) T_{0} (r_{b i}) {\hat{T}}_{m}^{†} (r_{b i}) \\ + \frac{1}{π} | μ_{0} | T_{0} (r_{b i}) \sum_{j = 1}^{N_{μ} / 2} \sum_{k = 1}^{N_{φ j}} w_{μ j}^{-} w_{φ j k}^{-} | μ_{j}^{-} | {\hat{I}}_{d}^{†} (r_{b i}, μ_{j}^{-}, φ_{k}^{-}) \\ + \frac{1}{π} μ_{m} (r_{b i}) {\hat{T}}_{m}^{†} (r_{b i}) \sum_{j = 1}^{N_{μ} / 2} \sum_{k = 1}^{N_{φ j}} w_{μ j}^{-} w_{φ j k}^{-} | μ_{j}^{-} | I_{d} (r_{b i}, μ_{j}^{-}, φ_{k}^{-}) \\ + \frac{1}{π} \sum_{j = 1}^{N_{μ} / 2} \sum_{k = 1}^{N_{φ j}} \sum_{p = 1}^{N_{μ} / 2} \sum_{q = 1}^{N_{φ p}} w_{μ j}^{-} w_{μ p}^{-} w_{φ j k}^{-} w_{φ p q}^{-} | μ_{j}^{-} | | μ_{p}^{-} | \\ \times {\hat{I}}_{d}^{†} (r_{b i}, μ_{j}^{-}, φ_{k}^{-}) I_{d} (r_{b i}, μ_{p}^{-}, φ_{q}^{-}) . \end{matrix}

(A44)

Appendix C

Referring to Figure 2, the first-order interpolation basis functions for a rectangular cuboid element

L_{\bar{i}}

\bar{i} = 1, \dots 8

at a point

r = (x, y, z)

are given by

\begin{matrix} L_{1} (ξ, η, ς) & = n_{-} (ξ) n_{-} (η) n_{-} (ς), \end{matrix}

(A45)

\begin{matrix} L_{2} (ξ, η, ς) & = n_{+} (ξ) n_{-} (η) n_{-} (ς), \end{matrix}

(A46)

\begin{matrix} L_{3} (ξ, η, ς) & = n_{-} (ξ) n_{+} (η) n_{-} (ς), \end{matrix}

(A47)

\begin{matrix} L_{4} (ξ, η, ς) & = n_{+} (ξ) n_{+} (η) n_{-} (ς), \end{matrix}

(A48)

\begin{matrix} L_{5} (ξ, η, ς) & = n_{-} (ξ) n_{-} (η) n_{+} (ς), \end{matrix}

(A49)

\begin{matrix} L_{6} (ξ, η, ς) & = n_{+} (ξ) n_{-} (η) n_{+} (ς), \end{matrix}

(A50)

\begin{matrix} L_{7} (ξ, η, ς) & = n_{-} (ξ) n_{+} (η) n_{+} (ς), \end{matrix}

(A51)

\begin{matrix} L_{8} (ξ, η, ς) & = n_{-} (ξ) n_{+} (η) n_{+} (ς), \end{matrix}

(A52)

where

n_{\pm} (ξ) = \frac{1}{2} (1 \pm ξ)

(A53)

are the one-dimensional linear interpolation functions, and

(ξ, η, ς)

are the normalized coordinates of the point in a local coordinate system attached to the center of the cell

ρ_{c} = (x_{c}, y_{c}, z_{c})

, i.e.,

ξ = \frac{x - x_{c}}{l_{x}}, η = \frac{y - y_{c}}{l_{y}}, ς = \frac{z - z_{c}}{l_{z}} .

(A54)

Appendix D

We consider the case of a detector with a fixed position

r_{D}

, that measures the radiance in K directions

Ω_{m k}^{0}

k = 1, \dots K

. For the measurement direction

Ω_{m k}^{0}

, let

S_{tm} (Ω_{m k}^{0})

be the detector footprint centered at

r_{t k}^{0}

. We aim to compute the state vector

ξ

by minimizing the sum of squared residuals

R (ξ) = \frac{1}{2} \sum_{k} {[I (ξ, Ω_{m k}^{0}) - I_{mes} (Ω_{m k}^{0})]}^{2},

(A55)

where

I (ξ, Ω_{m k}^{0})

and

I_{mes} (Ω_{m k}^{0})

are the simulated and measured signals of the detector, respectively. In the steepest-descent method characterized by a linear convergence rate, the objective function is approximated by a linear model and the search direction is taken as the negative of the gradient

g = \nabla R

. By taking account of Equation (15), we deduce that the pth component of the gradient vector is given by

\begin{matrix} g_{p} & = \frac{\partial R}{\partial ξ_{p}} \\ = \sum_{k} [I (ξ, Ω_{m k}^{0}) - I_{mes} (Ω_{m k}^{0})] \frac{\partial I}{\partial ξ_{p}} (ξ, Ω_{m k}^{0}) \\ = \int_{S_{t}} \sum_{k} [I (ξ, Ω_{m k}^{0}) - I_{mes} (Ω_{m k}^{0})] \\ \times h_{k} (r_{t}) \frac{\partial I}{\partial ξ_{p}} (r_{t}, Ω_{m k}^{0}) d S_{t}, \end{matrix}

(A56)

where

h_{k} (r_{t}) = \{\begin{matrix} \frac{1}{A_{tm k}}, r_{t} \in S_{tm} (Ω_{m k}^{0}) \\ 0, otherwise \end{matrix},

(A57)

A_{tm k} = R_{D k}^{2} Ω_{FOV} / μ_{m k}^{0}

and

R_{D k} = | r_{D} - r_{t k}^{0} |

. Further on, from Equation (25), we find that the adjoint source function reads as

\begin{matrix} Q^{†} (r, Ω) & = \sum_{k} F_{0 k}^{†} (r_{t}) δ (z - L_{z}) δ (Ω - Ω_{m k}^{0}), \end{matrix}

(A58)

where

F_{0 k}^{†} (r_{t}) = [I (ξ, Ω_{m k}^{0}) - I_{mes} (Ω_{m k}^{0})] h_{k} (r_{t}) .

(A59)

Thus, in each inversion step, we first solve the forward problem and compute the residual

I (ξ, Ω_{m k}^{0}) - I_{mes} (Ω_{m k}^{0})

, and then solve the adjoint problem with the adjoint source function (A58) and compute each component of the gradient of the residual function

g_{p}

by means of Equation (26). In this approach, the pseudo-forward direct beam is a superposition of direct beams corresponding to each

Ω_{m k}^{0}

. Consequently, if the angular domain covered by the measurements is large, it may happens that the spatial discontinuity of the pseudo-forward direct beam is not present.

References

Ingmann, P.; Veihelmann, B.; Langen, J.; Lamarre, D.; Stark, H.; Bazalgette Courrèges-Lacoste, G. Requirements for the GMES Atmosphere Service and ESA’s implementation concept: Sentinels-4/-5 and -5p. Remote Sens. Environ. Sentin. Missions Oppor. Sci. 2012, 120, 58–69. [Google Scholar] [CrossRef]
Kokhanovsky, A. The influence of horizontal inhomogeneity on radiative characteristics of clouds: An asymptotic case study. IEEE Trans. Geosci. Remote Sens. 2003, 41, 817–825. [Google Scholar] [CrossRef]
Marshak, A.; Platnick, S.; Várnai, T.; Wen, G.; Cahalan, R.F. Impact of three-dimensional radiative effects on satellite retrievals of cloud droplet sizes. J. Geophys. Res. 2006, 111. [Google Scholar] [CrossRef]
Efremenko, D.S.; Schüssler, O.; Doicu, A.; Loyola, D. A stochastic cloud model for cloud and ozone retrievals from UV measurements. J. Quant. Spectrosc. Radiat. Transf. 2016, 184, 167–179. [Google Scholar] [CrossRef]
Zhuravleva, T.B.; Nasrtdinov, I.M.; Russkova, T.V. Influence of 3D cloud effects on spatial-angular characteristics of the reflected solar radiation field. Atmos. Ocean. Opt. 2017, 30, 103–110. [Google Scholar] [CrossRef]
Efremenko, D.S.; Doicu, A.; Loyola, D.; Trautmann, T. Fast Stochastic Radiative Transfer Models for Trace Gas and Cloud Property Retrievals Under Cloudy Conditions. In Springer Series in Light Scattering; Springer International Publishing: Basel, Switzerland, 2018; pp. 231–277. [Google Scholar]
Evans, K. The spherical harmonic discrete ordinate method for three-dimensional atmospheric radiative transfer. J. Atmos. Sci. 1998, 55, 429–446. [Google Scholar] [CrossRef]
Pincus, R.; Evans, K. Computational cost and accuracy in calculating three-dimensional radiative transfer: Results for new implementations of Monte Carlo and SHDOM. J. Atmos. Sci. 2009, 66, 3131–3146. [Google Scholar] [CrossRef]
Spherical Harmonic Discrete Ordinate Method (SHDOM) for Atmospheric Radiative Transfer. Available online: http://coloradolinux.com/shdom/ (accessed on 20 April 2019).
Doicu, A.; Efremenko, D.; Trautmann, T. A multi-dimensional vector spherical harmonics discrete ordinate method for atmospheric radiative transfer. J. Quant. Spectrosc. Radiat. Transf. 2013, 118, 121–131. [Google Scholar] [CrossRef]
Emde, C.; Barlakas, V.; Cornet, C.; Evans, F.; Korkin, S.; Ota, Y.; Labonnote, L.C.; Lyapustin, A.; Macke, A.; Mayer, B.; et al. IPRT polarized radiative transfer model intercomparison project—Phase A. J. Quant. Spectrosc. Radiat. Transf. 2015, 164, 8–36. [Google Scholar] [CrossRef]
Spurr, R.; Kurosu, T.; Chance, K. A linearized discrete ordinate radiative transfer model for atmospheric remote-sensing retrieval. J. Quant. Spectrosc. Radiat. Transf. 2001, 68, 689–735. [Google Scholar] [CrossRef]
Spurr, R. LIDORT and VLIDORT. Linearized pseudo-spherical scalar and vector discrete ordinate radiative transfer models for use in remote sensing retrieval problems. In Light Scattering Reviews; Kokhanovsky, A., Ed.; Springer: Berlin/Heidelberg, Germany, 2008; Volume 3, pp. 229–275. [Google Scholar]
Marchuk, G. Equation for the value of information from weather satellite and formulation of inverse problems. Cosm. Res. 1964, 2, 394–409. [Google Scholar]
Marchuk, G.I. Adjoint Equations and Analysis of Complex Systems; Springer: Cham, Switzerland, 1995. [Google Scholar]
Landgraf, J.; Hasekamp, O.; Box, M.; Trautmann, T. A linearized radiative transfer model for ozone profile retrieval using the analytical forward-adjoint perturbation theory approach. J. Geophys. Res. Atmos. 2001, 106, 27291–27305. [Google Scholar] [CrossRef] [Green Version]
Ustinov, E.A. Adjoint sensitivity analysis of radiative transfer equation: Temperature and gas mixing ratio weighting functions for remote sensing of scattering atmospheres in thermal IR. J. Quant. Spectrosc. Radiat. Transf. 2001, 68, 195–211. [Google Scholar] [CrossRef]
Box, M. Radiative perturbation theory: A review. Env. Model. Softw. 2002, 17, 95–106. [Google Scholar] [CrossRef]
Walter, H.; Landgraf, J.; Hasekamp, O. Linearization of a pseudo-spherical vector radiative transfer model. J. Quant. Spectrosc. Radiat. Transf. 2004, 85, 251–283. [Google Scholar] [CrossRef]
Ustinov, E. Atmospheric weighting functions and surface partial derivatives for remote sensing of scattering planetary atmospheres in thermal spectral region: General adjoint approach. J. Quant. Spectrosc. Radiat. Transf. 2005, 92, 351–371. [Google Scholar] [CrossRef]
Rozanov, V.; Rozanov, A. Relationship between different approaches to derive weighting functions related to atmospheric remote sensing problems. J. Quant. Spectrosc. Radiat. Transf. 2007, 105, 217–242. [Google Scholar] [CrossRef]
Martin, W.; Cairns, B.; Bal, G. Adjoint methods for adjusting three-dimensional atmosphere and surface properties to fit multi-angle/multi-pixel polarimetric measurements. J. Quant. Spectrosc. Radiat. Transf. 2014, 144, 68–85. [Google Scholar] [CrossRef] [Green Version]
Martin, W.G.; Hasekamp, O.P. A demonstration of adjoint methods for multi-dimensional remote sensing of the atmosphere and surface. J. Quant. Spectrosc. Radiat. Transf. 2018, 204, 215–231. [Google Scholar] [CrossRef]
Wiscombe, W. The delta-M method: Rapid yet accurate radiative flux calculations for strongly asymmetric phase functions. J. Atmos. Sci. 1977, 34, 1408–1422. [Google Scholar] [CrossRef]
Nakajima, T.; Tanaka, M. Algorithms for radiative intensity calculations in moderately thick atmos using a truncation approximation. J. Quant. Spectrosc. Radiat. Transf. 1988, 40, 51–69. [Google Scholar] [CrossRef]
Henyey, L.C.; Greenstein, J.L. Diffuse radiation in the Galaxy. Astrophys. J. 1941, 93, 70. [Google Scholar] [CrossRef]
Mie, G. Beitraege zur Optik trueber Medien, speziell kolloidaler Metalloesungen. Ann. Phys. 1908, 330, 377–445. [Google Scholar] [CrossRef]
Bodhaine, B.A.; Wood, N.B.; Dutton, E.G.; Slusser, J.R. On rayleigh optical depth calculations. J. Atmos. Ocean. Technol. 1999, 16, 1854–1861. [Google Scholar] [CrossRef]
Anderson, G.; Clough, S.; Kneizys, F.; Chetwynd, J.; Shettle, E. AFGL Atmospheric Constituent Profiles (0–120 km); AFGL-TR-86-0110; Air Force Geophysics Laboratory: Hanscom Air Force Base, MA, USA, 1986.
Alexandrov, M.; Marshak, A.; Ackerman, A. Cellular Statistical Models of Broken Cloud Fields. Part I: Theory. J. Atmos. Sci. 2010, 67, 2125–2151. [Google Scholar] [CrossRef]
García, V.M.; Sasi, S.; Efremenko, D.S.; Doicu, A.; Loyola, D. Linearized radiative transfer models for retrieval of cloud parameters from EPIC/DSCOVR measurements. J. Quant. Spectrosc. Radiat. Transf. 2018, 213, 241–251. [Google Scholar] [CrossRef] [Green Version]
Platt, U.; Stutz, J. Differential Optical Absorption Spectroscopy: Principles and Applications; Springer: Berlin/Heidelberg, Germany, 2008. [Google Scholar]
Doicu, A.; Efremenko, D.; Trautmann, T. An analysis of the short-characteristic method for the spherical harmonic discrete ordinate method (SHDOM). J. Quant. Spectrosc. Radiat. Transf. 2013, 119, 114–127. [Google Scholar] [CrossRef]

Figure 1. Illustration of the radiative transfer problem.

Figure 2. Rectangular cuboid element with lengths

2 l_{x}

2 l_{y}

, and

2 l_{z}

. C is the center point of the cell, and

\bar{j} = 1, \dots, 8

is the local index of a grid point.

Figure 2. Rectangular cuboid element with lengths

2 l_{x}

2 l_{y}

, and

2 l_{z}

. C is the center point of the cell, and

\bar{j} = 1, \dots, 8

is the local index of a grid point.

Figure 3. Extinction field in the plane

y = L_{y} / 2

Figure 3. Extinction field in the plane

y = L_{y} / 2

Figure 4. Partial derivatives

\partial I / \partial σ_{ext p}

at the base grid points (

x = x_{t}^{0}

y = y_{t}^{0}

z = \bar{p} Δ z

\bar{p} = 1, \dots, 11

, computed with the linearized forward-adjoint approach (upper panel) and the linearized forward approach (lower panel) for different values of the splitting accuracy

ε_{split}

. The detector zenith angle is

θ_{m}^{0} = 0^{\circ}

, and the solar azimuthal angle is

φ_{0} = 0^{\circ}

Figure 4. Partial derivatives

\partial I / \partial σ_{ext p}

at the base grid points (

x = x_{t}^{0}

y = y_{t}^{0}

z = \bar{p} Δ z

\bar{p} = 1, \dots, 11

, computed with the linearized forward-adjoint approach (upper panel) and the linearized forward approach (lower panel) for different values of the splitting accuracy

ε_{split}

. The detector zenith angle is

θ_{m}^{0} = 0^{\circ}

, and the solar azimuthal angle is

φ_{0} = 0^{\circ}

Figure 5. Adaptive grids in the plane

y = L_{y} / 2

corresponding to the linearized forward-adjoint approach with

θ_{m}^{0} = 0^{\circ}

(upper panel) and

θ_{m}^{0} = 45^{\circ}

(middle panel), and the linearized forward approach (lower panel). The solar azimuthal angle is

φ_{0} = 0^{\circ}

. The horizontal axis is the x-axis, while the vertical axis is the z-axis.

Figure 5. Adaptive grids in the plane

y = L_{y} / 2

corresponding to the linearized forward-adjoint approach with

θ_{m}^{0} = 0^{\circ}

(upper panel) and

θ_{m}^{0} = 45^{\circ}

(middle panel), and the linearized forward approach (lower panel). The solar azimuthal angle is

φ_{0} = 0^{\circ}

. The horizontal axis is the x-axis, while the vertical axis is the z-axis.

Figure 6. Partial derivatives

\partial I / \partial σ_{ext p}

at the base grid points (

x = x_{t}^{0}

y = y_{t}^{0}

z = \bar{p} Δ z

\bar{p} = 1, \dots, 11

, for different values of the detector zenith angle

θ_{m}^{0}

. The results correspond to the solar azimuthal angle

φ_{0} = 0^{\circ}

and are computed with the linearized forward-adjoint approach (FAA), the linearized forward approach (FA), and the finite-difference approach (FDA).

Figure 6. Partial derivatives

\partial I / \partial σ_{ext p}

at the base grid points (

x = x_{t}^{0}

y = y_{t}^{0}

z = \bar{p} Δ z

\bar{p} = 1, \dots, 11

, for different values of the detector zenith angle

θ_{m}^{0}

. The results correspond to the solar azimuthal angle

φ_{0} = 0^{\circ}

and are computed with the linearized forward-adjoint approach (FAA), the linearized forward approach (FA), and the finite-difference approach (FDA).

Figure 7. Partial derivatives

\partial I / \partial σ_{ext p}

at the base grid points (

x = \bar{p} Δ x

y = y_{t}^{0}

z = L_{z} / 2

\bar{p} = 1, \dots, 19

, for different values of the solar azimuthal angle

φ_{0}

. The results correspond to the detector zenith angle

θ_{m}^{0} = 45^{\circ}

and are computed with the linearized forward-adjoint approach (FAA), the linearized forward approach (FA), and the finite-difference approach (FDA).

Figure 7. Partial derivatives

\partial I / \partial σ_{ext p}

at the base grid points (

x = \bar{p} Δ x

y = y_{t}^{0}

z = L_{z} / 2

\bar{p} = 1, \dots, 19

, for different values of the solar azimuthal angle

φ_{0}

. The results correspond to the detector zenith angle

θ_{m}^{0} = 45^{\circ}

and are computed with the linearized forward-adjoint approach (FAA), the linearized forward approach (FA), and the finite-difference approach (FDA).

Figure 8. Number density (left) and absorption cross section (right) profiles for

{NO}_{2}

Figure 8. Number density (left) and absorption cross section (right) profiles for

{NO}_{2}

Table 1. Relative errors

ε_{\partial I}

and

ε_{I}

for different values of the splitting accuracy

ε_{split}

. The results correspond to the linearized forward-adjoint approach (FAA) and the linearized forward approach (FA).

Table 1. Relative errors

ε_{\partial I}

and

ε_{I}

for different values of the splitting accuracy

ε_{split}

. The results correspond to the linearized forward-adjoint approach (FAA) and the linearized forward approach (FA).

$ε_{split}$	FAA		FA
$ε_{split}$	$ε_{\partial I}$	$ε_{I}$	$ε_{\partial I}$	$ε_{I}$
$8 \times 10^{- 4}$	$1.06 \times 10^{- 2}$	$1.16 \times 10^{- 4}$	$3.00 \times 10^{- 2}$	$8.70 \times 10^{- 5}$
$1 \times 10^{- 3}$	$1.42 \times 10^{- 2}$	$1.49 \times 10^{- 4}$	$3.83 \times 10^{- 2}$	$2.07 \times 10^{- 4}$
$3 \times 10^{- 3}$	$4.22 \times 10^{- 2}$	$8.36 \times 10^{- 4}$	$6.22 \times 10^{- 2}$	$7.05 \times 10^{- 4}$
$5 \times 10^{- 3}$	$5.58 \times 10^{- 2}$	$1.04 \times 10^{- 3}$	$7.75 \times 10^{- 2}$	$1.42 \times 10^{- 3}$

Table 2. The first part of the table shows the relative errors

ε_{\partial I}

in the partial derivatives

\partial I / \partial σ_{ext p}

at the base grid points (

x = x_{t}^{0}

y = y_{t}^{0}

z = \bar{p} Δ z

\bar{p} = 1, \dots, 11

, for different values of the detector zenith angle

θ_{m}^{0}

, while the second part shows the relative errors

ε_{\partial I}

in the partial derivatives

\partial I / \partial σ_{ext p}

at the base grid points (

x = \bar{p} Δ x

y = y_{t}^{0}

z = L_{z} / 2

\bar{p} = 1, \dots, 19

, for different values of the solar azimuthal angle

φ_{0}

. The results are computed with the linearized forward-adjoint approach (FAA) and the linearized forward approach (FA).

Table 2. The first part of the table shows the relative errors

ε_{\partial I}

in the partial derivatives

\partial I / \partial σ_{ext p}

at the base grid points (

x = x_{t}^{0}

y = y_{t}^{0}

z = \bar{p} Δ z

\bar{p} = 1, \dots, 11

, for different values of the detector zenith angle

θ_{m}^{0}

, while the second part shows the relative errors

ε_{\partial I}

in the partial derivatives

\partial I / \partial σ_{ext p}

at the base grid points (

x = \bar{p} Δ x

y = y_{t}^{0}

z = L_{z} / 2

\bar{p} = 1, \dots, 19

, for different values of the solar azimuthal angle

φ_{0}

. The results are computed with the linearized forward-adjoint approach (FAA) and the linearized forward approach (FA).

$θ_{m}^{0}$	$φ_{0}$	$ε_{\partial I}$
$θ_{m}^{0}$	$φ_{0}$	$FAA$	$FA$
$0^{°}$	$0^{°}$	$2.50 \times 10^{- 2}$	$4.02 \times 10^{- 3}$
$15^{°}$	$0^{°}$	$2.23 \times 10^{- 2}$	$4.28 \times 10^{- 3}$
$30^{°}$	$0^{°}$	$2.42 \times 10^{- 2}$	$3.22 \times 10^{- 3}$
$45^{°}$	$0^{°}$	$2.05 \times 10^{- 2}$	$1.29 \times 10^{- 3}$
$45^{°}$	$0^{°}$	$7.35 \times 10^{- 3}$	$9.74 \times 10^{- 4}$
$45^{°}$	$30^{°}$	$1.64 \times 10^{- 2}$	$6.52 \times 10^{- 4}$
$45^{°}$	$60^{°}$	$1.56 \times 10^{- 2}$	$3.97 \times 10^{- 4}$
$45^{°}$	$90^{°}$	$1.21 \times 10^{- 2}$	$3.10 \times 10^{- 3}$

Table 3. Partial derivative

\partial I / \partial σ_{ext}^{\max}

, and the relative errors

ε_{I}

and

ε_{\partial I}

for different values of the detector zenith angle

θ_{m}^{0}

. The results correspond to the solar azimuthal angle

φ_{0} = 0^{\circ}

and are computed with the linearized forward-adjoint approach (FAA) and the linearized forward approach (FA).

Table 3. Partial derivative

\partial I / \partial σ_{ext}^{\max}

, and the relative errors

ε_{I}

and

ε_{\partial I}

for different values of the detector zenith angle

θ_{m}^{0}

. The results correspond to the solar azimuthal angle

φ_{0} = 0^{\circ}

and are computed with the linearized forward-adjoint approach (FAA) and the linearized forward approach (FA).

$θ_{m}^{0}$	$\partial I / \partial σ_{ext}^{\max}$	$ε_{I}$	$ε_{\partial I}$
$θ_{m}^{0}$	$\partial I / \partial σ_{ext}^{\max}$	$ε_{I}$	FAA	FA
$0^{\circ}$	$5.776 \times 10^{- 4}$	$4.20 \times 10^{- 4}$	$3.46 \times 10^{- 3}$	$2.03 \times 10^{- 3}$
$15^{\circ}$	$1.439 \times 10^{- 3}$	$5.53 \times 10^{- 4}$	$4.86 \times 10^{- 3}$	$2.76 \times 10^{- 4}$
$30^{\circ}$	$2.930 \times 10^{- 3}$	$2.30 \times 10^{- 5}$	$3.40 \times 10^{- 3}$	$2.53 \times 10^{- 4}$
$45^{\circ}$	$5.499 \times 10^{- 3}$	$8.19 \times 10^{- 5}$	$1.81 \times 10^{- 3}$	$1.52 \times 10^{- 4}$

Table 4. Partial derivative

\partial I / \partial X

, and the relative errors

ε_{I}

and

ε_{\partial I}

for different values of the detector zenith angle

θ_{m}^{0}

. The results correspond to the solar azimuthal angle

φ_{0} = 0^{°}

and are computed with the linearized forward-adjoint approach (FAA) and the linearized forward approach (FA).

Table 4. Partial derivative

\partial I / \partial X

, and the relative errors

ε_{I}

and

ε_{\partial I}

for different values of the detector zenith angle

θ_{m}^{0}

. The results correspond to the solar azimuthal angle

φ_{0} = 0^{°}

and are computed with the linearized forward-adjoint approach (FAA) and the linearized forward approach (FA).

$θ_{m}^{0}$	$\partial I / \partial X$	$ε_{I}$	$ε_{\partial I}$
$θ_{m}^{0}$	$\partial I / \partial X$	$ε_{I}$	FAA	FA
$0^{°}$	$- 1.184 \times 10^{- 3}$	$4.20 \times 10^{- 4}$	$1.65 \times 10^{- 2}$	$3.42 \times 10^{- 3}$
$15^{°}$	$- 1.284 \times 10^{- 3}$	$5.53 \times 10^{- 4}$	$2.27 \times 10^{- 2}$	$3.70 \times 10^{- 3}$
$30^{°}$	$- 1.428 \times 10^{- 3}$	$2.30 \times 10^{- 5}$	$2.26 \times 10^{- 2}$	$2.82 \times 10^{- 3}$
$45^{°}$	$- 2.529 \times 10^{- 3}$	$8.19 \times 10^{- 5}$	$2.20 \times 10^{- 2}$	$3.06 \times 10^{- 3}$

Table 5. CPU time in min:sec for the test problems considered in this study and corresponding to the linearized forward-adjoint approach (FAA), the linearized forward approach (FA), and the finite-difference approach (FDA). The detector zenith angle is

θ_{m}^{0} = 45^{°}

and the solar azimuthal angle is

φ_{0} = 0^{°}

θ_{m}^{0} = 45^{°}

and the solar azimuthal angle is

φ_{0} = 0^{°}

Partial Derivatives	$FAA$	$FA$	$FDA$
$\partial I / \partial σ_{ext p}$ (Figure 6)	1:16	27:12	138:19
$\partial I / \partial σ_{ext p}$ (Figure 7)	1:47	32:34	153:26
$\partial I / \partial σ_{ext}^{\max}$	0:32	2:21	2:24
$\partial I / \partial X$	0:31	2:20	2:20

Table 6. The first part of the table shows the number of adaptive cells

N_{cells}

and the final number of grid points

N_{points}

for computing the derivatives

\partial I / \partial σ_{ext p}

with the splitting accuracy

ε_{split} = 10^{- 3}

, while the second part of the table shows

N_{cells}

and

N_{points}

for computing the derivatives

\partial I / \partial σ_{ext}^{\max}

and

\partial I / \partial X

with the splitting accuracy

ε_{split} = 3 \times 10^{- 3}

. The results correspond to the linearized forward-adjoint approach (FAA) and the linearized forward approach (FA).

Table 6. The first part of the table shows the number of adaptive cells

N_{cells}

and the final number of grid points

N_{points}

for computing the derivatives

\partial I / \partial σ_{ext p}

with the splitting accuracy

ε_{split} = 10^{- 3}

, while the second part of the table shows

N_{cells}

and

N_{points}

for computing the derivatives

\partial I / \partial σ_{ext}^{\max}

and

\partial I / \partial X

with the splitting accuracy

ε_{split} = 3 \times 10^{- 3}

. The results correspond to the linearized forward-adjoint approach (FAA) and the linearized forward approach (FA).

$ε_{split}$	$θ_{m}^{0}$	FAA		FA
$ε_{split}$	$θ_{m}^{0}$	$N_{cells}$	$N_{points}$	$N_{cells}$	$N_{points}$
$10^{- 3}$	$0^{°}$	$66,652$	$40,727$	$293,492$	$162,279$
$10^{- 3}$	$15^{°}$	$81,174$	$50,381$	$293,492$	$162,279$
$10^{- 3}$	$30^{°}$	$95,152$	$59,595$	$293,492$	$162,279$
$10^{- 3}$	$45^{°}$	$115,186$	$73,332$	$293,492$	$162,279$
$3 \times 10^{- 3}$	$0^{°}$	$17,804$	$13,196$	$54,164$	$33,576$
$3 \times 10^{- 3}$	$15^{°}$	$20,094$	$14,942$	$54,164$	$33,576$
$3 \times 10^{- 3}$	$30^{°}$	$22,674$	$16,989$	$54,164$	$33,576$
$3 \times 10^{- 3}$	$45^{°}$	$25,806$	$19,607$	$54,164$	$33,576$

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Doicu, A.; Efremenko, D.S. Linearizations of the Spherical Harmonic Discrete Ordinate Method (SHDOM). Atmosphere 2019, 10, 292. https://doi.org/10.3390/atmos10060292

AMA Style

Doicu A, Efremenko DS. Linearizations of the Spherical Harmonic Discrete Ordinate Method (SHDOM). Atmosphere. 2019; 10(6):292. https://doi.org/10.3390/atmos10060292

Chicago/Turabian Style

Doicu, Adrian, and Dmitry S. Efremenko. 2019. "Linearizations of the Spherical Harmonic Discrete Ordinate Method (SHDOM)" Atmosphere 10, no. 6: 292. https://doi.org/10.3390/atmos10060292

APA Style

Doicu, A., & Efremenko, D. S. (2019). Linearizations of the Spherical Harmonic Discrete Ordinate Method (SHDOM). Atmosphere, 10(6), 292. https://doi.org/10.3390/atmos10060292

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Linearizations of the Spherical Harmonic Discrete Ordinate Method (SHDOM)

Abstract

1. Introduction

2. Problem Formulation

3. Linearized SHDOM

3.1. Linearized Forward Approach

3.2. Linearized Forward-Adjoint Approach

4. Numerical Simulations

4.1. Example 1

4.2. Example 2

5. Conclusions

Author Contributions

Funding

Acknowledgments

Conflicts of Interest

Appendix A

Appendix B

Appendix C

Appendix D

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI