# Estimates of remote sensing retrieval errors by the GRASP algorithm: application to ground-based observations, concept and validation

Milagros Herrera, Oleg Dubovik, Benjamin Torres, Tatyana Lapyonok, David

Fuertes, Anton Lopatin, Pavel Litvinov, Cheng Chen, Jose Antonio

Benavent-Oltra, Juan Bali, et al.

► To cite this version:

Milagros Herrera, Oleg Dubovik, Benjamin Torres, Tatyana Lapyonok, David Fuertes, et al. Esti- mates of remote sensing retrieval errors by the GRASP algorithm: application to ground-based obser- vations, concept and validation. Atmospheric Measurement Techniques, 2022, 15 (20), pp.6075-6126. 10.5194/amt-15-6075-2022 . hal-03830279

HAL Id: hal-03830279

https://hal.science/hal-03830279

Submitted on 28 Oct 2022

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

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

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

# Estimates of remote sensing retrieval errors by the GRASP algorithm: application to ground-based observations, concept and validation

Milagros E. Herrera<fnref n="1" />, Oleg Dubovik<fnref n="1" />, Benjamin Torres<fnref n="1" />, Tatyana Lapyonok<fnref n="1" />, David Fuertes<fnref n="2" />, Anton Lopatin<fnref n="2" />, Pavel Litvinov<fnref n="2" />, Cheng Chen<fnref n="2" />, Jose Antonio Benavent-Oltra<fnref n="3" />, Juan L. Bali<fnref n="4" />, and Pablo R. Ristori<fnref n="5" />

**Correspondence:** Milagros E. Herrera (milagros.herrera@univ-lille.fr)

Received: 22 March 2022 – Discussion started: 8 April 2022

Revised: 17 August 2022 – Accepted: 9 September 2022 – Published: 21 October 2022

**Abstract.** Understanding the uncertainties in the retrieval of aerosol and surface properties is very important for an adequate characterization of the processes that occur in the atmosphere. However, the reliable characterization of the error budget of the retrieval products is a very challenging aspect that currently remains not fully resolved in most remote sensing approaches. The level of uncertainties for the majority of the remote sensing products relies mostly on post-processing validations and intercomparisons with other data, while the dynamic errors are rarely provided. Therefore, implementations of fundamental approaches for generating dynamic retrieval errors and the evaluation of their practical efficiency remains of high importance.

This study describes and analyses the dynamic estimates of uncertainties in aerosol-retrieved properties by the GRASP (Generalized Retrieval of Atmosphere and Surface Properties) algorithm. The GRASP inversion algorithm, described by Dubovik et al. (2011, 2014, 2021), is designed based on the concept of statistical optimization and provides dynamic error estimates for all retrieved aerosol and surface properties. The approach takes into account the effect of both random and systematic uncertainties propagations. The algorithm provides error estimates both for directly retrieved parameters included in the retrieval state vector and for the characteristics derived from these parameters. For ex-

ample, in the case of the aerosol properties, GRASP directly retrieves the size distribution and the refractive index that are used afterwards to provide phase function, scattering, extinction, single scattering albedo, etc. Moreover, the GRASP algorithm provides full covariance matrices, i.e. not only variances of the retrieval errors but also correlations coefficients of these errors. The analysis of the correlation matrix structure can be very useful for identifying less than obvious retrieval tendencies. This appears to be a useful approach for optimizing observation schemes and retrieval set-ups.

In this study, we analyse the efficiency of the GRASP error estimation approach for applications to ground-based observations by a sun/sky photometer and lidar. Specifically, diverse aspects of the error generations and their evaluations are discussed and illustrated. The studies rely on a series of comprehensive sensitivity tests when simulated sun/sky photometer measurements and lidar data are perturbed by random and systematic errors and inverted. Then, the results of the retrievals and their error estimations are analysed and evaluated. The tests are conducted for different observations of diverse aerosol types, including biomass burning, urban, dust and their mixtures. The study considers observations of AErosol RObotic NETwork (AERONET) sun/sky photometer measurements at 440, 675, 870 and 1020 nm and multiwavelength elastic lidar measurements at 355, 532 and

<footnote n="1">¹Laboratoire d'Optique Atmosphérique, CNRS/Université de Lille, 59655 Villeneuve d'Ascq, France</footnote>

<footnote n="2">²GRASP SAS, Remote Sensing Developments, 59260 Lezennes, France</footnote>

<footnote n="3">³Department of Electrical, Electronic, Automatic Control Engineering and Applied Physics, Escuela Técnica Superior de Ingeniería y Diseño Industrial, Universidad Politécnica de Madrid, 28012 Madrid, Spain</footnote>

<footnote n="4">⁴National Scientific and Technical Research Council (CONICET), C1425FQB CABA, Buenos Aires, Argentina</footnote>

<footnote n="5">⁵CEILAP-UNIDEF (MINDEF – CONICET) – CITEDEF, B1603ALO Buenos Aires, Argentina</footnote>
1064 nm. The sun/sky photometer data are inverted alone or
together with lidar data.

The analysis shows overall successful retrievals and er-
ror estimations for different aerosol characteristics, includ-
ing aerosol size distribution, complex refractive index, sin-
gle scattering albedo, lidar ratios, aerosol vertical profiles,
etc. Also, the main observed tendencies in the error dynamic
agree with known retrieval experience. For example, the
main accuracy limitations for retrievals of all aerosol types
relate to the situations with low optical depth. Also, in sit-
uations with multicomponent aerosol mixtures, the reliable
characterization of each component is possible only in lim-
ited situations, for example, from radiometric data obtained
for low solar zenith angle observations or from a combina-
tion of radiometric and lidar data. At the same time, the total
optical properties of aerosol mixtures are always retrieved
satisfactorily.

In addition, the study includes an analysis of the detailed
structure of the correlation matrices for the retrieval errors in
mono- and multicomponent aerosols. The conducted analysis
of error correlation appears to be a useful approach for opti-
mizing observation schemes and retrieval set-ups. The appli-
cation of the approach to real data is provided.

# 1 Introduction

Remote sensing is one major tool for monitoring atmosphere
and surface properties at large scales. These observations
have a nondestructive character and allow for dynamic lo-
cal, regional or global monitoring of the ambient atmosphere.
Correspondingly, diverse remote sensing observations are
employed for routine observations and characterization of
the Earth atmosphere. One of the key challenges in imple-
menting remote sensing is the development of the retrieval al-
gorithms. While remote sensing retrievals have substantially
evolved during the last few decades, a significant need for
further advancing various aspects of the retrieval algorithms
remains. One of the most challenging and important, while
underdeveloped, aspects is the evaluation of the errors in the
retrieval products. For example, the review by Sayer et al.
(2020) emphasizes that, for most aerosol satellite retrievals,
the quality of the retrieval uncertainty estimates has not been
routinely assessed.

Here we discuss and analyse the approach implemented in
the GRASP (Generalized Retrieval of Atmosphere and Sur-
face Properties) retrieval algorithm. The GRASP concept is
based on the statistical optimization fitting designed for the
retrieval of detailed aerosol properties from diverse obser-
vations (Dubovik et al., 2011, 2014, 2021). This algorithm
uses statistical estimates of random error propagation and
provides the dynamic error estimates for both retrieved pa-
rameters (such as the size distribution, refractive index, etc.)
and characteristics derived from those parameters (total scat-

tering, extinction, single scattering albedo, etc.). Specifically,
in this work, we discuss and analyse the error estimates of
the aerosol properties by GRASP for aerosol retrievals from
ground-based observations by sun/sky radiometers and li-
dars.

One of the most visible data sets of ground-based ra-
diometric observations is provided by AERONET (AErosol
RObotic NETwork; Holben et al., 1998), a network of
more than 500 operational sites distributed over the world.
AERONET provides column-integrated aerosol properties of
the different sites distributed over the world. AERONET pro-
vides column-integrated aerosol properties originally pro-
vided by the Nakajima et al. (1996) algorithm and later by
the Dubovik and King (2000) retrieval. Several studies access
the AERONET retrieval errors. First, Dubovik et al. (2000)
provided a rather comprehensive analysis of retrieval uncer-
tainties caused by both random measurement errors and sys-
tematic errors originating from potential biases in the mea-
surements and imperfections in the modelling aerosol prop-
erties. This analysis was revisited by Torres et al. (2017),
whose studies overall confirmed most of the uncertainty ten-
dencies revealed by Dubovik et al. (2000). Recently, Sinyuk
et al. (2020) published a concept for uncertainty in aerosol
retrievals that have been adapted in the version 3 aerosol op-
erational product of AERONET (Giles et al., 2019). In the
frame of this concept, the uncertainties are estimated using
the spread of retrieved parameters generated by 27 distinct
combinations of retrieval implemented with perturbed input
data (aerosol optical depth (AOD), sky radiances, solar spec-
tral irradiance and surface reflectance). A somewhat similar
concept for the estimating error was employed earlier in the
LiRIC (Lidar and Radiometer Inversion Code) approach for
the synergetic processing of co-located lidar and AERONET
sun/sky photometer observations (Chaikovsky et al., 2016).
LiRIC provided some uncertainty, which was simulated us-
ing a series of retrievals with the perturbed input data. A
large number of factors affect the retrievals, whose variation
is complex and nonlinear. Indeed, modelling all the factors
and circumstances that can affect the retrieval in all situa-
tions is theoretically impossible, and practically challenging,
within the limited series of perturbed runs, especially for sit-
uations when a large number of parameters is retrieved. In
these situations, the error propagation approaches based on
statistical estimation theory and described in numerous text-
books (e.g. Edie et al., 1971; Fourgeaud and Fuchs, 1967;
Rodgers, 2000) provide asymptotically comprehensive esti-
mates for random retrieval errors. At the same time, it should
be noted that both the result of perturbation tests and sta-
tistical estimates of propagated error relying on the forward
model employed may not fully represent the inaccuracies re-
lated to the limitations of this model. Some additional evalu-
ations and considerations are always desirable for accessing
the adequacy of the chosen forward model and its potential
limitations.
The GRASP is a state-of-the-art, highly versatile inversion algorithm of a new generation that can be applied to a variety of remote sensing and laboratory observations. For example, GRASP has been applied for several satellite instruments (Dubovik et al., 2019, 2021; Li et al., 2019; Chen et al., 2018, 2019, 2020; Puthukkudy et al., 2020) and demonstrated the capability to provide accurate information about detailed properties of aerosol and underlying surface reflectance. The GRASP has been successfully used for retrieving aerosol properties from polar nephelometer data (Espinosa et al., 2017, 2019; Schuster et al., 2019). The applications of the GRASP algorithm for retrieving aerosol properties from sun/sky photometer and radiometer observations, either alone or in combination with lidar observations, have been employed and discussed in numerous studies (e.g. Torres et al., 2017) as the synergy of sun/sky photometer plus lidar (Lopatin et al., 2013, 2021; Torres et al., 2017; Benavent-Oltra et al., 2017, 2019, 2021; Hu et al., 2019; Tsekeri et al., 2017; Román et al., 2018; Herreras et al., 2019; Titos et al., 2019). The successful application of GRASP has been illustrated for the interpretation of terrestrial observations with a celestial camera (Román et al., 2017). Recent studies by Torres et al. (2017) and Torres and Fuertes (2021) demonstrated the high potential of the GRASP retrieval concept for inverting only direct sun photometric observations.

The evaluation of the retrieval accuracy in all those studies was realized by comparing and validating the retrieval results with independent reliable data. It should be noted, however, that practically none of these studies discuss the retrieval error estimates, while the formalism of the error estimate has been realized in GRASP for a while. As a result, the validity of the retrieval error estimates provided by the GRASP approach remains unattended and unverified. Therefore, in order to address this issue, the current work proposes a discussion of the main aspects of the GRASP error generation and attempts to provide an evaluation of the retrieval error estimates provided by GRASP. The study is focused on the considerations of aerosol retrieval from sun/sky photometers and lidar ground-based observations.

## 2 Modelling of error estimates in the GRASP algorithm

As mentioned in the introduction, in this work we make use of the GRASP (Generalized Retrieval of Atmosphere and Surface Properties) algorithm. It is a rigorous, versatile and open-source algorithm capable of providing information of aerosol properties from the measurements of different instruments and dynamic error estimates (Dubovik et al., 2011, 2014, 2021). It is a flexible, generalized algorithm that relies on two independent modules, i.e. the forward model and the numerical inversion. The forward model contains the full description of the physical model, including various interactions of electromagnetic solar radiations, such as

aerosol scattering, surface reflectance and gaseous absorption. The multiple scattering interactions in the atmosphere are accounted for by solving the vector of radiative transfer equation. Thus, the GRASP forward model is capable of simulating diverse measurements in the laboratory and atmospheric remote sensing, including passive and active observations from space and the ground. On the other hand, numerical inversion is not directly related to any physical problem and realizes the formal inversion of the measurements using a statistical estimation approach. Specifically, GRASP employs the multiterm least squares method (LSM) that allows for a flexible utilization of multiple a priori constraints. This approach is very convenient for designing diverse remote sensing retrievals, as discussed in detail by Dubovik et al. (2021).

The retrieval error estimates in GRASP are calculated by modelling the propagation of measurement errors based on a statistical estimation approach. In addition, the formulation used for estimating errors accounts for some contribution of the systematic errors that could originate from biases in the measurement or some modifications implemented in the algorithm for improving the retrieval convergence of nonlinear solutions. Below, the descriptions of the overall concept and specific key implementations of the error estimation in GRASP are provided.

### 2.1 The numerical inversion based on a statistical optimization concept

The multiterm LSM employed in GRASP searches for the solution using statistically optimized fitting under multiple a priori constraints (Dubovik, 2004; Dubovik et al., 2011, 2021). It considers both measurements and a priori data in a similar manner by considering them as being data from different and independent data sources, as follows:

$$f_k^* = f_k(a) + \Delta f_k^*, \quad (1)$$

where $f_k(a)$ is the physical forward model, $a$ is the vector of unknown parameters, $\Delta f_k^*$ represents the uncertainty associated with the measurement $f_k^*$ and $k$ denotes different data sets that are not correlated and may have different levels of uncertainties described by different covariance matrices $C_k$. Such an explicit differentiation of the input data makes the retrieval more transparent because it clearly identifies the different data sets used. Correspondingly, the joint probability density function (PDF) of independent data sets $f_1^*, f_2^*, \dots, f_K^*$ can be obtained by the simple multiplication of the PDFs of data from all $K$ sources, as follows:

$$ \begin{aligned} P(f(a)|f^*) &= P(f_1(a), \dots, f_K(a)|f_1^*, \dots, f_K^*) \\ &= \prod_{k=1}^{K} P(f_k(a)|f_k^*). \end{aligned} \quad (2) $$

It can be noted that Eq. (1) does not assume any relation between the forward models $f_k(a)$, i.e. forward mod-
els $f_k(a)$ can be the same or different. In the frame of the LSM approach, i.e. under the assumptions of a normal PDF of the error $\Delta f_k^*$, the solution of the Eq. (1) corresponds to the minimum of the following functional:

$$2\hat{\Psi}(a) = \frac{1}{2} \sum_{k=1}^{K} (f_k(a) - f_k^*)^T C_k^{-1} (f_k(a) - f_k^*) = \min. \quad (3)$$

It can be noted that, if any of correlation coefficients are 1 or -1, the covariance matrix is degenerated with det($C_k$) = 0.

For the general case of nonlinear functions $f_k(a)$, the solution of Eq. (3) is sought iteratively, as follows:

$$a^{p+1} = a^p - \Delta a^p, \quad (4)$$

where $\Delta a^p$ is the solution that can be found by solving the system of so-called normal equations. This is as follows:

$$\left(\sum_{k=1}^{K} \mathbf{K}_{k,p}^{T} \mathbf{C}_{k}^{-1} \mathbf{K}_{k,p}\right) \Delta \mathbf{a}^{p} = \sum_{k=1}^{K} \mathbf{K}_{k,p}^{T} \mathbf{C}_{k}^{-1} \Delta \mathbf{f}_{k}^{p}, \quad (5)$$

where $\Delta f^p = f(a^p) - f^*$, and $\mathbf{K}_p$ is a Jacobean matrix at the pth iteration of the functions $f_k(a)$ in the vicinity of $a^p$ with the elements $\{\mathbf{K}_{k,p}\}_{j,i} = \frac{\partial f_{k,j}(a)}{\partial a_i}\bigg|_{a=a^p}$.

The asymptotic limit of the minimized quadratic form, for most applications, can be written as follows:

$$2\Psi(a) = \min \rightarrow \sum_{k=1}^{K} N_k - n, \quad (6)$$

where $N_k$ is a number of measurements (inputs) in the kth data set, and n is the number of retrieved unknowns (parameters).

It should be noted that the LSM solution of Eq.(3) corresponds to the minimum of the quadratic form $\hat{\Psi}(a)$ and does not depend on the value of this minimum. Considering this fact, in a practical application it is often convenient to renormalize the minimized quadratic $\hat{\Psi}(a)$, and in situations when only one data set is inverted, it is convenient to use a weighting matrix $\mathbf{W} = \mathbf{C}/\varepsilon_1^2$ and minimize the quadratic form $2\Psi'(a) = 2\varepsilon_1^2\Psi(a)$. In such an approach, one does not need to know the exact value of the variance $\varepsilon_1^2$. Moreover, $\varepsilon_1^2$ can be estimated from asymptotic LSM expectations provided by Eq. (9), as shown below.

In the frame of a multiterm approach, the use of weighting matrices additionally allows for making the contribution of different data sources more explicit. Indeed, using the weighting matrices $\mathbf{W}_k$ instead of covariance matrices $\mathbf{C}_k$ in Eq. (5) can be written as follows:

$$\left(\sum_{k=1}^{K} \gamma_k \mathbf{K}_k^T \mathbf{W}_k^{-1} \mathbf{K}_k\right) \Delta \mathbf{a}^p = \sum_{k=1}^{K} \gamma_k \mathbf{K}_k^T \mathbf{W}_k^{-1} \Delta \mathbf{f}_k^p. \quad (7)$$

In this formulation, the relative contributions of the data from different data sources are scaled by the corresponding Lagrange parameters, $\gamma_i$, defined as follows:

$$\mathbf{W}_i = \frac{1}{\varepsilon_i^2} \mathbf{C}_i \quad \text{and} \quad \gamma_i = \frac{\varepsilon_i^2}{\varepsilon_i^2}, \quad (8)$$

where $\varepsilon_i^2$ is the first diagonal element of $\mathbf{C}_i$, i.e. $\varepsilon_i^2 = \{\mathbf{C}_i\}_{11}$, and $\gamma_i$ is the ratio of the variances of scattered radiances and variances of the corresponding data set. Evidently, $\gamma_1 = 1$, as discussed by Dubovik and King (2000); Dubovik (2004); Dubovik et al. (2011) and others. This renormalization strategy is especially convenient for a multiterm LSM approach once some of data sets correspond to a priori information. In addition, the renormalized definition of the minimized quadratic function (or residual) is $\Psi'(a) = \varepsilon_1^2\Psi(a)$, and the measurement error variance $\varepsilon_1^2$ can be estimated from the residual of the fit. Indeed, once the weighting matrices are used in the solution, and taking into account that the values of minimized quadratic form of $\Psi(\hat{a})$ follow a $\chi^2$ distribution, Eq. (7) minimizes the quadratic form with the limit depending on $\varepsilon_1^2$, as follows:

$$2\Psi'(a) = 2\varepsilon_1^2 \Psi(a) = \min \rightarrow \varepsilon_1^2 \left(\sum_{k=1}^{K} N_k - n\right) \\ \text{and} \quad \hat{\varepsilon}_1^2 \approx \frac{2\Psi'(a^p)}{\sum_{k=1,...,K} N_k - n}. \quad (9)$$

Therefore, $\varepsilon_1^2$ can be estimated from the minimum value of the residual $\Psi(\hat{a})$.

### A priori constraints in a multiterm LSM approach and in GRASP algorithm

As discussed in detail by Dubovik et al. (2021), the multiterm LSM concept has been proposed as being a methodologically convenient approach for integrating different types of a priori constraints in remote sensing applications (Dubovik, 2004; Dubovik and King, 2000; Dubovik et al., 1995, 2000, 2008, 2011). In the multiterm LSM, a priori estimates are considered to be equivalent to the measurements, i.e. characterized by their PDF and treated equivalently to the actual measurements. In this regard, Eqs. (1)–(7) do not show any distinction between different $f_k(a)$.

At the same time, in practice, there are always two different types of data sets, i.e. measurements and the a priori constraint on the unknowns $a$. Therefore, the vector of the measurement $(f^*)^T = (f_1^*, f_2^*, \dots, f_K^*)^T$ can be written as follows:

$$(f^*)^T = (f_1^*, f_2^*, \dots, f_K^*, f_1^a, f_2^a, \dots, f_N^a)^T, \quad (10)$$

where $f_i^* = f_i^*(a)$ represent K directly measured characteristics, and $f_i^a = f_i^*(a)$ represent N a priori known characteristics of unknowns $a$. Correspondingly, the right-hand side of
Eq. (2) can be formally split into the following two groups:

$$P(f(a)|f^*) = \prod_{k=1}^{K} P(f_k(a)|f_k^*) \prod_{n=1}^{N} P(f_n^a(a)|f_n^{a,*}). \quad (11)$$

Therefore, Eq. (7) can also be formally arranged to identify the contribution of measurements and a priori terms, as follows:

$$\left( \sum_{k=1}^{K} \gamma_k \mathbf{K}_{k,p}^T \mathbf{W}_k^{-1} \mathbf{K}_{k,p} + \sum_{n=1}^{N} \gamma_{a,n} \mathbf{K}_{a,n,p}^T \mathbf{W}_{a,n}^{-1} \mathbf{K}_{a,n,p} \right) \Delta \mathbf{a}^p = \sum_{k=1}^{K} \gamma_k \mathbf{K}_{k,p}^T \mathbf{W}_k^{-1} \Delta \mathbf{f}_k^{a,*} + \sum_{n=1}^{N} \gamma_{a,n} \mathbf{K}_{a,n,p}^T \mathbf{W}_{a,n}^{-1} \Delta \mathbf{f}_n^{a,*}, \quad (12)$$

where the two groups of the terms on left and right parts of the equation represent the contributions of the set of K measured characteristics $f_k(a)$ and the set of N a priori $f_n^a(a)$ characteristics, and the Lagrangian parameters are defined as follows:

$$\gamma_k = \frac{\varepsilon_{k=1}^2}{\varepsilon_k^2} \quad \text{and} \quad \gamma_{a,n} = \frac{\varepsilon_{n=1}^2}{\varepsilon_{a,n}^2}. \quad (13)$$

As discussed by Dubovik (2004) and Dubovik et al. (2021), the multiterm approach is a simple rearranging of the base LSM formulation, while the resulting Eq. (7) provides a solid basis for unifying many known formulas of the constrained inversion in a single formalism and is practical, convenient and efficient for developing remote sensing algorithms using diverse complementary observations and a priori constraints.

While the multiterm LSM concept allows flexible utilizations of nearly arbitrary a priori constraints, the GRASP algorithm is fully adapted for using the most popular and physically transparent a priori constraints, such as direct a priori estimates of unknowns a and smoothness constraints in situations when the unknown vector a or any group of unknowns included in this vector represent continuous smooth function. For example, if vector a represents an aerosol size distribution that is known to be rather smooth, the system given by Eq. (1) can be explicitly written as follows:

$$\left\{
\begin{aligned}
f_{k=1}^{*} &= f_{k=1}^{*}(a) + \Delta f_{k=1}^{*} \\
f_{n=1}^{a,*} &= f_{n=1}^{a,*}(a) + \Delta f_{n=1}^{a,*} \\
f_{n=2}^{a,*} &= f_{n=2}^{a,*}(a) + \Delta f_{n=2}^{a,*}
\end{aligned}
\right.
=
\left\{
\begin{aligned}
f_{1}^{*} &= f_{k=1}^{*}(a) + \Delta f_{1}^{*} \\
f_{1}^{a,*} &= f_{1}^{a,*}(a) + \Delta f_{1}^{a,*} \\
f_{2}^{a,*} &= f_{2}^{a,*}(a) + \Delta f_{2}^{a,*}
\end{aligned}
\right.
= \left\{
\begin{aligned}
f_{1}^{*} &= f_{k=1}^{*}(a) + \Delta f_{1}^{*} \\
a^{*} &= a + \Delta a^{*} \\
0^{*} &= G_{m}^{a} + \Delta g^{*}.
\end{aligned}
\right.
\quad (14)$$

The a priori constraints defined by the second line $a^* = a + \Delta a^*$ represent the most common constraints of the solution by direct a priori estimates of unknowns $a^*$, where $\Delta a^*$ are the uncertainties in the estimates $a^*$ and are generally considered to be unbiased random errors within the covariance matrix $C_{a^*}$. These constraints can be easily included in Eq. (12) by defining the $K_a = 1$ unity matrix, i.e. $K_a^T W_a^{-1} K_a = W_a^{-1}$ and $K_a^T W_a^{-1} f_1^{a,*} = W_a^{-1} a$. The utilization of a priori estimates $a^*$ was introduced in the pioneering studies by Twomey (1963) and later evolved and discussed in detail in the Rodgers (2000) textbook on inversion. The third term represents another common type of a priori constraint known as the smoothness constraints that limit the variability in retrieved functions by using a priori knowledge about limitations on derivatives of those functions. For example, a priori knowledge limits high-frequency variations in continuous functions $v(x)$, such as the aerosol size distribution. In GRASP, the smoothness constraints are related to a priori known limited values of the derivatives, i.e. with their mth derivative deviations from zero, as follows:

$$\frac{\partial^m v(x)}{\partial x^m} \approx 0. \quad (15)$$

For the vector of unknowns $a = (a_1, a_2, \dots, a_n)^T$ that contain discrete elements describing the continuous function $v(x)$, the knowledge on the smoothness of function $v(x)$ can be defined using a vector-matrix linear system (e.g. see Dubovik et al., 2021). $0^* = G_m^a + \Delta g^*$, where $G_m$ is the Jacobean matrix of the matrix of the mth derivatives $G_m^a$. In practice, these are often approximated by matrices of the mth finite difference estimated in point $a$. The errors $\Delta g^*$ reflect the uncertainty in the knowledge of the deviations of $y(x)$ from the assumed constant ($m=1$), straight line ($m=2$), parabola ($m=3$) and so on. Under the assumption that the $\Delta g^*$ have a normal distribution, with mean zero and covariance matrix $C_g$, these constraints can be easily included in Eq. 12 by defining $K_{a,2} = G_m$ and $f_2^{a,*} = 0^*$, i.e. $K_a^T W_a^{-1} K_a = G_m^T W_{\Delta g}^{-1} G_m^T$ and $K_2^T W_2^{-1} \Delta f_2^{a,*} = G_m^T W_{\Delta g}^{-1} (a^p - 0^*) = G_m^T W_{\Delta g}^{-1} a^p$. Utilization of such smoothness constraints was suggested by one of the first formulations of constrained inversion by Phillips (1962) and was also considered in an article by Tikhonov (1963) and his later studies.

Thus, for a case where only direct a priori estimates and smoothness constraints are used, Eq. (12) can be explicitly written via weighting matrices as follows:

$$\left( \mathbf{K}^T \mathbf{W}_f^{-1} \mathbf{K} + \gamma_a \mathbf{W}_a^{-1} + \gamma_g \mathbf{\Omega}_m \right) \Delta \mathbf{a}^p = \mathbf{K}^T \mathbf{W}_f^{-1} \Delta \mathbf{f}^p + \gamma_a \mathbf{W}_a^{-1} (\mathbf{a}^p - \mathbf{a}^*) + \gamma_g \mathbf{\Omega}_m \mathbf{a}^p, \quad (16)$$

where $\mathbf{\Omega}_m$ denotes the smoothness matrix. This is defined as follows:

$$\mathbf{G}_m^T \mathbf{W}_{\Delta g}^{-1} \mathbf{G}_m^T = \mathbf{\Omega}_m, \quad (17)$$

and the explicit formulation of $\mathbf{\Omega}_m$ can be found in the paper by Dubovik et al. (2011). Equation (14) generalizes
the commonly used base equations of constrained inversion by Phillips (1962), Twomey (1975, 1977), Tikhonov (1963) and Rodgers (1976, 1990, 2000). It should be noted that Eq. (14) is written for the simplest situation, when the vector **a** represents only one continuous function $v(x)$, while, in many GRASP applications, the vector of unknowns includes several components $\mathbf{a}^T = (\mathbf{a}_{\text{sd}}^T, \mathbf{a}_{n(\lambda)}^T, \mathbf{a}_{k(\lambda)}^T, \mathbf{a}_h^T, \dots)^T$, where each component is relevant to continuous functions representing such physical characteristics as the aerosol particle size distribution ($\mathbf{a}_{\text{sd}}$), spectral dependence of real ($\mathbf{a}_{n(\lambda)}$) and complex ($\mathbf{a}_{k(\lambda)}$) parts of refractive index, vertical distribution ($\mathbf{a}_h$), etc. Each of those characteristics is continuous function and therefore, in retrieval, the smoothness constraints can be applied on each of corresponding component of the vector of unknowns. Evidently, direct a priori constraints can be applied to each single element of the vector **a**, while, from a practical viewpoint, separating and outlining the contribution of a priori estimates for each component, e.g. $(\mathbf{a}^*)^T = ((\mathbf{a}_{\text{sd}}^*)^T, (\mathbf{a}_{n(\lambda)}^*)^T, (\mathbf{a}_{k(\lambda)}^*)^T, (\mathbf{a}_h^*)^T, \dots)^T$. Similarly, the inverted measurements may come from different sources and therefore have different levels of accuracy and different weighting matrices. As a result, in practice, all the first, second and third terms in Eq. (14) may have many quite different components, and therefore, the actual formulation of the solution can be significantly more complex. Some of the explicit equations can be found in the paper by Dubovik et al. (2011).

The realization of the inversion in GRASP, in principle, is based on the general Eq. (7), while for convenience there is a logical separation, as indicated in Eq. (12), into actual measurements and a priori constraints. For each measurement data set $f_k^*$, two types of errors can be set, i.e. the relative or absolute, and the magnitude of the errors is defined by the standard deviation and a weighting matrix $\mathbf{W}_i$. The standard deviation is used inside the code to calculate the corresponding Lagrange parameters $\gamma_i$. The weighting matrix $\mathbf{W}_i$ is assumed as being the unity matrix by default, while it can also be set as diagonal, with different values at the diagonal, and in a more general way with non-zero non-diagonal values too. For applying a priori constraints, as discussed above, there are two main possibilities, namely using direct a priori constraints or applying smoothness constraints for the parameters that define continuous functions.

The direct a priori estimates $a_i^*$ for each of value $a_i$ in the vector of unknowns $\mathbf{a} = (a_1, a_2, \dots, a_n)^T$ can be provided with the corresponding Lagrangian parameters $\gamma_{a_i}$. There is also a possibility of assuming a vector $\mathbf{a}^*$ of a priori estimates for all the retrieved parameters or for selected groups (e.g. parameters describing size distribution) with the common Lagrangian parameter $\gamma_a$. In this case, the weighting matrix $\mathbf{W}_a$ is also provided that is assumed to be the unity matrix by default, or it can be set as diagonal, with different values at the diagonal, or in more general way with non-zero non-diagonal values.

The smoothness a priori constraints can be applied for each group of parameters describing a continuous function (e.g. $\mathbf{a}_{\text{sd}}^T, \mathbf{a}_{n(\lambda)}^T, \mathbf{a}_{k(\lambda)}^T, \mathbf{a}_h^T, \dots$) by defining the order $m$ of limited derivatives ($m=0$ is a constant; $m=1$ is a straight line; $m=2$ is a parabola, etc.), and the strength of the applied a priori smoothness constraints is defined by Lagrangian parameters $\gamma_n$. The smoothness matrix $\Omega_m$ is defined as in Eq. (15), where the weighting matrix $\mathbf{W}_{\Delta g}$ is the unity matrix by default and can be set as diagonal with different values on the diagonal in case the retrieved continuous function has a different level of variability for different ordinates.

It should be noted that GRASP considers two types of a priori constraints, namely the single-pixel a priori constraints for the retrieved parameters that correspond to simultaneous and co-located observations, and the multipixel constraints that limit the variability for unknowns in different groups of similar parameters when several such groups of unknowns are retrieved simultaneously from coordinated but not fully coincident or not fully co-located observations (see details in Dubovik et al., 2011, 2021). In the current paper, only single-pixel constraints are used.

## 2.2 Nonlinear inversion in GRASP using the Levenberg-Marquardt optimization

Since most of atmospheric remote sensing applications are strongly nonlinear, the Levenberg-Marquardt optimization (Press et al., 1992; Ortega and Rheinboldt, 1970) is realized to optimize the convergence of GRASP solutions. Specifically, as described by Dubovik et al. (2021), in GRASP it is assumed that the correction of the solution at the $p$th iteration $\Delta \mathbf{a}^p$ should be limited, especially at the initial iterations when the linearization error is the largest. For such cases, in GRASP, for the determination of $\Delta \mathbf{a}^p$ in the iterative procedure, an additional constraint on the correction $\Delta \mathbf{a}^p$ is added at each iteration, as follows:

$$ \Delta \mathbf{a}^{p,*} = \mathbf{0}^* + \Delta \mathbf{a}. \qquad (18) $$

Correspondingly, using this additional requirement, an additional term will be introduced in Eq. (7), as follows:

$$ \left( \sum_{k=1}^{K} \gamma_k \mathbf{K}_k^T (\mathbf{W}_k)^{-1} \mathbf{K}_k + \mathbf{D}_{\Delta a}^p \right) \Delta \mathbf{a}^p \\ = \sum_{k=1}^{K} \gamma_k \mathbf{K}_k^T (\mathbf{W}_k)^{-1} \Delta f_k^p, \qquad (19) $$

where the matrix $\mathbf{D}_{\Delta a}$ is a diagonal matrix. This has the following elements:

$$ \mathbf{D}_{\Delta a ii} = \gamma \Delta a_i = \frac{\varepsilon_1^2}{\varepsilon_{\Delta a_i}^2}. \qquad (20) $$

The variance $\varepsilon_{\Delta a_i}$ can be determined, for example, assuming that whole known range of each parameter $a_i$ variability should be covered by $3\varepsilon_{\Delta a_i}$, i.e. $a_{i,max} - a_{i,min} \approx 3\varepsilon_{\Delta a_i}$.
Also, following the common Levenberg-Marquardt procedure, the impact of the correction $\Delta a^p$ is always scaled by a factor $t_p$ in Eq. (4), as follows:

$$a^{p+1} = a^p - t_p \Delta a^p, \quad (21)$$

where $t_p$ is in the range $0 < t_p \le 1$. It is selected empirically to provide convergence, by decreasing $t_p = t_p/2$, until the decrease in the residual $\Psi'(a^p) \le \Psi'(a^{p-1})$ is achieved (see Dubovik et al., 2011).

Thus, in the case of nonlinear $f_k(a)$ and/or $f_i^a(a)$, the inversion in GRASP includes Levenberg-Marquardt-like optimizations and is implemented in the frame of Eqs. (4) and (5). While this optimization certainly helps to achieve a successful convergence of the solution in practice, it should also be considered to be one of possible sources of uncertainties, as pointed out by Dubovik et al. (2021), and will be discussed below.

## 2.3 Error propagation estimates in GRASP

Estimations of the retrieval errors in GRASP are based on LSM equations expressed for the case of multiterm solutions written via weighting matrixes (Dubovik et al., 2021). Both the contribution of random and systematic error components are estimated as follows:

$$C_{\hat{a}} = C_{\Delta\hat{a}_{\text{(ran)}}} + (\hat{a}_{\text{bias}})(\hat{a}_{\text{bias}})^T, \quad (22)$$

where, in the following,

$$C_{\Delta\hat{a}_{\text{(ran)}}} = \langle \Delta\hat{a}_{\text{(ran)}} \Delta\hat{a}_{\text{(ran)}}^T \rangle \\ \approx \left( \sum_{k=1}^{K} \gamma_k \mathbf{K}_k^T \mathbf{W}_k^{-1} \mathbf{K}_k \right)^{-1} \hat{\varepsilon}_1^2 \quad (23)$$

$$\hat{a}_{\text{bias}} = \left( \sum_{k=1}^{K} \gamma_k (\mathbf{K}_k^T \mathbf{W}_k^{-1} \mathbf{K}_k) \right)^{-1} \\ \left( \sum_{k=1}^{K} \gamma_k (\mathbf{K}_k^T \mathbf{W}_k^{-1} \mathbf{b}_k^*) \right), \quad (24)$$

where $\mathbf{b}_k^*$ denotes the bias vector in the $k$th data set $f_k$, and $\hat{\varepsilon}_1^2$ is estimated from the resulting misfit of the data using Eq. (9).

The estimation of not only the random retrieval error but also the error retrieval bias $a_{\text{bias}}$ is important for the adequate evaluation of retrieval uncertainty, especially in the case when multiple a priori constraints are used. For example, for the case of the retrieval given by Eq. (14), $C_{\Delta\hat{a}_{\text{(ran)}}}$ is expressed as follows:

$$C_{\Delta\hat{a}_{\text{(ran)}}} \approx (\mathbf{K}^T \mathbf{W}^{-1} \mathbf{K} + \gamma_a \mathbf{W}_a^{-1} + \gamma_g \mathbf{\Omega}_m)^{-1} \hat{\varepsilon}_1^2. \quad (25)$$

A rather obvious tendency can be seen from the analysis of this equation because the higher the contributions of the

second and the third terms, the smaller the random errors are, i.e. the stronger a priori constraints used, the lower the random errors will be in the retrieval. However, in practice, a priori constraints can be unintentionally inadequate and therefore introduce some systematic uncertainties, i.e. biases. In principle, there is no guaranteed approach for detecting those biases unless a comprehensive analysis and validation of the retrievals have been done. Nonetheless, some biases can manifest themselves via misfit of measurements $\Delta f_k^{\text{bias}} = f_k(a^{\text{solution}}) - f_k^*$ or the misfit of a priori constraints. For example, for Eq. (14), the bias can be introduced by a priori estimate $a_{\text{bias}}^* = a^{\text{solution}} - a^*$ or unsmoothed features in the retrieved solution, i.e. $a_{\text{bias}}^{\text{smooth}} = \mathbf{\Omega}_m a^{\text{solution}} \neq \mathbf{0}$. Correspondingly, the bias for single-pixel retrieval is estimated as follows:

$$\hat{a}_{\text{bias}} \approx (\mathbf{K}^T \mathbf{W}^{-1} \mathbf{K} + \gamma_a \mathbf{W}_a^{-1} + \gamma_g \mathbf{\Omega}_m)^{-1} (\mathbf{K}^T \mathbf{W}^{-1} \Delta f_{\text{bias}} + \gamma_a \mathbf{W}_a^{-1} a_{\text{bias}}^* + \gamma_g \mathbf{\Omega}_m a_{\text{bias}}^{\text{smooth}}). \quad (26)$$

In this equation, the contribution of a priori estimates to bias is probably the most significant in many applications, since it is never possible to have fully accurate a priori values (widely used in optimum estimation approaches) for constraining. In a similar way, the a priori biases are estimated in the case when multipixel a priori constraints are used.

The Levenberg-Marquardt optimization of the convergence, discussed in Sect. 2.2, may also introduce bias. Indeed, this optimization makes the iterations converge from a given initial guess to fit the data, even if the basic linear system is singular. Therefore, once the Levenberg-Marquardt optimization is used, there is an evident dependence on the initial guess that can bias the solution. In order to take this into account, Eqs. (23) and (24) are modified as follows:

$$C_{\Delta\hat{a}_{\text{(ran)}}} \approx \left( \sum_{k=1}^{K} \gamma_k \mathbf{K}_k^T \mathbf{W}_k^{-1} \mathbf{K}_k + \mathbf{D}_{\Delta a}^p \right)^{-1} \hat{\varepsilon}_1^2, \quad (27)$$

and

$$\hat{a}_{\text{bias}} = \left( \sum_{k=1}^{K} \gamma_k \mathbf{K}_k^T \mathbf{W}_k^{-1} \mathbf{K}_k + \mathbf{D}_{\Delta a}^p \right)^{-1} \left( \sum_{k=1}^{K} \gamma_k \mathbf{K}_k^T \mathbf{W}_k^{-1} \mathbf{b}_k^* + \mathbf{D}_{\Delta a}^p (a^{\text{solution}} - a^{p=0}) \right). \quad (28)$$

Note that, from fundamental viewpoint, the a priori information is used in order to make solution unique, and if it is fully and adequately added, no dependence on the initial guess should be observed. At the same time, in practice, such a dependence often appears to some extent, especially in cases when the state vector includes a large number of unknowns. Moreover, if the retrieval is not optimally set, such a dependence can be rather significant, while unnoticed, because the retrieval continues to converge to local minima
once the Levenberg-Marquardt optimization of the retrieval convergence is used. Therefore, in order to account for such an effect, the Levenberg-Marquardt contribution was added into the formalism for accounting possible biases. According to our evaluation, this term is nearly negligible if there is no dependence on the initial guess, while it increases if such a dependence appears.

Equations (27 and 28) allow one to obtain the error estimates for the retrieved parameters. That is, for example, when the configuration is from sun photometer and lidar measurements, then the expected retrieved parameters are $dV/d\ln(r)$, and the real and imaginary part of refractive index, sphericity fraction and aerosol volume concentration are vertically distributed.

Also, in practice, the users may not directly need the retrieved parameters $\hat{a}$ but rather their functions $m(\hat{a})$ that can be calculated from the retrieved parameters. For example, GRASP retrieves the parameters of aerosol microphysics (particle sizes, refractive indices, etc.), but users need aerosol optical depth (AOD). For such a situation, GRASP provides a set of such diverse indirect characteristics with the possibilities of providing the uncertainties that are calculated as follows:

$$ 
\begin{aligned}
C_{\Delta\hat{m}} & \approx M(C_{\Delta\hat{a}_{\text{ran}}} + \hat{a}_{\text{bias}}\hat{a}_{\text{bias}}^T)M^T \\
& = MC_{\Delta\hat{a}_{\text{ran}}}M^T + M\hat{a}_{\text{bias}}(M\hat{a}_{\text{bias}})^T \\
& = C_{\Delta\hat{m}_{\text{ran}}} + \hat{m}_{\text{bias}}\hat{m}_{\text{bias}}^T,
\end{aligned} 
\quad (29) $$

where $M$ is the matrix of first derivatives $M_{ji} = \frac{\partial m_j}{\partial a_i}$, $a^{\text{solution}}$.

Finally, the effect of biases in the measurements on the solution bias $\hat{a}_{\text{bias}}$ is accounted for in Eq. (26), based on the assumption that the presence of biases is manifested in the non-zero misfits $\Delta f_k^{\text{bias}}$. Indeed, this is true in many cases when systematic errors are present in the inverted measurements or the accurate fit of inverted data cannot be achieved (e.g. see the illustrations provided by numerical sensitivity tests for AERONET retrievals by Dubovik et al., 2000). At the same time, there are many situations in which the biases in the measurements may not significantly affect the residual (Eq. 9) and the misfits $\Delta f_k^{\text{bias}}$. For example, the retrievals of aerosol single scattering albedo (SSA) from AERONET ground-based measurements are highly sensitive to the calibration biases in the direct sun measurements, while the fitting of these direct measurements is always quite accurate (see the discussion by Dubovik et al., 2000). The effects of such measurement biases can be estimated by implementing proxy numerical tests applied to the measurements perturbed by possible biases. For example, the recent approach for evaluation retrieval errors in AERONET operational products is estimated using a series of ~27 numerical proxy inversion tests with the sets of perturbations in both the input measurements and auxiliary input parameters (Sinyuk et al., 2020). A similar strategy can be used for the evaluation of the potential effects of undetected biases. Specifically, the bias term $(\hat{a}_{\text{bias}})(\hat{a}_{\text{bias}})^T$ in Eq. (22) can be estimated as follows:

$$ (\hat{a}_{\text{bias}})(\hat{a}_{\text{bias}})^T \rightarrow \langle (\hat{a}_{\text{bias}})(\hat{a}_{\text{bias}})^T \rangle_{\text{bias proxy set}} \quad (30) $$

where the values of the retrieval biases are estimated as being an average effect from a preselected set of possible biases in measurements and auxiliary inputs. Therefore, if we assume a positive and negative bias in the equation for the systematic component, then the contribution to Eq. (26) can be written as follows:

$$ \hat{a}_{\text{bias}} \approx (\mathbf{K}^T \mathbf{W}^{-1} \mathbf{K} + \gamma_a \mathbf{W}_a^{-1} + \gamma_g \mathbf{\Omega}_m)^{-1} (\mathbf{K}^T \mathbf{W}^{-1} \mathbf{b}_{f_{\text{bias}}} + \gamma_a \mathbf{W}_a^{-1} \mathbf{a}_{\text{bias}}^* + \gamma_g \mathbf{\Omega}_m \mathbf{a}_{\text{bias}}^{\text{smooth}}), \quad (31) $$

where the vectors $\mathbf{b}_{f_{\text{bias}}}$ represent the new bias related to the measurement.

In addition, in this work we also study the structure of the covariance matrix for different aerosols and configurations. Apparently, such a matrix provides interesting information about the error estimates (focusing on the diagonal elements) and the relation between the retrieval parameters (from the covariance values, i.e. non-diagonal elements). The representation of the covariance matrix for the parameters has the following structure:

$$ \text{Cov}(\mathbf{a}) = \begin{pmatrix} \sigma_1^2 & \sigma_1\sigma_2\rho_{12} & \sigma_1\sigma_3\rho_{13} & \cdots \\ \sigma_2\sigma_1\rho_{21} & \sigma_2^2 & \sigma_2\sigma_3\rho_{23} & \cdots \\ \sigma_3\sigma_1\rho_{31} & \sigma_3\sigma_2\rho_{32} & \sigma_3^2 & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{pmatrix}, \quad (32) $$

where, in the diagonal, one finds the variance of each element, and the non-diagonal elements represent the covariance of each retrieved element $a_i$ with the others. The variances, i.e. diagonal elements, are always used for estimating retrieval errors and providing the error bars. The non-diagonal elements are rarely considered, while they provide the very interesting and less-than-obvious information about error correlations.

In order to study the error correlation structure of the error, the following correlation matrix will be considered in this work that can be obtained from the covariance matrix (Eq. 32), as follows:

$$ \text{Corr}(\mathbf{a}) = \begin{pmatrix} 1 & \rho_{12} & \rho_{13} & \cdots \\ \rho_{21} & 1 & \rho_{23} & \cdots \\ \rho_{31} & \rho_{32} & 1 & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{pmatrix}, \quad (33) $$

where each diagonal element corresponds to the correlation with itself which is equal to 1, and the non-diagonal elements are the correlations related to each parameter that can vary between -1 and 1.
## 3 Methodology of error analysis

The calculation of retrieval estimates in GRASP is based on rigorous formulations of statistical estimations described above. At the same time, a practical evaluation of the developed error formalism and possible tuning is desirable for a comprehensive evaluation of the approach and gaining full confidence in the practical efficiency of the approach. In this regard, one can probably state that the error estimate always tends to be less accurate than the retrievals themselves. Indeed, in remote sensing, the retrieval relies on the formalism of the electromagnetic light interaction theory that is fundamentally very accurate and well established, while the factors contributing to the uncertainties can be very diverse, not fully formalized and often not even fully understood. For example, the forward model is nonlinear, while the error propagations are usually (and in this work specifically) estimated in linear approximations, as commented previously, and the retrieval can be affected by not fully predicted biases in the measurements or by impercipient of aerosols or surface models (in our applications). Therefore, an important part of establishing error estimates is their evaluation and validation.

As mentioned above, in this study, we attempt to evaluate the GRASP estimations based on an extensive series of the numerical tests with added random noise that covers a wide range of practical situations. Moreover, we complement this study, assuming different biases, in order to see how the error estimates are represented in the cases with both random noise and bias. This section describes the design of the numerical experiment, including the following:

- the instruments and retrieval scenarios used,

- the description of the overall experiment,

- the assumed atmospheric properties and covered distinct specific situations of interest, and

- the assumptions made for generated random errors, the considered retrieved parameter, the considered error characteristics and so on.

As mentioned earlier, this study evaluates the GRASP error estimates produced for the aerosol properties retrieved from ground-based observations. The details of the used observations and considered aerosol retrieval scenarios are provided in the next sections.

### 3.1 Observations and aerosol retrieval approaches considered

The analysis is focused on two widely known, and probably the most popular, retrieval scenarios used for deriving detailed aerosol optical properties.

i. Retrieval of columnar properties of aerosol from the measurements by ground-based sun/sky-scanning radiometers alone.

ii. Simultaneous retrieval of both columnar aerosol properties and their vertical distribution from the combined observations by sun/sky-scanning radiometers and multiwavelength lidar.

#### 3.1.1 Aerosol retrieval from sun/sky radiometers alone

All the tests and analyses in this study include the spectral observations by the ground-based sun/sky-scanning radiometers. These radiometers were used for more than 30 years by the worldwide AERONET project (https://aeronet.gsfc.nasa.gov, last access: 17 August 2022; Holben et al., 1998) that unities a federation of a large number of ground-based remote sensing aerosol international networks. At present, AERONET observations are widely recognized as being a benchmark validation data set for satellite aerosol retrieval and as source of some unique information about detailed aerosol properties that are used in diverse studies on monitoring and predicting regional and global pollution evolution and climate change. The standard set of AERONET observations includes spectral direct sun measurements and spectral measurements of angular sky radiance obtained from sky scans by the radiometer. The direct sun observations provide, with rather straightforward processing, spectral aerosol optical depth (AOD) that is itself highly valuable for satellite product validation and diverse aerosol studies (e.g. Eck et al., 1999). The combination of spectral AOD and sky radiances at four wavelengths, 440, 675, 870 and 1020 nm (see Table 1), are used for the retrieval of detailed aerosol size distributions (the aerosol concentration in 22 logarithmically equidistant-sized bins in the range from 0.05 to 15 µm) together with spectral dependence complex index of refraction (Dubovik and King, 2000). The retrieval also provides aerosol absorption characterized by single scattering albedo (SSA), the parameters of fine- and coarse-mode size distribution and other diverse detailed properties of columnar aerosol (e.g. see Dubovik et al., 2002b). In addition, based on the concept developed by Dubovik et al. (2002a); Dubovik et al. (2006), AERONET retrieval considers aerosol as being a mixture of two components, spherical and nonspherical, and provides a fraction of spherical particles as an additional parameter. The nonspherical fraction is modelled as a mixture of randomly oriented spheroids using fixed-axis ratio distribution, equal to the one retrieved by Dubovik et al. (2006), by inverting the full-phase matrices of the K-feldspar dust sample measured in the laboratory by Volten et al. (2001). A rather complete description of this retrieval concept is also provided in the paper by Dubovik et al. (2011). The set of the aerosol parameters retrieved in the AERONET standard operational processing is shown in Table 1, where the number of the retrieved parameters is 31 (22 size bins, with four values for the real refractive index (RRI) and imaginary refractive index (IRI) at each wavelength and one value related to the sphericity fraction). In addition, the parameters obtained for the case of mixed aerosol properties are also shown, in
which the properties for each mode, fine and coarse, are re-
trieved. For this last case, the number of the retrieved pa-
rameters is 43 (25 size bins, four values for RRI and IRI for
each mode, i.e. fine and coarse, and the sphericity fraction).
The AERONET operational retrievals were mainly provided
for solar almucantar geometry, while the AERONET has re-
cently also provide retrievals for more complex hybrid ob-
servational geometries (Sinyuk et al., 2020). These studies
are focused on retrievals in the solar almucantar only, while
the GRASP algorithm provides the error calculations for any
geometry. Thus, in this study, the direct sun measurements
and sky radiances at four different wavelengths, 440, 675,
870 and 1020 nm, for both are used in the inversion tests.
These sky radiances measurements are measured in the solar
almucantar (fixed-view zenith angle equal to the solar zenith
angle, SZA) with a varying azimuth angle ranging from ±3.5
to ±180° (Table 1).

The detailed aerosol properties in the total atmospheric
column provided by the AERONET inversion of sun/sky-
scanning radiometers has been widely recognized as being
rather unique, reliable data. For example, AERONET re-
trievals provided the first reliable data about aerosol spectral
absorption and other detailed aerosol optical characteristics
(e.g. see Dubovik et al., 2002b; Giles et al., 2012, and oth-
ers). These detailed data are of vital importance for evaluat-
ing the impact of aerosol on such important aspects as a cli-
mate change and diverse pollution effects and can be reliable
as they are estimated nearly uniquely from remote sensing
observations (Kaufman et al., 2002). Therefore, this retrieved
aerosol information has been proven to be very useful for the
assessment of climate change dynamics in the Intergovern-
mental Panel on Climate Change (IPCC) reports (Boucher et
al., 2013; IPCC, 2021) and other high-profile analyses.

As mentioned in earlier sections, the evaluations of the ac-
curacy of retrieved aerosol parameters mainly relied on ex-
tensive sensitivity studies by Dubovik et al. (2000). The re-
sults were used for providing quality assurance criteria and
expected accuracy estimation (see Dubovik et al., 2002b;
Holben et al., 2006). Sinyuk et al. (2020) recently pre-
sented the approach to estimate retrieval uncertainties used
in AERONET version 3 data. The approach estimates the er-
ror using the variability in retrieved parameters generated by
27 perturbations in both input measurements and auxiliary
input parameters. In comparison with these previous efforts,
this study evaluates the dynamic error estimates generated
by the GRASP approach for each retrieval based on the mea-
surement error propagation and bias estimations. This study
estimates the complete covariance matrices of retrieval er-
rors. In addition, the analysis of retrieval errors and their cor-
relations conducted here also aims at demonstrating the value
of the obtained estimates for understanding the retrieval error
tendencies and optimizing the retrieval approaches.

3.1.2 Aerosol retrieval from a combination of sun/sky radiometers and lidar data

The inversion of co-located observations by sun/sky ra-
diometers and lidar is another popular retrieval approach in
the aerosol community. Indeed, radiometer direct sun and
multiangular polarimetric observations of diffuse Sun radi-
ation transmitted through the atmosphere have a significant
sensitivity to the atmospheric aerosol amount, its particles
size, shape and morphology; however, they have practically
no sensitivity to the vertical variability in aerosols. The li-
dar observations, on the other hand, provide the information
about the vertical distribution of aerosol, while their sensi-
tivity to other aerosol properties is more limited compared
to radiometer observations. Therefore, the information from
co-located photometric measurements and lidar systems is
complementary and always desirable for the enhanced char-
acterization of aerosol properties. This complementarity is
well recognized by the research community, and a large num-
ber of joint observational sites with both radiometer and li-
dar observations have been established in last decade. In
these regards, the European ACTRIS (Aerosols, Clouds and
Trace gases Research Infrastructure Network) infrastructure
(https://www.actris.eu, last access: 12 August 2022) is one
of the good examples of networks emphasizing the acquisi-
tion of diverse complementary observations at each site. All
ACTRIS observational supersites possess both sun/sky radio-
metric and complex multiwavelength lidar systems.

GRASP retrieval has been successfully adapted by Lopatin
et al. (2013) for processing such combined observations in an
algorithm initially known as Generalized Aerosol Retrieval
from Radiometer and Lidar Combined data (GARRLiC).
This algorithm has been used in numerous studies (e.g.
Granados-Muñoz et al., 2014; Granados-Muñoz et al., 2016;
Tsekeri et al., 2017; Benavent-Oltra et al., 2017, 2019, 2021,
and others) and was also adapted for the operational process-
ing of lidar/radiometer observations in the frame of ACTRIS
infrastructure. Later, the capabilities of GARRLiC/GRASP
were significantly extended by Lopatin et al. (2021) for pro-
cessing diverse, vertically resolved observations alone in a
diverse combination with radiometric observations.

In these studies, we consider the aerosol retrieval from
the base GARRLiC/GRASP input data set that includes
AERONET sun/sky-scanning observations in the solar almuc-
cantar at four wavelengths and the lidar backscattering atten-
uation profile at three wavelengths at 355, 532 and 1064 nm
(see Table 1). Thus, for the synergy of sun/sky radiometers
and lidar measurements, we considered the same sun/sky
radiometer input data combined with the correlative range-
corrected signal (RCS) values at 355, 532 and 1064 nm. The
lidar signal provided in GRASP as input data is normalized
at 60 log-spaced bins at different heights, as in Lopatin et al.
(2013, 2021), giving minimum and maximum heights. It is
because all lidars provide observations within a certain dis-
tance range, which varies from instrument to instrument, and
<table><caption>**Table 1.** Summary of general input data and the set of parameters retrieved by the GRASP algorithm used in this work for two configurations, i.e. sun/sky radiometer only and sun/sky radiometer plus lidar.</caption><thead><tr><th>Sun/sky radiometer only</th><th>Sun/sky radiometer plus lidar</th></tr></thead><tbody><tr><td>Input</td><td>Input</td></tr><tr><td>- AOD*</td><td>- Sun/sky radiometer data</td></tr><tr><td>- Calibrated radiances*</td><td>- AOD*</td></tr><tr><td></td><td>- Calibrated radiances*</td></tr><tr><td></td><td>- Lidar data</td></tr><tr><td></td><td>- Range-corrected profiles (RCS**)</td></tr><tr><td>* At 440, 675, 870 and 1020 nm</td><td>normalized at 60 log-spaced bins at different heights<br>** At 355, 532 and 1064 nm</td></tr><tr><td>Retrieved aerosol properties in the total atmospheric column**</td><td>Retrieved aerosol properties (column integrated and vertical distributions**)</td></tr><tr><td>- dV(r<sub>i</sub>)/dlnr<sub>i</sub></td><td>- dV(r<sub>i</sub>)/dlnr<sub>i</sub> (in total atmospheric column)</td></tr><tr><td>- C<sub>sph</sub></td><td>- C<sub>sph</sub> (in total atmospheric column)</td></tr><tr><td>- n(λ<sub>i</sub>)</td><td>- n(λ<sub>i</sub>) (in total atmospheric column)</td></tr><tr><td>- k(λ<sub>i</sub>)</td><td>- k(λ<sub>i</sub>) (in total atmospheric column)</td></tr><tr><td></td><td>- C<sub>V</sub>(h) (vertical distribution)</td></tr><tr><td>Retrievals provided for total, fine and coarse modes.</td><td>Retrievals provided for total, fine and coarse modes.</td></tr></tbody></table>

it is limited by emitter/receiver field of view overlap in the lower part and by the signal-to-noise ratio in the upper part. The GRASP inversion of these data derives, in addition to columnar aerosol properties provided from radiometer only inversion, the vertical profile of aerosol concentration. Moreover, the aerosol can be considered to be an external mixture of two aerosol components (fine and coarse). In such a case, all retrieved parameters are provided for both aerosol modes, as shown in the Table 1, where the number of the retrieved parameters in this case is 174 (120 values of the aerosol vertical concentration for fine and coarse modes, 25 size bins, seven values of RRI and IRI at each wavelength and for each mode, fine and coarse, and the sphericity fraction).

## 3.2 Analysis of the structure of different error parameters

As was already mentioned, GRASP has a capability to provide the full covariance matrix of the retrieval errors, and this study aimed to evaluate and illustrate the efficiency of these estimated covariance matrices. At the same time, retrieval error evaluations in most of the practical applications rely on the consideration of mainly diagonal elements of the covariance matrices, while non-diagonal elements of covariance matrices are much less common. Indeed, in spite of the fact that non-diagonal elements of covariance matrices provide valuable and interesting information about retrieval error correlations, these non-diagonal elements are not often available in practice, and the analysis of error correlations requires more sophisticated considerations compared to straightforward analysis diagonal elements only, and therefore, it is less

popular. Considering these aspects, in the present study, as a first step, we make a more detailed and extensive analysis of the error variances, and then, as a second step, we illustrate the usefulness of obtained non-diagonal elements.

The performance of the GRASP error variances estimated, provided by Eqs. (25)–(27), is studied using a series of numerical tests. Figure 1 illustrates the general scheme of the organization of these tests.

First, as shown in Fig. 1, the parameters $a_{assumed}$ for the assumed detailed aerosol properties ($dV(r_i)/d\ln r_i$, $n(\lambda_i)$, $k(\lambda_i)$, $C_{sph}$ and $C_V(h)$ in the case of lidar) are used to obtain the synthetic observations using the GRASP forward model. These synthetic observations include the spectral AOD, sky radiances and range-corrected signals (RCSs) of lidar. These data are used then in the inversion tests where the aerosol parameters and their errors are estimated from these synthetic observations using the GRASP algorithm. In order to study the effects of the different uncertainties, both random and systematic errors are added to the synthetic measurements before the inversion, and then the retrieved parameters $\hat{a}_{retr}$ are compared with $a_{assumed}$ (see Fig. 1). Therefore, from the retrieved parameters, the retrieval errors provided by the GRASP algorithm and actual retrieval errors can be compared. Thus, these actual errors are calculated comparing $a_{assumed}$ and $\hat{a}_{retr}$ as follows:

$$ 
\begin{cases}
\Delta \hat{a}_{\text{abs}} = \hat{a}_{\text{retr}} - a_{\text{assumed}} \\
\Delta \hat{a}_{\text{rel}} = \frac{\Delta \hat{a}_{\text{abs}}}{a_{\text{assumed}}} \cdot 100\%,
\end{cases}
\quad (34) $$

<footnote>* Azimuth angles, for sky radiances in the almucantar geometry, relative to the Sun (in °): 3.0, 3.5, 4.0, 5.0, 6.0, 7.0, 8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 110.0, 120.0, 140.0, 160.0 and 180.</footnote>

<footnote>** Number of retrieved parameters is different for the different situations, and it is specifically described for each case in the text.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** General scheme for the validation of the error estimates.</figcaption></figure>

where $\hat{a}_{retr}$ is the retrieved parameter by the GRASP algorithm, and $a_{assumed}$ is the parameter assumed in the input data for the generation of the synthetic observation. Equation (34) is used for each retrieved parameter, including the size distribution value at each size bin, the values of complex refractive index at each wavelength, the values of aerosol vertical profile at each altitude and the values of spherical particle fraction. We also implemented the evaluations of the errors for aerosol SSA and other parameters that are not part of the directly retrieved parameter while it is a function of the retrieved parameters, and it is estimated based on $\hat{a}_{retr}$. Thus, the retrieval error variances estimated by GRASP can be compared with the calculated actual retrieval errors. It should be noted here that we have always verified that the errors in the retrieval realized by GRASP from the error-free synthetic data (i.e. with no error specifically added) are negligibly small.

The GRASP-generated variances of the retrieval errors are evaluated in the presence of random errors and analysed using a series on the numerical tests conducted for a statistically representative set of random error realizations. These results are then summarized for the whole series of the tests by figures and tables. The tests with added systematic errors are discussed for most of the separate systematic error types, while some overall summaries are also provided.

As was mentioned before, in addition to the standard deviation, the non-diagonal elements of covariance matrices provide additional important insight about the retrieval quality. This additional information mainly relates to non-zero correlation coefficients. Therefore, in order to illustrate the correlation

structure, in this work we also analysed the correlation matrix that contains the covariance matrix elements normalized by the respective variances, as shown by Eq. (33). Our studies are not attempting to evaluate correlation matrices provided by the GRASP algorithm, since this would require the efforts exceeding the scope of this paper. Instead, we try to provide several demonstrations of how the structure of the correlation matrix may help to understand several interesting observations in the existing retrieval experience.

### 3.3 Aerosols models and realizations used in the tests

The synthetic tests were performed for several preselected realizations of aerosol in the atmosphere. These realizations were selected based on extensive experience with aerosol retrieval from sun/sky radiometer data and their combination with co-located lidar data. It is expected that the selected aerosol realization scenarios are representative of the majority of distinct actual observations of atmospheric aerosols.

Two main observational scenarios are considered for different total aerosol optical depth at 440 nm.

i. Single aerosol, such as biomass burning (BB), urban and dust for different aerosol loads $\tau(440) = 0.3, 0.6$ and $0.9$.

ii. The mixture of dust with BB and with urban (BB-Dust and Urban-Dust) is given. For each mixture, we have selected nine different scenarios that correspond to three different aerosol loads, $\tau(440) = 0.2, 0.5$ and $1.0$, where the different cases of the partition between
the fine and coarse mode were as follows: $τ_f/τ_c = 4.0$, $τ_f/τ_c = 1$, and $τ_f/τ_c = 0.25$.

The single-aerosol and aerosol mixture observational scenarios are used in the generation of synthetic tests with sun/sky-photometer-only observations. By considering both a single-aerosol and two aerosol types, in this work we evaluated how the accuracy of the retrieved data evolves once a larger number of parameters are derived from the same information content. In contrast, the retrieval based on the synergy between lidar and sun/sky photometers aimed for the retrieval of the properties of two fine- and coarse-mode aerosol components; therefore, the numerical tests for this type of the retrieval rely on the mixed aerosol observation scenario. At the same time, the error estimation is also checked in the case when the joint radiometer and lidar observations of single aerosol are analysed. The aerosol properties description used for the synthetic cases can be founded in Table 1 of Torres et al. (2017). They were modelled using the climatology of aerosol retrievals from AERONET observations described by Dubovik et al. (2002a). The dynamic climatological model from Mongu (Zambia) was used for BB aerosol, the model from the Goddard Space Flight Center (GSFC, Maryland, USA) for urban aerosol and the model from the Solar Village (Riyadh, Saudi Arabia) for dust aerosol. The real refractive index (RRI) and imaginary refractive index (IRI) for $λ = 355, 532$ and $1064$ nm (lidar measurements) were obtained by the extrapolation of the values from Dubovik et al. (2002a), as was suggested by Torres et al. (2017). All scenarios were simulated assuming a solar zenith angle (SZA) equal to 75°.

The retrieval settings were used similar to those that conventionally used in the retrieval of aerosol from AERONET sun/sky radiometer observations by Dubovik and King (2000) and from combined observations by sun/sky radiometer observations and lidar by Lopatin et al. (2013, 2021). Specifically, in the retrievals from sun/sky-radiometer-only observations, the size distribution (SD) was simulated using 22 logarithmically equidistant size bins between 0.05 and 15 µm. In the retrieval of the aerosol mixture from combined observations by sun/sky radiometers, the size distribution is modelled using 10 logarithmically equidistant bins between 0.05 and 0.58 µm for the fine mode and 15 logarithmically equidistant bins between 0.33 and 15 µm for the coarse mode. A similar approach was employed in the retrieval from sun/sky-radiometer-only observations when bi-component aerosol model was retrieved.

As mentioned above, in the case of the joint processing of AERONET radiometer and lidar data, we considered the approach developed earlier by Lopatin et al. (2013). In frame of this approach, the aerosol is modelled as an external mixture of two components. These components are characterized by height-independent microphysical properties, including the size distribution (represented by several size bins) and spectrally dependent complex refractive index. Moreover, each

component is described by the detailed vertical profile of the volume concentration. Therefore, the retrieval provides height-independent columnar properties of each component (size distribution and complex refractive index) and two profiles of fine- and coarse-mode volume concentrations. It is expected that this model is sufficient to adequately describe both radiometric and lidar observations.

# 4 Test results

Several tests were realized to evaluate the error estimates reliability and usefulness in the presence of both random and systematic uncertainties for aerosol retrievals from the observations of sun/sky radiometers alone and in combination with lidar. In this section, the results for the following two scenarios are presented: (i) a simpler case when only one type of aerosol is present and (ii) a more complex case in which two distinct types of aerosol are present at the same time. Moreover, we estimate the correlation matrices for both scenarios and illustrate their usefulness for understanding retrieval error tendencies, thus optimizing the retrieval approach.

## 4.1 Random error analysis

In a series of these tests of all inverted the synthetic measurements, we added random noise with a standard deviation of $ε_{Δτ}(λ) = 0.01$ for AOD, $ε_{ΔI/τ}(λ) = 5\%$ for radiances in order to model realistic uncertainties in AERONET observations (Holben et al., 1998; Eck et al., 1999; Dubovik et al., 2000; Sinyuk et al., 2020), $ε_{355} = 0.2$, $ε_{532} = 0.15$ and $ε_{1064} = 0.1$ for lidar attenuation measurements that vary with the altitude, as explained by Lopatin et al. (2013, 2021).

### 4.1.1 Retrieval of the single-aerosol component from radiometer measurements

This section describes the evaluation of the error estimates, assuming the presence of only one type of aerosol, i.e. BB, urban or dust. As mentioned before, the retrieval aerosol properties under the assumption of the presence of a single-aerosol type composed of homogeneous particles is a well-established approach for deriving detailed aerosol properties from ground-based observations by a sun/sky radiometer that is adapted by the operational AERONET retrievals by Dubovik and King (2000). The detailed error analysis of AERONET inversion aerosol product was provided by Dubovik et al. (2000), Torres et al. (2017) and by the recent study of Sinyuk et al. (2020) that described the uncertainty approach adapted in for AERONET version 3 retrieval products.

Figures 2–4 illustrate the error variances estimated by GRASP for all retrieved aerosol parameters in the selected synthetic tests for the observation of BB, urban and dust with different aerosol loads, i.e. $τ(440) = 0.3, 0.6$ and $0.9$. The displayed error bars for the standard deviation are calculated
from the diagonal elements of the covariance matrix (Eq. 32). Some tendencies can be seen from these illustrations. For example, the errors in the SD in the extremes (for the largest and smallest particles) are the biggest. This is an expected tendency, since these particles typically have a lower contribution to the measured signal (radiances and aerosol optical depths) compared to the particles of intermediate radius.

The retrievals improve and the errors decrease when the aerosol load increases, specially for IRI and SSA (absorption information). For BB and urban, the SSA error increases with the wavelength. On the other hand, SSA error decreases with the wavelength for dust. This is an expected behaviour, since the scattering efficiency is more pronounced at short wavelengths for small particles, while it is somewhat increasing with wavelength for large particles. Furthermore, as shown in Fig. 2a, the observed underestimation in the SD fine mode seems to be related to an overestimation in RRI.

To evaluate the error estimates in the presence of random noise, a set of the simulations for 300 different realizations of noise modelled using a random number generator has been analysed in this work. The results of such numerical tests conducted with a statistically representative set of random errors are summarized and illustrated, using box plots of the errors, as demonstrated in Fig. 5 for SSA(675) values. In the upper part of the figure, the box represents 50% of the data, with the whiskers representing 5th and 95th percentiles of the data, the solid line in the box plot representing the median, and the points are the mean values.

Figure 6 shows the distributions of the error estimates provided by GRASP (in red) and the calculated errors (in blue) for the cases when $τ(440) = 0.6$ and the following convergence criteria are satisfied: $Δτ ≤ 0.01\%$ and $ΔI/I ≤ 5\%$. It can be seen that, overall, the error estimates provided by GRASP capture the actual error tendencies, with some overestimation of their values, quite well. Thus, the retrieved errors can be considered to be the upper estimates of actual errors. This observed general overestimation can be, at least partially, explained by the fact that the error estimates by Eqs. (23) and (24) rely on a linear approximation. In this respect, it is known from practice that the nonlinear effects often lead to some saturation, while that cannot be captured by linear estimates. Some interesting tendencies can be appreciated in the obtained illustrations. For example, errors in SSA increase with the wavelength for BB and urban and decrease for dust.

On the other hand, the RRI errors seem to be similar at the different wavelengths. This is likely related to the fact that, spectrally, RRI retrievals rely on rather strong smoothness constraints on spectral variability in RRI (e.g. Dubovik and King, 2000). In contrast to the RRI, a large variability in the calculated errors is observed in the distribution of the errors for the IRI. Indeed, in order to capture the possible real spectral variability in IRI as being that of dust (e.g. see Dubovik et al., 2002a), the IRI is retrieved under milder smoothness constraints on the spectral variability (see Dubovik and King,

2000). Some of the aforementioned and other tendencies in the retrieval errors will be further discussed and evaluated in Sect. 4.3.1, which deals with error correlation matrices.

Table 2 summarizes the evaluation of the error estimates represented in the box plots. It provides the mean values for each parameter (RRI, IRI and SSA) at different wavelengths. These values correspond to the situation with a solar zenith angle equal to 75°. The obtained estimates compare reasonably with the corresponding values provided by in Table 4 in the paper by Dubovik et al. (2006). Specifically, the RRI error at 440 nm provided by GRASP for BB is 0.079 (0.04), where the values in parenthesis are from Dubovik et al. (2000), and for urban it is 0.056. The IRI error at 440 nm is 24% (30%) for BB, 54.1% for urban and 24.4% (50%) for dust. The values of the SSA errors are 0.028 (0.03) for BB, 0.013 for urban and 0.014 (0.03) for dust. At the same time, the RRI error at 440 nm provided by GRASP for dust is 0.201 (0.04) and is quite different, although Dubovik et al. (2000) considered only spherical particles. Moreover, the error estimates for SSA are consistent with the U27 estimates provided by Sinyuk et al. (2020). For example, at 440 nm for AOD = 0.6, the corresponding value for GSFC is 0.017, while GRASP provides values of 0.013; for Mongu, it is 0.023, while GRASP provides the error as 0.028.

### 4.1.2 Retrieval of mixed aerosol properties from measurements of a radiometer only

As already mentioned, most conventional aerosol retrievals from ground-based radiometer measurements (e.g. Dubovik and King, 2000; Nakajima et al., 2020) assume that aerosol is represented by homogeneous polydisperse particles with the size-independent refractive index. At the same time, this condition is not always correct in reality. Moreover, it is likely somewhat incorrect in the majority of the cases. Dubovik et al. (2000) showed that in some cases the retrieval assuming homogeneous particles would provide an effective index of refraction that allows the reproduction of the scattering properties of mixed aerosol rather adequately. Nonetheless, the assumption of homogeneous particles is often questioned and revisited (Xu et al., 2015); therefore, considerations of aerosol inhomogeneity are also included in the present study. In this regard, while the retrieval of the multicomponent aerosol is not a part of the standard AERONET inversion, the GRASP algorithm allows the retrieval of several aerosol components from diverse remote sensing observations, including the case of aerosol retrieval from radiometer measurements only.

At the same time, since the retrieval of multicomponent aerosol from radiometers only is not often used and not employed for operational retrievals, the tests in this section are limited to several illustrations only, and no statistical evaluation is performed. The illustrations are produced for the observations of a mixture of Urban-Dust and BB-Dust (see Sect. 3.3) for three cases of total $τ(440) = 0.2, 0.5$ and $1.0$.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Aerosol properties retrieved from simulated sun photometer data, with random noise added for BB aerosols for τ(440) = 0.3, 0.6 and 0.9 (left to right). The solid lines indicate the simulated properties (SD, RRI, IRI and SSA), and the dashed lines are the retrieved parameters. The shaded areas indicate error estimated by GRASP algorithm.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Aerosol properties retrieved from simulated sun photometer data, with random noise added for urban aerosol for τ (440) = 0.3, 0.6 and 0.9 (left to right). The solid lines indicate the simulated properties (SD, RRI, IRI and SSA), and the dashed lines are the retrieved parameters. The shaded areas indicate error estimated by GRASP algorithm.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Aerosol properties retrieved from simulated sun photometer data, with random noise added for dust aerosol for $\tau(440) = 0.3$, 0.6 and 0.9 (left to right). The solid lines indicate the simulated properties (SD, RRI, IRI and SSA), and the dashed lines are the retrieved parameters. The shaded areas indicate error estimated by GRASP algorithm.</figcaption></figure>

<table><caption>**Table 2.** Errors provided by GRASP for the RRI, IRI and SSA are represented by the mean values of each box plot for their respective wavelength. Absolute errors are given for RRI and SSA and relative errors for IRI. Mean values of actual errors are provided in parenthesis.</caption><thead><tr><th rowspan="2"></th><th colspan="3">BB</th><th colspan="3">Urban</th><th colspan="3">Dust</th></tr><tr><th>RRI</th><th>IRI (%)</th><th>SSA</th><th>RRI</th><th>IRI (%)</th><th>SSA</th><th>RRI</th><th>IRI (%)</th><th>SSA</th></tr></thead><tbody><tr><td>440</td><td>0.079<br>(0.014)</td><td>24.0<br>(6.66)</td><td>0.028<br>(0.005)</td><td>0.056<br>(0.016)</td><td>54.1<br>(17.6)</td><td>0.013<br>(0.004)</td><td>0.201<br>(0.03)</td><td>24.4<br>(11.52)</td><td>0.014<br>(0.006)</td></tr><tr><td>675</td><td>0.082<br>(0.018)</td><td>10.9<br>(7.8)</td><td>0.033<br>(0.005)</td><td>0.053<br>(0.015)</td><td>33.3<br>(19.8)</td><td>0.014<br>(0.005)</td><td>0.17<br>(0.013)</td><td>24.2<br>(11.36)</td><td>0.008<br>(0.004)</td></tr><tr><td>870</td><td>0.084<br>(0.019)</td><td>11.29<br>(11.61)</td><td>0.042<br>(0.009)</td><td>0.051<br>(0.015)</td><td>31.55<br>(28.57)</td><td>0.017<br>(0.009)</td><td>0.17<br>(0.011)</td><td>28.2<br>(14.3)</td><td>0.008<br>(0.003)</td></tr><tr><td>1020</td><td>0.081<br>(0.017)</td><td>13.17<br>(13.52)</td><td>0.043<br>(0.014)</td><td>0.052<br>(0.015)</td><td>35.1<br>(27.3)</td><td>0.021<br>(0.011)</td><td>0.167<br>(0.011)</td><td>30.4<br>(11.7)</td><td>0.007<br>(0.002)</td></tr></tbody></table>

In particular, we illustrate the case for $\tau(440) = 1.0$, with $\tau_f = 0.8$ and $\tau_c = 0.2$, $\tau_f = \tau_c = 0.5$ and $\tau_f = 0.2$, and $\tau_c = 0.8$, anticipating more potential for the adequate retrieval of multicomponent aerosol since the effect of AOD errors decreases for higher AOD. The analysis is focused on the possibilities of the differentiation between the properties of fine- and coarse-mode aerosol parameters such as complex refractive indices, size distributions, single scattering albedo.

Figures 7 and 8 illustrate the results of bicomponent retrievals and their error estimates from observations of the sun/sky radiometers of mixed aerosol. In the same figures, we also show a magnified plot for the effective RRI and IRI and the total SSA with their errors. Several retrieval tendencies are evident from the figures. For example, in the presence of one mode dominating in optical thickness, the retrievals and error estimates of the dominating component are more accurate. For example, in Fig. 7c for $\tau_f = 0.2$ and $\tau_c = 0.8$, the re-
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Comparison of the variance SSA(675) values estimated by the GRASP algorithm with actual errors obtained for extensive tests with randomly added modelled errors. In the upper panel, the box represents 50% of the data, with the whiskers representing 5th and 95th percentiles of the data, and the solid line in the box plot representing the median value.</figcaption></figure>

trievals of the coarse-mode properties are more accurate. An opposite behaviour can be seen in Fig. 7a for $τ_f = 0.8$ and $τ_c = 0.2$, when the predominance is in the fine mode. The clear trend can be observed in spectral dynamic of the error values for SSA because the error increases with the wavelengths in the fine mode and decreases for the coarse mode.

The most obvious difficulties in the separation of modes are evident when the properties of each mode are not very different. For example, such a situation can be seen for IRI of the Urban–Dust mixture (Fig. 7) and for RRI of the BB–Dust mixture (Fig. 8). In such situations, the error variances of each parameter are large and likely correlated (more details are provided in the discussion of covariance matrices). However, it is very important to note that, while the discrimination of some parameters of each component separately is not evident, most of the total and effective properties (magnified plots) can be estimated rather accurately.

It should be noticed that the retrieval of the multicomponent aerosol from radiometric observations was added here to illustrate error transformation tendencies that are not evident. Indeed, the retrieval of multicomponent aerosol from the AERONET observation is very uncertain, as was shown, for example, by Dubovik et al. (2000). At the same time, some sensitivity to the presence of multicomponent aerosols exists, and the attempts of multicomponent retrievals distinguishing the two components are often encouraged by the aerosol community, especially for favourable situations when two aerosol components have comparable influence on AERONET measurements (the situation used in our tests). At the same time, such cases are very suitable for demonstrating that, if constraints are not sufficient, the errors can

be unacceptably high and strongly correlated. At the same time, one can see that when some estimates are highly uncertain and strongly correlated, they still can be used for accurate estimation of their functions. For example, we showed that a property such as the total SSA of mixed aerosol can be rather accurately obtained from the retrieved SSA of fine and coarse modes. In the next section, we will also illustrate the improvements in multicomponent retrieval and error reduction when extra information is included, such as lidar measurements. In addition, the retrieval of the multicomponent aerosol used often becomes robust when co-located AERONET data are inverted together with data from co-located lidar (Lopatin et al., 2013, 2021).

### 4.1.3 Retrieval of mixed aerosol properties from the measurements of radiometers in combination with lidar

The GRASP aerosol retrieval from combined sun/sky radiometers and lidar observations were always designed for the retrieval of bicomponent aerosol (Lopatin et al., 2013), and the approach is employed for operational processing in the frame of ACTRIS activities. Therefore, the evaluation of the random error effect in the aerosol retrieval from radiometer and lidar observations of aerosol mixtures includes both the analysis of the selected illustrations and the statistically representative series of numerical tests with random errors. The considered synthetic data include synthetic observations produced for the same examples of aerosol mixtures (Urban–Dust and BB–Dust) as used in Sect. 4.1.2.

Figures 9 to 12 illustrate the retrievals and their error estimates obtained for aerosol properties of both fine- and coarse-mode aerosols. The good agreement of the actual retrieved parameters (solid lines) with the assumed values (dashed lines) can be seen for all cases. From a comparison of Figs. 7–8 with Figs. 9–10, it is easy to see that the retrieval error estimate is lower when lidar data are also used. The improvements (compared to the results from radiometer only retrievals) are especially evident in the separation of the retrieved aerosols properties, especially when the contribution of the aerosol load is lower (this will be shown in Sect. 4.2 when bias is also assumed). At the same time, the total properties are accurately estimated in both retrieval scenarios.

Figure 11 shows the lidar ratio of the fine mode, coarse mode and total aerosol for the three aforementioned cases. In general, good agreements of retrieved and assumed values are obtained, especially for the total lidar ratio (LR). However, there are some discrepancies at short wavelengths for fine-mode lidar ratios.

Figure 12 illustrates the retrieval of aerosol vertical profile for each case. The agreement between the retrieved and assumed values of the vertical profiles is good, mainly for the coarse mode at the altitudes where it has maximum values and dominates. At the altitudes where there is a superposition of aerosol layers with a comparable presence of both
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Comparison of estimated and actual error distributions for spectrally dependent aerosol parameters retrieved from sun/sky-photometer-simulated measurements (a case with τ(440) = 0.6). The distributions were obtained using 300 realizations of added random errors. The median values of the errors are shown with a line in the box plot, along with the 25th–75th percentiles indicated by a box and the 5th–95th percentiles indicated using whiskers. The mean values are represented by the black dot. The red colour shows the error estimates provided by GRASP, and the blue shows the calculated actual errors (Eq. 34).</figcaption></figure>

aerosols, the retrieval struggles to discriminate the contribu-
tion of both modes, and a clear overestimation of the fine
mode (and, consequently, an underestimation of the coarse
mode) can be seen.

In order to evaluate the error estimates in the presence of
random errors, a set of simulations, adding 300 realization of
random noise values, is analysed. Figures 13 to 19 show the
comparisons of all retrieved aerosol parameters separately
for fine and coarse aerosol modes. In addition, the retrieval of
total SSA and LR are shown. The case for total τ(440) = 1.0
is shown more extensively, as in the previous section, due to
the interest in the retrieval of the situation with higher aerosol
loads. The main result that can be gained from illustrations
is that the GRASP error estimates are typically higher than
actual errors; this same result was obtained for the retrieval
of only sun/sky radiometer data.

Figures 13 and 14 illustrate the comparisons of distribu-
tions of the GRASP error estimates and actual errors for RRI

and IRI of fine and coarse aerosol modes for situations when
mixtures of Urban–Dust and BB–Dust are observed. It can
be seen that the accuracy of the refractive index retrievals for
each mode depends strongly on the contribution of the mode
to the signal, as was observed by Lopatin et al. (2013). For
example, if we analyse the performance of the fine mode,
the higher the contribution of the fine optical thickness, the
better the accuracy in the retrievals of fine-mode aerosol pa-
rameters.

Figures 15 and 16 show the situation for the error distribu-
tion for SSA of fine mode, coarse mode and total. Similar to
that observed in earlier tests for the retrieval of fine-mode pa-
rameters, the errors increase with the wavelength, while for
the retrieval of the coarse-mode parameters, the errors de-
crease with the wavelengths. Also, the results show that the
errors in the total SSA are rather small, even if the SSA of
fine and coarse modes are quite high.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Aerosol properties retrieved from simulated sun photometer data with random noise added for a mixture of Urban-Dust aerosols. The solid lines indicate the simulated properties (SD, RRI, IRI and SSA), and the dashed lines are the retrieved parameters. The shaded areas indicate the error estimated by the GRASP algorithm. The magnified plots represent the effective refractive index and total SSA.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Aerosol properties retrieved from simulated sun photometer data with random noise added for a mixture of BB-Dust aerosols. The solid lines indicate the simulated properties (SD, RRI, IRI and SSA), and the dashed lines are the retrieved parameters. The shaded areas indicate the error estimated by the GRASP algorithm. The magnified plots represent the effective refractive index and total SSA.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Aerosol properties retrieved from simulated sun photometer and lidar data with random noise added for a mixture of Urban–Dust aerosols. The solid lines indicate the simulated properties (SD, RRI, IRI and SSA), and the dashed lines are the retrieved parameters. The shaded areas indicate the error estimated by the GRASP algorithm. The magnified plots represent the effective refractive index and total SSA.</figcaption></figure>

The error evaluation for LR is represented in Figs. 17 and 18. In most of the cases, we see good agreements between the error estimations and actual error. The only exception is the errors in the LR of the fine mode at short wavelengths, where the actual errors are higher than the errors provided by GRASP. This tendency seems to be anticorrelated with the results found for the coarse-mode LR error estimates, where the GRASP error estimates are notably higher than the actual values.

The results illustrated by the figures are summarized in Tables 3 to 5. These tables show the mean values of the GRASP error estimates for the cases when the total $τ(440) = 1.0$ and $τ_f = τ_c = 0.5$, i.e. when there is no predominance of either one of the modes. The values are provided for the aerosol parameter considered at different wavelengths, both for the simulation of Urban–Dust and BB–Dust observations calculated for a case of the SZA $= 75°$.

For retrieval errors in fine mode in the case of urban aerosol parameters, the mean values for RRI are around 0.05, and the values do not present much variability with the wavelength. For the retrieval of IRI, the mean values of the GRASP retrieval errors are at the level of around 73%, showing a pronounced underestimation with respect to the actual error at short wavelengths. With respect to SSA errors provided by GRASP, a clear tendency is observed because the error increases with the wavelengths from 0.024 to 0.061.

Finally, the mean values of LR errors provided by GRASP decrease with the wavelength between 15% and 10%, with notable underestimations with respect to the actual errors at short wavelengths. In the case of the retrieval of fine-mode BB parameters, the mean values for RRI errors provided by GRASP are around 0.05. Some underestimations with respect to the actual errors are observed at short wavelengths. The mean values for IRI errors are around 60%, and the errors for SSA show a clear tendency to increase with the wavelengths between 0.04 to 0.09. Mean values of LR errors provided by GRASP decrease with the wavelength between 18% and 14%.

The mean values of the error estimates provided by GRASP for dust present a good agreement in the case of both mixtures. In general, the mean values for RRI error estimates vary between 0.07 and 0.09, and they do not present much variability with the wavelength, while smaller values of errors are seen for the Urban–Dust mixture case. The mean values of IRI error estimates are around 50%, while for the BB–Dust mixture, we observe some underestimations of the actual errors by GRASP calculations. The errors in SSA show a clear tendency that decreases with the wavelengths between 0.04 and 0.009. The mean values for the LR retrievals increase with the wavelength from 37% to 60%, with bigger errors observed for the BB–Dust mixture.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Aerosol properties retrieved from simulated sun photometer and lidar data with random noise added for a mixture of BB-Dust aerosols. The solid lines indicate the simulated properties (SD, RRI, IRI and SSA), and the dashed lines are the retrieved parameters. The shaded areas indicate the error estimated by the GRASP algorithm. The magnified plots represent the effective refractive index and total SSA.</figcaption></figure>

Once again, it is important to note that the errors in the parameters characterizing the total aerosol are generally accurately estimated. For both cases of Urban-Dust and BB-Dust mixtures, the mean values of the total SSA error estimates vary between 0.02 and 0.009, and the mean values of the total LR error estimates range from 23% to 55%.

Figure 19 shows the relative errors in AVP retrievals for fine and coarse aerosol modes for the Urban-Dust and BB-Dust aerosols mixture. The errors estimated by GRASP are a bit higher than the errors obtained by the simulations of random errors; correspondingly, the GRASP errors can be safely used as the upper estimates of actual retrieval uncertainties. Table 5 summarizes the evaluation of the errors estimates for all the scenarios discussed above. The GRASP estimates of the retrieval errors for both mixtures are between 50%–70% for the fine mode and 50%–57% for the coarse mode.

Finally, a lower sensitivity to the retrieval of fine-mode properties can be observed as a clear tendency in the evaluation of the retrieval errors for the cases when mixed aerosols are analysed. In particular, quite high errors were obtained for the complex refractive index. Then, these errors consequently propagate to the errors in other optical properties, such as the SSA of fine mode, as was found in the earlier study by Lopatin et al. (2013).

## 4.2 Analysis of the retrieval in presence of the systematic uncertainties

In Sect. 4.1, the evaluation and validation of the errors in the different aerosol properties, considering propagation of the random noise from measurements into retrieval, was presented. The analysis confirmed the rather satisfactory performance of the approach adapted in GRASP for the estimation of retrieval errors in the presence of random noise. This section discusses the approach for estimating the possible contributions of the systematic errors in the retrieval uncertainties. In principle, each retrieval methodology assumes that there are no systematic uncertainties in the measurements or in the forward model used. If any systematic bias is identified, it is corrected in the measurements or in their interpretation. However, in practice, the systematic uncertainties may remain unidentified and make a significant contribution to the retrieval uncertainties.

As mentioned above, in Eq. (24), the apparent misfit was used as an indicator of bias; however, in real situations, not all biases can be seen from the misfit. Thus, in this section, the results are presented considering a possible solution to this problem. Therefore, the contribution of potential bias is commonly included in the estimation of the retrieval errors (e.g. see Dubovik et al., 2000; Sinyuk et al., 2020). Using similar logic, an extra term was added in the present method-
## Urban and Dust

$$ \tau_f = 0.8 \quad \tau_c = 0.2 $$

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

$$ \tau_f = 0.5 \quad \tau_c = 0.5 $$

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

$$ \tau_f = 0.2 \quad \tau_c = 0.8 $$

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

## BB and Dust

$$ \tau_f = 0.8 \quad \tau_c = 0.2 $$

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

$$ \tau_f = 0.5 \quad \tau_c = 0.5 $$

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

$$ \tau_f = 0.2 \quad \tau_c = 0.8 $$

<figure><img src="image_18.png" /><figcaption>**Figure 11.** The aerosol lidar ratio (LR) retrieved from simulated sun photometer and lidar data with random noise added for a mixture Urban–Dust aerosols (above) and BB–Dust (below). The solid lines indicate the simulated properties (SD, RRI, IRI and SSA), and the dashed lines are the retrieved parameters. The shaded areas indicate the error estimated by the GRASP algorithm. The magnified plots represent the results for total LR retrievals.</figcaption></figure>

ogy, Eq. (31), and then that accounts for the propagation of possible bias from the measurements. The propagation accounted for the linear approximation in similar manner to the systematic term in Eq. (28) that accounts for bias from misfit. Thus, this section analyses the potential effect of realistic biases and their overall importance for reliable estimations of the retrieval errors in practice.

The potential effect of the systematic errors is analysed in series of numerical tests with possible assumed systematic errors. Following previous studies by Dubovik et al. (2000); Torres et al. (2017); Sinyuk et al. (2020) in ground-based photometric and radiometric data, we consider two types of potential main biases in measured AOD and sky radiances. These biases could originate from the miscalibration of direct sun or diffuse sky sensors (Eck et al., 1999). The biases are assumed to be wavelength independent, and since spectral systematic deviations are easier to identify in a direct analysis of the observation, they are likely to be manifested in the misfit and may compensate for each other's influence on the retrievals. Specifically, two possible levels of biases, nominal and maximum, are considered to be as follows:

i. In AOD, there is a nominal bias of ±0.01 and maximum bias of ±0.02.

ii. In radiances, there is a nominal bias of ±3% and maximum bias of ±5%.

To evaluate the effects of biases, the above values were added to the synthetic direct measurements of AOD and sky radiance by an AERONET-like ground-based radiometer. These data were inverted by the GRASP code, and the retrieved values of aerosol parameters were compared to the values assumed in the synthetic simulation as a truth. In addition, the deviations in the retrieved values from the true ones are compared to the errors estimates generated by GRASP based on Eq. (31) using the known values of the added biases. In a similar manner, the influence of the potential systematic errors in aerosol retrieval from the combined observations of a ground-based radiometer and lidar was analysed. In these tests, the biases in the lidar attenuation measurements were assumed for each wavelength following the studies by Lopatin et al. (2013, 2021), with $ε_{355} = ±0.2$, $ε_{532} = ±0.15$ and $ε_{1064} = ±0.1$. It should be noted that the conducted synthetic tests not only allow the verification of the accuracy of the systematic error estimates by GRASP but also allow us to analyse the effects of biases on the retrievals for different retrieval scenarios in diverse situations.
<figure><img src="image_19.png" /><figcaption>Figure 12. The aerosol vertical profiles (AVPs) retrieved from simulated sun photometer and lidar data with random noise added for a mixture of Urban–Dust aerosols (a, b, c) and BB–Dust (d, e, f). The solid lines indicate the simulated properties (AVP), and the dashed lines are the retrieved parameters. The shaded areas indicate the error estimated by the GRASP algorithm.</figcaption></figure>

### 4.2.1 Effects of measurement bias in the retrieval of single-aerosol components from radiometer measurements

In this section, the study is focused on the analysis of the effects of the biases and on estimating contribution of systematic errors in the retrievals of aerosol from ground-based observations by a radiometer. In a similar manner to the analysis of random errors, we first considered the observations dominated by two types of aerosols, i.e. BB and dust. The effect of measurement biases is expected to be manifested in the situations with low and moderate aerosol loading; therefore, the analysis is focused on the scenarios with AOD(440) = 0.1, 0.3 and 0.6.

The following two situations were considered:

i. when a single bias in AOD or radiances is present, and

ii. when the biases can be present in both AOD and radiances simultaneously. In this case, the different combinations of positive and negative biases in AOD and radiances are considered.

The estimations of the errors introduced by the biases were calculated as follows:

$$ \sigma_{\text{bias}}^2 = \sigma_{\text{lm}}^2 + \sigma_{\text{misfit}}^2 + \frac{1}{N} \sum_{k=1}^{N} \sigma_k^2, \qquad (35) $$

where $\sigma_{\text{lm}}^2$ corresponds to contributions from systematic errors introduced by the Levenberg–Marquardt procedure, and $\sigma_{\text{misfit}}^2$ are the errors manifested by the misfit estimated by Eq. (28), and each $\sigma_k^2$ is the contribution adding + bias and – bias in the measurements.

Figures 20 to 23 illustrate the results of the analysis for the different retrieved properties for the situation with a bias of
<figure><img src="image_20.png" /><figcaption>**Figure 13.** Comparison of estimated and actual error distributions for spectrally dependent aerosol parameters retrieved from measurements by simulated sun/sky photometer values and lidar for a mixture of Urban–Dust aerosol. The distributions were obtained using 300 realizations of added random errors. The median values of the errors are shown by a line in the box plot, along with the 25th–75th percentiles indicated by a box, and 5th–95th percentiles are indicated using whiskers. The mean values are represented by the black dot. The red colour shows the error estimates provided by GRASP, and the blue colour shows the calculated actual errors (Eq. 34).</figcaption></figure>

±0.01 and ±0.02 in AOD only. These results show the spe-
cific effects from AOD bias. The figures have two blocks; on
the left are the retrievals with added positive bias in AOD,
and on the right, retrievals with negative bias are illustrated.
In both cases, the error bars represent the systematic com-

ponent adding the positive or negative bias, respectively. In
all the figures, the solid lines show the assumed value of the
parameters in the simulation, the dotted lines show the re-
trieved values, and the magnitudes of the estimated bias are
shown by the shaded areas. It should be noted that, for the
<figure><img src="image_21.png" /><figcaption>**Figure 14.** Comparison of estimated and actual error distributions for spectrally dependent aerosol parameters retrieved from measurements by simulated sun/sky photometer values and lidar for a mixture of BB-Dust aerosol. The distributions were obtained using 300 realizations of added random errors. The median values of the errors are shown by a line in the box plot, along with the 25th–75th percentiles indicated by a box, and 5th–95th percentiles are indicated using whiskers. The mean values are represented by the black dot. The red colour shows the error estimates provided by GRASP, and the blue colour shows the calculated actual errors (Eq. 34).</figcaption></figure>

case of BB with AOD(440) = 0.1, the results with negative bias are not shown. This is because the AOD for BB decreases very strongly with the wavelength, and for the case of AOD(440) = 0.1, the AOD at 1020 nm is ~ 0.01.

The figures show different and clear tendencies, which are in agreement with general expectations, and with the tenden-

cies already observed in previous studies by Dubovik et al. (2000) and Torres et al. (2014). For example, it can be seen that bias in AOD most strongly affects the estimate of the parameters characterizing the aerosol absorption, such as the imaginary part of the refractive index and single scattering absorption. This is an anticorrelation, where the positive bias
<figure><img src="image_22.png" /><figcaption>**Figure 15.** Comparison of estimated and actual error distributions for spectrally dependent aerosol parameters retrieved from measurements by simulated sun/sky photometer values and lidar for a mixture of Urban–Dust aerosol. The distributions were obtained using 300 realizations of added random errors. The median values of the errors are shown by a line in the box plot, along with the 25th–75th percentiles indicated by a box, and 5th–95th percentiles are indicated using whiskers. The mean values are represented by the black dot. The red colour shows the error estimates provided by GRASP, and the blue colour shows the calculated actual errors (Eq. 34).</figcaption></figure>

results in the overestimation of absorption (higher RRI and lower SSA) and the negative in the underestimation of absorption (lower RRI and higher SSA), respectively. The result was expected, since radiance values in this first experience do not vary. Thus, if we keep the scattering component (which is derived from radiances) but we enlarge the extinction component (by enlarging the AOD) necessary, then the retrieval understands that the absorption should be larger (imaginary part of the refractive index). Conversely, if we reduce the value of extinction, then the retrieval would reduce the value of absorption. Also, the strongest effect is observed for optically thin situations when a small absolute error in the optical thickness becomes comparable with the magnitude of aerosol optical thickness. This is especially clear for BB observations, where AOD(1020) is always rather small as also discussed earlier by Dubovik et al. (2000). For the observa-

tions of dust aerosol, the effect of biases in AOD are significantly smaller than for BB. It can be explained since dust has a small value of the Ångström exponent and therefore larger values of AOD at longer wavelengths. For the retrieval of the size distribution, the bias in AOD has a rather minor effect, though we found a general overestimation for positive bias values and an underestimation for negative values.

Overall, the estimated systematic error agrees well with the actual manifestations of the bias in the retrieval. The quantitative estimations are also quite convincing and shown in Figs. 20 and 22 for biases of ±0.01. In some cases, some underestimations of the bias effects can be observed. For example, the largest differences are identified for the case of the higher value of bias (±0.02), as shown in Figs. 21 and 23, while a significant increase in the systematic component of the retrieved error is also well captured by the error estimates.
<figure><img src="image_23.png" /><figcaption>**Figure 16.** Comparison of estimated and actual error distributions for spectrally dependent aerosol parameters retrieved from measurements by simulated sun/sky photometer values and lidar for a mixture of BB-Dust aerosols. The distributions were obtained using 300 realizations of added random errors. The median values of the errors are shown by a line in the box plot, along with the 25th-75th percentiles indicated by a box, and 5th-95th percentiles are indicated using whiskers. The mean values are represented by the black dot. The red colour shows the error estimates provided by GRASP, and the blue colour shows the calculated actual errors (Eq. 34).</figcaption></figure>

It can be seen that, among all considered aerosol parameters, the main differences between the bias effects and the obtained error estimates are observed for the real part of the refractive index (RRI). In these cases, the bias is not fully covered by the systematic component of the retrieved error results. Similarly, an apparent underestimation of the RRI errors was also seen by Sinyuk et al. (2020), who attributed these underestimations to different factors such as, for example, the effect of not accounting for the pointing bias. In the present simulations, there is no pointing bias considered, and the discrepancy is likely coming from the fact that Eqs. (24) and (31) rely on the derivatives estimated in the vicinity of the solution and based on linear approximation. Indeed, the dependence of both AOD and radiances scattered by aerosol is very complex and nonlinear. Therefore, with both taking derivatives in the vicinity of obtained solution instead of the

vicinity of the true values and the nonlinear character of AOD and radiances may explain the differences. At the same time, it is important to note that, as can be seen from the analysis of the random component of the RRI error (Sect. 4.1), the random error effect is likely to dominate over the effect of AOD bias, and therefore, the estimation of the total error (described below in this section) seems to allow us to make an objective and complete observation on this parameter.

Figures 24 to 27 show the effects of the biases in the radiances of the two magnitudes of ±3% and ±5% for the observations of BB and dust. In general, it can be seen from the results that, in both cases of BB and dust, the retrievals are less affected by bias in the radiances than by the biases in the AOD, even when the bias in the radiances is ±5%. A similar tendency was also reported in the studies by Dubovik et al. (2000), Torres et al. (2014) and Sinyuk et al. (2020).
<figure><img src="image_24.png" /><figcaption>**Figure 17.** Comparison of estimated and actual error distributions for spectrally dependent aerosol parameters retrieved from measurements by simulated sun/sky photometer values and lidar for a mixture of Urban–Dust aerosols. The distributions were obtained using 300 realizations of added random errors. The median values of the errors are shown by a line in the box plot, along with the 25th–75th percentiles indicated by a box, and 5th–95th percentiles are indicated using whiskers. The mean values are represented by the black dot. The red colour shows the error estimates provided by GRASP, and the blue colour shows the calculated actual errors (Eq. 34).</figcaption></figure>

At the same time, it should be noted that the present analysis is focused on the measurement configuration corresponding to solar almucantar when the SZA is 75° (see Table 1) and when the measurements include are taken in a wide range of scattering angles. In this respect, Dubovik et al. (2000) showed that the effect of the sky radiance bias increases when the range of observed scattering angles is limited, i.e. in almucantar observation corresponding to a SZA of less than 60°. Moreover, according to recent tendencies in observational practices, the use of such measurements is limited, and most analyses are focused on observational scenarios with a sufficient range of observed scattering angles. For example, AERONET start to establish so-called hybrid observational scenario during high SZA times (Giles et al., 2019). With regard to the performance of the error estimation, the effect of

the bias in the sky radiances seems to be well captured by the GRASP error estimates.

It should be also noted that we observe an anticorrelation between the radiance bias and the retrieval of the imaginary part of the refractive index. This effect is opposite to the one observed in the case of AOD and with significantly smaller differences. Thus, when the biases are positive (+3% and +5%), there is a decrease in the imaginary part of the refractive index. The fact that the value of AOD remains the same and that there is an increase in the value of the scattering is interpreted by the code as a decrease in the aerosol absorption. Conversely, the negative bias in the radiances produces an increase in the imaginary part of the refractive index which can be explained by the same reason.

Figure 28 shows the results of the analysis of the situation when the systematic biases present in both AOD and radi-
<figure><img src="image_25.png" /><figcaption>**Figure 18.** Comparison of estimated and actual error distributions for spectrally dependent aerosol parameters retrieved from measurements by simulated sun/sky photometer values and lidar for a mixture of BB-Dust aerosols. The distributions were obtained using 300 realizations of added random errors. The median values of the errors are shown by a line in the box plot, along with the 25th-75th percentiles indicated by a box, and 5th-95th percentiles are indicated using whiskers. The mean values are represented by the black dot. The red colour shows the error estimates provided by GRASP, and the blue colour shows the calculated actual errors (Eq. 34).</figcaption></figure>

ances are simultaneously assumed. The results for BB are on the left and for dust on the right. The illustrations are shown for the specific situation with two positive biases, namely +0.01 in AOD and +5% in radiances. It should be noted that the tests were produced for both situations with the biases of the same and opposite signs. This case with the biases of the same sign showed the most interesting results with the strongest manifestation of bias effects, and therefore, they are presented here. In the situation with biases of opposite signs, the effects on the retrievals are rather minor due to the internal compensations of the influences of the biases. Additionally, the misfit of the observations is more pronounced, which helps to identify the issues and account for the biases in the error estimates. Also, the analysis here is focused on the simultaneous biases of moderate values (±0.01 in AOD and ±5% in radiances), since the appearance of simultane-

ous biases of the highest bias values (i.e. ±0.02 in the AOD and ±5% in radiancies) lead to very strong effects in the retrievals. Those situations can be easily seen and screened out by quality filters (e.g. by the high value of the misfit). Also, it is quite unlikely to have such strong systematic errors in practical observations as those by AERONET.

As can be seen from Fig. 28, the biggest differences and highest bias values in the retrieval are found for low AOD (0.1). As seen earlier for this situation, the errors for the RRI remain notably underestimated. As already mentioned, the situation is expected to be improved once the effects of both random and systematic errors are considered.

Figure 29 illustrates such a situation for the retrieval of BB and dust for the three different aerosol loads (0.1, 0.3 and 0.6) when the total error estimate includes both random
<table><caption>**Table 3.** The mean values of RRI, IRI, SSA and LR retrieval errors estimated by GRASP for the synthetic test for a mixture of the Urban–Dust aerosol mixture. The mean values represent the distributions obtained using 300 realizations of the added random errors for the situation with total τ (440) = 1.0, with τ<sub>f</sub> = τ<sub>c</sub> = 0.5 and SZA = 75°. The absolute errors are provided for RRI and SSA and relative errors for IRI and LR. Mean values for the actual errors are provided in parentheses.</caption><thead><tr><th></th><th colspan="7">Urban–Dust</th></tr><tr><th></th><th>355 nm</th><th>440 nm</th><th>532 nm</th><th>675 nm</th><th>870 nm</th><th>1020 nm</th><th>1064 nm</th></tr></thead><tbody><tr><th>RRI<sub>f</sub></th><td>0.050<br>(0.033)</td><td>0.049<br>(0.030)</td><td>0.048<br>(0.029)</td><td>0.045<br>(0.031)</td><td>0.045<br>(0.036)</td><td>0.045<br>(0.036)</td><td>0.045<br>(0.036)</td></tr><tr><th>RRI<sub>c</sub></th><td>0.076<br>(0.028)</td><td>0.073<br>(0.028)</td><td>0.071<br>(0.028)</td><td>0.080<br>(0.021)</td><td>0.085<br>(0.017)</td><td>0.086<br>(0.015)</td><td>0.086<br>(0.015)</td></tr><tr><th>IRI<sub>f</sub> (%)</th><td>73.31<br>(103.03)</td><td>71.18<br>(103.4)</td><td>70.67<br>(94.19)</td><td>70.33<br>(81.46)</td><td>71.34<br>(76.62)</td><td>72.18<br>(76.78)</td><td>72.26<br>(76.79)</td></tr><tr><th>IRI<sub>c</sub> (%)</th><td>50.59<br>(36.17)</td><td>45.00<br>(26.45)</td><td>44.58<br>(22.32)</td><td>45.12<br>(18.50)</td><td>49.01<br>(18.23)</td><td>51.22<br>(17.92)</td><td>51.51<br>(17.89)</td></tr><tr><th>SSA<sub>f</sub></th><td>0.024<br>(0.017)</td><td>0.026<br>(0.019)</td><td>0.029<br>(0.020)</td><td>0.033<br>(0.022)</td><td>0.045<br>(0.027)</td><td>0.057<br>(0.034)</td><td>0.061<br>(0.037)</td></tr><tr><th>SSA<sub>c</sub></th><td>0.039<br>(0.042)</td><td>0.031<br>(0.025)</td><td>0.024<br>(0.014)</td><td>0.015<br>(0.007)</td><td>0.011<br>(0.004)</td><td>0.009<br>(0.003)</td><td>0.009<br>(0.003)</td></tr><tr><th>SSA<sub>T</sub></th><td>0.017<br>(0.009)</td><td>0.015<br>(0.004)</td><td>0.014<br>(0.004)</td><td>0.010<br>(0.005)</td><td>0.008<br>(0.004)</td><td>0.008<br>(0.003)</td><td>0.008<br>(0.003)</td></tr><tr><th>LR<sub>f</sub> (%)</th><td>14.93<br>(40.45)</td><td>12.11<br>(29.39)</td><td>10.41<br>(15.37)</td><td>8.88<br>(7.58)</td><td>8.92<br>(8.44)</td><td>9.79<br>(8.96)</td><td>9.96<br>(8.63)</td></tr><tr><th>LR<sub>c</sub> (%)</th><td>37.04<br>(17.42)</td><td>35.71<br>(12.01)</td><td>35.23<br>(9.48)</td><td>42.39<br>(11.81)</td><td>47.77<br>(14.90)</td><td>48.99<br>(15.17)</td><td>49.32<br>(14.99)</td></tr><tr><th>LR<sub>T</sub> (%)</th><td>23.31<br>(9.01)</td><td>26.46<br>(7.48)</td><td>28.21<br>(6.85)</td><td>35.61<br>(9.71)</td><td>41.45<br>(13.3)</td><td>43.36<br>(13.9)</td><td>43.97<br>(13.9)</td></tr></tbody></table>

and systematic components as follows:

$$ \sigma_{\text{tot}} = \sqrt{\sigma_{\text{ran}}^2 + \sigma_{\text{bias}}^2}, \qquad (36) $$

where $\sigma_{\text{bias}}^2$ is calculated as indicated in Eq. (35).

It can be seen that the total error estimates capture the deviations for all parameters in the presence of random and systematic noise. These results confirm that the estimations using Eq. (31), based on the additional assumptions of potential presence of bias in the measurements, improve the results of the error estimated compared to the approach discussed in Sect. 4.1, when the effects of biased were taken into account only based on the value of the observation misfit. The observed tendencies in the effects of biases on the retrieval are consistent with all the results previously described in earlier studies. The obtained results are expected to be representative of most practical situations, while some additional tests and analysis could certainly be useful. Therefore, in the examples presented below, and for the real cases analysed, the total error will be used as described in Eq. (36). This means that the representation of the error will take into account the contribution of the random and systematic component. This last component contains the contribution of Levenberg–

Marquardt and the misfit and the measurements in which the contributions of ± bias added in the measurements are considered (Eq. 35). These values of the assumed biases in our applications are consistent with AERONET as mentioned previously, i.e. ±0.01 in AOD and ±5% in radiances.

## 4.2.2 Effects of measurement bias in the retrieval of mixed aerosol properties from measurements of radiometers only

The present section tries to understand how the bias is affected when inhomogeneous aerosol are observed. The example is not commonly considered in practical application, e.g. in AERONET operational processing. At the same time, since GRASP can consider this type of bicomponent inversion that is fundamentally of high interest, we are analysing this situation in the presence of biases. In Sect. 4.2.1, we have shown different examples, considering the bias in each measurement separately, and we have also illustrated the complete example with presence of both random errors and bias in all the measurements. Here we directly illustrate the results, considering the presence of both random and bias in
<table><caption>**Table 4.** The mean values of RRI, IRI, SSA and LR retrieval errors estimated by GRASP for the synthetic test for a mixture of the BB-Dust aerosol mixture. The mean values represent the distributions obtained using 300 realizations of the added random errors for the situation with total τ (440) = 1.0, with τ<sub>f</sub> = τ<sub>c</sub> = 0.5 and SZA = 75°. The absolute errors are provided for RRI and SSA and relative errors for IRI and LR. Mean values for the actual errors are provided in parentheses.</caption><thead><tr><th></th><th colspan="7">BB-Dust</th></tr><tr><th></th><th>355 nm</th><th>440 nm</th><th>532 nm</th><th>675 nm</th><th>870 nm</th><th>1020 nm</th><th>1064 nm</th></tr></thead><tbody><tr><th>RRI<sub>f</sub></th><td>0.052<br>(0.064)</td><td>0.050<br>(0.067)</td><td>0.048<br>(0.061)</td><td>0.044<br>(0.052)</td><td>0.045<br>(0.046)</td><td>0.046<br>(0.047)</td><td>0.046<br>(0.047)</td></tr><tr><th>RRI<sub>c</sub></th><td>0.084<br>(0.019)</td><td>0.083<br>(0.019)</td><td>0.079<br>(0.019)</td><td>0.085<br>(0.017)</td><td>0.089<br>(0.016)</td><td>0.091<br>(0.014)</td><td>0.091<br>(0.014)</td></tr><tr><th>IRI<sub>f</sub> (%)</th><td>60.51<br>(35.65)</td><td>57.61<br>(35.40)</td><td>57.91<br>(37.56)</td><td>58.17<br>(41.63)</td><td>60.33<br>(43.74)</td><td>61.48<br>(44.14)</td><td>61.59<br>(44.13)</td></tr><tr><th>IRI<sub>c</sub> (%)</th><td>48.75<br>(33.75)</td><td>44.57<br>(41.39)</td><td>43.81<br>(48.20)</td><td>43.89<br>(60.36)</td><td>47.43<br>(62.61)</td><td>49.56<br>(55.55)</td><td>49.87<br>(55.44)</td></tr><tr><th>SSA<sub>f</sub></th><td>0.041<br>(0.031)</td><td>0.044<br>(0.031)</td><td>0.048<br>(0.033)</td><td>0.052<br>(0.040)</td><td>0.067<br>(0.050)</td><td>0.085<br>(0.057)</td><td>0.089<br>(0.059)</td></tr><tr><th>SSA<sub>c</sub></th><td>0.049<br>(0.032)</td><td>0.040<br>(0.031)</td><td>0.032<br>(0.027)</td><td>0.021<br>(0.019)</td><td>0.015<br>(0.013)</td><td>0.013<br>(0.010)</td><td>0.013<br>(0.009)</td></tr><tr><th>SSA<sub>T</sub></th><td>0.022<br>(0.011)</td><td>0.016<br>(0.004)</td><td>0.016<br>(0.005)</td><td>0.011<br>(0.004)</td><td>0.010<br>(0.003)</td><td>0.009<br>(0.003)</td><td>0.009<br>(0.003)</td></tr><tr><th>LR<sub>f</sub> (%)</th><td>17.99<br>(33.43)</td><td>15.17<br>(16.46)</td><td>12.53<br>(11.29)</td><td>10.74<br>(18.05)</td><td>11.75<br>(21.76)</td><td>13.48<br>(21.43)</td><td>13.82<br>(21.30)</td></tr><tr><th>LR<sub>c</sub> (%)</th><td>49.38<br>(17.24)</td><td>48.73<br>(13.65)</td><td>47.91<br>(11.60)</td><td>54.82<br>(15.04)</td><td>58.57<br>(19.94)</td><td>59.54<br>(22.16)</td><td>59.76<br>(22.23)</td></tr><tr><th>LR<sub>T</sub> (%)</th><td>30.30<br>(8.83)</td><td>34.11<br>(9.03)</td><td>36.09<br>(10.52)</td><td>44.35<br>(15.65)</td><td>50.82<br>(20.38)</td><td>53.51<br>(22.13)</td><td>54.18<br>(22.14)</td></tr></tbody></table>

<table><caption>**Table 5.** The mean values of aerosol vertical profile (AVP) retrieval errors estimated by GRASP for the synthetic test for a mixture of the Urban-Dust aerosol mixture. The mean values represent the distributions obtained using 300 realizations of the added random errors for the situation with total τ (440) = 1.0, with τ<sub>f</sub> = τ<sub>c</sub> = 0.5 and SZA = 75°. The shown relative errors for AVP (1 km) are represented by the mean values for three layers. Mean values for the actual errors are provided in parentheses.</caption><thead><tr><th></th><th></th><th>Up to 1.5 km</th><th>1.5-3.5 km</th><th>Above 3.5 km</th></tr></thead><tbody><tr><td rowspan="2">Urban-Dust</td><td>AVP<sub>f</sub> (%)</td><td>61.91<br>(25.79)</td><td>66.90<br>(27.73)</td><td>59.95<br>(22.06)</td></tr><tr><td>AVP<sub>c</sub> (%)</td><td>57.25<br>(14.43)</td><td>54.91<br>(4.11)</td><td>56.34<br>(23.51)</td></tr><tr><td rowspan="2">BB-Dust</td><td>AVP<sub>f</sub> (%)</td><td>61.46<br>(25.80)</td><td>68.50<br>(30.53)</td><td>59.17<br>(22.28)</td></tr><tr><td>AVP<sub>c</sub> (%)</td><td>55.74<br>(14.36)</td><td>53.90<br>(4.08)</td><td>55.30<br>(22.91)</td></tr></tbody></table>

all the measurements, since this is complete situation that is closest to the most practical situation.

Different tests were performed for this study. In particular, we focus on the case of BB-Dust, since Sect. 4.2.1 has already demonstrated the bias affects when each type of aerosol is observed separately. The effects of each bias sepa-

rately were analysed, while the corresponding illustrations are not shown since the results presented similar tendencies to those previously discussed and are observed for each type of aerosol separately in the last section. Figure 30 illustrates the examples of BB-Dust when τ (440) = 1.0, for different aerosol loads (τ<sub>f</sub> = 0.8 and τ<sub>c</sub> = 0.2, τ<sub>f</sub> = τ<sub>c</sub> = 0.5 and
<figure><img src="image_26.png" /><figcaption>**Figure 19.** Comparison of estimated and actual error distributions for AVP retrieved from measurements by simulated sun/sky photometer values and lidar for a mixture of Urban–Dust aerosols (a, b) and BB–Dust (c, d). The distributions were obtained using 300 realizations of added random errors. The mean values are represented by the black dot, and the median values of the errors are shown by a line in the box plot, along with the 25th–75th percentiles indicated by a box, and 5th–95th percentiles are indicated using whiskers. The mean values are represented by the black dot. The red colour shows the error estimates provided by GRASP, and the blue colour shows the calculated actual errors (Eq. 34).</figcaption></figure>

$\tau_f = 0.2$ and $\tau_c = 0.8$), assuming bias and random noise. The shaded areas represent the estimated total errors, as shown in Eq. (36). An important observation is that the error estimates for all retrieved and derived parameters characterize the actual errors well. As can be gained from the figure, the retrieval of the properties of the minor component appears to be the most challenging. As a matter of fact, the biggest errors in the retrieval are observed for the fine-mode properties, particularly in the case of $\tau_f = 0.2$ and $\tau_c = 0.8$. The largest discrepancies between estimated errors and the actual ones are observed in this situation. On the other hand, the properties of coarse mode are well represented in almost all cases, showing a good accuracy compared to the properties of the fine mode, even in the most challenging cases with the smallest presence of the coarse mode. This can probably be explained by the fact that desert dust AOD has rather moderate spectral changes.

On the other hand, in this section, some illustrations for the lidar ratio are also provided in order to demonstrate how the retrievals and the error estimates are affected by the bias in the measurements. Figure 30 illustrates the lidar ratios in this situation of mixed aerosol. The retrieval results and estima-

tion of the LR errors are rather satisfactory, with exceptions of low AOD cases, mainly in the case where the fine mode has only a very minor presence (of $\tau_f = 0.2$ and $\tau_c = 0.8$). Also, it should be emphasized that the errors estimated for total SSA and refractive index (RI) are rather adequate, while in Sect. 4.1, where only random errors were considered, the results showed some apparent underestimations.

### 4.2.3 Effects of measurement bias in the retrieval of mixed aerosol properties from observations of radiometers in combination with lidar

This section considers the same example as in Sect. 4.2.2 and analyses the effects of measurement biases in the synergy retrieval using coincident measurements from the sun/sky photometer and lidar measurements. At the same time, the results are presented for the most practical situation when both random and biased values are present in the measurements and accounted in the error estimates.

Figure 31 shows the results for the example of BB–Dust described in the previous section for $\tau(440) = 1.0$ and assumes the presence of both bias and random noise in all the
<figure><img src="image_27.png" /><figcaption>Figure 20. Aerosol properties retrieved from simulated sun/sky photometer data with assumed bias in AOD-simulated data for BB aerosol for $\tau(440) = 0.1, 0.3$ and $0.6$ (left to right). Retrievals after adding the positive bias +0.01 are represented in the block on the left and negative bias -0.01 in the block on the right. The case for $\tau(440) = 0.1$ is not represented on the right side, since the AOD at 1020 nm is ~0.01. The solid lines are the simulated properties (SD, RRI, IRI and SSA), and the dashed lines are the retrieved parameters. The shaded area indicates the systematic errors estimated by the GRASP algorithm.</figcaption></figure>

synthetic measurements for AOD, radiometers and lidar. As can be seen, the results of these different tests illustrate the positive influence of using radiometer and lidar synergy. The error estimates seem to be rather accurate too. For example, the most notable enhancement is in the lidar ratio accuracy, especially when the fine mode is the smallest, i.e. for the case with the following: $\tau_f = 0.2$ and $\tau_c = 0.8$. This behaviour was also seen by Lopatin et al. (2013), who explained that these were expected results, since the lidar ratio has a high sensitivity to the lidar signal. Nevertheless, we can see some improvements in the retrieval of complex refractive index using both lidar and photometer data.

With regard to the accuracy of the error estimation, in Sect. 4.1, we have illustrated the retrieval of error estimates for LR and showed some apparent underestimation when only random errors were considered. Figure 31 illustrates an important improvement in the estimation of the errors, once both random noise and bias are considered and Eq. (31) is used for accounting the effect of the systematic component.

Thus, using the synergy of both instruments can provide more accurate retrievals of LR, and the error can be estimated rather accurately using the developed methodology for both aerosol components, even for an aerosol mode with a lower presence. Figure 32 illustrates the retrievals of the vertical

aerosol profile in all three cases. The results show similar tendencies to Sect. 4.1.

## 4.3 Illustration and description of the correlation matrices

The values of non-diagonal elements of covariance provide important and interesting information about the retrieved parameters. For example, if the values $\rho_{ii'} \neq 0$ are close to $\pm 1$, the similitude of the influences of the parameters $a_i$ and $a_{i'}$ on the inverted measurements $f^*$ may explain the large variances in the retrieval error for these parameters. Also, knowledge about $\rho_{ii'} \neq 0$ is highly useful for the situation when several parameters from a set of simultaneously retrieved parameters $a_i$ need be jointly used in the applications. This can be easily seen from Eq. (29). For example, let us consider the estimates of two parameters, $a_1$ and $a_2$, which have errors, $\Delta a_1$ and $\Delta a_2$, characterized by covariance matrix, as follows:

$$C_{\Delta} = \begin{pmatrix} \sigma_1^2 & \sigma_1\sigma_2\rho_{12} \\ \sigma_1\sigma_2\rho_{12} & \sigma_2^2 \end{pmatrix}, \qquad (37)$$

where $a$ is a vector defined as $a = (a_1, a_2)^T$. Correspondingly, if in the application one needs to use the characteristic $m$ that is a liner function of $m = K_1a_1 + K_2a_2$, then the vari-
<figure><img src="image_28.png" /><figcaption>**Figure 21.** Aerosol properties retrieved from simulated sun/sky photometer data with assumed bias in AOD-simulated data for BB aerosol for $τ(440) = 0.1$, 0.3 and 0.6 (left to right). Retrievals after adding the positive bias +0.02 are represented in the block on the left and negative bias -0.02 in the block on the right. The case for $τ(440) = 0.1$ is not represented on the right side, since the AOD at 1020 nm is ~ 0.01. The solid lines are the simulated properties (SD, RRI, IRI and SSA), and the dashed lines are the retrieved parameters. The shaded area indicates the systematic errors estimated by the GRASP algorithm.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Figure 22.** Aerosol properties retrieved from simulated sun/sky photometer data with assumed bias in AOD-simulated data for dust aerosol for $τ(440) = 0.1$, 0.3 and 0.6 (left to right). Retrievals after adding the positive bias +0.01 are represented in the block on the left and negative bias -0.01 in the block on the right. The solid lines are the simulated properties (SD, RRI, IRI and SSA), and the dashed lines are the retrieved parameters. The shaded area indicates the systematic errors estimated by the GRASP algorithm.</figcaption></figure>
<figure><img src="image_30.png" /><figcaption>**Figure 23.** Aerosol properties retrieved from simulated sun/sky photometer data with assumed bias in AOD simulated data for dust aerosol for $τ(440) = 0.1, 0.3$ and $0.6$ (left to right). Retrievals after adding positive bias +0.02 are represented in the block on the left and negative bias -0.02 in the block on the right. The solid lines are the simulated properties (SD, RRI, IRI and SSA), and the dashed lines are the retrieved parameters. The shaded area indicates the systematic errors estimated by the GRASP algorithm.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>**Figure 24.** Aerosol properties retrieved from simulated sun/sky photometer data with assumed bias in AOD-simulated data for BB aerosol for $τ(440) = 0.1, 0.3$ and $0.6$ (left to right). Retrievals after adding the positive bias +0.03 are represented in the block on the left and negative bias -0.03 in the block on the right. The solid lines are the simulated properties (SD, RRI, IRI and SSA), and the dashed lines are the retrieved parameters. The shaded area indicates the systematic errors estimated by the GRASP algorithm.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>**Figure 25.** Aerosol properties retrieved from simulated sun/sky photometer data with assumed bias in AOD-simulated data for BB aerosol for $τ(440) = 0.1, 0.3$ and $0.6$ (left to right). Retrievals after adding the positive bias +0.05 are represented in the block on the left and negative bias -0.05 in the block on the right. The solid lines are the simulated properties (SD, RRI, IRI and SSA), and the dashed lines are the retrieved parameters. The shaded area indicates the systematic errors estimated by the GRASP algorithm.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>**Figure 26.** Aerosol properties retrieved from simulated sun/sky photometer data with assumed bias in AOD-simulated data for dust aerosol for $τ(440) = 0.1, 0.3$ and $0.6$ (left to right). Retrievals after adding the positive bias +0.03 are represented in the block on the left and negative bias -0.03 in the block on the right. The solid lines are the simulated properties (SD, RRI, IRI and SSA), and the dashed lines are the retrieved parameters. The shaded area indicates the systematic errors estimated by the GRASP algorithm.</figcaption></figure>
<figure><img src="image_34.png" /><figcaption>**Figure 27.** Aerosol properties retrieved from simulated sun/sky photometer data with assumed bias in AOD-simulated data for dust aerosol for $τ(440) = 0.1, 0.3$ and $0.6$ (left to right). Retrievals after adding the positive bias +0.05 are represented in the block on the left and negative bias -0.05 in the block on the right. The solid lines are the simulated properties (SD, RRI, IRI and SSA), and the dashed lines are the retrieved parameters. The shaded area indicates the systematic errors estimated by the GRASP algorithm.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>**Figure 28.** Aerosol properties retrieved from simulated sun/sky photometer data with assumed bias in AOD- and radiance-simulated data for BB (left) and dust (right) aerosol for $τ(440) = 0.1, 0.3$ and $0.6$ (left to right). Retrievals after adding positive bias +0.01 in AOD and +5% in radiances are represented in both cases. The solid lines are the simulated properties (SD, RRI, IRI and SSA), and the dashed lines are the retrieved parameters. The shaded area indicates the systematic error estimated by the GRASP algorithm.</figcaption></figure>
<figure><img src="image_36.png" /><figcaption>**Figure 29.** Aerosol properties retrieved from simulated sun/sky photometer data with assumed bias in AOD- and radiance-simulated data for BB (left) and dust (right) aerosol for $τ(440) = 0.1, 0.3$ and $0.6$ (left to right). Retrievals after adding positive bias +0.01 in AOD and +5% in radiances are represented in both cases. The solid lines are the simulated properties (SD, RRI, IRI and SSA), and the dashed lines are the retrieved parameters. The shaded area indicates the systematic error estimated by the GRASP algorithm.</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>**Figure 30.** Aerosol properties retrieved from simulated sun/sky photometer data with assumed random noise and bias in AOD- and radiance-simulated data for BB–Dust for $τ(440) = 1.0$. Retrievals after adding positive bias +0.01 in AOD and +5% in radiances are represented in both cases. The solid lines are the simulated properties (SD, RRI, IRI and SSA), and the dashed lines are the retrieved parameters. The shaded area indicates the systematic error estimated by the GRASP algorithm.</figcaption></figure>
<figure><img src="image_38.png" /><figcaption>**Figure 31.** Aerosol properties retrieved from simulated sun/sky photometer and lidar data with assumed random noise and bias in AOD-, radiance- and lidar-simulated data for BB-Dust for $\tau(440) = 1.0$. Retrievals after adding positive bias +0.01 in AOD, +5% in RAD and +0.2, +0.15 and +0.1 for each lidar wavelength at 355 nm, 532 nm and 1064 nm, respectively. The solid lines are the simulated properties (SD, RRI, IRI, SSA and LR), and the dashed lines are the retrieved parameters. The shaded area indicates the total errors estimated by the GRASP algorithm.</figcaption></figure>

$\tau_f = 0.8 \quad \tau_c = 0.2$

$\tau_f = \tau_c = 0.5$

$\tau_f = 0.2 \quad \tau_c = 0.8$

<figure><img src="image_39.png" /><figcaption>**Figure 32.** Aerosol properties retrieved from simulated sun/sky photometer and lidar data with assumed random noise and bias in AOD-, radiance- and lidar-simulated data for BB-Dust for $\tau(440) = 1.0$. Retrievals after adding positive bias +0.01 in AOD, +5% in RAD and +0.2, +0.15 and +0.1 for each lidar wavelength at 355 nm, 532 nm and 1064 nm, respectively. The solid lines are the simulated properties (SD, RRI, IRI, SSA and LR), and the dashed lines are the retrieved parameters. The shaded area indicates the total errors estimated by the GRASP algorithm.</figcaption></figure>
ance $\sigma_m^2$ can be obtained from Eq. (29) as follows:

$$ 
\begin{aligned}
\sigma_m^2 &= \mathbf{KC}_{\Delta a}\mathbf{K}^T \\
&= (K_1 \quad K_2) \begin{pmatrix} \sigma_1^2 & \sigma_1\sigma_2\rho_{12} \\ \sigma_1\sigma_2\rho_{12} & \sigma_2^2 \end{pmatrix} \begin{pmatrix} K_1 \\ K_2 \end{pmatrix} \\
&= K_1^2\sigma_1^2 + K_2^2\sigma_2^2 + 2K_1K_2\sigma_1\sigma_2\rho_{12}.
\end{aligned} 
\quad (38) $$

From this equation, the importance of the correlation coefficient $\rho_{12}$ is quite evident. Specifically, if $\rho_{12} = 0$, then the variance $\sigma_m^2$ is just a simple sum $K_1^2\sigma_1^2 + K_2^2\sigma_2^2$. Therefore, the error propagation from $\Delta a_i$ to $\Delta m$ is straightforward, and only the values of the sensitivities $K_i^2$ determine the contribution of $\Delta a_i$ (decreasing or increasing) to $\Delta m$.

When $\rho_{12} \neq 0$, the situation is more complex. However, the estimation of the main tendencies can be simplified in some cases. For instance, in the following:

$$ 
\begin{aligned}
& \text{if } \sigma_1^2 = \sigma_2^2 \text{ and } K_1 = K_2, \\
& \text{then } \sigma_m^2 = 2K_1^2\sigma_1^2(1 + \rho_{12}),
\end{aligned} 
\quad (39) $$

or

$$ 
\begin{aligned}
& \text{if } \sigma_1^2 = \sigma_2^2 \text{ and } K_1 = -K_2, \\
& \text{then } \sigma_m^2 = 2K_1^2\sigma_1^2(1 - \rho_{12}).
\end{aligned} 
\quad (40) $$

From these equations, it can be seen that if the correlation coefficient $\rho_{12} \to 1$ or $\rho_{12} \to (-1)$, then, depending on the case in Eqs. (39) and (40), $\sigma_m^2$ can be close to zero or up to $4K_1^2\sigma_1^2$. Therefore, knowledge about the non-zero correlation coefficients $\rho_{ii'} \neq 0$ is very important for understanding how the error is propagated to derived parameters obtained from the primary set a.

In practical cases, when the derived parameter m is a function of a large number of parameters $a_i$, the contributions to $\sigma_m^2$ increasingly become very complex with the increase in the number of involved parameters $a_i$. Therefore, unfortunately, the general qualitative analysis, similar to the one demonstrated by Eqs. (38)–(40), becomes very difficult and often practically impossible. Nonetheless, as will be shown below, the visualization of the correlation matrices in Eq. (33) can be very useful for an analysis of the retrieval tendencies.

### 4.3.1 Retrieval of single-aerosol components from radiometer measurements

Figure 33 shows the correlation matrices of random retrieval errors for BB (spherical particles) and dust (nonspherical particles) for a conventional AERONET-like inversion. The first 22 parameters ($22 \times 22$) represent the SD. This is followed by two blocks of $4 \times 4$. These two blocks are related to the RRI and IRI for four wavelengths (440 to 1020 nm). The last parameter is the sphericity fraction ($1 \times 1$). The colours represent the values of the correlation coefficients, where the red

colour denotes positive correlations, and the blue colour indicates negative correlations. The density of the colours indicates the values of the correlation coefficients changing from zero (the white colour) to dense red or blue colours corresponding to values of 1 and $-1$ accordingly.

The correlation for biomass burning case is shown in Fig. 33a. As can be seen, in general, size bin retrievals have a rather moderate correlation between them, though large positive correlations between the retrieval errors in neighbours can be observed. This is more evident for size bins at the smallest and largest particle sizes. This indicates that the size distribution values for those sizes have a tendency to be overestimated or underestimated together, which can mostly be explained by the use of the typical smoothness constraints imposed on the size distributions. The errors in RRI and IRI are negatively correlated with the SD parameters. The correlations seem especially pronounced between RRI and the values of the size distribution for the fine mode. Correspondingly, the overestimations of the size distribution values may tend to be accompanied by an underestimation of RRI, and vice versa. The errors in the fraction of spherical particle seem to show a positive correlation with the SD retrieval errors. This correlation is more evident when there is a fine-mode-dominated aerosol, since the scattering of fine-mode particles has a quite similar shape for spheres and spheroids. Therefore, when there is fine-mode domination it is more difficult to differentiate between spheres or spheroids. The positive sign of the correlation can be explained by the fact that extinction cross sections for the equivalent radii are a bit higher for spheroids. Thus, a higher percentage of spheres can be optically compensated by an increase in the volume concentration without a big impact in the total residual. The fraction of the spherical particle shows a negative strong correlation with the errors in the refractive index. Strong positive correlations can be seen between the spectral values of RRI. The positive correlations are present but lower between spectral values of IRI. As already noticed, this likely relates to the use of the rather strong smoothness constraint of the spectral variability in RRI and a weaker constraint of the spectral variability in IRI in the retrieval (Dubovik and King, 2000). The essential positive correlation can also be noticed between errors in RRI and IRI.

Figure 33b shows the correlation matrix for the retrieval of dust aerosol. The structure of the correlation for SD exhibits some differences compared to the BB case. Specifically, the positive correlations between the neighbouring size bins for the smallest and largest particle sizes are even more pronounced. Also, somewhat stronger negative correlations can be seen in the intermediate-sized ones. The strong negative correlation between RRI and SD retrieval errors remains only between concentrations of very small particles and the values of RRI at shortest wavelengths. The notable positive correlation is present only between the spectral values of RRI at the shortest wavelengths and between the spectral values of RRI at the longest wavelengths. At the same time, overall,
<figure><img src="image_40.png" /><figcaption>**Figure 33.** Correlation matrices of the estimated errors for aerosol retrieval from sun/sky radiometer observations (a) for biomass burning aerosols and (b) for desert dust using the GRASP algorithm. The values close to 1 or -1 mean stronger correlations between the properties, positive or negative, respectively.</figcaption></figure>

the correlation of retrieval errors in both RRI and IRI, between themselves and with other parameters, decreases compared to the case of BB. The errors in the fraction of the spherical particle for the dust case correlates much less with the errors in the other parameter compared to BB case. This can be explained by the fact that the light scattering of large particles is significantly more sensitive to the deviation in aerosol particles from spheres compared to spheroids than the light scattering of fine fraction particles (Dubovik et al., 2006). Therefore, when coarse particles dominate, the discrimination between spheres and spheroids becomes more evident.

Thus, the analysis of the correlation matrices itself provide very useful insight that helps to understand and interpret retrieval results. For example, such artefacts as the appearance of tails (unrealistically high concentrations) at extremes of the size distribution have been noticed and widely discussed (Dubovik et al., 2002a; Dubovik et al., 2006; Torres et al., 2014). Such retrieval artefacts as the underestimation of RRI accompanied by an overestimation in the size distribution of very fine particles has been widely discussed in studies by Dubovik et al. (2000); Dubovik et al. (2002b, a). These artefacts were strongly reduced by accounting for the particle nonsphericity of desert dust particles (Dubovik et al., 2002b; Dubovik et al., 2006), but nonetheless the less-pronounced appearance of such artefacts remains in AERONET-like retrievals (Torres et al., 2014, 2017). These artefacts are clearly related to the above-observed presence of the strong negative correlation between the values of RRI and the size distribution of very fine particles. It should be noted that the presence of high correlations is an indication that adding information

about one of the correlated parameters should improve the retrieval, not only for the constrained parameter itself but also for the parameters that strongly correlated with this parameter. For example, the addition of polarimetric observations to the traditional set of AERONET observations results in a clear improvement in the retrieval of RRI and the size distribution of very fine particles (Li et al., 2009; Fedarenka et al., 2016). Indeed, the degree of linear polarization is known to be very sensitive to the amount and especially the RRI of fine particles (Dubovik et al., 2006). This is why the addition of polarimetric observations helps to reduce the correlations between the errors in RRI and size distribution of fine particles and helps the overall improvement in the retrieval accuracy of these particles.

### 4.3.2 Retrieval of mixed aerosol components from radiometer measurements only

In this section, the correlation matrix for bicomponent aerosol retrieved from the synthetic observations of sun/sky radiometers of two aerosol mixtures is illustrated in Fig. 34 for BB-Dust and Urban-Dust. The structure of this matrix consists of 25 parameters related to the SD that are separated into two blocks, with 10 parameters for fine mode and 15 parameters for coarse mode. The following four blocks of 4 × 4 are related to the RRI and IRI of fine and coarse modes at four wavelengths. These blocks are followed by a single value of the sphericity fraction.

The area of the correlation matrix that contains SD, RRI, IRI and the sphericity fraction is quite similar to the correlation matrix obtained for aerosol AERONET-like retrieval. The main difference is the separation into two modes, since
<figure><img src="image_41.png" /><figcaption>**Figure 34.** Correlation matrices of the estimated errors for mixed aerosol retrieval from sun/sky radiometer observations (a) for the mixture of biomass burning and dust aerosols and (b) for the mixture of urban and dust aerosols using the GRASP algorithm. The values close to 1 or -1 mean stronger correlations between the properties, positive or negative, respectively.</figcaption></figure>

strong negative correlations can be observed between the corresponding parameters of the fine and coarse mode. For example, strong negative correlations can be observed between IRI fine and coarse mode. These correlations mean that overestimating the amount or absorption of one aerosol mode is likely compensated by underestimation of the amount or absorption of another aerosol mode. Another interesting anticorrelation can be observed for the last three bins of SD fine mode and the first three bins of SD coarse mode. Actually, both volume distributions have these three bins in common. This overlapping zone is never easy to properly separate for the code, but at the same time, it coincides with a local minimum value of most size distributions found in the real retrievals.

### 4.3.3 Retrieval of mixed aerosol properties from measurements of radiometers in combination with lidar

Figure 35 shows the correlation matrix for the case when the bicomponent aerosol retrieved from synthetic observations for the sun/sky photometer and lidar of aerosol mixtures is given. In the figure, the different blocks are identified. The first 25 parameters represent the SD that are separated into two blocks, with 10 parameters for fine mode and 15 parameters for coarse mode. The following four blocks of 7 × 7 are related to the RRI and IRI of fine and coarse modes at seven wavelengths. These blocks are followed by a single value of the sphericity fraction. The two last and largest blocks of

<figure><img src="image_42.png" /><figcaption>**Figure 35.** Correlation matrices of the estimated errors for aerosol retrieval from joint sun/sky radiometer and lidar observations for a mixture of urban and desert dust using the GRASP algorithm. The values close to 1 or -1 mean stronger correlations between the properties, positive or negative, respectively.</figcaption></figure>

60 × 60 parameters each correspond to the AVP values of the two modes given at 60 different altitudes.

The area of the correlation matrix that contains SD, RRI, IRI and the sphericity fraction is quite similar to the correlation matrix described in the previous section, considering
that the aerosol mixture retrievals are from AERONET-like
observations only.

As expected, the block of the correlations of AVP retrieval
shows strong negative correlations between the errors in the
retrieved parameters of fine and coarse mode. Thus, an over-
estimation of one mode is highly correlated with an under-
estimation of another AVP mode. Furthermore, strong posi-
tive correlations can be observed between AVP values corre-
sponding to the same fine or coarse mode, i.e. the AVP val-
ues of each mode have tendency to be simultaneously over-
estimated or underestimated. This is related to the limited
sensitivity of the used lidar data for distinguishing the contri-
butions of different modes and also to the use of smoothness
constraints on vertical variations in AVP of each fraction. On
the other hand, a near-zero correlation can be seen for AVP
parameters at altitudes with a significant presence of one or
both aerosol components. Correspondingly, there is a high
sensitivity of both lidar and radiometer observations to the
aerosol parameters at those altitudes. It should be noted that
all the above-discussed retrievals suggested from the anal-
ysis of the correlation matrices were actually observed in
the retrievals from real data, as discussed by Lopatin et al.
(2013, 2021).

## 4.4 Illustration of the GRASP error estimates with real observations

This section illustrates the GRASP error estimates perfor-
mance for the retrieval from real data. For that purpose, the
lidar and sun/sky photometer measurements collected at the
Aeroparque (34°33'51" S, 58°25'02" W) and Villa Martelli
(34°33'21" S, 58°30'23" W) stations in Buenos Aires, Ar-
gentina, have been used. These instruments are part of
the LALINET (Latin America Lidar Network; Guerrero-
Rascado et al., 2016) and AERONET networks. Both sites
are located in an industrialized city dominated by continen-
tal and urban/industrial aerosols and are affected during win-
ter and spring by biomass burning from the north and centre
of the country and neighbouring countries (mainly Brazil).
The Aeroparque station is located at the Jorge Newbery In-
ternational Airport within the city limits. This station does
not have a co-located sun photometer, but its location is 7 km
from the Villa Martelli station where the sun photometer
is installed. On the other hand, the Villa Martelli station is
found in the limits of Buenos Aires in a highly populated
and industrialized area.

Observations from two different biomass burning events in
Argentina were selected for the illustrations. Specifically, 3 d
were chosen with different aerosol loads and SZA. The lidar
range-corrected signals (RCSs) corresponding to each event
are shown in Fig. 36. They have been calculated from the
lidar signal, with background and dark current corrections
and multiplied by the height squared. In addition, the back-
 trajectories calculated from the HYSPLIT (Hybrid Single-
 Particle Lagrangian Integrated Trajectory; Stein et al., 2015;

Rolph et al., 2017) models are presented in order to confirm
where the air masses come from (Fig. 37).

The two first cases selected correspond to an important
event of biomass burning that occurred in the bordering
countries to the north of Argentina in August 2014, partic-
ularly in the south of Brazil and Paraguay. It was detected
in Buenos Aires between 19 and 23 August. For illustration
purposes, the measurements corresponding to 19 August are
used, which present a low aerosol load at 440 nm (~ 0.11)
and SZA > 50°. Figure 36a indicated the presence of sev-
eral layers of aerosols up to 1.5 km. The lidar measurements
on 19 August were taken between 12:15 and 12:45 UTC,
and the AERONET measurements correspond to 11:25 UTC.
Figure 37 shows the HYSPLIT back trajectories that vali-
date the source of the air masses. The measurements corre-
sponding to 22 August are shown in Fig. 36b, where sev-
eral layers of aerosol up to 3 km are observed. For this day,
the aerosol load increases (AOD at 440 nm ~ 0.31), and the
SZA is < 50°. The inversion was realized with the average li-
dar data between 15:00 and 15:20 UTC, and the AERONET
measurements were considered at 16:59 UTC. The satellite
image corresponding to 22 August (Fig. 37a) shows the pres-
ence of aerosols that extend from the north of Argentina to-
wards the centre, passing through the province of Buenos
Aires. Moreover, MODIS hot spots are detected in the satel-
lite image. The source of the air masses can be validated from
the HYSPLIT back trajectories (Fig. 37c). The last selected
case corresponds to the biomass burning event on 25 Septem-
ber 2017, which occurred in the north of Argentina and in
the bordering countries (Fig. 36c). In this work, lidar mea-
surements from the Aeroparque station between 19:20 and
20:10 UTC and the AERONET measurements correspond-
ing to 19:20 UTC, whose AOD value at 440 nm is 0.57 and
SZA > 50°, were used.

Figure 38 illustrates the retrieved columnar properties for
each day obtained by GRASP from a combination of ra-
diometer and lidar data and the comparison with the cor-
responding standard AERONET retrievals. The results pro-
vided by GRASP are represented in solid lines, with blue for
the fine mode and green the coarse mode. Shaded areas rep-
resent the error provided by GRASP for each retrieved and
derived property. Magnified plots represent the effective re-
fractive index and total SSA for GRASP (black solid line)
and AERONET (black dashed line). From the illustrations,
one can see that almost all the GRASP-retrieved properties in
the three cases present a good agreement with the AERONET
retrievals.

The error tendencies for SD that can be seen from Fig. 38
agree with those identified above in the present study and
with results presented in previous sections and also with the
results in some other works by Dubovik et al. (2000) and
Lopatin et al. (2013). For example, the retrieval errors clearly
increase at the extremes of SD. Moreover, one clear and
known tendency can be mentioned. The size distribution shift
towards higher radii in the three cases could be explained by
<figure><img src="image_43.png" /><figcaption>**Figure 36.** RCSs at 1064 nm in arbitrary units from Villa Martelli, Argentina, on 19 August 2014 (a) and 22 August 2014 (b) and from the Aeroparque station, Argentina, on 25 September 2017 (c). The two red lines indicate the analysed lidar data interval. The black solid line indicates the sun photometer measurements.</figcaption></figure>

the use of lidar data in the inversions that provide additional information at scattering angles of 180° (Lopatin et al., 2013; Bovchaliuk et al., 2016; Benavent-Oltra et al., 2017). As seen in Fig. 38, these deviations, in almost all cases, are included in the error (shaded areas).

The errors in RRI, IRI and SSA were retrieved for each mode separately by GRASP, and they are significantly higher than the error for RRI, IRI and SSA of the total components. The effective RRI and IRI and the total SSA obtained by GRASP are in the middle of the retrieved values for fine and coarse mode separately. The total values shown in the magnified plots agree well with the RRI, IRI and SSA provided by AERONET.

On the other hand, the case corresponding to 25 September has an AOD > 0.4 and SZA > 50°, allowing us to have the uncertainty in the SSA provided by AERONET. Thus, in Fig. 38, for this particular case, we can observe the comparison of the uncertainty in SSA from AERONET and the

SSA error provided by GRASP. Note that it is an advantage of GRASP that it provides the errors for each parameter in all the situations. Regarding the values of SSA and the tendencies of their variability, the results show that SSA is representative of biomass burning cases, namely that the values of SSA decreasing with the wavelength agree with AERONET climatology given by Dubovik et al. (2002a). As expected, based on the results of previous studies by Dubovik et al. (2000), Lopatin et al. (2013) and Tsekeri et al. (2017), the best agreements are obtained as the aerosol load increases. More specifically, we observe that the estimated errors in total SSA in the two first cases increase. This can be associated with unfavourable configurations of the observation. On 19 August the AOD at 440 nm is ~ 0.11 (SZA is > 50°) and on 22 August the AOD at 440 nm is ~ 0.31 (SZA is < 50°). Thus, the measurements in situations with a low amount of aerosol and with small SZA may not contain enough information to adequately retrieve the SSA (Dubovik et al., 2000;
<figure><img src="image_44.png" /><figcaption>**Figure 37.** (a) Satellite image with hot spots corresponding to 22 August 2014 (NASA Worldview), (b) air mass back-trajectories for the Villa Martelli measurement site on 19 August 2014 and (c) air mass back trajectories for the Villa Martelli measurement site on 22 August 2014.</figcaption></figure>

Lopatin et al., 2013; Torres et al., 2014). Therefore, the case of 25 September corresponds to the most favourable situation for realizing reliable aerosol retrieval, since the AOD at 440 nm value is > 0.4, and SZA is > 50°. Indeed, the GRASP and AERONET retrievals have the best agreement for this day.

Figure 39 shows the retrieved vertical distributions of fine and coarse modes. The vertical structure of the aerosols of different types is clearly discriminated and shows good agreements with the back-trajectory analysis for each day. Furthermore, the error estimates show good agreements with the previous results provided in the last sections for the simulated cases.

Moreover, Fig. 40 shows the retrieved LR and their error estimates using the GRASP algorithm. The magnified plots show the total LR provided by GRASP (black solid line) and by AERONET (black dashed line). The associated errors are represented in shaded areas for GRASP. As can be seen, it is only possible to compare the values with AERONET retrieval for Fig. 40, which corresponds to higher aerosol load and SZA > 50° on 25 September 2017.

Thus, the retrieved parameters and error estimates from the GRASP application to the real data and their comparisons to the AERONET retrieval results showed an encouraging agreement between the columnar properties of aerosol. At the same time, GRASP provide the error estimates for the retrieved properties in both the fine and coarse mode and also for the total components. Moreover, GRASP has also the advantage of providing the dynamic error estimates in all configurations. As seen above, AERONET error estimates are only provided in some particular situations when the AOD at 440 nm is > 0.4, and SZA > 50°.

## 5 Conclusions

In this work, we reviewed the approach realized in the GRASP algorithm for estimating errors in the parameters retrieved from remote sensing observations. The employed approach relies on the rigorously realized concept of statistical estimations and tends to account for the propagation of both random and systematic errors. Then we evaluated the performance of the GRASP error estimates for aerosol parameters retrieved from ground-based observations. We considered AERONET-like retrievals from observations by sun/sky-scanning radiometers and GRASP synergy aerosol retrievals from joint observations by radiometers and lidar. GRASP generates the full covariance matrices that are expected to be used for generating error bars for retrieved parameters and provides an interesting insight for understanding retrieval tendencies. Therefore, we studied the quantitative reliability of the obtained covariance diagonal elements and analysed the structure of correlation coefficient of covariance matrices.

The performance of the GRASP estimates of error variances in the presence of random errors was evaluated in a series of numerical tests and illustrated the capabilities of the GRASP algorithm to provide rigorous estimates of dynamic retrieval errors. In the frame of these tests, the synthetic proxy observations perturbed by 300 random-noise-generated realizations were inverted using the GRASP algorithm. Then, the retrieved parameters were compared to those used for the generation of the synthetic data, and the obtained error estimates were compared with actual deviations in the retrieved parameters from assumed values. This analysis was realized for the synthetic observations for three different
<figure><img src="image_45.png" /><figcaption>**Figure 38.** Comparison of columnar properties retrieved by GRASP from a combination of sun/sky photometer and lidar data and data retrieved conventionally by AERONET. SD, RRI, IRI and SSA retrieved by GRASP are shown in solid lines of blue (fine mode) and green (coarse mode). The shaded area in colours of blue and green represents the total error provided by GRASP, and black shaded areas are the uncertainties provided by AERONET. Magnified panels show the RI effective and total SSA provided by GRASP (black solid line) and AERONET (black dashed line). Their associated errors are represented in the grey shaded area for GRASP and with error bars for AERONET.</figcaption></figure>

types of aerosols and for the mixture of them. Observations of dust were modelled using AERONET retrieval climatology at the Solar Village (Riyadh, Saudi Arabia) site. The AERONET retrieval climatologies from the African savanna (Zambia) and the GSFC (Maryland, USA) were used to simulate urban and BB aerosol observations, respectively. The Urban–Dust aerosols and BB–Dust mixture were considered for modelling the properties of mixed aerosols. For each observed aerosol type or mixture thereof, different aerosol loads were tested. First, we modelled observations of aerosols of only one type aerosol for $τ(440) = 0.3, 0.6$ and $0.9$. For aerosol mixtures, we also considered scenarios with different aerosol loads, while we present most of the illustrations of the situations with a large aerosol load ($τ(440) = 1.0$, combining different aerosol loads for each mode $τ_f = 0.2$ and $τ_c = 0.8$, $τ_f = τ_c = 0.5$ and $τ_f = 0.8$, and $τ_c = 0.2$). The

data were perturbed by random noise before applying the retrieval algorithm. For all these simulations, the SZA at 75° was used.

The tests evaluated the situations when only radiometer data were inverted, and then radiometer data were inverted jointly with coincident lidar data. Two GRASP retrieval set-ups were tested, where (i) the retrieval assumes that aerosol is composed of homogeneous particles, and parameters of only one aerosol component are retrieved, and (ii) then the aerosol is assumed to be an external mixture of two aerosol components, and the parameters of each component are retrieved separately. In the case when the lidar data were used, the vertical profiles of the concentration were also retrieved for each aerosol component. The illustrations for aerosol retrieval from all simulated data sets were provided, using both approaches. However, the full statistical analysis
<figure><img src="image_46.png" /><figcaption>**Figure 39.** Retrieved aerosol vertical profiles (AVPs) by GRASP from a combination of sun/sky photometer and lidar data. The blue solid line represents the AVP fine mode and green the AVP coarse mode. The shaded areas correspond to the total error provided by GRASP.</figcaption></figure>

<figure><img src="image_47.png" /><figcaption>**Figure 40.** Retrieved LR by GRASP from a combination of sun/sky photometer and lidar data. The blue solid line represents the LR fine mode and green the LR coarse mode. Magnified plots show, with black solid lines, the LR provided by GRASP and, with dashed lines, the LR provided by AERONET. Their GRASP-associated errors are represented in shaded areas.</figcaption></figure>

provided only the two conventional retrieval scenarios, i.e. AERONET-like, single-component retrieval from radiometric observations and GRASP bicomponent aerosol retrieval from combinations of radiometer and lidar data.

The results of the tests showed that the complete set of aerosol parameters for each aerosol component can be robustly derived with acceptable accuracy in almost all considered situations. The retrieval of bicomponent aerosol was evaluated using radiometer-only simulated measurements and then adding lidar observations. These tests allowed us to observe that, by using the synergy of two instruments, there are some improvements in the retrieval of the aerosol properties of each component of the observed aerosol mixture and in the estimations of the retrieval errors. The test for selected cases with different presences of different aerosol components ($\tau_f = 0.2$, $\tau_c = 0.8$ and $\tau_f = 0.8$, $\tau_c = 0.2$) showed that optical properties of the dominant mode can be retrieved in a significantly more accurate manner than expected. It is interesting to note that in all situations that use only radiometer data or add lidar-simulated measurements, properties such as total SSA and effective refractive index can be retrieved rather accurately, even in cases where the retrieval of the properties of each mode separately is questionable.

The results of the statistical tests with randomly generated noise showed that the GRASP error estimates, in most cases, are comparable to or exceed the actual errors by 20% to 30% and therefore can be safely used for assuring uncertainties in the actual retrieval products. In addition, the observation of typical error values was summarized for different situations and retrieval scenarios. Namely, the study confirmed that the detailed properties of aerosol mixtures can be rather reliably retrieved from a combination of radiometer and lidar data, provided that there is a sufficient amount of both aerosol components. For example, for the case when the total $\tau(440) = 1.0$ with a comparable presence of both components $\tau_f(440) = \tau_c(440) = 0.5$, and SZA is 75°, the mean values for RRI errors are ~0.05 for BB and urban and vary between 0.07 to 0.09 for dust, and IRI errors are around 60% for BB, 75% for urban and 50% for dust. SSA errors vary between 0.024 and 0.061 for urban, 0.042 and 0.086 for BB and 0.04 and 0.009 for dust, showing a clear tendency to increase with the wavelengths for values of fine mode and decrease for coarse mode. However, even for this case, the separation of the LR values for both modes showed high uncertainties at short wavelengths, in particular for the fine mode, while the values of the total LR errors were found to be real-
sonable in the range of 20% to 55%. The relative error estimates for AVP for both aerosol mixture (Urban–Dust and BB–Dust) cases were ranging for fine mode between 50% and 70% and for coarse dust mode between 50% and 55%.

The effects of the systematic errors in the retrievals were also analysed in a series of limited dedicated numerical tests. The results of the tests were used for adjusting the GRASP retrieval estimates to the potential effects of the systematic errors. The results show enhancements of total error estimates from the assumption of bias in the equation of systematic components.

In addition to the evaluation of the error bar estimates and the effects of systematic errors, in this paper we illustrated and discussed the correlation structures of the error covariance matrices for all main considered retrieval scenarios. The results showed that analysis of the correlation structure can be very useful for understanding the observed retrieval tendencies and optimizing retrieval. For example, for conventional AERONET-like aerosol retrievals from radiometer data only, the strong negative correlation between errors is in the real part of the refractive index and size distribution values for small sizes. This agrees well with the tendency commonly observed in actual retrieval when the underestimations of the real part are coincident with the overestimation of the fine-mode size distribution. Also, the presence of a high positive correlation between the errors in the size distribution for extreme sizes and between the errors in the refractive index at different wavelengths agrees well with the known possibilities of possible overestimations of aerosol concentrations for very small or very large particles and joint overestimations/underestimations of the refractive value at different wavelengths. For bicomponent retrievals, strong negative correlations can be observed between nearly all corresponding parameters of the fine and coarse mode. This means, for example, that the overestimation of the amount or absorption of one aerosol mode is likely compensated by underestimation of the amount or absorption of another aerosol mode. The decrease in some of these correlations was observed when inverted radiometer data were inverted simultaneously with the lidar data. The high positive correlations were seen for the errors in the vertical profile of the fine and coarse concentrations, with the exception of the values for the altitudes where one or both of the aerosol modes had substantial loads. These and other less obvious, but quite interesting, correlation structures and tendencies can be identified using the analysis of the correlation matrix structure. Thus, the availability and analysis of not only the error variances but also the correlation patterns appear to be a useful and promising approach for optimizing observation schemes and retrieval set-ups.

Finally, the utilization of GRASP for deriving detailed aerosol properties and estimations of their errors was demonstrated for the coincident lidar and sun photometer observations from Buenos Aires, Argentina. The GRASP retrievals and the error estimates of the columnar aerosol properties

were shown to be fully adequate in a comparative analysis with the aerosol products available from AERONET operational retrievals. The retrieval of the vertical profiles of fine and coarse aerosol modes showed results consistent with the expectation and the predictions of back-trajectory analysis.

Thus, the results presented in this work show promising potential for the utilization of GRASP-retrieved dynamic error estimates for the detailed retrieved aerosol parameters from measurements of ground-based radiometers and lidars, considering different geometries and the presence of diverse aerosol loads. These studies are expected to be completed in future by a more extensive analysis of the error estimates for such detailed parameters as vertical profiles of SSA and LR.

**Code availability.** GRASP is an open-source software that is available, upon registration, from https://access-request.grasp-cloud.com/service/gitlab (last access: 14 August 2022; GRASP Open repository, 2022).

**Data availability.** The data are available upon request to the corresponding author. AERONET data used for this study can be downloaded from NASA website https://aeronet.gsfc.nasa.gov/ (last access: 14 August 2022; NASA, 2022). LALINET data are available upon request.

**Author contributions.** OD contributed to the development of the overall algorithm methodology, research planning and article writing. BT contributed to the result discussions and the article writing. BT, TL, DF, AL, PL, CC and MEH contributed to the algorithm designs, development and software support. PL, CC, JABO and JLB contributed to the result discussions and polishing the article text. PR contributed to the data preparation and analysis. MEH prepared the paper, including the co-author contributions, and performed the main part of the presented research.

**Competing interests.** At least one of the (co-)authors is a member of the editorial board of *Atmospheric Measurement Techniques*. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.

**Disclaimer.** Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Acknowledgements.** The authors acknowledge the support from the H2020 Marie Skłodowska-Curie RISE Actions (grant no. 778349), and the authors acknowledge the support from the European Metrology Program for Innovation and Research (EMPIR) within the joint research project of EMPIR 19ENV04 MAPP "Metrology for aerosol optical properties". The EMPIR is jointly funded by the EMPIR participating countries within EURAMET and the Euro-
pean Union. Tatyana Lapyonok has been supported by the Chem-
ical and Physical Properties of the Atmosphere Project funded
by the French National Research Agency through the Programme
d'Investissement d'Avenir (contract no. ANR-11-LABX-0 0 05-
01), the Regional Council, Hauts-de-France, the European Funds
for Regional Economic Development (grant no. URF/1/2180-01-
01), and the Combined Radiative and Air Quality Effects of An-
thropogenic project. The authors also acknowledge AERONET and
LALINET, for the scientific and technical support. We acknowledge
the use of imagery from the NASA Worldview application (https:
//worldview.earthdata.nasa.gov, last access: 22 July 2022), part of
the NASA Earth Observing System Data and Information System
(EOSDIS). The authors also would like to acknowledge the NOAA
Air Resources Laboratory (ARL), for the provision of the HYS-
PLIT transport and dispersion model (https://www ready.noaa.gov,
last access: 22 July 2022) used in this publication.

*Financial support.* This research has been supported by the Eu-
ropean Metrology Programme for Innovation and Research (grant
no. 19ENV04 MAPP) and the French National Research Agency
(grant no. ANR-II-LABX-0005-01).

*Review statement.* This paper was edited by Daniel Perez-Ramirez and reviewed by Feng Xu and one anonymous referee.

References

Benavent-Oltra, J. A., Román, R., Granados-Muñoz, M. J., Pérez-Ramírez, D., Ortiz-Amezcua, P., Denjean, C., Lopatin, A., Lyamani, H., Torres, B., Guerrero-Rascado, J. L., Fuertes, D., Dubovik, O., Chaikovsky, A., Olmo, F. J., Mallet, M., and Alados-Arboledas, L.: Comparative assessment of GRASP algorithm for a dust event over Granada (Spain) during ChArMEx-ADRIMED 2013 campaign, Atmos. Meas. Tech., 10, 4439–4457, https://doi.org/10.5194/amt-10-4439-2017, 2017.

Benavent-Oltra, J. A., Román, R., Casquero-Vera, J. A., Pérez-Ramírez, D., Lyamani, H., Ortiz-Amezcua, P., Bedoya-Velásquez, A. E., de Arruda Moreira, G., Barreto, Á., Lopatin, A., Fuertes, D., Herrera, M., Torres, B., Dubovik, O., Guerrero-Rascado, J. L., Goloub, P., Olmo-Reyes, F. J., and Alados-Arboledas, L.: Different strategies to retrieve aerosol properties at night-time with the GRASP algorithm, Atmos. Chem. Phys., 19, 14149–14171, https://doi.org/10.5194/acp-19-14149-2019, 2019.

Benavent-Oltra, J. A., Casquero-Vera, J. A., Román, R., Lyamani, H., Pérez-Ramírez, D., Granados-Muñoz, M. J., Herrera, M., Cazorla, A., Titos, G., Ortiz-Amezcua, P., Bedoya-Velásquez, A. E., de Arruda Moreira, G., Pérez, N., Alastuey, A., Dubovik, O., Guerrero-Rascado, J. L., Olmo-Reyes, F. J., and Alados-Arboledas, L.: Overview of the SLOPE I and II campaigns: aerosol properties retrieved with lidar and sun–sky photometer measurements, Atmos. Chem. Phys., 21, 9269–9287, https://doi.org/10.5194/acp-21-9269-2021, 2021.

Boucher, O., Randall, D., Artaxo, P., Bretherton, C., Feingold, G., Forster, P., Kerminen, V.-M., Kondo, Y., Liao, H., Lohmann, U., Rasch, P., Satheesh, S., Sherwood, S., Stevens, B., Zhang, X.,

Qin, D., Plattner, G., Tignor, M., Allen, S., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P.: Clouds and aerosols, in: Climate change 2013: the physical science basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, 571–657, ISBN 9781107661820, 2013.

Bovchaliuk, V., Goloub, P., Podvin, T., Veselovskii, I., Tanre, D., Chaikovsky, A., Dubovik, O., Mortier, A., Lopatin, A., Korenskiy, M., and Victori, S.: Comparison of aerosol properties retrieved using GARRLiC, LIRIC, and Raman algorithms applied to multi-wavelength lidar and sun/sky-photometer data, Atmos. Meas. Tech., 9, 3391–3405, https://doi.org/10.5194/amt-9-3391-2016, 2016.

Chaikovsky, A., Dubovik, O., Holben, B., Bril, A., Goloub, P., Tanré, D., Pappalardo, G., Wandinger, U., Chaikovskaya, L., Denisov, S., Grudo, J., Lopatin, A., Karol, Y., Lapyonok, T., Amiridis, V., Ansmann, A., Apituley, A., Allados-Arboledas, L., Binietoglou, I., Boselli, A., D'Amico, G., Freudenthaler, V., Giles, D., Granados-Muñoz, M. J., Kokkalis, P., Nicolae, D., Oshchepkov, S., Papayannis, A., Perrone, M. R., Pietruczuk, A.,rocadenbosch, F., Sicard, M., Slutsker, I., Talianu, C., De Tomasi, F., Tsekeri, A., Wagner, J., and Wang, X.: Lidar-Radiometer Inversion Code (LIRIC) for the retrieval of vertical aerosol properties from combined lidar/radiometer data: development and distribution in EARLINET, Atmos. Meas. Tech., 9, 1181–1205, https://doi.org/10.5194/amt-9-1181-2016, 2016.

Chen, C., Dubovik, O., Henze, D. K., Lapyonok, T., Chin, M., Ducos, F., Litvinov, P., Huang, X., and Li, L.: Retrieval of desert dust and carbonaceous aerosol emissions over Africa from POLDER/PARASOL products generated by the GRASP algorithm, Atmos. Chem. Phys., 18, 12551–12580, https://doi.org/10.5194/acp-18-12551-2018, 2018.

Chen, C., Dubovik, O., Henze, D. K., Chin, M., Lapyonok, T., Schuster, G. L., Ducos, F., Fuertes, D., Litvinov, P., Li, L., Lopatin, A., Hu, Q., and Torres, B.: Constraining global aerosol emissions using POLDER/PARASOL satellite remote sensing observations, Atmos. Chem. Phys., 19, 14585–14606, https://doi.org/10.5194/acp-19-14585-2019, 2019.

Chen, C., Dubovik, O., Fuertes, D., Litvinov, P., Lapyonok, T., Lopatin, A., Ducos, F., Derimian, Y., Herman, M., Tanré, D., Rémer, L. A., Lyapustin, A., Sayer, A. M., Levy, R. C., Hsu, N. C., Descloitres, J., Li, L., Torres, B., Karol, Y., Herrera, M., Herreras, M., Aspetsberger, M., Wanzenboeck, M., Bindreiter, L., Marth, D., Hangler, A., and Federspiel, C.: Validation of GRASP algorithm product from POLDER/PARASOL data and assessment of multi-angular polarimetry potential for aerosol monitoring, Earth Syst. Sci. Data, 12, 3573–3620, https://doi.org/10.5194/essd-12-3573-2020, 2020.

Dubovik, O.: Optimization of numerical inversion in photopolarimetric remote sensing, in: Photopolarimetry in remote sensing, Springer, 65–106, https://doi.org/10.1007/1-4020-2368-5_3, 2004.

Dubovik, O. and King, M. D.: A flexible inversion algorithm for retrieval of aerosol optical properties from Sun and sky radiance measurements, J. Geophys. Res.-Atmos, 105, 20673–20696, https://doi.org/10.1029/2000jd900282, 2000.

Dubovik, O., Lapyonok, T., and Oshchepkov, S.: Improved technique for data inversion: optical sizing of
multicomponent aerosols, Appl. Optics, 34, 8422-8436,
https://doi.org/10.1364/AO.34.008422, 1995.

Dubovik, O., Smirnov, A., Holben, B., King, M., Kaufman, Y., Eck, T., and Slutsker, I.: Accuracy assessments of aerosol optical properties retrieved from Aerosol Robotic Network (AERONET) Sun and sky radiance measurements, J. Geophys. Res.-Atmos., 105, 9791-9806, https://doi.org/10.1029/2000jd900040, 2000.

Dubovik, O., Holben, B., Eck, T. F., Smirnov, A., Kaufman, Y. J., King, M. D., Tanré, D., and Slutsker, I.: Variability of absorption and optical properties of key aerosol types observed in worldwide locations, J. Atmos. Sci., 59, 590-608, https://doi.org/10.1175/1520-0469(2002)059<0590:voaaop>2.0.co;2, 2002a.

Dubovik, O., Holben, B., Lapyonok, T., Sinyuk, A., Mishchenko, M., Yang, P., and Slutsker, I.: Non-spherical aerosol retrieval method employing light scattering by spheroids, Geophys. Res. Lett., 29, 541-544, https://doi.org/10.1029/2001GL014506, 2002b.

Dubovik, O., Sinyuk, A., Lapyonok, T., Holben, B. N., Mishchenko, M., Yang, P., Eck, T. F., Volten, H., Munoz, O., Veihelmann, B., Van der Zande, W. J., Leon, J., Sorokin, M., and Slutsker, I.: Application of spheroid models to account for aerosol particle nonsphericity in remote sensing of desert dust, J. Geophys. Res.-Atmos., 111, D11208, https://doi.org/10.1029/2005JD006619, 2006.

Dubovik, O., Lapyonok, T., Kaufman, Y. J., Chin, M., Ginoux, P., Kahn, R. A., and Sinyuk, A.: Retrieving global aerosol sources from satellites using inverse modeling, Atmos. Chem. Phys., 8, 209-250, https://doi.org/10.5194/acp-8-209-2008, 2008.

Dubovik, O., Herman, M., Holdak, A., Lapyonok, T., Tanré, D., Deuzé, J. L., Ducos, F., Sinyuk, A., and Lopatin, A.: Statistically optimized inversion algorithm for enhanced retrieval of aerosol properties from spectral multi-angle polarimetric satellite observations, Atmos. Meas. Tech., 4, 975-1018, https://doi.org/10.5194/amt-4-975-2011, 2011.

Dubovik, O., Lapyonok, T., Litvinov, P., Herman, M., Fuertes, D., Ducos, F., Lopatin, A., Chaikovsky, A., Torres, B., Derimian, Y., Huang, X., Aspetsberger, M., and Federspiel, C.: GRASP: a versatile algorithm for characterizing the atmosphere, SPIE Newsroom, 25, https://doi.org/10.1117/2.1201408.005558, 2014.

Dubovik, O., Li, Z., Mishchenko, M. I., Tanre, D., Karol, Y., Bojkov, B., Cairns, B., Diner, D. J., Espinosa, W. R., Goloub, P., Gu, X., Hasekamp, O., Hong, J., Hou, W., Knobelspiesse, K. D., Landgraf, J., Li, L., Litvinov, P., Liu, Y., Lopatin, A., Marbach, T., Maring, H., Martins, V., Meijer, Y., Milinevsky, G., Mukai, S., Parol, F., Qiao, Y., Remer, L., Rietjens, J., Sano, I., Stammes, P., Stamnes, S., Sun, X., Tabary, P., Travis, L. D., Waquet, F., Xu, F., Yan, C., and Yin, D.: Polarimetric remote sensing of atmospheric aerosols: Instruments, methodologies, results, and perspectives, J. Quant. Spectrosc. Ra., 224, 474-511, https://doi.org/10.1016/j.jqsrt.2018.11.024, 2019.

Dubovik, O., Fuertes, D., Litvinov, P., Lopatin, A., Lapyonok, T., Doubovik, I., Xu, F., Ducos, F., Chen, C., Torres, B., Derimian, Y., Li, L., Herreras-Giralda, M., Herrera, M., Karol, Y., Matar, C., Schuster, G. L., Espinosa, R., Puthukkudy, A., Li, Z., Fischer, J., Preusker, R., Cuesta, J., Kreuter, A., Cede, A., Aspetsberger, M., Marth, D., Bindreiter, L., Hangler, A., Lanzinger, V., Holter, C., and Federspiel, C.: A Comprehensive Description of Multi-Term LSM for Applying Multiple a Priori Con-

straints in Problems of Atmospheric Remote Sensing: GRASP Algorithm, Concept, and Applications, Front. Remote Sens., 2, 706851, https://doi.org/10.3389/frsen.2021.706851, 2021.

Eck, T. F., Holben, B., Reid, J., Dubovik, O., Smirnov, A., O'neill, N., Slutsker, I., and Kinne, S.: Wavelength dependence of the optical depth of biomass burning, urban, and desert dust aerosols, J. Geophys. Res.-Atmos., 104, 31333-31349, 1999.

Edie, W. T., Drijard, D., James, F. E., Roos, M., and Sadoulet, B.: Statistical methods in experimental physics, New York North-Holland, 155 pp., ISBN 9780720402391, 1971.

Espinosa, W. R., Remer, L. A., Dubovik, O., Ziemba, L., Beyersdorf, A., Orozco, D., Schuster, G., Lapyonok, T., Fuertes, D., and Martins, J. V.: Retrievals of aerosol optical and microphysical properties from Imaging Polar Nephometer scattering measurements, Atmos. Meas. Tech., 10, 811-824, https://doi.org/10.5194/amt-10-811-2017, 2017.

Espinosa, W. R., Martins, J. V., Remer, L. A., Dubovik, O., Lapyonok, T., Fuertes, D., Puthukkudy, A., Orozco, D., Ziemba, L., Thornhill, K. L., Levy, R.: Retrievals of aerosol size distribution, spherical fraction, and complex refractive index from airborne in situ angular light scattering and absorption measurements, J. Geophys. Res.-Atmos., 124, 7997-8024, https://doi.org/10.1029/2018JD030009, 2019.

Fedarenka, A., Dubovik, O., Goloub, P., Li, Z., Lapyonok, T., Litvinov, P., Barel, L., Gonzalez, L., Podvin, T., and Crozel, D.: Utilization of AERONET polarimetric measurements for improving retrieval of aerosol microphysics: GSFC, Beijing and Dakar data analysis, J. Quant. Spectrosc. Ra., 179, 72-97, https://doi.org/10.1016/j.jqsrt.2016.03.021, 2016.

Fourgeaud, C. and Fuchs, A.: Statistique, 2nd edition, Dunod, Paris, 1967.

Giles, D. M., Holben, B. N., Eck, T. F., Sinyuk, A., Smirnov, A., Slutsker, I., Dickerson, R. R., Thompson, A. M., and Schafer, J. S.: An analysis of AERONET aerosol absorption properties and classifications representative of aerosol source regions, J. Geophys. Res., 117, D17203, https://doi.org/10.1029/2012JD018127, 2012.

Giles, D. M., Sinyuk, A., Sorokin, M. G., Schafer, J. S., Smirnov, A., Slutsker, I., Eck, T. F., Holben, B. N., Lewis, J. R., Campbell, J. R., Welton, E. J., Korkin, S. V., and Lyapustin, A. I.: Advances in the Aerosol Robotic Network (AERONET) Version 3 database - automated near-real-time quality control algorithm with improved cloud screening for Sun photometer aerosol optical depth (AOD) measurements, Atmos. Meas. Tech., 12, 169-209, https://doi.org/10.5194/amt-12-169-2019, 2019.

Granados-Muñoz, M., Guerrero-Rascado, J., Bravo-Aranda, J., Navas-Guzmán, F., Valenzuela, A., Lyamani, H., Chaikovsky, A., Wandinger, U., Ansmann, A., Dubovik, O., Grudo, J. O., and Alados-Arboledas, L.: Retrieving aerosol microphysical properties by Lidar-Radiometer Inversion Code (LIRIC) for different aerosol types, J. Geophys. Res.-Atmos., 119, 4836-4858, https://doi.org/10.1002/2013JD021116, 2014.

Granados-Muñoz, M. J., Bravo-Aranda, J. A., Baumgardner, D., Guerrero-Rascado, J. L., Pérez-Ramírez, D., Navas-Guzmán, F., Veselovskii, I., Lyamani, H., Valenzuela, A., Olmo, F. J., Titos, G., Andrey, J., Chaikovsky, A., Dubovik, O., Gil-Ojeda, M., and Alados-Arboledas, L.: A comparative study of aerosol microphysical properties retrieved from ground-based remote sensing and aircraft in situ measurements during a Saharan dust event,