Open AccessArticle

A Proof-of-Concept Algorithm for the Retrieval of Total Column Amount of Trace Gases in a Multi-Dimensional Atmosphere

Adrian Doicu

^*,

Dmitry S. Efremenko

and

Thomas Trautmann

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

Author to whom correspondence should be addressed.

Remote Sens. 2021, 13(2), 270; https://doi.org/10.3390/rs13020270

Submission received: 10 December 2020 / Revised: 4 January 2021 / Accepted: 7 January 2021 / Published: 14 January 2021

(This article belongs to the Special Issue Current Advances in Radiative Transfer Modeling for Satellite Optical Remote Sensing Applications)

Download

Browse Figures

Figure 1
Indicator function <math display="inline"><semantics> <mrow> <mi>f</mi> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>,</mo> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>)</mo> </mrow> </semantics></math> with <math display="inline"><semantics> <mrow> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>=</mo> <mi>i</mi> <mo>Δ</mo> <mi>x</mi> </mrow> </semantics></math> and <math display="inline"><semantics> <mrow> <msub> <mi>y</mi> <mi>j</mi> </msub> <mo>=</mo> <mi>j</mi> <mo>Δ</mo> <mi>y</mi> </mrow> </semantics></math> for the 8 cloudy scenes (cases (a–h)); i and j refer to the pixel index in the xy plane, while each pixel corresponds to a 500 m × 500 m area. "> Figure 2
Slices of the indicator function <math display="inline"><semantics> <mrow> <mi>f</mi> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> </semantics></math> with <math display="inline"><semantics> <mrow> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>=</mo> <mi>i</mi> <mo>Δ</mo> <mi>x</mi> </mrow> </semantics></math> at <math display="inline"><semantics> <mrow> <mi>y</mi> <mo>=</mo> <mn>7</mn> </mrow> </semantics></math> km for the cloudy scenes shown in <a href="#remotesensing-13-00270-f001" class="html-fig">Figure 1</a> (a–d,f,h). "> Figure 3
A priori vertical profile of <math display="inline"><semantics> <msub> <mi>NO</mi> <mn>2</mn> </msub> </semantics></math> volume mixing ratio. "> Figure 4
Ring correction spectrum included in the retrieval. "> Figure 5
Residuals versus the iteration index for Differential Radiance Model with Internal closure (DRMI) (left) and Differential Radiance Model with External closure (DRME) (right). "> Figure 6
Correction spectra due to three-dimensional effects. ">

Versions Notes

Abstract

An algorithm for the retrieval of total column amount of trace gases in a multi-dimensional atmosphere is designed. The algorithm uses (i) certain differential radiance models with internal and external closures as inversion models, (ii) the iteratively regularized Gauss–Newton method as a regularization tool, and (iii) the spherical harmonics discrete ordinate method (SHDOM) as linearized radiative transfer model. For efficiency reasons, SHDOM is equipped with a spectral acceleration approach that combines the correlated k-distribution method with the principal component analysis. The algorithm is used to retrieve the total column amount of nitrogen for two- and three-dimensional cloudy scenes. Although for three-dimensional geometries, the computational time is high, the main concepts of the algorithm are correct and the retrieval results are accurate.

Keywords:

NO2 retrieval; SHDOM; multi-dimensional radiative trasnfer

1. Introduction

The retrieval of trace gas products from UV/VIS spectrometers is strongly affected by the presence of clouds. In general, there are three different cloud contributions: (i) the albedo effect, associated with the enhancement of reflectivity for cloudy scenes compared to cloud-free scenes; (ii) the shielding effect when a part of the trace gas column below the cloud is hidden; and (iii) the increase in absorption, related to multiple scattering inside clouds. The albedo and in-cloud absorption effects increase the detectability of trace gases at and above the cloud-top (in case of low clouds), while the shielding effect normally results in an underestimation of the trace gas column [1].

The majority of retrieval algorithms are based on one-dimensional radiative transfer modeling and use the so-called independent pixel approximation (IPA). The IPA for an atmosphere containing partial cloudiness means to compute separately the radiances for completely cloudy and clear skies, and then to express the partially cloudy radiance as a weighted sum of the separate radiances; the weighting factor being provided by the intensity weighted cloud fraction. The clouds within each pixel are assumed to be plane-parallel and homogeneous in both horizontal and vertical directions. The main advantage of the IPA is its computational efficiency, as it requires the solution of only two independent one-dimensional radiative transfer problems. The disadvantage is that for cloudy scenes of small horizontal extent, the errors due to the three-dimensional effects may be very significant [2,3,4]. This is the case of the new generation of atmospheric composition UV–VIS–NIR sensors, such as Sentinel-5P, Sentinel-4, and Sentinel-5. In [5] the importance of three-dimensional effects on computations of air mass factors has been shown.

To analyze the impact of cloud inhomogeneity on the retrieval, multi-dimensional radiative transfer models are used. In particular, for a specific cloudy scene, a spectral signal is simulated by a multi-dimensional radiative transfer model and then included in a one-dimensional retrieval algorithm to derive for example, the total column amount of a trace gas. By comparing the retrieved total column with the true value, a bias due to cloud effects can be deduced, and possible correction strategies can be explored.

In this paper, we go one step further and present the main concepts of an algorithm for retrieving the total column amount of trace gases in a multi-dimensional atmosphere. To achieve this goal we combine inverse models and regularization tools from inversion theory with a fast linearized version of the spherical harmonic discrete ordinate method [6]. These theoretical results are summarized in Section 2, while in Section 3, the accuracy and efficiency of the algorithm, used to retrieve the total column amount of nitrogen dioxide (

N O_{2}

) for two- and three-dimensional cloudy scenes, are numerically investigated. We conclude our analysis with some recommendations for operational retrievals.

2. Retrieval Algorithm

For a multi-dimensional geometry, the spectral signal of an instrument that measures the radiance at the top of the atmosphere in direction

Ω_{m}

at wavelength

λ

in the spectral interval

[λ_{\min}, λ_{\max}]

, is given by

\begin{matrix} I (λ) & = \int_{λ_{\min} - s / 2}^{λ_{\max} + s / 2} g (λ - λ^{'}) I_{FOV} (λ^{'}) d λ^{'}, \end{matrix}

(1)

where, under the assumption that the distance from the top of the atmosphere to the instrument is large,

I_{FOV} (λ^{'}) = \int_{S_{t}} h (r_{t}) I (λ^{'}, r_{t}, Ω_{m}) d S_{t},

(2)

is the signal integrated over the field of view (FOV) of the instrument at wavelength

λ^{'}

(see, e.g., in [7] for a general description),

g (λ)

the slit function (defined on all

R

), s the spectral bandwidth,

h (r_{t})

the characteristic function that takes the value

1 / A_{tm}

for

r_{t} \in S_{tm}

and 0 otherwise,

S_{tm}

the footprint of the instrument on the top face

S_{t}

, and

A_{tm}

the area of the instrument footprint. Essentially, Equations (1) and (2) state that the mean radiance across the wavelength and the area

S_{t}

at the top of the atmosphere is measured. The top of atmosphere radiance at the exiting point

r_{t}

along the characteristic

Ω_{m}

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

, is computed by using the integral form of the radiative transfer equation

\begin{matrix} I (λ, r_{t}, Ω_{m}) & = I (λ, r_{b}, Ω_{m}) e^{- \int_{0}^{s_{0}} σ_{ext} (λ, r_{b} + s Ω_{m}) d s} \\ + \int_{0}^{s_{0}} σ_{ext} (λ, r_{b} + s Ω_{m}) J (λ, r_{b} + s Ω_{m}, Ω_{m}) \\ \times [e^{- \int_{s}^{s_{0}} σ_{ext} (λ, r_{b} + s^{'} Ω_{m}) d s^{'}}] d s, \end{matrix}

(3)

where

I (λ, r_{b}, Ω_{m})

is the radiance at the surface point

r_{b}

J (λ, r, Ω_{m})

the source function, and

s_{0} = | r_{t} - r_{b} |

the distance between the entering and exiting points, respectively. Equation (3) can be derived by integrating the radiative transfer equation [8,9]. The first term corresponds to the upwelling radiance at the bottom of the atmosphere and attenuated exponentially as it propagates through the atmosphere. The second term refers to the gain due to sources. The source function J contains only the solar term and the multiple scattering term [10].

In the case of an atmosphere consisting of gas molecules and a cloud, we assume that (i) the optical coefficients of the gas molecules depend on the altitude level and the wavelength, while (ii) the optical coefficients of the cloud depend on the spatial coordinates but not on the wavelength. The second assumption is justified by the fact that a narrow spectral interval is considered in the retrieval. According to the additivity property of optical cross sections [11], in the case of

N_{g}

gases, the extinction coefficient is computed as

σ_{ext} (λ, r) = σ_{ext}^{cloud} (r) + σ_{sct}^{mol} (λ, z) + \sum_{g = 1}^{N_{g}} σ_{abs g}^{gas} (λ, z),

(4)

where

σ_{ext}^{cloud} (r)

is the extinction coefficient in the cloud (methods for computing

σ_{ext}^{cloud} (r)

are summarized in [12]),

σ_{sct}^{mol} (λ, z)

the molecular scattering coefficient due to Rayleigh scattering [13],

σ_{abs g}^{gas} (λ, z)

the absorption coefficient of gas g, and

(x, y, z)

the Cartesian coordinates of point

r

. Considering a discretization of the atmosphere in

N_{z}

levels, i.e.,

{z_{j}}_{j = 1}^{N_{z}}

, the absorption coefficient on level

z_{j}

is given by

σ_{abs g}^{gas} (λ, z_{j}) = C_{abs g} (λ, z_{j}) n_{g} (z_{j}),

(5)

where

C_{abs g} (λ, z_{j})

and

n_{g} (z_{j})

are the absorption cross section and the number density of gas j on level

z_{j}

, respectively. Under the assumption that on a layer j bounded below by

z_{j}

and above by

z_{j + 1}

, the number density

n_{g} (z)

varies linearly with respect to z, the partial and total column of gas g are defined, respectively, by

x_{g, j} = \int_{z_{j}}^{z_{j + 1}} n_{g} (z) d z \approx \frac{Δ z_{j}}{2} [n_{g} (z_{j}) + n_{g} (z_{j + 1})]

(6)

and

X_{g} = \int_{z_{1}}^{z_{N_{z}}} n_{g} (z) d z = \sum_{j = 1}^{N_{z} - 1} x_{g, j},

(7)

where

Δ z_{j} = z_{j + 1} - z_{j}

The retrieval algorithm is based on the assumption that cloud information are available from co-located imagers, as, for example, the Moderate Resolution Imaging Spectroradiometer (MODIS) on board the Terra and Aqua satellites and the Visible Infrared Imaging Radiometer Suite (VIIRS) on board the Suomi National Polar-Orbiting Partnership spacecraft. The main cloud products delivered by MODIS/VIIRS are the cloud mask (indicator function of the cloud field), cloud optical thickness, cloud-top height (cloud-top pressure), and effective radius of the size distribution. In the MODIS/VIIRS retrieval algorithm, the cloud geometrical thickness is derived under the assumption that the clouds are homogeneous in the vertical direction, or equivalently, that the liquid water content does not change with the altitude. Therefore, it is appropriate to use stochastic cloud models (in particular statistical cloud generators) which operate with the same assumption. For this reason, we will use stochastic cloud models relying on the vertical homogeneity assumption, and in particular, the two-dimensional broken cloud model of Alexandrov et al. [14].

In general, the design of a retrieval algorithm requires the specification of

an inverse model,
a solution method for solving the inverse problem, and
a linearized multi-dimensional radiative transfer model for computing the forward model and the Jacobian matrix at each iteration step.

2.1. Inverse Models

The retrieval of the total columns

X_{g}

N_{g}

gases can be performed by using two differential radiance models (see details in [15,16,17] and references therein):

the Differential Radiance Model with Internal closure (DRMI), in which the measured and simulated differential spectral signals are fitted, and
the Differential Radiance Model with External closure (DRME), in which the measured differential spectral signal and a simulated spectral signal with its smooth component extracted are fitted.

Actually, in DRMI we solve the nonlinear equation

\begin{matrix} R_{mes}^{δ} (λ_{k}) & = R_{sim} (λ_{k}, X) + \sum_{j = 1}^{N_{s}} b_{j} S_{j} (λ_{k}, X_{a}), k = 1, \dots, N_{λ}, \end{matrix}

(8)

where

R_{mes}^{δ} (λ_{k}) = ln I_{mes}^{δ} (λ_{k}) - P_{mes} (λ_{k}, c_{mes})

(9)

and

R_{sim} (λ_{k}, X) = ln I_{sim} (λ_{k}, X) - P_{sim} (λ_{k}, c_{sim} (X))

(10)

are the measured and simulated differential spectral signals, respectively, while in DRME we solve the nonlinear equation

\begin{matrix} R_{mes}^{δ} (λ_{k}) & = ln I_{sim} (λ_{k}, X) + \sum_{j = 1}^{N_{s}} b_{j} S_{j} (λ_{k}, X_{a}) - P (λ_{k}, c), k = 1, \dots, N_{λ} . \end{matrix}

(11)

In Equations (8)–(11),

X = [X_{1}, \dots, X_{N_{g}}]

is the vector comprising the total columns of trace gases;

X_{a}

an a priori estimate of

X

;

I_{mes}^{δ} (λ_{k})

the (noisy) spectral signal measured by the instrument at wavelength

λ_{k}

, where

k = 1, \dots N_{λ}

N_{λ}

the number of measurement wavelengths;

I_{sim} (λ_{k}, X)

the simulated spectral signal computed by a multi-dimensional radiative transfer model;

S_{j} (λ_{k}, X_{a})

with

j = 1, \dots, N_{s}

, the correction spectra including, for example, the polarization correction spectrum and the Ring spectrum,

N_{s}

the number of correction spectra; and

b_{j}

with

j = 1, \dots, N_{s}

, the (wavelength independent) amplitudes of the correction spectra. The polynomials

P_{mes} (λ, c_{mes})

P_{sim} (λ, c_{sim} (X))

, and

P (λ, c)

with degree

N_{p} - 1

account for the low-order spectral structure due to the scattering by clouds and aerosols, while the differential spectral signals

R_{sim}

and

R_{mes}^{δ}

are mostly due to the gaseous absorption processes (see Ref. [18] for details). Specifically, in DRMI, the coefficients

c_{mes} = [c_{mes 0}, \dots, c_{mes N_{p} - 1}] and c_{sim} (X) = [c_{sim 0} (X), \dots, c_{sim N_{p} - 1} (X)]

of the smoothing polynomials

P_{mes} (λ, c_{mes})

and

P_{sim} (λ, c_{sim} (X))

, respectively, are computed as the solutions of the least-squares problems

c_{mes} = arg min_{c} \sum_{k = 1}^{N_{λ}} {[ln I_{mes}^{δ} (λ_{k}) - P_{mes} (λ_{k}, c)]}^{2}

(12)

and

c_{sim} (X) = arg min_{c} \sum_{k = 1}^{N_{λ}} {[ln I_{sim} (λ_{k}, X) - P_{sim} (λ_{k}, c)]}^{2},

(13)

respectively. Thus, these coefficients, which are uniquely determined by

ln I_{mes}^{δ} (λ_{k})

and

ln I_{sim} (λ_{k}, X)

, are not a part of the retrieval, while in DRMI, the coefficients

c = [c_{0}, \dots, c_{N_{p} - 1}]

of the smoothing polynomial

P (λ, c)

are included in the retrieval.

Similar inversion models can be used for tropospheric column retrieval. These models together with the underlying retrieval algorithms are described in Appendix A.

To speed up the calculation, two approximate inverse models can be designed.

In DRME, we assume the linearized model [19]

$\begin{matrix} ln I_{sim} (λ_{k}, X) & = ln I_{sim} (λ_{k}, X_{a}) + \sum_{g = 1}^{N_{g}} (X_{g} - X_{a g}) \frac{\partial ln I_{sim}}{\partial X_{g}} (λ_{k}, X_{a}) + δ R_{p} (λ_{k}), \end{matrix}$

(14)

where the second term in the right-hand side of the equation describes the variations of $ln I_{sim}$ with respect to the total column amounts of trace gases $X_{g}$ with $g = 1, \dots, N_{g}$ , and the third term $δ R_{p} (λ_{k})$ comprises the contributions of all other atmospheric parameters. Approximating $δ R_{p} (λ_{k})$ by a smoothing polynomial, i.e., $δ R_{p} (λ_{k}) \approx P_{p} (λ_{k}, c_{p})$ , and putting $P (λ_{k}, c) - P_{p} (λ_{k}, c_{p}) \to P (λ_{k}, c)$ , we are led to the linear equation

$\begin{matrix} R_{mes}^{δ} (λ_{k}) & = ln I_{sim} (λ_{k}, X_{a}) + \sum_{g = 1}^{N_{g}} (X_{g} - X_{a g}) \frac{\partial ln I_{sim}}{\partial X_{g}} (λ_{k}, X_{a}) \\ + \sum_{j = 1}^{N_{s}} b_{j} S_{j} (λ_{k}, X_{a}) - P (λ_{k}, c), k = 1, \dots, N_{λ}, \end{matrix}$

(15)

which is equivalent with the Differential Optical Absorption Spectroscopy (DOAS) equation [18]. This equivalence is proved in Appendix B.
We adopt an approximate inverse model in which we compute the spectral signal by an one-dimensional radiative transfer model, and (similar to the Ring and polarization correction spectra) introduce a correction spectrum $S_{3 D} (λ_{k}, X_{a})$ that accounts on three-dimensional effects. To define the correction spectrum, we take into account that, for example, in DRME, we can write

$\begin{matrix} ln I_{sim}^{3 D} (λ_{k}, X) - P (λ_{k}, c) \\ = ln I_{sim}^{1 D} (λ_{k}, X) + ln [1 + δ I_{sim}^{3 D} (λ_{k}, X)] - P (λ_{k}, c) \\ \approx ln I_{sim}^{1 D} (λ_{k}, X) + δ I_{sim}^{3 D} (λ_{k}, X) - P (λ_{k}, c) \\ = ln I_{sim}^{1 D} (λ_{k}, X) + [δ I_{sim}^{3 D} (λ_{k}, X) - P_{3 D} (λ_{k}, c_{3 D} (X))] - P_{1 D} (λ_{k}, c_{1 D}), \end{matrix}$

(16)

where $I_{sim}^{3 D} (λ_{k}, X)$ and $I_{sim}^{1 D} (λ_{k}, X)$ , are the spectral signals computed by a three- and an one-dimensional radiative transfer model, respectively, and

$δ I_{sim}^{3 D} (λ_{k}, X) = \frac{I_{sim}^{3 D} (λ_{k}, X) - I_{sim}^{1 D} (λ_{k}, X)}{I_{sim}^{1 D} (λ_{k}, X)},$

(17)

is a relative correction spectrum. In deriving Equation (16), we used the first-order Taylor approximation $ln (1 + δ I_{sim}^{3 D}) \approx δ I_{sim}^{3 D}$ for $δ I_{sim}^{3 D} ≪ 1$ , and the decomposition $P (λ_{k}, c) = P_{3 D} (λ_{k}, c_{3 D} (X) + P_{1 D} (λ_{k}, c_{1 D})$ , in which, the coefficients $c_{3 D} (X)$ are uniquely determined by $δ I_{sim}^{3 D} (λ_{k}, X)$ (see below), while the coefficients $c_{1 D}$ are unknown (as the coefficients $c$ are). In view of Equation (16), we assume that in the retrieval, the differential spectral term

$δ I_{sim}^{3 D} (λ_{k}, X) - P_{3 D} (λ_{k}, c_{3 D} (X))$

can be described by $b_{3 D} S_{3 D} (λ_{k}, X_{a})$ , where

$S_{3 D} (λ_{k}, X_{a}) = δ I_{sim}^{3 D} (λ_{k}, X_{a}) - P_{3 D} (λ_{k}, c_{3 D} (X_{a}))$

(18)

is the (differential) correction spectrum due to three-dimensional effects, $P_{3 D} (λ_{k}, c_{3 D} (X_{a}))$ the associated smoothing polynomial, and $b_{3 D}$ the amplitude of the correction spectrum (to be included in the retrieval). From Equation (18), we see that the correction spectrum $S_{3 D} (λ_{k}, X_{a})$ reproduces the differential structure of the relative correction spectrum $δ I_{sim}^{3 D} (λ_{k}, X_{a})$ .

2.2. Solution Method

The nonlinear Equations (8) and (11) can be written in compact form as

y^{δ} = F (x),

(19)

where

F : R^{n} \to R^{m}

is the forward model,

x \in R^{n}

the state vector, and

y^{δ} \in R^{m}

the noisy data vector. In DRMI and DRME, the state vectors are

x = {[X, b]}^{T}

and

x = {[X, b, c]}^{T}

, respectively, where

b = [b_{1}, \dots, b_{N_{s}}]

is the vector comprising the amplitudes of the correction spectra

b_{j}

j = 1, \dots, N_{s}

. Furthermore, the number of data points is

m = N_{λ}

, and the measured differential spectral signals affected by noise

R_{mes}^{δ} (λ_{k})

for

k = 1, \dots N_{λ}

are the components of the noisy data vector

y^{δ}

Because the nonlinear Equation (19) is usually ill-posed, a regularization method, as, for example, the iteratively regularized Gauss–Newton method, is used to compute a solution with physical meaning [15,20,21]. In this approach, at the iteration step k, we consider a linearization of

F (x)

around the current iterate

x_{k}^{δ}

and solve the linearized equation by means of Tikhonov regularization with the penalty term

α_{k} | | L (x - x_{a}) {| |}^{2}

, where

α_{k}

is the regularization parameter at the iteration step k,

L

the regularization matrix, and

x_{a}

the a priori state vector, the best beforehand estimate of the solution. Thus, the new iterate minimizing the function

F_{l k} (x) = {∥y_{k}^{δ} - K_{k} (x - x_{a})∥}^{2} + α_{k} {∥L (x - x_{a})∥}^{2},

(20)

is given by

x_{k + 1}^{δ} = x_{a} + K_{k}^{†} y_{k}^{δ},

(21)

where

K_{k}^{†} = {(K_{k}^{T} K_{k} + α_{k} L^{T} L)}^{- 1} K_{k}^{T}

(22)

is the regularized generalized inverse,

K_{k} = K (x_{k}^{δ}) = \partial F (x_{k}^{δ}) / \partial x

the Jacobian matrix evaluated at

x_{k}^{δ}

, and

y_{k}^{δ} = y^{δ} - F (x_{k}^{δ}) + K_{k} (x_{k}^{δ} - x_{a})

(23)

the linearized data vector at the iteration step k. The following peculiarities of the iteratively regularized Gauss–Newton method deserve to be mentioned [15].

In contrast to the method of Tikhonov regularization [22], the regularization parameter is not constant during the iterative process. Instead, the regularization parameters $α_{k}$ are the terms of a decreasing (geometric) sequence, i.e., $α_{k} = q α_{k - 1}$ with $q < 1$ . In this way, the amount of regularization is gradually decreased during the iterative process.
For iterative regularization methods, the number of iteration steps k plays the role of the regularization parameter, and the iterative process is stopped after an appropriate number of steps $k^{☆}$ in order to avoid an uncontrolled expansion of the errors in the data. The stopping rule used in this study is the discrepancy principle [23], according to which, the iterative process is terminated after $k^{☆}$ steps such that

$∥r_{k^{☆}}^{δ}∥ \leq τ Δ < ∥r_{k}^{δ}∥, 0 \leq k < k^{☆},$

(24)

where $r_{k}^{δ} = y^{δ} - F (x_{k}^{δ})$ is the residual vector at $x_{k}^{δ}$ , $τ > 1$ a control parameter, and $Δ$ the noise level (an upper bound for the noise in the data). Because in practice the noise level cannot be a priori estimated, we adopt a practical approach based on the observation that the residual $| | r_{k}^{δ} | |$ decreases during the iterative process and attains a plateau at approximately $Δ$ . Thus, if the nonlinear residuals $| | r_{k}^{δ} | |$ converge to $| | r_{\infty}^{δ} | |$ within a prescribed tolerance, we use the estimate $Δ = | | r_{\infty}^{δ} | |$ .
The regularization matrix is chosen as a diagonal matrix, that is, the penalty term $α_{k} | | L (x - x_{a}) {| |}^{2}$ is taken as

$α_{k} [\sum_{g = 1}^{N_{g}} \frac{w_{g}}{X_{a g}} {(X_{g} - X_{a g})}^{2} + \sum_{j = 1}^{N_{s}} \frac{w_{b j}}{b_{a j}} {(b_{j} - b_{a j})}^{2}]$

(25)

for DRMI and

$α_{k} [\sum_{g = 1}^{N_{g}} \frac{w_{g}}{X_{a g}} {(X_{g} - X_{a g})}^{2} + \sum_{j = 1}^{N_{s}} \frac{w_{b j}}{b_{a j}} {(b_{j} - b_{a j})}^{2} + \sum_{p = 0}^{N_{p} - 1} \frac{w_{c p}}{c_{a p}} {(c_{p} - c_{a p})}^{2}]$

(26)

for DRME. Here, the scalars ${w_{g}}_{g = 1}^{N_{g}}$ , ${w_{b j}}_{j = 1}^{N_{s}}$ , and ${w_{c p}}_{p = 0}^{N_{p} - 1}$ give the weight of each component of the state vector into the regularization matrix.

As compared to the method of Tikhonov regularization, the iteratively regularized Gauss–Newton method is less sensitive to overestimations of the regularization parameter, but requires more iteration steps to achieve convergence.

We conclude this section by mentioning that in the case of DRME, one iteration step of the iteratively regularized Gauss–Newton method corresponds to the linear problem described by Equations (11) and (14).

2.3. Linearized Radiative Transfer Model

In order to solve the nonlinear Equation (19) by means of the iteratively Gauss–Newton method, the forward model

F (x)

and the Jacobian matrix

K (x)

have to be calculated at each iteration step. From Equation (3) we infer that the computation of

F (x)

requires the computation of the radiance

I (λ, r_{b}, Ω)

and the source function

J (λ, r, Ω)

at all surface and grid points, respectively, while the computation of

K (x)

requires the computation of the partial derivatives of the extinction field

\partial σ_{ext} (\cdot, r) / \partial X_{g}

, the surface radiance

\partial I (\cdot, r_{b}, \cdot) / \partial X_{g}

, and the source function

\partial J (\cdot, r, \cdot) / \partial X_{g}

For derivative calculations, we employ the same assumption as in the standard DOAS model, that is, the (discrete) number density profile

{n_{g, j}}_{j = 1}^{N_{z}}

, where

n_{g, j} = n_{g} (z_{j})

, is assumed to be a scaled version of an a priori profile

{n_{a g, j}}_{j = 1}^{N_{z}}

. Consequently, if

s_{g}

is the scale factor of gas g, i.e.,

n_{g, j} = s_{g} n_{a g, j}

, for all

j = 1, \dots, N_{z}

, the following equalities,

s_{g} = \frac{n_{g, j}}{n_{a g, j}} = \frac{σ_{abs g, j}^{gas} (λ)}{σ_{abs a g, j}^{gas} (λ)} = \frac{x_{g, j}}{x_{a g, j}} = \frac{X_{g}}{X_{a g}}

(27)

with

σ_{abs g, j}^{gas} (λ) = σ_{abs g}^{gas} (λ, z_{j})

, are valid. Thus, for any value of the total column

X_{g}

, the absorption coefficient is updated as

σ_{abs g, j}^{gas} (λ) = \frac{σ_{abs a g, j}^{gas} (λ)}{X_{a g}} X_{g},

(28)

and the corresponding top of atmosphere radiance

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

is determined by means of a radiative transfer calculation. In general, for computing the derivative of a spectral quantity

F (λ)

with respect to the total column

X_{g}

, we regard

F (λ)

as a function of

σ_{abs g, j}^{gas} (λ)

, i.e.,

F = F (σ_{abs g, 1}^{gas}, \dots, σ_{abs g, N_{z}}^{gas})

, and assume that the partial derivatives

\partial F / \partial σ_{abs g, j}^{gas}

j = 1, \dots, N_{z}

, are known. Taking into account the one-to-one correspondence (28), we are led to the computational formula

\frac{\partial F}{\partial X_{g}} (λ) = \sum_{j = 1}^{N_{z}} \frac{\partial F}{\partial σ_{abs g, j}^{gas}} (λ) \frac{\partial σ_{abs g, j}^{gas}}{\partial X_{g}} (λ) = \sum_{j = 1}^{N_{z}} \frac{σ_{abs a g, j}^{gas} (λ)}{X_{a g}} \frac{\partial F}{\partial σ_{abs g, j}^{gas}} (λ) .

(29)

Thus, the radiative transfer model should be linearized with respect to the absorption coefficient, that is, the linearized radiative transfer model should deliver the partial derivatives

\partial F / \partial σ_{abs g, j}^{gas}

for

j = 1, \dots, N_{z}

and

g = 1, \dots, N_{g}

The multi-dimensional radiative transfer model used in our retrieval algorithm is the spherical harmonic discrete ordinate method (SHDOM) [6,24]. Linearizations of SHDOM by means of a forward and a forward-adjoint approach have been discussed in [25]. These can be summarized as follows.

The linearized forward approach relies on an analytical computation of the derivatives. The method is accurate and has the advantage that no assumptions rather than those of the forward model have to be made. However, the method is time consuming and memory demanding when the number of parameters to be retrieved is large. 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.
The linearized forward-adjoint approach relies on the application of the adjoint radiative transfer theory. The method requires less storage for derivatives calculation, is 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.

Because in our application the number of retrieved quantities is relatively small, we decided in the favor of the linearized forward approach. Specifically, the partial derivatives are computed in an analytical manner by using the chain rule in a sort of backward procedure starting with the output radiance and following the chain of dependencies back to the inputs to the model.

In [26], a simplified approach, leading to a considerable reduction of the computation time, was proposed. The idea is to neglect in Equation (3) the contributions of the partial derivatives of the surface radiance

\partial I (\cdot, r_{b}, \cdot) / \partial X_{g}

and the source function

\partial J (\cdot, r, \cdot) / \partial X_{g}

. In other words, when computing the partial derivative of the top of atmosphere radiance

\partial I (\cdot, r_{t}, \cdot) / \partial X_{g}

, only the contributions of the partial derivatives of the extinction field

\partial σ_{ext} (\cdot, r) / \partial X_{g}

are taken into account.

To speed up the computations, the linearized SHDOM is equipped with a spectral acceleration approach that combines the correlated k-distribution method with dimensionality reduction techniques. This approach was applied for computing the spectral signal in [27], and can be extended to derivative calculations as follows.

Correlated k-distribution method. Consider a discretization of the spectral interval $[λ_{\min} - s / 2, λ_{\max} + s / 2]$ into a set ${{\bar{λ}}_{k}}_{k = 1}^{\bar{W}}$ of $\bar{W}$ equally spaced wavelengths with the discretization step $Δ λ$ , and assume that the transmission within a spectral interval depends only on the distribution of the gas absorption coefficient $σ_{abs}^{gas} (λ)$ within the spectral interval [28]. Letting $F = F (σ_{abs k}^{gas})$ be the cumulative density function of $σ_{abs}^{gas} (λ)$ in the spectral interval $[{\bar{λ}}_{k} - Δ λ / 2, {\bar{λ}}_{k} + Δ λ / 2]$ , $σ_{abs k}^{gas} (F)$ the inverse distribution function, and ${F_{l}, ϖ_{l}}_{l = 1}^{N_{q}}$ a set of $N_{q}$ quadrature points and weights in the interval $[0, 1]$ , the spectral signal (1) and its partial derivative with respect to the total column are computed as

$\begin{matrix} I (λ) & = \sum_{w = 1}^{W} ω_{w} g (λ - λ_{w}) I_{FOV} (σ_{abs}^{gas} (λ_{w})), \end{matrix}$

(30)

$\begin{matrix} \frac{\partial I}{\partial X_{g}} (λ) & = \sum_{w = 1}^{W} ω_{w} g (λ - λ_{w}) \frac{\partial I_{FOV}}{\partial X_{g}} (σ_{abs}^{gas} (λ_{w})), \end{matrix}$

(31)

where $λ_{w} = {\bar{λ}}_{k}$ , $ω_{w} = Δ λ ϖ_{l}$ , and $σ_{abs}^{gas} (λ_{w}) = σ_{abs k}^{gas} (F_{l})$ for $w = l + (k - 1) N_{q}$ , $k = 1, \dots, \bar{W}$ , $l = 1, \dots, N_{q}$ , and $W = \bar{W} N_{q} .$ Thus, in the framework of the correlated k-distribution method, W monchromatic radiative transfer calculations are required for computing $I_{FOV} (λ_{w})$ and $\partial I_{FOV} (λ_{w}) / \partial X_{g}$ , and so, $I (λ)$ and $\partial I (λ) / \partial X_{g}$ . A further acceleration can be achieved when $I_{FOV} (λ_{w})$ and $\partial I_{FOV} (λ_{w}) / \partial X_{g}$ are computed by using dimensionality reduction techniques, as for example, the principal component analysis.
Principal component analysis. At wavelength $λ_{w}$ , the integrated signal $I_{FOV} (λ_{w})$ is related to the integrated signal calculated by a simplified (approximate) radiative transfer model $I_{FOV}^{s} (λ_{w})$ through the relation

$ln \frac{I_{FOV} (λ_{w})}{I_{FOV}^{s} (λ_{w})} = f (λ_{w}),$

(32)

i.e.,

$I_{FOV} (λ_{w}) = I_{FOV}^{s} (λ_{w}) e^{f (λ_{w})} .$

(33)

The correction factor $f (λ_{w})$ in Equation (32) is actually the quantity that is calculated by means of the principal component analysis. To summarize this approach, we assume that for each wavelength $λ_{w}$ , the spectral variability of the optical parameters can be described by a vector $x_{w} \in R^{N}$ , defined by

$x_{w}^{T} = [ln σ_{abs 1}^{gas} (λ_{w}), . . ., ln σ_{abs L}^{gas} (λ_{w}), ln σ_{sct 1}^{mol} (λ_{w}), . . ., ln σ_{sct L}^{mol} (λ_{w}), ln F_{0} (λ_{w})],$

(34)

where $σ_{abs k}^{gas}$ and $σ_{sct k}^{mol}$ are the optical coefficients in the kth level, $N = 2 N_{z} + 1$ , and $N_{z}$ is the number of altitude levels. Denoting by $\bar{x} = (1 / W) \sum_{w = 1}^{W} x_{w}$ the sample mean of the data, the aim is to find an M-dimensional subspace ( $M < N$ ) spanned by a set of linear independent vectors ${a_{k}}_{k = 1}^{M}$ , such that the centered data $x_{w} - \bar{x}$ belong to this subspace, i.e., $x_{w} \approx \bar{x} + \sum_{k = 1}^{M} y_{w k} a_{k} = \bar{x} + A y_{w}$ yielding $y_{w} = A^{†} (x_{w} - \bar{x})$ for $w = 1, \dots, W$ , where $A = {[a_{k}]}_{k = 1}^{M} \in R^{N \times M}$ , $A^{†} = {(A^{T} A)}^{- 1} A^{T} \in R^{M \times N}$ is the pseudoinverse of $A$ , and $y_{w k}$ is the kth component of the vector of parameters $y_{w} \in R^{M}$ . In the principal component analysis, the original N-dimensional data $X = {[x_{w}]}_{w = 1}^{W} \in R^{N \times W}$ are projected on the M-dimensional subspace spanned by the dominant singular vectors of the data covariance matrix $C_{x} = (1 / W) X X^{T} \in R^{N \times N}$ , that is, with $σ_{k}$ and $u_{k}$ being the kth singular value and vector of the matrix $C_{x}$ , we choose $A = U_{M} Σ_{M}^{1 / 2}$ and $A^{†} = Σ_{M}^{- 1 / 2} U_{M}^{T}$ , where $U_{M} = {[u_{k}]}_{k = 1}^{M} \in R^{N \times M}$ and $Σ_{M} = diag {[σ_{k}]}_{k = 1}^{M} \in R^{M \times M}$ . Furthermore, approximating the correction factor $f (x_{w})$ by a second-order Taylor expansion around $\bar{x}$ , and the gradient and the Hessian of f by central differences, we are led to the computational formula

$\begin{matrix} f (x_{w}) & \approx f (\bar{x}) + \frac{1}{2} \sum_{k = 1}^{M} [f (\bar{x} + a_{k}) - f (\bar{x} - a_{k})] y_{w k} \\ + \frac{1}{2} \sum_{k = 1}^{M} [f (\bar{x} + a_{k}) - 2 f (\bar{x}) + f (\bar{x} - a_{k})] y_{w k}^{2} . \end{matrix}$

(35)

To compute the derivative of the integrated signal with respect to the total column $X_{g}$ , we may use the principal component analysis to calculate the derivative correction factor $f_{X_{g}} (λ_{w})$ , defined by [29]

$f_{X_{g}} (λ_{w}) = ln [\frac{\partial I_{FOV}}{\partial X_{g}} (λ_{w}) / \frac{\partial I_{FOV}^{s}}{\partial X_{g}} (λ_{w})],$

(36)

in which case, we have

$\frac{\partial I_{FOV}}{\partial X_{g}} (λ_{w}) = \frac{\partial I_{FOV}^{s}}{\partial X_{g}} (λ_{w}) e^{f_{X_{g}} (λ_{w})} .$

(37)

Alternatively, taking the derivative of Equation (33) with respect to $X_{g}$ , i.e.,

$\frac{\partial I_{FOV}}{\partial X_{g}} (λ_{w}) = \frac{\partial I_{FOV}^{s}}{\partial X_{g}} (λ_{w}) e^{f (λ_{w})} + I_{FOV} (λ_{w}) \frac{d f}{d X_{g}} (λ_{w}),$

(38)

we may use the principal component analysis to calculate the derivative correction factor $f_{X_{g}} (λ_{w})$ , defined by

$f_{X_{g}} (λ_{w}) = \frac{d f}{d X_{g}} (λ_{w}) = \frac{1}{I_{FOV} (λ_{w})} \frac{\partial I_{FOV}}{\partial X_{g}} (λ_{w}) - \frac{1}{I_{FOV}^{s} (λ_{w})} \frac{\partial I_{FOV}^{s}}{\partial X_{g}} (λ_{w}) .$

(39)

3. Numerical Simulations

In this section, we analyze the accuracy and efficiency of the algorithms for retrieving total column amount of

{NO}_{2}

under cloudy conditions. The cloudy scenes considered in our simulations are similar to those described in [27]. Specifically, the geometrical and optical parameters are chosen as follows.

The domain of analysis is a rectangular prism of lengths $L_{x} = L_{y} = 15$ km and $L_{z} = 50$ km. The discretization steps along the horizontal directions are $Δ x = Δ y = 0.5$ km. Along the vertical direction, the atmosphere between 0 and 50 km is discretized with a step of 0.5 km between 0 and 3 km, 0.1 km between 3 and 4 km, 0.5 km between 4 and 10 km, 1.0 km between 10 and 14 km, 2 km between 14 and 30 km, and 5 km between 30 and 50 km.
A homogeneous cloud is placed between 3 and 4 km. The cloud extinction field is given by $σ_{ext}^{cloud} (x, y) = σ_{0} f (x, y)$ , where $σ_{0} = 6 {km}^{- 1}$ and $f (x, y)$ is the indicator function (note that $f (x, y)$ takes the values 1 and 0 inside and outside the cloud, respectively). The cloud phase function is a Henyey–Greenstein phase function [30] with the asymmetry parameter $g = 0.8$ , and the cloud single-scattering albedo is $ω_{cloud} =$ 0.99. Eight cloudy scenes are generated by a two-dimensional broken cloud model [14] with a cloud fraction of about $0.4 .$ The extinction field $σ_{ext}^{cloud} (x, y)$ is smoothed at the boundary of a cloudy region in order to avoid abrupt changes in the horizontal plane. The indicator functions corresponding to the eight cloudy scenes are illustrated in Figure 1a–h. For two-dimensional geometries, slices at $y = 7$ km are selected from the cloud fields 1, 2, 3, 5, 6, and 8. The corresponding indicator functions are shown in Figure 2. Note that the slices corresponding to the cloudy scenes 4 and 7 are similar to the slice corresponding to the cloudy scene 2, and are therefore omitted.
The number of discrete zenith and azimuth angles are $N_{μ} = 16$ and $N_{φ} = 2 N_{μ}$ , respectively, the solar and instrument zenith angles are $θ_{0} = 30^{\circ}$ and $θ_{m}^{0} = 45^{\circ}$ , respectively, and the relative azimuth angle is $Δ φ = 0$ .
A Lambertian reflecting surface with the surface albedo $A = 0.2$ is considered.
The footprint of the detector is a square of length $2 a = 10 Δ x$ centered at $x_{0} = y_{0} = L_{x} / 2$ , and $z_{0} = L_{z}$ , and a wavelength-dependent slit function corresponding to the TROPOMI instrument is assumed.
In addition to the scattering and absorption by the cloud, molecular Rayleigh scattering and the absorption by ${NO}_{2}$ , ozone ( $O_{3}$ ), oxygen dimer ( $O_{4}$ ), and water vapor ( $H_{2} O$ ) are considered. The measurement spectral grid roughly resembles the TROPOMI’s spectral resolution and consists of 119 spectral points between 425 nm and 450 nm.

The process of generating synthetic measurement spectra is organized as follows.

For a clean scenario, we use the a priori partial column profile of ${NO}_{2}$ [17] illustrated in Figure 3 to generate the true (exact) partial column profile. In this regard, denoting the a priori partial columns of gas g by $x_{a g, j}$ , we choose the true partial columns as $x_{g, j}^{†} = s_{g} x_{a g, j}$ , $j = 1, \dots, N_{z}$ with $s_{{NO}_{2}} = 2.0$ and $s_{O_{3}} = s_{O_{4}} = s_{H_{2} O} = 1.2$ . The true total column $X_{g}^{†}$ of gas g is then computed as $X_{g}^{†} = \sum_{j = 1}^{N_{z}} x_{g, j}^{†} = s_{g} X_{a g}$ ; thus, $X_{{NO}_{2}}^{†} = 2 X_{a {NO}_{2}}$ .
For $X^{†} = [X_{{NO}_{2}}^{†}, X_{O_{3}}^{†}, X_{O_{4}}^{†}, X_{H_{2} O}^{†}]$ , we generate the simulated spectral signal $I_{sim} (λ_{k}, X^{†})$ by means of SHDOM.
For cubic smoothing polynomials, i.e., $N_{p} = 4$ , we determine the coefficients $c_{sim} (X^{†})$ of the polynomial $P_{sim} (λ, c_{sim} (X^{†}))$ as the solutions of the least-squares problem (13).
We compute the noisy spectral signal as $I_{sim}^{δ} (λ_{k}, X^{†}) = I_{sim} (λ_{k}, X^{†}) + δ_{k}$ , where the measurement errors $δ_{k}$ are assumed to be independent Gaussian random variables with zero mean and standard deviation $σ_{k} = I_{sim} (λ_{k}, X^{†}) / SNR$ , where SNR is the signal-to-noise ratio. It should be pointed out that in view of the approximation

$\begin{matrix} ln I_{sim}^{δ} (λ_{k}, X^{†}) \approx ln I_{sim} (λ_{k}, X^{†}) + \frac{δ_{k}}{I_{sim} (λ_{k}, X^{†})}, \end{matrix}$

the error in $ln I_{sim}^{δ} (λ_{k}, X^{†})$ is $δ_{ln k} = δ_{k} / I_{sim} (λ_{k}, X^{†})$ yielding $σ_{ln k} = 1 / SNR$ for all k. In other words, the measurement error vector is white noise with the covariance matrix $(1 / {SNR}^{2}) I_{m}$ , where $I_{m}$ is the identity matrix. Because in our simulations we are mainly interested in multi-dimensional effects, we take $SNR = 10^{4}$ , that is, we assume an almost perfect instrument and neglect the forward model errors.
We include the Ring correction spectrum $S_{R} (λ_{k}, X_{a})$ illustrated in Figure 4 in the retrieval, and choose the a priori and true Ring amplitudes as $b_{aR} = 5 \times 10^{- 2}$ and $b_{R}^{†} = 2 b_{aR}$ , respectively. Note that the inelastic scattering is described by a first-order Rayleigh scattering model, i.e., by applying a first-order iteration scheme to the one-dimensional radiative transfer equation for inelastic scattering [31].
For DRMI, we compute the measured differential spectral signal as

$R_{mes}^{δ} (λ_{k}) = [ln I_{sim}^{δ} (λ_{k}, X^{†}) - P_{sim} (λ, c_{sim} (X^{†}))] + b_{R}^{†} S_{R} (λ_{k}, X_{a}),$

while for DRME, we choose $c^{†} = 0.5 c_{sim} (X^{†})$ and compute the measured differential spectral signal as

$R_{mes}^{δ} (λ_{k}) = ln I_{sim}^{δ} (λ_{k}, X^{†}) + b_{R}^{†} S_{R} (λ_{k}, X_{a}) - P (λ_{k}, c^{†}) .$

SHDOM is run with a solution accuracy of

10^{- 4}

and by using

an adaptive grid with a splitting accuracy of $10^{- 4}$ ,
the principal component analysis with $M = 4$ , and the derivative correction factor $f_{X_{g}} (λ_{w})$ as in Equation (39),
periodic boundary conditions,
the delta-M approximation [32], and the untruncated phase function single-scattering solution (TMS correction of Nakajima and Tanaka [33].

Some parameters characterizing the iteratively regularized Gauss–Newton algorithm are chosen as follows.

The regularization parameters $α_{k}$ are the terms of a geometric sequence with ratio $q = 0.1$ and initial value $α_{0} = 10^{- 3}$ . Thus, at the first iteration step, the regularization parameter is $α_{1} = 10^{- 4}$ .
The weighting factors specifying the contribution of each component of the state vector into the regularization matrix are $w_{{NO}_{2}} = 1.0$ , $w_{O_{3}} = w_{O_{4}} = w_{H_{2} O} = 10^{2}$ , $w_{R} = 10^{- 3}$ , and $w_{c p} = 1.0$ for all $p = 0, \dots, N_{p} - 1$ . By this choice, the total columns of the auxiliary gases are stronger constrained to the a priori than the total column of ${NO}_{2}$ and the Ring correction spectrum.
The control parameter $τ$ in the discrepancy principle Equation (24) is $τ = 1.2$ .

The simulations were performed on a computer Intel(R) Core(TM) i5-3340M CPU @ 2.70 GHz with 7858 MB RAM.

3.1. Test Example 1

In the first test example we consider a two-dimensional geometry and perform the retrieval by using the DRMI and DRME inverse models. It may in general be noted that

the inverse problem is severely ill-posed (for the DRMI model, the condition number of the Jacobian matrix at the initial guess is $κ_{DRMI} = 2.92 \times 10^{7}$ ) and
there is a strong correlation between the ${NO}_{2}$ and Ring effect signatures (when the Ring correction spectrum is not included in the retrieval, the condition number is $κ_{DRMI} = 2.16 \times 10^{4}$ ; thus, $κ_{DRMI}$ decreases by three order of magnitude).

The histories of the residuals with respect to the iteration index are shown in Figure 5. The following conclusions can be drawn.

At the initial guess, the residual corresponding to DRMI is much smaller than that corresponding to DRME. This occurs because the discrepancies between the differential spectra are usually small.
In DRMI, the residual decreases very fast at the first iteration step and then more steadily, while in DRME, the residual gradually decreases.

The relative errors and the computation times for DRMI and DRME are given in the first two columns of Table 1 and Table 2, respectively. Because the SNR is high, in both cases, the relative errors are smaller than

0.01 %

. On average, the computation times for DRMI and DRME are of about 40 and 35 CPU minutes, respectively. The higher computation time for DRMI is due to the fact that, according to Equation (13), the partial derivative of

c_{sim} (X)

with respect to

X

needs also to be computed.

In the second two columns of Table 1, we list the relative errors and the computation times for DRMI when the partial derivative of the top of atmosphere radiance

\partial I (\cdot, r_{t}, \cdot) / \partial X_{g}

is computed by taking into account only the contributions of the partial derivatives of the extinction field

\partial σ_{ext} (\cdot, r) / \partial X_{g}

. The reduction of the computation time by about 25% does not justify the large relative errors. Thus, the simplified approach suggested in [26] for the retrieval of cloud extinction field, seems not to work in the case of total column trace gas retrieval.

In the second two columns of Table 2, we list the relative errors and the computation times for DRME when the iterative process is stopped after one iteration, that is, when the linear model (14) is considered. A reasonable accuracy can be obtained when the regularization parameter

α

is optimally chosen. Some hints regarding the selection of the optimal value of the regularization parameter

α_{opt}

can be found by analyzing the history of the residual curve in Figure 5. From this curve, we infer that for

α = 10^{- 4}

the residual is of about

10^{- 5}

, for

α = 10^{- 5}

of about

10^{- 6}

, and for

α = 10^{- 6}

of about

10^{- 8}

. Therefore, by choosing

α

in the range

10^{- 6}, \dots, 10^{- 5}

we expect good reconstruction errors. However, as it can be seen in Table 2,

α_{opt}

depends on the cloud field. Moreover, for some cloudy scenes,

α_{opt}

should be precisely determined (for example, in the case of the cloudy scene 8, the relative errors are

8.30 \times 10^{- 2}

1.87 \times 10^{- 2}

, and

9.02 \times 10^{- 2}

when

α

takes the values

1 \times 10^{- 5}

2 \times 10^{- 5}

, and

3 \times 10^{- 5}

, respectively). In this regard, it should be mentioned that although less inefficient, the iteratively regularized Gauss–Newton method is the best option, i.e., when the initial value of regularization parameter

α_{0}

is chosen sufficiently large, accurate results will always be found in comparatively small number of iteration steps (which depends on the initial value

α_{0}

and the ratio q).

3.2. Test Example 2

For the cloudy scenes considered in Figure 1, we illustrate, in Table 3, the relative errors and the computation times for DRME with several and one iteration steps. The solution accuracies are comparable to those corresponding to two-dimensional geometries. However, the computation time is extremely high; on average, it is of about 14 h and 15 min for DRME with several iteration steps, and 2 h and 30 min for DRME with one iteration step.

By taking a closer look at the cloudy scenes depicted in Figure 1, we see that the footprint of the detector is almost cloud-free. In this regard, we may try to retrieve the total column amount by using an one-dimensional radiative transfer model for a cloud-free atmosphere. The results given in the first two columns of Table 4 and Table 5 show that in this case, the relative error is about 15%. In these simulations, the spectral signal measured by the instrument is computed by a three-dimensional radiative transfer model, while in the retrieval, the simulated spectral signal is computed by an one-dimensional radiative transfer model. The relative errors can be approximately halved when the differential correction spectrum due to three-dimensional effects (18) is included in the inverse model and its amplitude is retrieved. This is shown in the second two columns of Table 4 and Table 5. Note that the computation time given in these tables does not take into account the time for computing the correction spectrum by a three-dimensional radiative transfer model. As shown in [27], this computation time is of about 1 h and 20 min, so that roughly speaking, this is the time for retrieving the total column amount by an one-dimensional inverse model that accounts on three-dimensional effects. The correction spectra included in the inverse model are shown in Figure 6, and it is apparent that they depend on the cloud field (some similitude can be observed for the cloudy scenes 2 and 7).

4. Discussion

The main concepts of an algorithm for the retrieval of the total column amount of trace gases in a multi-dimensional atmosphere have been presented. The algorithm uses (i) as inverse models, certain differential radiance models with internal and external closures; (ii) as regularization tool, the iteratively regularized Gauss-Newton method; and (iii) as linearized radiative transfer model, the spherical harmonics discrete ordinate method with a spectral acceleration approach that combines the correlated k-distribution method with the principal component analysis.

The algorithm has been applied to the retrieval of the total column amount of

{NO}_{2}

under cloudy conditions. For two-dimensional geometries, we come to the following findings.

The differential radiance models with internal and external closures yield accurate results with reasonable computation times (of about 35–40 min).
An inverse model based on an approximate computation of the partial derivative leads to a reduction of the computation time by about 25%, but to large relative errors.
Provided that the regularization parameter is optimally chosen, reasonable accurate results with a computation time of about 6 min can be obtained when the iterative process, corresponding to the differential radiance model with external closure, is stopped after one iteration step. The fact that the optimal value of the regularization parameter depends on the cloudy scene, makes it more difficult to apply this one-step retrieval algorithm, or equivalently, DOAS-type models.

For three-dimensional geometries, we come to the following results.

Although accurate, the retrievals based on the differential radiance model with external closures are inefficient; the computation time is of about 14 h and 15 min for a full-step retrieval algorithm, and 2 h and 30 min for an one-step retrieval algorithm. The one-step retrieval algorithm is less accurate than in the case of two-dimensional geometries, but the results are still acceptable.
The application of a fast one-dimensional radiative transfer model to retrieve the ${NO}_{2}$ column amount for a three-dimensional cloudy scene leads to relative errors up to 15%. These errors can be reduced when a differential correction spectrum due to three-dimensional effects is included in the retrieval.

The main conclusion of our analysis is that although for three-dimensional geometries the computational time is high, the main concepts of the algorithm are correct and the retrieval results are accurate.

For operational retrievals, we recommend

to construct a database for the spectral signal $I_{sim} (λ_{k}, X_{a})$ and its derivatives $\partial ln I_{sim} (λ_{k}, X_{a}) / \partial X_{g}$ , and to use these spectra in an one-step, three-dimensional retrieval algorithm, or
to construct a database for the differential correction spectra accounting on three-dimensional effects $S_{3 D} (λ_{k}, X_{a})$ , and to use these spectra in a full-step, one-dimensional retrieval algorithm.

The databases should be built up for different cloud scenarios, surface albedo, and solar and satellite geometries. The main problem which has to be solved is the design of a classification algorithm for broken clouds. One option is to use a nearest neighbor classifier based on the two-dimensional principal component analysis. Another option is to use a classification algorithm in which only the distributions of the cloud fields in the two-dimensional domains around the azimuth planes of the solar and viewing directions are of interest. The width of these domains corresponds to the instrument footprint. Note that these domains are also relevant in the tilted independent pixel approximation [34]: the direct solar beam is computed at all grid points which belong to the vertical columns inside the instrument footprint and for all characteristics within the entrance (solar direction) domain, while the measured radiance is computed by integrating the source function along all characteristics within the exit (viewing direction) domain. This will be the topic of a forthcoming paper.

Author Contributions

Conceptualization, A.D.; Software, A.D. and D.S.E.; formal analysis, D.S.E. and T.T.; writing—original draft preparation, A.D.; writing—review and editing, D.S.E. and T.T. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Appendix A

In this appendix we describe two inversion models for tropospheric column retrieval.

Assuming that the troposphere extends from level 1 to

N_{t}

, and the stratosphere from level

N_{t}

N_{z}

, the tropospheric and stratospheric total columns are given, respectively, by

X_{t g} = \sum_{j = 1}^{N_{t} - 1} x_{g, j}

(A1)

and

X_{s g} = \sum_{j = N_{t}}^{N_{z} - 1} x_{g, j} .

(A2)

Obviously, we have

X_{g} = X_{t g} + X_{s g}

. The retrieval of the tropospheric column of a gas

\bar{g}

(specifically,

{NO}_{2}

) is performed under the assumption that we have some a priori knowledge about the stratospheric column. In other words, we suppose that

X_{s \bar{g}}

can be approximated by

X_{s \bar{g}} \approx X_{s \bar{g}}^{☆}

, where

X_{s \bar{g}}^{☆}

is delivered for example, by the reference sector method. In this method, the stratospheric

{NO}_{2}

columns is estimated from measurements over the remote regions, under the assumptions of low longitudinal variations of stratospheric

{NO}_{2}

, and negligible tropospheric

{NO}_{2}

[35,36,37]. Denoting by

X_{- \bar{g}}

the set of all total columns excepting

X_{\bar{g}}

, i.e.,

X = {X_{\bar{g}}} \cup X_{- \bar{g}}

, a tropospheric column retrieval algorithm can be summarized as follows [17].

Solve the nonlinear Equation (8) of the DRMI inversion model for $x = {[X, b]}^{T}$ , or the nonlinear Equation (11) of the DRME inversion model for $x = {[X, b, c]}^{T}$ .
With $X_{s \bar{g}}^{☆}$ delivered by the reference sector method, and $X_{- \bar{g}}$ and $b$ determined at Step 1, solve the nonlinear equation

$\begin{matrix} R_{mes}^{δ} (λ_{k}) & = R_{sim} (λ_{k}, X_{t \bar{g}}, X_{s \bar{g}}^{☆}, X_{- \bar{g}}) + \sum_{j = 1}^{N_{s}} b_{j} S_{j} (λ_{k}, X_{a}), k = 1, \dots, N_{λ}, \end{matrix}$

(A3)

of the DRMI inversion model for $x = [X_{t \bar{g}}]$ , or the nonlinear equation

$\begin{matrix} R_{mes}^{δ} (λ_{k}) & = ln I_{sim} (λ_{k}, X_{t \bar{g}}, X_{s \bar{g}}^{☆}, X_{- \bar{g}}) + \sum_{j = 1}^{N_{s}} b_{j} S_{j} (λ_{k}, X_{a}) - P (λ_{k}, c), k = 1, \dots, N_{λ}, \end{matrix}$

(A4)

of the DRME inversion model for $x = {[X_{t \bar{g}}, c]}^{T}$ .

Appendix B

In this appendix we present a simplified derivation of the standard DOAS equation. For a more rigorous treatment, we refer to the work in [19].

We start with the linearized model (14) written as

\begin{matrix} ln I_{sim} (λ_{k}, X) & = \sum_{g = 1}^{N_{g}} X_{g} W_{g} (λ_{k}, X_{a}) + ln I_{sim} (λ_{k}, X_{a}) - \sum_{g = 1}^{N_{g}} X_{a g} W_{g} (λ_{k}, X_{a}) + P_{p} (λ_{k}, c_{p}), \end{matrix}

(A5)

where

W_{g} (λ_{k}, X_{a}) = \frac{\partial ln I_{sim}}{\partial X_{g}} (λ_{k}, X_{a}) .

(A6)

Further, we consider the first-order Taylor expansion

\begin{matrix} ln I_{sim} (λ_{k}, X) & = ln I_{sim} (λ_{k}, X_{a}) + \sum_{g = 1}^{N_{g}} (X_{g} - X_{a g}) \frac{\partial ln I_{sim}}{\partial X_{g}} (λ_{k}, X_{a}) + ε_{lin} (λ_{k}, X - X_{a}) \end{matrix}

(A7)

X = 0

, that is,

\begin{matrix} ln I_{sim} (λ_{k}, 0) & = ln I_{sim} (λ_{k}, X_{a}) - \sum_{g = 1}^{N_{g}} X_{a g} W_{g} (λ_{k}, X_{a}) + ε_{lin} (λ_{k}, - X_{a}), \end{matrix}

(A8)

where

I_{sim} (λ_{k}, 0)

corresponds to an atmosphere without gaseous absorption and

ε_{lin} (λ_{k}, - X_{a})

is the linearization error. Defining

ln {\bar{I}}_{sim} (λ_{k}, 0) = ln I_{sim} (λ_{k}, X_{a}) - \sum_{g = 1}^{N_{g}} X_{a g} W_{g} (λ_{k}, X_{a}),

(A9)

which yields

ln I_{sim} (λ_{k}, 0) = ln {\bar{I}}_{sim} (λ_{k}, 0) + ε_{lin} (λ_{k}, - X_{a}),

(A10)

we express Equation (A8) as

ln I_{sim} (λ_{k}, X) = \sum_{g = 1}^{N_{g}} X_{g} W_{g} (λ_{k}, X_{a}) + ln {\bar{I}}_{sim} (λ_{k}, 0) + P_{p} (λ_{k}, c_{p}) .

(A11)

In general,

ln I_{sim} (λ_{k}, 0)

can be approximated by a polynomial, but

ln {\bar{I}}_{sim} (λ_{k}, 0)

not. However, if

ε_{lin} (λ_{k}, - X_{a})

is negligible, then we have

ln I_{sim} (λ_{k}, 0) \approx ln {\bar{I}}_{sim} (λ_{k}, 0)

; thus,

ln {\bar{I}}_{sim} (λ_{k}, 0)

can be approximated by a polynomial, which can be included in

P_{p} (λ_{k}, c_{p})

. Consequently, we obtain

ln I_{sim} (λ_{k}, X) = \sum_{g = 1}^{N_{g}} X_{g} W_{g} (λ_{k}, X_{a}) + P_{p} (λ_{k}, c_{p}),

(A12)

so that by putting

P (λ_{k}, c) - P_{p} (λ_{k}, c_{p}) \to P (λ_{k}, c)

we are led to the following representation for the DRME model (11),

\begin{matrix} R_{mes}^{δ} (λ_{k}) & = \sum_{g = 1}^{N_{g}} X_{g} W_{g} (λ_{k}, X_{a}) + \sum_{j = 1}^{N_{s}} b_{j} S_{j} (λ_{k}, X_{a}) - P (λ_{k}, c) . \end{matrix}

(A13)

In view of the computational formula (cf. Equation (29))

\begin{matrix} W_{g} (λ_{k}, X_{a}) & = \frac{\partial ln I_{sim}}{\partial X_{g}} (λ_{k}, X_{a}) = \frac{1}{X_{a g}} \sum_{j = 1}^{N_{z}} C_{abs g} (λ_{k}, z_{j}) n_{a g, j} \frac{\partial ln I_{sim}}{\partial σ_{abs g, j}^{gas}} (λ_{k}, X_{a}) \end{matrix}

(A14)

we define

{\bar{C}}_{abs g} (λ_{k}) = \frac{\sum_{j = 1}^{N_{z}} C_{abs g} (λ_{k}, z_{j}) n_{a g, j} \frac{\partial ln I_{sim}}{\partial σ_{abs g, j}^{gas}} (λ_{k}, X_{a})}{\sum_{j = 1}^{N_{z}} n_{a g, j} \frac{\partial ln I_{sim}}{\partial σ_{abs g, j}^{gas}} (λ_{k}, X_{a})},

(A15)

and rewrite Equation (A13) as

\begin{matrix} R_{mes}^{δ} (λ_{k}) & = - \sum_{g = 1}^{N_{g}} A_{g} (λ_{k}, X_{a}) {\bar{C}}_{abs g} (λ_{k}) X_{g} + \sum_{j = 1}^{N_{s}} b_{j} S_{j} (λ_{k}, X_{a}) - P (λ_{k}, c), \end{matrix}

(A16)

where

A_{g} (λ_{k}, X_{a}) = - \frac{1}{{\bar{C}}_{abs g} (λ_{k})} W_{g} (λ_{k}, X_{a})

(A17)

is the air mass factor. Equation (A16) represents the standard DOAS equation. Note that if the computation of the partial derivatives in Equation (A15) is an expensive process, we may define

{\bar{C}}_{abs g} (λ_{k})

by the relation

{\bar{C}}_{abs g} (λ_{k}) = \frac{\sum_{j = 1}^{N_{z}} C_{abs g} (λ_{k}, z_{j}) n_{a g, j}}{\sum_{j = 1}^{N_{z}} n_{a g, j}} .

(A18)

References

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; Kokhanovsky, A., Ed.; Springer International Publishing: Cham, Switzerland, 2018; pp. 231–277. [Google Scholar] [CrossRef]
Belov, V.V.; Kirnos, I.V.; Tarasenkov, M.V. Estimation of the influence of cloudiness on the Earth observation from space through a gap in a cloudy field. In Proceedings of the 21st International Symposium Atmospheric and Ocean Optics: Atmospheric Physics, SPIE, Tomsk, Russia, 22–26 June 2015; Romanovskii, O.A., Ed.; [Google Scholar] [CrossRef]
Tarasenkov, M.V.; Kirnos, I.V.; Belov, V.V. Observation of the Earth’s surface from the space through a gap in a cloud field. Atmos. Ocean. Opt. 2017, 30, 39–43. [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]
Schwaerzel, M.; Emde, C.; Brunner, D.; Morales, R.; Wagner, T.; Berne, A.; Buchmann, B.; Kuhlmann, G. Three-dimensional radiative transfer effects on airborne and ground-based trace gas remote sensing. Atmos. Meas. Tech. 2020, 13, 4277–4293. [Google Scholar] [CrossRef]
Evans, K. The spherical harmonic discrete ordinate method for three-dimensional atmospheric radiative transfer. J. Atmos. Sci. 1998, 55, 429–446. [Google Scholar] [CrossRef]
Cracknell, A.; Hayes, L. Introduction to Remote Sensing, 2nd ed.; CRC Press: New York, NY, USA, 2007. [Google Scholar]
Chandrasekhar, S. Radiative Transfer. In Dover Books on Intermediate and Advanced Mathematics; Dover Publications: New York, New York, USA, 1960. [Google Scholar]
Sobolev, V.V. Light Scattering in Planetary Atmospheres; Pergamon Press: New York, USA, 1972. [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]
Bohren, C.; Huffman, D. Absorption and Scattering of Light by Small Particles; Wiley: New York, NY, USA, 1998. [Google Scholar] [CrossRef] [Green Version]
Kokhanovsky, A. Cloud Optics; Springer: Dordrecht, The Netherlands, 2006. [Google Scholar] [CrossRef]
Bodhaine, B.; Wood, N.; Dutton, E.; Slusser, J. On Rayleigh optical depth calculations. J. Atmos. Ocean. Technol. 1999, 16, 1854–1861. [Google Scholar] [CrossRef]
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]
Doicu, A.; Trautmann, T.; Schreier, F. Numerical Regularization for Atmospheric Inverse Problems; Springer: Berlin/Heidelberg, Germany, 2010. [Google Scholar] [CrossRef] [Green Version]
Efremenko, D.; 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]
Liu, S. Inversion Models for the Retrieval of Total and Tropospheric NO2 Columns. Atmosphere 2019, 10, 607. [Google Scholar] [CrossRef] [Green Version]
Platt, U.; Stutz, J. Differential Optical Absorption Spectroscopy; Springer: Berlin/Heidelberg, Germany, 2008. [Google Scholar] [CrossRef] [Green Version]
Rozanov, V.V.; Rozanov, A.V. Differential optical absorption spectroscopy (DOAS) and air mass factor concept for a multiply scattering vertically inhomogeneous medium: Theoretical consideration. Atmos. Meas. Tech. 2010, 3, 751–780. [Google Scholar] [CrossRef] [Green Version]
Bakushinskii, A. The problem of the convergence of the iteratively regularized Gauss–Newton method. Comput. Math. Math. Phys. 1992, 32, 1353–1359. [Google Scholar]
Xu, J.; Rao, L.; Schreier, F.; Efremenko, D.; Doicu, A.; Trautmann, T. Insight into Construction of Tikhonov-Type Regularization for Atmospheric Retrievals. Atmosphere 2020, 11, 1052. [Google Scholar] [CrossRef]
Tikhonov, A. On the solution of ill-posed problems and the method of regularization. Dokl. Akad. Nauk SSSR 1963, 151, 501–504. [Google Scholar]
Morozov, V. On the solution of functional equations by the method of regularization. Dokl. Akad. Nauk SSSR 1966, 167, 510–512. [Google Scholar]
Spherical Harmonic Discrete Ordinate Method (SHDOM) for Atmospheric Radiative Transfer. Available online: http://coloradolinux.com/shdom/ (accessed on 3 January 2020).
Doicu, A.; Efremenko, D. Linearizations of the Spherical Harmonic Discrete Ordinate Method (SHDOM). Atmosphere 2019, 10, 292. [Google Scholar] [CrossRef] [Green Version]
Levis, A.; Schechner, Y.; Aides, A.; Davis, A. Airborne Three-Dimensional Cloud Tomography. In Proceedings of the 2015 IEEE International Conference on Computer Vision (ICCV), Santiago, Chile, 7–13 December 2015. [Google Scholar] [CrossRef]
Doicu, A.; Efremenko, D.; Trautmann, T. A Spectral Acceleration Approach for the Spherical Harmonics Discrete Ordinate Method. Remote Sens. 2020, 12, 3703. [Google Scholar] [CrossRef]
Ambartzumyan, V. The effect of the absorption lines on the radiative equilibrium of the outer layers of the stars. Publ. Obs. Astron. Univ. Leningrad 1936, 6, 7–18. [Google Scholar]
Efremenko, D.; Doicu, A.; Loyola, D.; Trautmann, T. Optical property dimensionality reduction techniques for accelerated radiative transfer performance: Application to remote sensing total ozone retrievals. J. Quant. Spectrosc. Radiat. Transfer. 2014, 133, 128–135. [Google Scholar] [CrossRef]
Henyey, L.C.; Greenstein, J.L. Diffuse radiation in the Galaxy. Astrophys. J. 1941, 93, 70. [Google Scholar] [CrossRef]
Rozanov, V.; Vountas, M. Radiative transfer equation accounting for rotational Raman scattering and its solution by the discrete-ordinates method. J. Quant. Spectrosc. Radiat. Transf. 2014, 133, 603–618. [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] [2.0.CO;2" target='_blank'>Green Version]
Nakajima, T.; Tanaka, M. Algorithms for radiative intensity calculations in moderately thick atmos using a truncation approximation. J. Quant. Spectrosc. Radiat. Transfer. 1988, 40, 51–69. [Google Scholar] [CrossRef]
Várnai, T.; Davies, R. Effects of Cloud Heterogeneities on Shortwave Radiation: Comparison of Cloud-Top Variability and Internal Heterogeneity. J. Atmos. Sci. 1999, 56, 4206–4224. [Google Scholar] [CrossRef]
Richter, A.; Burrows, J. Tropospheric NO₂ from GOME measurements. Adv. Space Res. 2002, 29, 1673–1683. [Google Scholar] [CrossRef]
Bucsela, E.J.; Krotkov, N.A.; Celarier, E.A.; Lamsal, L.N.; Swartz, W.H.; Bhartia, P.K.; Boersma, K.F.; Veefkind, J.P.; Gleason, J.F.; Pickering, K.E. A new stratospheric and tropospheric NO₂ retrieval algorithm for nadir-viewing satellite instruments: Applications to OMI. Atmos. Meas. Tech. 2013, 6, 2607–2626. [Google Scholar] [CrossRef] [Green Version]
Beirle, S.; Hörmann, C.; Jöckel, P.; Liu, S.; de Vries, M.; Pozzer, A.; Sihler, H.; Valks, P.; Wagner, T. The STRatospheric Estimation Algorithm from Mainz (STREAM): Estimating stratospheric NO₂ from nadir-viewing satellites by weighted convolution. Atmos. Meas. Tech. 2016, 9, 2753–2779. [Google Scholar] [CrossRef] [Green Version]

Figure 1. Indicator function

f (x_{i}, y_{j})

with

x_{i} = i Δ x

and

y_{j} = j Δ y

for the 8 cloudy scenes (cases (a–h)); i and j refer to the pixel index in the xy plane, while each pixel corresponds to a 500 m × 500 m area.

Figure 1. Indicator function

f (x_{i}, y_{j})

with

x_{i} = i Δ x

and

y_{j} = j Δ y

for the 8 cloudy scenes (cases (a–h)); i and j refer to the pixel index in the xy plane, while each pixel corresponds to a 500 m × 500 m area.

Figure 2. Slices of the indicator function

f (x_{i})

with

x_{i} = i Δ x

y = 7

km for the cloudy scenes shown in Figure 1 (a–d,f,h).

Figure 2. Slices of the indicator function

f (x_{i})

with

x_{i} = i Δ x

y = 7

km for the cloudy scenes shown in Figure 1 (a–d,f,h).

Figure 3. A priori vertical profile of

{NO}_{2}

volume mixing ratio.

Figure 3. A priori vertical profile of

{NO}_{2}

volume mixing ratio.

Figure 4. Ring correction spectrum included in the retrieval.

Figure 5. Residuals versus the iteration index for Differential Radiance Model with Internal closure (DRMI) (left) and Differential Radiance Model with External closure (DRME) (right).

Figure 6. Correction spectra due to three-dimensional effects.

Table 1. Relative errors and the computation times (CPU) in minutes: seconds for DRMI. The results correspond to a two-dimensional geometry.

Cloudy Scene	DRMI		DRMI Approx. Derivatives
Cloudy Scene	Rel. Errors	CPU	Rel. Errors	CPU
1	$1.66 \times 10^{- 5}$	39:51	$1.18 \times 10^{- 1}$	32:24
2	$1.10 \times 10^{- 5}$	40:12	$1.50 \times 10^{- 1}$	32:48
3	$9.22 \times 10^{- 6}$	39:14	$9.79 \times 10^{- 2}$	31:06
5	$1.08 \times 10^{- 5}$	38:54	$1.37 \times 10^{- 1}$	30:44
6	$1.03 \times 10^{- 5}$	39:32	$8.77 \times 10^{- 2}$	31:16
8	$1.94 \times 10^{- 5}$	40:06	$1.51 \times 10^{- 1}$	32:33

Table 2. Relative errors and the computation times (CPU) in minutes:seconds for DRME. The numbers in parenthesis are the values of the regularization parameter

α

. The results correspond to a two-dimensional geometry.

Table 2. Relative errors and the computation times (CPU) in minutes:seconds for DRME. The numbers in parenthesis are the values of the regularization parameter

α

. The results correspond to a two-dimensional geometry.

Cloudy Scene	DRME		DRME (One Iteration)
Cloudy Scene	Rel. Errors	CPU	Rel. Errors	CPU
1	$6.94 \times 10^{- 5}$	35:10	$2.84 \times 10^{- 2} (3.0 \times 10^{- 6})$	5:33
2	$3.23 \times 10^{- 5}$	36:47	$3.74 \times 10^{- 3} (1.0 \times 10^{- 6})$	5:54
3	$1.23 \times 10^{- 5}$	34:57	$7.74 \times 10^{- 3} (1.0 \times 10^{- 6})$	5:12
5	$4.84 \times 10^{- 5}$	34:23	$6.73 \times 10^{- 3} (1.0 \times 10^{- 6})$	5:01
6	$1.08 \times 10^{- 5}$	35:06	$2.60 \times 10^{- 2} (5.0 \times 10^{- 6})$	5:28
8	$2.31 \times 10^{- 5}$	36:23	$1.87 \times 10^{- 2} (2.0 \times 10^{- 5})$	5:48

Table 3. The same as in Table 2 but for a three-dimensional geometry.

Cloudy Scene	DRME		DRME (One Iteration)
Cloudy Scene	Rel. Errors	CPU	Rel. Errors	CPU
1	$4.74 \times 10^{- 4}$	14:15:31	$4.05 \times 10^{- 2} (4.0 \times 10^{- 5})$	2:34:03
2	$3.04 \times 10^{- 4}$	14:21:43	$1.04 \times 10^{- 2} (1.0 \times 10^{- 6})$	2:36:42
3	$1.35 \times 10^{- 4}$	14:15:57	$3.01 \times 10^{- 2} (2.0 \times 10^{- 6})$	2:32:26
5	$3.63 \times 10^{- 4}$	14:15:23	$2.38 \times 10^{- 2} (1.0 \times 10^{- 6})$	2:32:11
6	$1.46 \times 10^{- 4}$	14:14:36	$3.56 \times 10^{- 2} (5.0 \times 10^{- 6})$	2:33:38
8	$1.88 \times 10^{- 4}$	14:20:13	$3.18 \times 10^{- 2} (4.0 \times 10^{- 5})$	2:34:16

Table 4. Relative errors and the computation times (CPU) in minutes:seconds for DRMI. The measured and simulated spectral signals are computed by a three- and an one-dimensional radiative transfer model, respectively.

Cloudy Scene	DRMI		DRMI (Correction. Spectrum)
Cloudy Scene	Rel. Errors	CPU	Rel. Errors	CPU
1	$1.21 \times 10^{- 1}$	2:21	$7.59 \times 10^{- 2}$	2:42
2	$5.18 \times 10^{- 2}$	2:32	$2.76 \times 10^{- 2}$	2:58
3	$7.03 \times 10^{- 2}$	2:14	$3.74 \times 10^{- 2}$	2:36
4	$1.51 \times 10^{- 1}$	2:04	$8.40 \times 10^{- 2}$	2:24
5	$6.08 \times 10^{- 2}$	2:27	$3.76 \times 10^{- 2}$	2:47
6	$9.82 \times 10^{- 2}$	2:25	$5.50 \times 10^{- 2}$	2:45
7	$5.41 \times 10^{- 2}$	2:25	$3.26 \times 10^{- 2}$	2:46
8	$7.15 \times 10^{- 2}$	2:30	$4.16 \times 10^{- 2}$	2:50

Table 5. The same as in Table 4 but for DRME.

Cloudy Scene	DRME		DRME (Correction. Spectrum)
Cloudy Scene	Rel. Errors	CPU	Rel. Errors	CPU
1	$1.21 \times 10^{- 1}$	2:10	$8.60 \times 10^{- 2}$	2:16
2	$5.19 \times 10^{- 2}$	2:12	$2.70 \times 10^{- 2}$	2:18
3	$6.99 \times 10^{- 2}$	2:05	$3.70 \times 10^{- 2}$	2:11
4	$1.53 \times 10^{- 1}$	1:52	$8.45 \times 10^{- 2}$	2:01
5	$6.12 \times 10^{- 2}$	2:16	$3.82 \times 10^{- 2}$	2:22
6	$8.82 \times 10^{- 2}$	2:14	$5.53 \times 10^{- 2}$	2:18
7	$5.37 \times 10^{- 2}$	2:13	$3.21 \times 10^{- 2}$	2:19
8	$7.11 \times 10^{- 2}$	2:21	$4.03 \times 10^{- 2}$	2:17

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2021 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.; Trautmann, T. A Proof-of-Concept Algorithm for the Retrieval of Total Column Amount of Trace Gases in a Multi-Dimensional Atmosphere. Remote Sens. 2021, 13, 270. https://doi.org/10.3390/rs13020270

AMA Style

Doicu A, Efremenko DS, Trautmann T. A Proof-of-Concept Algorithm for the Retrieval of Total Column Amount of Trace Gases in a Multi-Dimensional Atmosphere. Remote Sensing. 2021; 13(2):270. https://doi.org/10.3390/rs13020270

Chicago/Turabian Style

Doicu, Adrian, Dmitry S. Efremenko, and Thomas Trautmann. 2021. "A Proof-of-Concept Algorithm for the Retrieval of Total Column Amount of Trace Gases in a Multi-Dimensional Atmosphere" Remote Sensing 13, no. 2: 270. https://doi.org/10.3390/rs13020270

APA Style

Doicu, A., Efremenko, D. S., & Trautmann, T. (2021). A Proof-of-Concept Algorithm for the Retrieval of Total Column Amount of Trace Gases in a Multi-Dimensional Atmosphere. Remote Sensing, 13(2), 270. https://doi.org/10.3390/rs13020270

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

Article Menu

A Proof-of-Concept Algorithm for the Retrieval of Total Column Amount of Trace Gases in a Multi-Dimensional Atmosphere

Abstract

1. Introduction

2. Retrieval Algorithm

2.1. Inverse Models

2.2. Solution Method

2.3. Linearized Radiative Transfer Model

3. Numerical Simulations

3.1. Test Example 1

3.2. Test Example 2

4. Discussion

Author Contributions

Funding

Conflicts of Interest

Appendix A

Appendix B

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI