# A Comprehensive Description of Multi-Term LSM for Applying Multiple a Priori Constraints in Problems of Atmospheric Remote Sensing: GRASP Algorithm, Concept, and Applications

O. Dubovik, David Fuertes, Pavel Litvinov, Anton Lopatin, Tatyana Lapyonok, Ivan Doubovik, Feng Xu, Fabrice Ducos, Cheng Chen, Benjamin Torres, et al.

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

## To cite this version:

O. Dubovik, David Fuertes, Pavel Litvinov, Anton Lopatin, Tatyana Lapyonok, et al.. A Comprehensive Description of Multi-Term LSM for Applying Multiple a Priori Constraints in Problems of Atmospheric Remote Sensing: GRASP Algorithm, Concept, and Applications. Frontiers in Remote Sensing, 2021, 2, pp.706851. 10.3389/frsen.2021.706851. insu-03435399

HAL Id: insu-03435399

https://insu.hal.science/insu-03435399

Submitted on 19 Nov 2021

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Check for
updates

# A Comprehensive Description of Multi-Term LSM for Applying Multiple a Priori Constraints in Problems of Atmospheric Remote Sensing: GRASP Algorithm, Concept, and Applications

Oleg Dubovik<fnref n="1" />*, David Fuertes², Pavel Litvinov², Anton Lopatin², Tatyana Lapyonok<fnref n="1" />, Ivan Doubovik<fnref n="1" />,², Feng Xu³, Fabrice Ducos<fnref n="1" />, Cheng Chen², Benjamin Torres<fnref n="1" />, Yevgeny Derimian<fnref n="1" />, Lei Li⁴, Marcos Herreras-Giralda², Milagros Herrera<fnref n="1" />, Yana Karol², Christian Matar², Gregory L. Schuster⁵, Reed Espinosa⁶, Anin Puthukkudy⁷, Zhengqiang Li⁸, Juergen Fischer⁹, Rene Preusker⁹, Juan Cuesta¹⁰, Axel Kreuter¹¹, Alexander Cede⁶,¹¹, Michael Aspetsberger¹², Daniel Marth¹², Lukas Bindreiter¹², Andreas Hangler¹², Verena Lanzinger¹², Christoph Holter¹² and Christian Federspiel¹²

**OPEN ACCESS**

**Edited by:**

Yingying Ma,
Wuhan University, China

**Reviewed by:**

Evgueni Kassianov,
Pacific Northwest National Laboratory
(DOE), United States
Alexander Kokhanovsky,
Telespazio Belgium, Germany

*Correspondence:

Oleg Dubovik
oleg.dubovik@univ-lille.fr

**Specialty section:**
This article was submitted to *Satellite Missions*, a section of the journal *Frontiers in Remote Sensing*

Received: 08 May 2021

Accepted: 09 August 2021

Published: 19 October 2021

**Citation:**
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 GL, 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 (2021) A Comprehensive Description of Multi-Term LSM for Applying Multiple a Priori Constraints in Problems of Atmospheric Remote Sensing: GRASP Algorithm, Concept, and Applications. *Front. Remote Sens.* 2:706851
doi: 10.3389/frsen.2021.706851

Advanced inversion Multi-term approach utilizing multiple a priori constraints is proposed. The approach is used as a base for the first unified algorithm GRASP that is applicable to diverse remote sensing observations and retrieving a variety of atmospheric properties. The utilization of GRASP for diverse remote sensing observations is demonstrated.

**Keywords:** inversion method, multiple a priori constraints, atmospheric remote sensing, GRASP algorithm, atmospheric aerosol

We describe an approach called the Multi-term Least Square Method (LSM) that has been used to develop complex aerosol inversion algorithms for a number of years and applied to retrievals of laboratory and ground-based measurements. Theoretically, it was shown how to unite the advantages of a variety of approaches and to provide transparency and flexibility in development of practically efficient retrievals. From a practical viewpoint, this approach provides a methodology for using multiple a priori constraints to atmospheric problems where rather different groups of parameters should be retrieved simultaneously. For example, Dubovik and King (J. Geophys. Res., 2000, 105, 673–696) used multi-term LSM for designing the algorithm that retrieves aerosol size distribution and spectrally dependent complex index of refraction from Sun/sky-radiometer ground-based observations. Furthermore, the significant potential of the multi-term LSM approach was demonstrated with the development of the GRASP (Generalized Retrieval of Aerosol and Surface Properties) algorithm. The GRASP algorithm is based on several generalization principles with the idea to develop a scientifically rigorous and versatile algorithm. It has significantly extended capabilities and areas of applicability and can be applied to diverse remote sensing observations. This paper also illustrates the practical applicability of GRASP and, therefore the multi-term LSM, in diverse situations.

<footnote n="1">¹Univ. Lille, CNRS, UMR 8518 - LOA - Laboratoire d'Optique Atmosphérique, Lille, France, ²GRASP-SAS, Villeneuve d'Ascq, France, ³School of Meteorology, The University of Oklahoma, Norman, OK, United States, ⁴State Key Laboratory of Severe Weather (LASW) and Key Laboratory of Atmospheric Chemistry (LAC), Chinese Academy of Meteorological Sciences, Beijing, China, ⁵NASA Langley Research Center, Hampton, VA, United States, ⁶NASA Goddard Space Flight Center, Greenbelt, MD, United States, ⁷Physics Department, University of Maryland Baltimore County, Baltimore, MD, United States, ⁸Aerospace Information Research Institute, CAS, Beijing, China, ⁹Institute for Space Science, Free University of Berlin, Berlin, Germany, ¹⁰Univ Paris Est Creteil and Université de Paris, CNRS, LISA, Créteil, France, ¹¹LuftBlick OG, Austria Instead of LuftBlick GmbH, Innsbruck, Austria, ¹²Cloudflight GmbH, Linz, Austria</footnote>
GRASP has two main independent modules. The first module is a numerical inversion that includes general mathematical operations not related to a particular physical nature of the inverted data. Numerical inversion is implemented as a statistically optimized fitting of observations following the multi-term LSM strategy. The presentation of the GRASP numerical inversion provides a profound description of the main methodological aspects used for establishing a multi-term LSM approach that is aimed at applying multiple a priori constraints in the retrieval. The foundation of this approach uses the fundamental frameworks of the Method of Maximum Likelihood (MML) and LSM statistical estimation concepts. We discuss the asymptotical optimal properties of MML and LSM estimates in order to emphasize the importance of statistical estimation methods in remote sensing.

We also compare the multi-term LSM with other established statistical optimization approaches of numerical inversion that are constrained by a priori information, such as Bayesian concepts and the Optimal Estimation approach (Rodgers, Inverse Methods for Atmospheric Sounding: Theory and Practice, 2000). In addition, we discuss several other methodological inversion optimization aspects, including accounting for non-negativity of measured and retrieved characteristics, optimizing inversion in non-linear cases, accounting for measurement redundancy, estimating contributions of random measurement and systematic errors on the retrieval uncertainties, etc. All of these aspects are complimentary to multi-term LSM that need to be fully addressed in the development of practically efficient inversion algorithms. These methodological developments have already been applied to several remote sensing algorithms, such as the algorithm by Dubovik and King (J. Geophys. Res., 2000, 105, 673–696) that has been employed for more than two decades for operational processing of observations from AERONET ground-based Sun-sky radiometers, and the algorithm by Dubovik et al. (J. Geophys. Res., 2006, 111, D11208) developed for retrieval of aerosol particle refractive index together with size and shape distributions from full phase matrix measurements.

The second main module of GRASP is the forward model. Similarly to the numerical inversion module, it has been developed in a universal way for simulating various atmospheric remote sensing observations with high accuracy. As a result, GRASP is a highly versatile algorithm that can be applied to diverse passive and active satellite and ground-based atmospheric observations and is inherently designed for synergetic retrievals when different observations are inverted simultaneously. Depending on the input data, GRASP can retrieve detailed columnar and vertical aerosol properties and surface reflectance. Diverse approaches for modeling aerosol and surface properties together with different a priori constraints can be used in GRASP retrievals.

Thus, the GRASP package can be considered as a platform for developing, testing, and refining novel retrieval concepts and their utilization in operational processing environment. In addition, GRASP is designed as a practically efficient, transparent, and accessible community open source algorithm that can be used as an advanced tool for verifying different

retrieval concepts and realizing those concepts in high-performance operational software. At present, GRASP is being adapted to reprocess the observations provided by several satellite instruments, ground-based networks, single instruments, aircraft and for various synergetic data sets that combine coordinated passive and active remote sensing observations. For example, GRASP has been adapted for operational reprocessing of observations from POLDER-1, -2, -3, MERIS, and AATSR/Envisat satellite missions. There are developments of operational retrievals of aerosol and surface properties from OLCI/Sentinel-3, Sentinel-5P, 3MI/Metop-SG, and Sentinel-4 geostationary observations. GRASP has also been used for synergetic aerosol retrievals by inverting a combination of ground-based lidar and radiometer data. GRASP has also been used for interpretation of airborne and laboratory nephelometer measurements, and it has been used for developing new approaches to derive diverse aerosol properties from ground-based direct sun and diffuse sky radiation measurements by radiometers and sky cameras.

GRASP algorithm and software organization are described in detail, and an overview of GRASP applications and data products is provided.

# 1 INTRODUCTION

The multi-term LSM (Least Square Method) strategy has been proposed and advocated for more than two decades in series of studies (Dubovik et al., 1995; Dubovik and King, 2000; Dubovik 2004; Dubovik et al., 2006; Dubovik et al., 2008; Dubovik et al., 2011; King and Dubovik, 2013, etc.). These realized algorithmic developments showed the multi-term LSM as fruitful methodology both for understanding the advantages of well-known approaches of constrained inversion and for constructing new practically efficient retrieval methods relying elaborated a priori constraints. The above algorithmic studies and implementations have also advanced many complimentary aspects addressing various issues and needs of successful implementation of the proposed retrieval methodology in practical applications. Specifically, the questions of implementing multi-term LSM strategy in both linear and non-linear cases, accounting for non-negativity nature of measured and retrieved characteristics, dealing with redundancy of observations, assuring consistency of multiple constraints, estimating retrieval errors, similarities and differences with other methodologies, etc. were discussed and analyzed in depth. It was shown that most of standard approaches for constraining inversion solution could be derived and optimized in the frame of multi-term LSM strategy. Moreover, this strategy is not limited to any particular application and advocated as quite universal concept for designing new practical procedures for advanced constrained inversion techniques. As a practical result of these multi-year methodological efforts a highly versatile algorithm GRASP (Generalized Retrieval of Aerosol and Surface Properties) for atmospheric remote sensing has been proposed and implemented. The description of GRASP algorithm structure
<figure><img src="image_4.png" /></figure>

and details of implementations will be used in this paper for illustrating the practical utilization of multi-term LSM approach. The GRASP development was pursuing the idea for creating the algorithm that works with diverse remote sensing applications, is convenient for application to diverse synergetic processing of the complimentary observations and well adapted for continuing evolution once used by the community. **Figure 1** illustrates the concept of the algorithm. The algorithm is a useful practical tool for interpretation of actual remote sensing data and illustrates the usefulness of utilizing the multi-term LSM in diverse applications.

The GRASP development was initiated by Dubovik et al. (2011) in the frame of the efforts on designing satellite algorithm of new generation for improving aerosol retrieval from POLDER-3/PARASOL imager over land surfaces where the high surface reflectance dwarfs the signal from aerosols. It derives simultaneously a large number of parameters characterizing aerosol and surface properties. GRASP also relies on the heritage of earlier efforts on developments for retrieving aerosol properties from the ground-based observations by AERONET network radiometers (Dubovik et al., 2000; Dubovik and King, 2000; Dubovik, 2004; King and Dubovik, 2013) and from laboratory measurements of aerosol phase matrices (Dubovik et al., 1995; Dubovik et al., 2006, etc.). All above algorithms had significant similarities and even common blocks used for mathematical inversion and observation modeling. Therefore, it has become quite logical and appealing from diverse considerations to merge these different retrievals in one single algorithm based on unified principles. Having such strategy is helpful for benefiting from

historical positive experience in new retrieval developments and for developing successive strategy for evaluating and improving these unified principles and therefore for advancing the retrieval strategy in general. The overall idea of such algorithm and many specific features are aimed to address the existing challenges of modern remote sensing data interpretation (Dubovik et al., 2021). Thus, GRASP unified algorithm and software was developed as a tool applicable to many different observations that is alternative to development of several retrieval codes dedicated to interpretation of different remote sensing observations. GRASP has been developed as freely assessable software from a platform for GRASP open source code (https://www.grasp-open.com/), introductory video is available at (http://www.youtube.com/watch?v=PcDeqwDF15A). Moreover, recently GRASP-CLOUD services have been offered to remote sensing community described in details by Aspetsberger et al. (2019).

GRASP is based on several generalization principles with the idea to develop scientifically rigorous, versatile, practically efficient, transparent and accessible algorithm. Two main modules of GRASP numerical inversion and forward model are independent, can be changed separately and even fully replaced if needed for a particular application. Indeed, numerical inversion includes general mathematical operations and forward model module is developed for simulating various atmospheric remote sensing observations. Both modules are described below in details. As it will be shown many specific features in both modules can be adapted to a particular application. Following the multi-term LSM strategy the numerical inversion was implemented as a statistically optimized fitting of observations that unites the advantages of
a variety of known approaches (Dubovik, 2004). For example, in aerosol retrievals, different smoothness constraints were applied simultaneously on aerosol size distributions and spectral dependencies of aerosol refractive index (see Dubovik and King, 2000; King and Dubovik, 2013). In addition, the multi-term LSM has been used for formulation of multi-pixel retrieval that implements simultaneous optimized inversion of a large group of independent observation (see Dubovik et al., 2011). This inversion scheme improves retrieval consistency by using known limitations on spatial and/or temporal variability of retrieved parameters. For example, in satellite retrieval, the horizontal pixel-to-pixel variations of aerosol and day-to-day variations of surface reflectance are enforced to be smooth by an additional set of a priori constraints. This concept is also promising for developing synergetic retrievals using observations that are not fully coincident in time or/and space as illustrated below in the paper.

The forward model simulates atmospheric radiation resulting from interaction of incident light with atmosphere and underlying surface. The aerosol model is assumed to be a mixture of several aerosol components with particles possibly having different compositions and shapes. A number of different concepts for parameterizing particle size distributions, refractive indices dependencies and shape mixtures can be used. A BRDM (Bi-directional Reflectance Distribution Matrix) that includes Bidirectional Reflectance Distribution Function (BRDF) and Bidirectional Polarization Distribution Function (BPDF) are used to account for the effects of surface reflectance, respectively (see GRASP Forward Model). The algorithm fully accounts for multiple interactions of scattered solar light with aerosol, gas and surface by means of solving radiative transfer equation. As one of the specific characteristics of GRASP retrieval, all calculations are done on-line without searching a pre-calculated look-up table. At the same time, several different strategies based on trade-offs by increasing retrieval speed by the price of reducing completeness of retrieval have been realized for adapting GRASP for processing large amount of observations in timely manner for example as required from Near-Real-Time (NRT) processing of satellite observations. Thus, the structure of the algorithm is convenient for adapting and testing alternative routines for handling aerosol, surface, gas or multiple scattering contributions. This flexibility is convenient for adapting the forward modeling strategy to a completeness of the input information and time performance requirements. This makes GRASP uniquely adaptable both to ground-based observations, where information content is high and a very complex retrieval can be realized, and to single-view satellites, where the amount of input information is limited, and the forward modeling needs to be simplified.

GRASP is a versatile algorithm that can be applied to various passive and active remote sensing observations made from satellites, ground, aircraft or in laboratory and also to diverse combinations of those measurements. Depending on the input data, GRASP can be configured to retrieve columnar and vertical aerosol properties and surface reflectance. Specifically, GRASP is being adapted for reprocessing complex POLDER observations (Dubovik et al., 2011; Dubovik et al., 2019; Li et al., 2019). Chen

et al. (2020) provide detailed description and discussion of currently available PODLER-3/GRASP aerosol products. Some demonstration of POLDER-3/GRASP retrievals are also documented in several other papers (Bovchaliuk, et al., 2013; Milinevsky et al., 2014; Kokhanovsky et al., 2015; Popp et al., 2016; Neukermans et al., 2018, Chen et al., 2018; Chen et al., 2019; Frouin et al., 2019, Remer et al., 2019, Sogacheva et al., 2020, Wei et al., 2020; Wei et al., 2021, etc.). It has been shown that using these advanced multi-angular polarimetric measurements GRASP can derive detailed aerosol information including aerosol absorption – a property that is very difficult to derive from satellite observations. This concept now has been adapted for research NRT retrieval from future 3MI/Metop-SG mission. At the same time, the application of GRASP was shown to be useful for deriving aerosol properties from single- and dual-view polar MERIS and AATSR/Envisat, OLCI/Sentinel-3, Sentinel 5P and geostationary Sentinel-4 observations, etc. Obviously, these observations contain less information and the number of parameters retrieved by GRASP is significantly lower than from polarimetric observation. Nonetheless, the base generalized approach of GRASP retrieval remains the same, i.e., from practical application point of view GRASP satellite retrievals are formally implemented in the same way for all observations:

1) The same assumptions are used everywhere, i.e., there are no location specific assumptions.

2) The same set of a priori constraints is applied globally.

3) The retrieval uses the same initial guess globally (e.g. for aerosol parameters);

4) Minimal or no post-processing used for the obtained results.

The outcome of such strategy will be discussed and illustrated below using the real results and; their validation.

GRASP can be equally applied to various ground-based observations. Since overall it inherits inversion strategy and many specific elements from the algorithm developed by (Dubovik and King, 2000; Dubovik et al., 2006, etc.) for inverting AERONET data (Holben et al., 1998) it can be naturally applied to the observations of ground-based sun/sky-radiometers. It should give very close results to AERONET operational code if applied in identical way as done in the network processing. At the same time, GRASP includes some additional features that can be of interest for obtaining some new results from ground-based radiometers. For example, GRASP can be applied to aerosol size information from direct Sun observations only (Torres et al., 2017), or it can implement simultaneous inversion of sequence of Sun and/or Sun/sky radiometric observations under a priori constraints on aerosol variability using multi-pixel approach. This will be demonstrated below. Schuster et al. (2019) used the code with Polarized Imaging Nephelometer measurements to mimic the AERONET retrieval algorithm and validated the results with simultaneous laboratory measurements of scattering, absorption, and particle size. Moreover, the GARRLiC (Generalized Aerosol Retrieval from Radiometer and Lidar Combined data) as a branch of GRASP was developed for inversion of coincident multi-spectral lidar and radiometer observations by Lopatin et al. (2013). Using this
synergy, GARRLiC/GRASP retrieves improved aerosol columnar properties together with details about vertical aerosol variability including profiles of fine and coarse aerosol mode concentrations and vertical profiles of $\omega_0(\lambda)$, as will be illustrated below. GARRLiC/GRASP has been used in a number of studies and employed in routine processing within Aerosol, Clouds and Trace Gases Research Infrastructure (ACTRIS) (Tsekeri, et al., 2013; Bovchaliuk, et al., 2016; Benavent-Oltra et al., 2017; Tsekeri, et al., 2017; Román et al., 2018; Benavent-Oltra et al., 2019; Hu et al., 2019; Herreras et al., 2019; Titos et al., 2019; Konsta, et al., 2021). The evolution of the GARRLiC/GRASP branch, its applications to more complex lidar systems and more complex synergy processing is discussed by Lopatin et al. (2021). There are ongoing developments for utilization of the GRASP for processing other ground-based observations as those by sun and lunar sky cameras (Román, et al., 2017; 2021). In addition, GRASP has been adapted and extensively used for obtaining the aerosol observations from airplane and laboratory measurements by polar-nephelometer (Espinosa et al., 2017; Espinosa W. R. et al., 2019; Schuster et al., 2019). Such measurements by polar-nephelometer of angular distribution of intensity and polarization, as well as, total aerosol scattering and absorption, that are not affected by multiple scattering effects, are very useful for validating used aerosol model and overall verification of the concept of aerosol retrieval from remote sensing.

It should be noted that the main potential and advantage of unified GRASP software is that there are many different options for implementing retrieval and for generating and displaying the outputs that are developed in core software and can be used in different applications. Specifically, different approaches for modeling atmospheric radiation can be used in retrieval. For example, modeling of scattering by aerosol polydispersions can be realized using different models for parameterizing aerosol size distributions. Also, the aerosol retrieval can be set to derive optical properties of aerosol particles such as complex refractive index along with the size information as was realized by Dubovik and King (2000). Alternatively, the aerosol can be modeled as an internal or external mixture of different chemical aerosol components and the retrieval could be aimed for deriving aerosol chemical composition from the remote sensing observations. Such approach has been integrated in GRASP for both satellite and ground-based measurements, as discussed by Li et al. (2019). Similarly, different surface reflectance models can be used in GRASP retrievals. Also, different possibilities for using a priori constraints are available. Namely, the retrievals can be constrained by direct a priori estimated of unknowns, by using a priori smoothness assumptions on variability of any retrieved functions (e.g., size distributions, spectral dependencies, etc.) or by using a priori smoothness assumptions on variability of retrieved parameters in time or in space if multi-pixel approach is used, or by using all or several of these constraints simultaneously. For all the retrieved parameters and functions derived from them GRASP can provide the error estimation by calculating the elements of the covariance matrix. Also, GRASP is designed to provide the values of broadband fluxes and radiative forcing using retrieved

parameters (see Derimian et al., 2016). In addition, for some situations the possibilities of estimating aerosol type, air quality (PM 2.5) and other characteristics of high practical and scientific interest have been introduced into GRASP.

The details of using different options in the retrieval and forward modeling, scientific and methodological concepts of setting up and utilization of GRASP in diverse applications are discussed the following sections in detail.

## 2 THEORY OF MULTI-TERM LSM INVERSION

The inversion module is one the key parts of the GRASP algorithm and software. It is implemented based on rather general principles that are not directly related with any specific application. This is why, this module can be used with diverse forward models that make the GRASP concept fundamentally generalizable. It is important to note that the inversion module is based on a rather extensive and positive heritage of inversion algorithm developments for atmospheric remote sensing (Dubovik et al., 1995; Dubovik et al., 2000; Dubovik and King, 2000; Dubovik, 2004; Dubovik et al., 2006; Dubovik et al., 2008; Dubovik et al., 2011, etc.). As a result, the GRASP inversion strategy and software module has a generalized character and, at the same time, it is well tuned for interpretation of ground-based and satellite observations.

The inversion is implemented as a statistically optimized fitting of observations following the Multi-term LSM strategy that unites the advantages of a variety of approaches; it provides transparency and flexibility in algorithm development inverting passive and/or active observations and deriving several groups of unknown parameters (e.g., Dubovik, 2004, etc.). This methodology is designed to use the underlying LSM concept for building statistically optimum solution for practical situations when different observations and a priori data should be processed and interpreted together. Overall, the approach is based on the fundamental Method of Maximum Likelihood (MML). If there is a vector $f^*$ composed by the measurements of physical characteristics $f_j(a)$ depending on the vector of parameters $a$ and this vector contains random measurement errors $\Delta f^*$, i.e.,

$$f^* = f(a) + \Delta f^*, \qquad (2.0.1)$$

then the vector of unknowns $a$ derived from Eq. 2.0.1 will inevitably contain some errors $\Delta a$. Here and everywhere below in the text, the vectors are denoted by lowercase italic letters in bold, while the matrices are denoted by the upper-case letters in bold. According to MML the solution $\hat{a}^{best}$ should result in modeled errors $\Delta \hat{f}^* = \hat{f}_j^* - \hat{f}_j(\hat{a})$, corresponding to the most probable error realization, i.e., to the maximum of Probability Density Function (PDF):

$$P(\Delta \hat{f}^*) = P(f(a) - f^*) = P(f(a)|f^*) = \max \qquad (2.0.2)$$

where $P(f(\hat{a})|f^*)$ is the PDF of unknown parameters $a$ conditioned on the measurements $f^*$ and called the Likelihood Function. The MML is one of the strategic principles of statistical
estimation that provides statistically the best solution in many senses (Edie et al., 1971). For example, under few rather general conditions (the derivatives of PDF should exist and be limited in the entire parameter space, etc.) the asymptotical error distribution (for infinite number of $\Delta f^*$ realizations) of MML estimates $\hat{a}$ have the smallest possible variances of $\Delta \hat{a}_i$. Most of statistical properties of the MML solution remain optimum for a limited number of observations (Edie et al., 1971). Under the assumptions of a normal PDF, i.e.,

$$P(f(\mathbf{a})|f^*) = \sim \exp\left(-\frac{1}{2}(f(\mathbf{a}) - f^*)^\mathrm{T}\mathbf{C}_f^{-1}(f(\mathbf{a}) - f^*)\right), \quad (2.0.3)$$

where $C_f$ is the covariance matrix of errors $\Delta f^*$. Based on MML, the maximum of PDF corresponds to the minimum of exponent of Eq. 2.0.3, i.e., the solution should satisfy to the so-called LSM condition:

$$2\hat{\Psi}(\mathbf{a}) = (f(\mathbf{a}) - f^*)^\mathrm{T}\mathbf{C}_f^{-1}(f(\mathbf{a}) - f^*) = \min. \quad (2.0.4)$$

The minimum of the multi-term quadratic form $\hat{\Psi}(\mathbf{a})$ corresponds to its extremum, i.e., to a point with zero gradient $\nabla\hat{\Psi}(\mathbf{a})$:

$$\nabla\hat{\Psi}(\mathbf{a}) = 0. \quad (2.0.5)$$

For the general case of non-linear functions $f_k(\mathbf{a})$ the solution can be found iteratively:

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

where $\Delta\mathbf{a}^p$ is the solution that can be found by solving the system of so-called normal equations, i.e.:

$$(\mathbf{K}_p^T \mathbf{C}_f^{-1} \mathbf{K}_p) \Delta\mathbf{a}^p = \mathbf{K}_p^T \mathbf{C}_f^{-1} \Delta f^p, \quad (2.0.7)$$

where $\Delta f^p = f(\mathbf{a}^p) - f^*$, and $\mathbf{K}_p$ is Jacobian matrix with the elements $\{\mathbf{K}_p\}_{ji} = \frac{\partial f_i(\mathbf{a})}{\partial a_j}|_{\mathbf{a}=\mathbf{a}^p}$. Then, the minimized quadratic form given by Eq. 2.0.4 has approaches a $\chi^2$ distribution with known the asymptotic limit (e.g., Fourgeaut and Fuchs 1967; Edie et al., 1971):

$$2\Psi(\mathbf{a}) = \min \rightarrow (N_f - N_a), \quad (2.0.8)$$

where $N_f$ - number of measurements and $N_a$ - number of unknowns, i.e. numbers of elements in the vectors $f^*$ and $\mathbf{a}$.

## 2.1 Multi-Term LSM Approach for Inverting Multi-Source Data

The multi-term LSM concept is adapted in GRASP for combining information from different data sources. This approach has been formulated for implementing the flexible constrained inversion as a statistically optimized inversion of multi-source data. The concept follows the developments of Dubovik et al. (1995), Dubovik and King (2000), Dubovik et al. (2000), Dubovik (2004), Dubovik et al. (2008), and Dubovik et al. (2011). Such an approach is naturally derived from Eq 2.0.7 and allows generalizing various inversion formulas into a single formalism. The approach considers the situation when

measurement vector $f^*$ in Eq. 2.0.1 is obtained from different and independent “sources”, i.e., $f^*$ and its covariance matrix has the following array structure:

$$f^* = \begin{pmatrix} f_1^* \\ f_2^* \\ \vdots \\ f_K^* \end{pmatrix} \quad \text{and} \quad \mathbf{C}_f = \begin{pmatrix} \mathbf{C}_1 & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{C}_2 & \mathbf{0} & \mathbf{0} \\ \vdots & \vdots & \vdots & \vdots \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{C}_K \end{pmatrix} \quad (2.1.1)$$

where $(f^*)^\mathrm{T} = (f_1^*, f_2^*, \dots, f_K^*)^\mathrm{T}$ is a column vector, $C_k$ is the covariance matrix of the k-th data set $f_k^*$ and index k denotes different data sets or “sources”. This assumes that the data from the same source have similar error structures that are independent of errors in the data from other sources. For example, direct Sun and diffuse sky radiances have different magnitudes and are measured by sensors with different sensitivities, i.e., their errors should be independent (due to the different sensors) and likely have different magnitudes. In this situation, it is convenient to consider Eq. 2.0.1 as a joint system:

$$f_k^* = f_k(\mathbf{a}) + \Delta f_k^*, \quad (k = 1, 2, \dots, K), \quad (2.1.2)$$

where the index k denotes different data sets. It should be noted here that from the formal viewpoint, the only difference between Eq. 2.1.2 and Eq. 2.0.1 is that Eq. 2.1.2 explicitly outlines an expectation of an array structure for the covariance matrix $C_f$. Such explicit demarcation of the input data makes the retrieval more transparent because the statistical optimization of the retrieval is driven by a covariance matrix of random errors. It should, be noted that Eq. 2.1.2 does not assume any relations between the different forward models $f_k(\mathbf{a})$, i.e. forward models $f_k(\mathbf{a})$ can be the same or different.

Thus, the joint 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:

$$P(f(\mathbf{a})|f^*) = P(f_1(\mathbf{a}), \dots, f_K(\mathbf{a})|f_1^*, \dots, f_K^*) = \prod_{k=1}^{K} P(f_k(\mathbf{a})|f_k^*). \quad (2.1.3)$$

Then, under the assumptions of a normal PDF, one can write

$$P(f(\mathbf{a})|f^*) = \prod_{k=1}^{K} P(f_k(\mathbf{a})|f_k^*) \sim \exp \left( -\frac{1}{2} \sum_{k=1}^{K} (f_k(\mathbf{a}) - f_k^*)^\mathrm{T} \mathbf{C}_k^{-1} (f_k(\mathbf{a}) - f_k^*) \right). \quad (2.1.4)$$

Thus, for multi-source data, the solution should correspond to minimum of K quadratic forms:

$$2\hat{\Psi}(\mathbf{a}) = 2 \sum_{k=1}^{K} \hat{\Psi}_k(\mathbf{a}) = \sum_{k=1}^{K} (f_k(\mathbf{a}) - f_k^*)^\mathrm{T} \mathbf{C}_k^{-1} (f_k(\mathbf{a}) - f_k^*) = \min, \quad (2.1.5)$$

and for the non-linear functions $f_k(\mathbf{a})$ the solution can be found iteratively as shown by Eq 2.0.6, 2.0.7, while the left and
right parts of normal system of Eq. 2.0.7 can be written as a sum of terms:

$$ \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 (2.1.6) $$

where $\mathbf{K}_{k,p}$ are Jacobian matrices at p-th iteration of the functions $\mathbf{f}_k(\mathbf{a})$ in the vicinity of $\mathbf{a}^p$. This is conventional LSM solution that can be easily derived for the linear $\mathbf{f}_k(\mathbf{a})$ and can also be used for the case when $\mathbf{f}_k(\mathbf{a})$ is non-linear (see, e.g., Dubovik, 2004 and discussion in Section 1.3.2 below). Similarly, the asymptotic limit of the minimized quadratic form, for most applications, can be written as:

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

where $N_K$ is a number of measurements (inputs) in the k-th data set and $n$ is the number of retrieved unknowns (parameters). Similarly, a priori data sets are included in the K data sets. They are statistically independent of the measurements, i.e. they have errors with a different level of accuracy uncorrelated with the remote-sensing errors.

It should be noted that in practical application it is often convenient to use weighting matrices $\mathbf{W}_k$ instead of covariance matrices $\mathbf{C}_k$ (e.g., Linnik 1962; Twomey 1977). When several data sets inverted together the use of weighting matrices and convenient transparency into approach because it makes it evident the relative contribution of the data from different data sources (e.g., Dubovik et al., 1995; Dubovik and King, 2000; Dubovik 2004, etc.). Specifically, using the weighting matrices Eq. 2.1.6 can be written as

$$ \left( \sum_{k=1}^{K} \gamma_k \mathbf{K}_{k,p}^{T} \mathbf{W}_k^{-1} \mathbf{K}_{k,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}^{p}. \quad (2.1.8) $$

Here the contribution of different terms in Eq. 2.1.1 are scaled by the corresponding Lagrange multipliers $\gamma_i$, defined as:

$$ \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 (2.1.9) $$

where $\varepsilon_i^2$ is the first diagonal element of $\mathbf{C}_i$, i.e., $\varepsilon_i^2 = \{\mathbf{C}_i\}_{11}$. Evidently, $\gamma_1 = 1$. As discussed by Dubovik and King (2000), Dubovik (2004), and Dubovik et al. (2011), etc., such renormalization is also convenient, because with the definition of the minimized quadratic function (or “residual”) given by Eq. 2.1.8, the measurement error $\varepsilon_1^2$ can be estimated from the residual of the fit. Indeed, once the weighting matrices used in the solution provided by Eq. 2.1.8 minimizes quadratic with the limit depending on $\varepsilon_1^2$:

$$ 2\Psi'(\mathbf{a}) = 2\varepsilon_1^2 \Psi(\mathbf{a}) = \min \rightarrow \varepsilon_1^2 \left( \sum_{k=1}^{K} N_k - n \right) \quad \text{and} \quad (2.1.10) $$

$$ \varepsilon_1^2 \approx \frac{2\Psi'(\mathbf{a}^p)}{\sum_{k=1,...,K} N_k - N_a}. \quad (2.1.11) $$

### 2.1.1 A Priori Constraints in Multi-Term LSM Approach

In many practical situations, the information from measurement is not sufficient. In such situation a solution of Eq. 2.0.7 can be non-unique and inclusion of additional a priori information (constraints) becomes necessary. Multi-term LSM concept has been proposed for the integration of different types of a priori constraints in remote sensing applications (Dubovik et al., 1995; Dubovik and King 2000; Dubovik et al., 2000; Dubovik 2004; Dubovik et al., 2008; Dubovik et al., 2011). In this approach a priori estimates are considered to be “equivalent” to the measurements in the sense that the a priori data are characterized by their PDF and could be treated equivalently to the actual measurements. In this regards Eq 2.1.6-Eq 2.1.9 do not show any formal distinction between different $\mathbf{f}_k(\mathbf{a})$. In practice there are always two different types of data sets: direct measurements and quantified a priori knowledge about unknowns $\mathbf{a}$. Therefore, the vector of the measurement $(\mathbf{f}^*)^T = (\mathbf{f}_1^*, \mathbf{f}_2^*, \dots, \mathbf{f}_K^*)^T$ can be formally written as:

$$ (\mathbf{f}^*)^T = (\mathbf{f}_1^*, \mathbf{f}_2^*, \dots, \mathbf{f}_K^*, \mathbf{f}_1^a, \mathbf{f}_2^a, \dots, \mathbf{f}_N^a)^T, \quad (2.1.12) $$

where $\mathbf{f}_i^* = \mathbf{f}_i^*(\mathbf{a})$ represent directly measured characteristics and $\mathbf{f}_i^a = \mathbf{f}_i^*(\mathbf{a})$ represent a priori known characteristics of unknowns $\mathbf{a}$. Correspondingly the right side of Eq. 2.1.3 can be formally split in two groups:

$$ P(\mathbf{f}(\mathbf{a})|\mathbf{f}^*) = \prod_{k=1}^{K} P(\mathbf{f}_k(\mathbf{a})|\mathbf{f}_k^*) \prod_{n=1}^{N} P(\mathbf{f}_n^a(\mathbf{a})|\mathbf{f}_n^{a,*}), \quad (2.1.13) $$

where the first group represents K sets of data obtained from actual measurements, the second group represent N sets of data $\mathbf{f}_n^{a,*}$ obtained from a priori knowledge.

Correspondingly, the resulting Eq. 2.1.6 can also be formally arranged to identify the contribution of measurements and a priori terms.

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

where two groups of the terms in left and right parts of the equation represent contributions of the set of K measured characteristics $\mathbf{f}_k(\mathbf{a})$ and the set of N a priori $\mathbf{f}_n^a(\mathbf{a})$ characteristics.

Thus, the above Multi-term approach is a simple rearranging of the base LSM formulation. At the same time, it provides the solid basis for unifying many known formulas of constrained inversion in a single formalism, as discussed by Dubovik (2004). In addition, it was shown to be practically convenient and efficient approach for developing remote sensing algorithms using diverse complimentary observations and a priori constraints. Specifically, the Multi-term approach has been successfully used in previous inversion algorithms for laboratory (Dubovik et al., 2006), ground-based (Dubovik and King, 2000; Dubovik et al., 2000; Dubovik et al., 2002), airborne (Gatebe et al., 2010), satellite observations (Dubovik et al., 2011) and in inverse modeling (Dubovik et al., 2008). All remote
sensing development have been integrated into a versatile GRASP algorithm (Dubovik et al., 2014). Many of the technical details adapted in the GRASP have been previously shown and discussed. Therefore, the Sections below outline only newly introduced elements and are focused on conceptual elements of inversion as well as on user-oriented explanations of the GRASP algorithm and software utilization.

## 2.1.2 A Priori Constraints Used in GRASP Algorithm

The multi-term LSM concept allows flexible utilizations of nearly arbitrary a priori constraints, i.e., knowledge about any a priori known function $f_n^a(a)$ of unknowns $a$. However, only a limited number of the most popular and physically transparent a priori constraints are presently used in GRASP. At the same time, the utilization of other a priori constraints that are not included in the proposed options can also be realized upon a request.

The base configuration of GRASP proposes two types of a priori limitations: “single-pixel” and “multi-pixel” constraints. The utilization the word “pixel” in this terminology originated from interpretation of satellite images that are composed of separate pixels. Correspondingly, “single-pixel” relates to the constraints on parameters within each pixel independently, while “multi-pixel” relates to the constraints on inter-pixel spatial or temporal variability of the retrieved parameters. This classification (with the introduction of “multi-pixel” constraints) is convenient not only for the interpretation of satellite images, but it is also convenient for many other applications (laboratory and ground-based measurements) where input data can be related in time or space.

### 2.1.2.1 “Single-Pixel” Constraints—A Priori Constraints for Coincident and Collocated Observations

This type of constraint includes a priori information about each retrieved parameter (or characteristic in each set of measurements) that can be used absolutely independently on the measurements or as a priori knowledge in other data sets. The most common example of this type of constraint is the application of direct a priori estimates of unknowns $a^*$, i.e.,

$$f_1^{a,*} = f_1^{a,*}(a) + \Delta f_1^{a,*} \rightarrow a^* = a + \Delta a, \quad (2.1.15)$$

where $\Delta a$ are the uncertainties of 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. 2.1.14 by defining: $K_{a,p} = 1$ - unity matrix; i.e., $K_{a,p}^T C_{a}^{-1} K_{a,p} = C_{a}^{-1}$ and $K_{a,p}^T C_{a}^{-1} f_1^{a,*} = C_{a}^{-1} a^*$. 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.

Another common type of a priori constraint are smoothness constraints that limit the variability of retrieved functions by using a priori knowledge about limitations on derivatives of those functions. For example, a priori knowledge limits high frequency variations of continuous functions $v(x)$, such as the aerosol size distribution, spectral dependence of the refractive index, parameters of surface reflectance, etc. From a formal point

of view, the smoothness constraints are related by limited values of the derivatives, i.e., with deviations of their m-th derivative deviations from zero:

$$\frac{\partial^m v(x)}{\partial x^m} \approx 0. \qquad (2.1.16)$$

For the vector of unknowns $a = (a_1, a_2, \dots, a_n)^T$ that contains the elements describing these continuous functions $v(x)$, the knowledge on the smoothness of the retrieved function can be defined using a vector-matrix linear system (e.g., see Dubovik 2004; Dubovik et al., 2011):

$$f_2^{a,*} = f_2^{a,*}(a) + \Delta f_2^{a,*} \rightarrow \mathbf{0}^* = G_m a + \Delta_g^*, \quad (2.1.17)$$

where $G_m$ is the Jacobian matrix of the matrix of the m-th derivatives, $\mathbf{0}^*$ is zero vector representing the fact that the a priori estimates of the corresponding derivatives equal to zeros. In practice, these are often approximated by matrices of the m-th finite difference estimated in point $a$. The errors $\Delta_g^*$ reflect the uncertainty in the knowledge of the deviations of $v(x)$ from the assumed constant ($m=1$), straight line ($m=2$), parabola ($m=3$), and so on.

Under the assumption that the $\Delta_g^*$ in Eq. 2.1.15 have a normal distribution with the unbiased covariance matrix $C_{\Delta m}$, these constraints can be easily included in Eq. 2.1.14 by defining: $K_{a,2,p} = G_m$, where $G_m$ is Jacobian of m-th derivatives, and $f_2^{a,*} = \mathbf{0}^*$. Correspondingly, the following can be written $K_{a,2,p}^T C_{a,2}^{-1} K_{a,2,p} = G_m^T C_{\Delta m}^{-1} G_m$ and $K_{2,p}^T C_{2}^{-1} \Delta f_2^{a,*} = G_m^T C_{\Delta m}^{-1} (a^p - \mathbf{0}^*) = G_m^T C_{\Delta m}^{-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 article by Tikhonov (1963) and Tikhonov's later studies.

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

$$\left( \sum_{k=1}^{K} \gamma_k K_{k,p}^T W_k^{-1} K_{k,p} + \gamma_a W_a^{-1} + \sum_{n=1}^{N} \gamma_n \Omega_n \right) \Delta a^p = \sum_{k=1}^{K_1} \gamma_k K_{k,p}^T W_k^{-1} \Delta f_k^p + \gamma_a W_a^{-1} (a^p - a^*) + \sum_{n=1}^{N} \gamma_n \Omega_n a^p, \quad (2.1.18a)$$

where $\Omega_n$ denotes the smoothness matrices defined as

$$G_m^T W_{2,m_n}^{-1} G_m = \Omega_{2,m_n} = \Omega_n \quad (2.1.18b)$$

Equation 2.1.18 include several smoothness constraints so that several physically independent functions can be retrieved simultaneously under different a priori smoothness constraints. For example, Dubovik and King (2000) retrieve aerosol size distribution and spectral dependence of refractive index.

Thus, Eq 2.1.18 generalize the commonly used equations of constrained inversion by Phillips (1962), Tikhonov (1963), Twomey (1975), Rodgers (1976), Twomey (1977), Rodgers (1990), and Rodgers (2000) and allows for significantly extended flexibility in practical use (Dubovik and King, 2000; Dubovik 2004, etc.)
The values of uncertainties in a priori constraints are considered to be directly comparable to uncertainties in the actual measurements and even defined in the GRASP software relatively to the measurement errors.

The definition of measurement uncertainty in GRASP code is the following:

1) Type of measurement errors:

- for each measurement data set $f_k^*$ two types of errors $\Delta f_k^*$ can be set: relative or absolute;

2) The magnitude of the errors:

- for each data set the magnitude of the errors is defined by the standard deviation $\epsilon_i$ and a weighting matrix $W_i$ (assumed as the unity matrix by default). The standard deviation is used inside of the code to calculate corresponding Lagrange multipliers $\gamma_i$.

Correspondingly the magnitudes of uncertainties in a priori estimates are defined in the same way as measurements. Specifically, the following “single-pixel” a priori constraints can be set in GRASP code.

“Single-pixel” a priori constraints:

1) Direct a priori estimates of unknown parameters $\mathbf{a} = (a_1, a_2, \dots, a_n)^T$:

- for each parameter $a_i$ an a priori value $a_i^*$ can be provided with the corresponding Lagrange multipliers $\gamma_i$;

- it is also possible to assume a vector $\mathbf{a}^*$ of a priori estimates for all the retrieved parameters or for selected groups (e.g. parameters describing size distribution) with common Lagrange multipliers $\gamma_a$ and weighting matrix $W_a$.

2) The vector $\mathbf{a}$ of unknowns includes groups of parameters describing specific retrieved characteristics $\mathbf{a} = (a_{sd}, a_n, a_{brdf}, a_h, \dots)^T$, where $a_{sd}, a_n, a_k, a_{brdf}, a_h$ are vectors containing only the parameters of size distribution, real and imaginary refractive indices and surface BRDF and aerosol vertical profile correspondingly. For each such group of parameters $a_{sd}, a_n, a_k, a_{brdf}, a_h$, etc. that describing describes a continuous function (e.g., $a_{sd}$) in the retrieved vector $\mathbf{a} = (a_{sd}, a_n, a_{brdf}, a_h, \dots)^T$ the smoothness constraints can be set by defining:

- the order $m$ of limited derivatives ($m = 0$ – constant; $m = 1$ – straight line; $m = 2$ – parabola, etc.)

- the strength of the applied a priori smoothness constraints is defined by Lagrange multipliers $\gamma_n$.

The effect on the solution of applying smoothness constraints using the limitations on the derivatives can be easily illustrated. Namely, in Eq. 2.1.18 one can use only the linear system representing the a priori constraints given by Eq. 2.1.17, i.e., $\mathbf{0}^* = \mathbf{G}_m \mathbf{a} + \Delta_g^*$. Correspondingly, the solution of such system is provided by iterations shown by Eq. 2.1.18, where

only one term is used; i.e., Eq 2.1.18 is reduced to $\Omega_n \Delta \mathbf{a}^p = \Omega_n \mathbf{a}^p$. Evidently, the solution of such system is non-unique because the matrix $\mathbf{G}_m$ has a smaller number of rows (n-m – number of finite differences that approximate the corresponding derivatives) than columns (n – number of unknowns) and therefore it is evident that the matrix $\Omega_n$ has zero determinant and the solution is non-unique, i.e., many $\mathbf{a}^p$ could provide the equality. Correspondingly, the result strongly depends on the initial guess, i.e., iterative solutions initialized with different initial guesses would converge to different solutions. At the same time, all the solutions should satisfy the smoothness equations, i.e., each solution should be a constant for $m = 1$, a straight line for $m = 2$ and a parabola for $m = 3$, while the exact placement of constant, straight line or parabola depends on the initial guess. Figure 2 illustrates this consideration by showing the solutions chosen by the iterative retrieval from the same initial guess using different a priori smoothness constraints.

It should be noted, that all above “single-pixel” and “multi-pixel” (discussed in the next sections) a priori constraints implemented in the software were already actively used in practical applications and, therefore should be considered as recommended. For example, in case of “single-pixel” a priori constraints, the use of a priori estimate is often useful for the parameters or characteristics when their absolute values known from climatology or other ancillary information and the retrieved values do not expect deviate from these values and overall sensitivity of observation is not very high to these parameters. For example, such constraints are used for complex refractive index in GRASP-AOD applications (see Ground-Based Passive Observations). In satellite retrievals discussed in Satellite and Airborne Passive Observations, climatological values are often used as direct a priori estimates for some parameters of land surface reflectance. “Single-Pixel” smoothness constraints are often used for retrieval continuous function, such as aerosol size distribution, spectral dependence of complex refractive index, spectral dependence of the parameters of land surface reflectance. At the same time, the strength of a priori constraints varies significantly for different parameters. For example, smoothness constraints are much stronger for such parameters as spectral dependence of aerosol complex refractive index, then for size distribution (e.g., see Dubovik and King, 2000). Similarly the first parameter in the models of BRDF land surface reflectance is expected to have strong variation while other parameters (second, third ...) describing land surface reflectance anisotropy are nearly spectrally constant (see Dubovik et al., 2011).

Finally, the practical libraries of GRASP-OPEN software provide large number of the examples of GRASP settings files that are useful examples showing how the a priori constraints in other users in realized applications.

### 2.1.2.2 “Multi-Pixel” Constraints—A Priori Constraints (2-D Case) for Coordinated but Non Coincident or/and Non-Collocated Observations

The GRASP “multi-pixel” fit is statistically optimized inversion that is implemented for several independent observations simultaneously (Dubovik et al., 2011), e.g., for a group of
<figure><img src="image_5.png" /><figcaption>**FIGURE 2** | Illustration of the effect of smoothness constraints limiting derivatives of different order on the solution. The illustrations show the solutions chosen for Eq. 2.1.17 alone (no other data inverted simultaneously) by the iterative retrieval from the same initial guess under different a priori smoothness constraint. The first, second and third derivatives are a priori limited in the (A), (B), (C) correspondingly.</figcaption></figure>

multiple satellite image pixels. Such simultaneous retrieval can be easily designed using Multi-term LSM concept defined in Eq. 2.1.1-Eq 2.1.14, as described in details in paper by Dubovik et al. (2011). In “multi-pixel” fitting the state vector is composed of the state vectors from several (n) pixels:

$$ \mathbf{a} = (\mathbf{a}_1, \mathbf{a}_2, \dots, \mathbf{a}_n)^T, \qquad (2.1.19) $$

where each component $\mathbf{a}_i$ represents a state vector for a set of observations (or a “pixel” for satellite images). The observations and a priori constraints for each “pixel” are defined in exactly the same way as the “single pixel” retrieval described in Section 2.1.2.1. At the same time, implementing “multi-pixel” fitting of the data allows one to additionally use inter-pixel constraints. In fact, a priori constraints about known limited inter-pixel variability of retrieved parameters can be realized by using a priori knowledge about limitations on derivatives on time or spatial variability of parameters retrieved from observations in different pixels. Indeed, observations of different pixels provide important information about the temporal and spatial characteristics of the retrieved parameters. For example, a satellite can observe the same pixel in time and several neighboring pixels simultaneously. In principle, the variability of each physical parameter $\mathbf{a}_i$ can be considered as a value of a continuous function $\mathbf{a}_i = \mathbf{a}_j(x, y, z, \dots, t)$. Therefore, the limitation on the variability of every parameter in time and space that can be used as additional constraints. Similar to Eq. 2.1.16, inter-pixel variability constraints are related with limited values of the derivatives, i.e., with their m-th derivative deviations from zero. At present, the time- and spatial-variation of each parameter in GRASP can be limited using the following a priori assumptions:

$$ \frac{\partial^m \mathbf{a}_i(x, y, z, t \dots)}{\partial x^m} \approx 0, \quad \frac{\partial^m \mathbf{a}_i(x, y, z, t \dots)}{\partial y^m} \approx 0, \\ \frac{\partial^m \mathbf{a}_i(x, y, z, t \dots)}{\partial z^m} \approx 0, \quad \text{and} \quad \frac{\partial^m \mathbf{a}_i(x, y, z, t \dots)}{\partial t^m} \approx 0, \qquad (2.1.20) $$

and can be presented in matrix as similar to Eq. 2.1.17 written for a single-pixel case as:

$$ \left\{ \begin{aligned} \mathbf{f}_2^{a,*} &= \mathbf{f}_2^{a,*}(\mathbf{a}) + \Delta\mathbf{f}_2^{a,*} \\ \mathbf{f}_3^{a,*} &= \mathbf{f}_3^{a,*}(\mathbf{a}) + \Delta\mathbf{f}_3^{a,*} \\ \mathbf{f}_4^{a,*} &= \mathbf{f}_4^{a,*}(\mathbf{a}) + \Delta\mathbf{f}_4^{a,*} \\ \mathbf{f}_5^{a,*} &= \mathbf{f}_5^{a,*}(\mathbf{a}) + \Delta\mathbf{f}_5^{a,*} \end{aligned} \right. \rightarrow \left\{ \begin{aligned} \mathbf{0}_x^* &= \mathbf{G}_{x,m_x}\mathbf{a} + \Delta\mathbf{0}_x^* \\ \mathbf{0}_y^* &= \mathbf{G}_{y,m_y}\mathbf{a} + \Delta\mathbf{0}_y^* \\ \mathbf{0}_z^* &= \mathbf{G}_{z,m_z}\mathbf{a} + \Delta\mathbf{0}_z^* \\ \mathbf{0}_t^* &= \mathbf{G}_{t,m_t}\mathbf{a} + \Delta\mathbf{0}_t^* \end{aligned} \right. \qquad (2.1.21) $$

where $m_x, m_y, m_z$ and $m_t$ are the orders of the derivatives used for limiting the time (t) and spatial (x, y, z) variation of retrieved parameters. These orders can be different for each dimension x, y, z and t.

Then, the solution for the multi-pixel fitting equivalent of Eq. 2.1.18 (written for a single-pixel case) can be presented as:

$$ \left( \begin{array}{c} \mathbf{A}_{1,p} \quad \mathbf{0} \quad \dots \quad \mathbf{0} \\ \mathbf{0} \quad \mathbf{A}_{2,p} \quad \dots \quad \mathbf{0} \\ \dots \quad \dots \quad \dots \quad \dots \quad + \mathbf{\Omega}_{\text{inter}} \end{array} \right) \Delta \mathbf{a}^p \\ = \left( \begin{array}{c} \nabla \Psi_1(\mathbf{a}_1^p) \\ \nabla \Psi_2(\mathbf{a}_2^p) \\ \dots \\ \nabla \Psi_N(\mathbf{a}_N^p) \end{array} \right) + \mathbf{\Omega}_{\text{inter}} \mathbf{a}^p, \qquad (2.1.22) $$

where $\mathbf{A}_{i,p}$ and $\nabla\Psi_i(\mathbf{a}_i^p)$ refer to the left and right parts of Eq. 2.1.8 defined for i-th single pixel, so that Eq. 2.1.8 can be denoted compactly as $\mathbf{A}_{i,p}\Delta\mathbf{a}_i^p = \nabla\Psi_i(\mathbf{a}_i^p)$. The matrix $\mathbf{\Omega}_{\text{inter}}$ of multi-pixel constraints is defined via smoothness matrices for the spatial and temporal variability of each retrieved parameter as:

$$ \mathbf{\Omega}_{\text{inter}} = \gamma_x \mathbf{\Omega}_{x,m_x} + \gamma_y \mathbf{\Omega}_{y,m_y} + \gamma_z \mathbf{\Omega}_{z,m_z} + \gamma_t \mathbf{\Omega}_{t,m_t}. \qquad (2.1.23) $$

The detailed description of multi-pixel constraints and their application is provided in the paper by Dubovik et al. (2011). Specifically, Dubovik et al. (2011) provide explicit expressions for $\mathbf{\Omega}_{x,m_x}$, $\mathbf{\Omega}_{y,m_y}$ and $\mathbf{\Omega}_{t,m_t}$, as well as corresponding Lagrange multipliers $\gamma_x$, $\gamma_y$ and $\gamma_t$. The equations are fully implemented in the GRASP code. The vertical multi-pixel constraints defined by $\mathbf{\Omega}_{z,m_z}$ and $\gamma_z$
that appear logically in general Eq 2.1.20-Eq 2.1.23 are not yet implemented in a similar way as is done for horizontal and temporal constraints, because in practice there are almost no fully independent vertical observations. Instead, as was discussed in a previous Section, the a priori limitations on vertical variability are included as a part of single-pixel a priori constraints that seem to address most needs in practical applications (Lopatin et al., 2013; Lopatin et al., 2021). The illustrations of a priori constraint on vertical variability are discussed in Multi-Term LSM: Retrieval Practice Using GRASP. It should be noted that the rather complex and explicit expressions provided by Dubovik et al. (2011) are of primary interest for methodological analysis. At the same, the implementation of vertical or any other multi-pixel constraints in practical algorithms can be efficiently realized in computer routines using conceptual Eq 2.1.20-Eq 2.1.23 without using explicit visualization of the final equations.

Thus, the following “multi-pixel” a priori constraints can be set in GRASP code.

“Multi-pixel” a priori constraints:

1) For each single parameter $a_i$ variability along coordinate x can be limited by applying a priori constraints on the derivatives of $a_i = a_i(x)$. For example, in satellite retrievals, the coordinate x relates with latitudinal variability. The smoothness constraints on $a_i = a_i(x)$ can be set by defining:

- order $m$ of limited derivatives [$m = 0 - a_i = a_i(x)$ is close to a constant; $m = 1 - a_i = a_i(x)$ is close to a straight line; $m = 2 - a_i = a_i(x)$ is close to a parabola, etc.]

- the strength of applied a priori smoothness constraints is defined by Lagrange multipliers $\gamma_{x,i}$.

2) For each single parameter $a_i$ a priori constraints on variability along coordinate y, i.e., $a_i = a_i(y)$ can be limited. For example, in satellite retrievals, the coordinate y relates with longitudinal variability. The smoothness constraints on $a_i = a_i(y)$ can be set by defining:

- order $m$ of limited derivatives [$m = 0 - a_i = a_i(y)$ is close to a constant; $m = 1 - a_i = a_i(y)$ is close to a straight line; $m = 2 - a_i = a_i(y)$ is close to a parabola, etc.]

- the strength of applied a priori smoothness constraints is defined by Lagrange multipliers $\gamma_{y,i}$.

3) For each single parameter $a_i$ a priori constraints on time variability, i.e., $a_i = a_i(t)$ can be limited. The smoothness constraints on $a_i = a_i(t)$ can be set by defining:

- order $m$ of limited derivatives [$m = 0 - a_i = a_i(t)$ is close to a constant; $m = 1 - a_i = a_i(t)$ is close to a straight line; $m = 2 - a_i = a_i(t)$ is close to a parabola, etc.]

- the strength of applied a priori smoothness constraints is defined by Lagrange multipliers $\gamma_{t,i}$.

For the parameters with multi-dimensional variability $a_i = a_i(x,y,t)$ the above constraints can be applied simultaneously.

For example, aerosol concentrations, values of refractive index, non-spherical fraction and any parameter used for the description of surface reflectance BRDF and BPDF can be considered as characteristics $a_i = a_i(x,y,t)$ that change horizontally and in time. In case of applying multi-pixel constraints on the retrieved multi-dimensional function $a_i = a_i(x,y,t)$, the smoothness constraints would force the projection of the solution at each coordinates x, y or t to take a from close to a constant, a straight line or a parabola depending which of corresponding partial derivatives $\frac{\partial^m a_i(x,y,t)}{\partial x^m}$, $\frac{\partial^m a_i(x,y,t)}{\partial y^m}$ or $\frac{\partial^m a_i(x,y,t)}{\partial t^m}$ with respect to x, y or t are assumed to be close to zero. Therefore limiting several different partial derivatives will force $a_i(x,y,t)$ to the corresponding 3D surface.

**Figure 3** illustrates the application of multi-pixel constraints in 2D. Similarly to the single-pixel case shown in **Figure 2**, one can use only reduced linear system Eq. 2.1.22 using terms representing multi-pixel a priori constraints. Evidently, the solution of such system is non-unique because the matrix in the left part of Eq. 2.1.22 is degenerated (it has smaller number of rows than columns), and the results will strongly depend on the initial guess, i.e., iterative solutions initialized with different initial guesses would converge to different solutions. Situations with a priori constraints limiting different partial derivatives are shown in **Figure 3**. One can see that the resulting solutions are always converge to smooth surfaces with projection on the axis x or y represented a constant, a straight line or a parabola depending which of corresponding partial derivatives in respect to coordinate x or y are limited.

In practical applications “multi-pixel” constraints were introduced and proven to be efficient in satellite applications. Specifically, in GRASP simultaneous retrievals of land surface and aerosol form satellite observations (see *Satellite and Airborne Passive Observations*, the “multi-pixel” smoothness constraints were applied to limit temporal variability for retrieved values of land reflectance and spatial variability for retrieved values of aerosol parameters (see Dubovik et al., 2011). Utilization of temporal “multi-pixel” constraints is also was shown for diverse synergy application (see discussion in *Ground-Based Active Observation and Synergy With Passive Observations and Synergetic Retrievals and Other Diverse GRASP Applications*). For example, Lopatin et al. (2021) used the “multi-pixel” smoothness constraints on temporal variability of aerosol parameters in simultaneous inversion of the different ground-based observations, that were obtained in the same place but not at the same time. In such approach, the “multi-pixel” smoothness constraints helped to significantly enrich the aerosol characterization.

It should be noted that, in principle, the choice of possible multi-pixel constraints is not limited to the set of constraint possibilities listed above or already realized in the software. For example, the a priori limitations can be expected for higher order or mixed partial derivatives:
<figure><img src="image_6.png" /><figcaption>FIGURE 3 | Illustration of the effect from multi-pixel smoothness constraints (2-D case). The illustrations show the solutions chosen for Eq. 2.1.2 alone (no other data inverted simultaneously) by the iterative retrieval from the same initial guess under different a priori smoothness constraint. Each panel corresponds to different set of a priori limited derivatives that shown by equation on the top of each panel.</figcaption></figure>

$$ \frac{\partial^{(k_1+k_2+k_1+k_2)} a_i(x, y, z, t)}{\partial x^{k_1} \partial y^{k_2} \partial z^{k_3} \partial t^{k_4}} \approx 0. \qquad (2.1.24) $$

Also, the derivatives may also have non-zero average trends which could also be useful practical constraints. Explicit derivation of the matrices $\Omega_{inter}$, as well as understanding and visual interpretation of a priori constraints on the partial derivatives shown by Eq. 2.1.24 can be quite challenging for mixed partial derivatives of high order in N-D cases. At the same time, handling these cases using computer routines in inversion algorithms is rather straightforward. In addition, as emphasized by Dubovik et al. (2011), the matrices $\Omega_{inter}$ are very sparse, therefore solving an apparently large matrix system of Eq. 2.1.22 can be efficiently achieved using standard numerical solvers developed for sparse systems.

In addition, Xu et al. (2019) proposed to integrate some elements of Principal Component (PC) analysis into the multi-pixel formalism. In this approach the vector of unknowns is represented via expansion using PC basis and the actual PC components and the coefficients are sought as a part of the solution. This approach allows reduction of the solution of dimensionality, decrease of

calculation time and other potentially interesting solution optimizations.

All above features are considered for integration in the future evolution of GRASP.

### 2.1.3 Comparison With Bayesian Strategy and Optimal Estimation Approach

As discussed above, in a contrast with original approaches of constrained inversion originated by Phillips (1962), Tikhonov (1963), Twomey (1975), and Twomey (1977) the Multi-term LSM is fully based on statistical estimation methodology that suggests use of MML and LSM as the most practically efficient approach. In this regard, Multi-term LSM approach is fundamentally similar to conventional Bayesian approach used in Optimal Estimation (OE) algorithms (Rodgers 2000). At the same time, Multi-term LSM has some distinct and practically useful differences with OE approach that will be discussed below. Specifically, OE is absolutely equal to a particular case of Multi-term LSM when only direct a priori estimates are used as constraints, meanwhile the Multi-term LSM approach allows rather straightforward application of much more refined a priori constraints.
In probability theory and statistics, Bayes theorem describes the probability of an event, based on a prior knowledge of conditions that might be related to the event (e.g., Idie et al., 1971). In a similar way, any retrieval can be considered as a modification of prior knowledge about unknowns by adding information from indirect measurements $f^*(a)$. Such vision can be expressed mathematically using MML formalism. Namely, in situation when available a priori information can be expressed as a priori estimates $a^*$ statistically independent of the measurements $f^*$, the a priori estimates $a^*$ can be described by Likelihood Function i.e., PDF function $P(a|a^*)$ that is independent of the Likelihood Function of the measurements $P(f(a)|f^*)$. Then, correspondingly the Likelihood Function written if both a priori estimates $a^*$ and measurements $f^*$ are available as a joint PDF $P(a|f^*, a^*)$ can be presented as the product of $P(a|a^*)$ and $P(f(a)|f^*)$:

$$P(a|f^*, a^*) = P(f(a)|f^*)P(a|a^*). \qquad (2.1.25)$$

Correspondingly, the MML for such a situation can be written as follows:

$$P(a|f^*, a^*) = P(f(a)|f^*)P(a|a^*) = \max. \qquad (2.1.26)$$

In most practical cases, it is reasonable to assume that there is some a priori knowledge $a^*$ about unknowns $a$ (even if only quite uncertain). Following such logic, Eq. 2.1.25 can be considered as a rather basic MML formulation as suggested by studies of Rodgers (2000) where it is also called the OE technique. In the simplest linear case, the OE solution corresponds to the solution of the following linear system:

$$(K^T C^{-1} K + C_a^{-1}) a = K^T C^{-1} f^* + C_a^{-1} a^*. \qquad (2.1.27)$$

This OE technique is rather popular for designing linear optimized estimators for constraint inversion in remote sensing applications and it is often considered as a universal approach for integrating a priori constraints in the retrieval. However, in this regard, Multi-term LSM can propose a somewhat different view on the definition of a priori constraints that can be more convenient in some practical situations. Indeed, the Multi-term LSM Eq. 2.1.13 can be formally considered as an equivalent of Bayesian or OE formulations. Namely, the first and second terms on the right-hand side of Eq. 2.1.13 become:

$$P(f(a)|f^*) = \prod_{k=1}^{K} P(f_k(a)|f_k^*) P(a|a^*), \qquad (2.1.28)$$

where PDF of $a^*$ is defined as

$$P(a|a^*) = \prod_{n=1}^{N} P(f_n^a(a)|f_n^{a,*}) \qquad (2.1.29)$$

For the cases where all PDFs are normal functions, the PDF $P(a|a^*)$ on the left side can be described only by the $a^*$ covariance matrix $C_a$, while the right side is a product of several PDFs, each of them given for a different function $f_n^a(a)$ and is defined by the mean vector $f_n^*$ and its covariance matrix $C_a$. Correspondingly, for many practical applications, using Eq. 2.1.26 directly and

finding the $a^*$ covariance matrix $C_a$ directly is not easy and even nearly infeasible. In this regard, a Multi-term LSM equation that considers an a priori PDF as a product of several different PDFs of a priori constraints is significantly more convenient for a number of practical situations. For example, Eq. 2.1.17 that was quite naturally defined above in the frame of Multi-term LSM can be obtained from optimum estimation Bayesian Eq. 2.1.25 by assuming the following $a^*$ and $C_a$:

$$a^* = \left( \gamma_a W_a^{-1} + \sum_{n=1}^{K_2} \gamma_n \Omega_n \right)^{-1} \gamma_a W_a^{-1} (a^p - a^*) \qquad (2.1.30)$$

and

$$C_a = \varepsilon_1^2 \left( \gamma_a W_a^{-1} + \sum_{n=1}^{K_2} \gamma_n \Omega_n \right)^{-1} \qquad (2.1.31)$$

respectively. Indeed, Eq. 2.1.28 are exactly the two key terms of Multi-term LSM expression introduced in the form of Eq. 2.1.13. Obviously, that such $a^*$ and $C_a$ can only be constructed using considerations of $P(a|a^*)$ and a product of many PDFs as shown in Eq 2.1.28, 2.1.29. Constructing such $a^*$ and $C_a$ could probably be useful in some situations, but generally is not necessary and often may even be impossible. For example, even the one of the simplest pioneering Phillips formulas written for the linear case:

$$(K^T C_f^{-1} K + \Omega_m) a = K^T C_f^{-1} f^*, \qquad (2.1.32)$$

Can be derived from OE Bayesian Eq. 2.1.27 assuming $a^*$ and $C_a$ as:

$$a^* = 0^* \quad \text{and} \quad C_a = \Omega_m^{-1}. \qquad (2.1.33)$$

However, $\Omega_m$ is a smoothness matrix (generally defined as $\Omega_m = G_m^T G_m$, e.g. Dubovik, 2004) that has det = 0, and therefore defining matrix $C_a$ explicitly as given in Eq 2.1.30, 2.1.31 is an ambiguous and unachievable task. Evidently, more complex Multi-term LSM based equations such as Eq 2.1.18, 2.1.22 can also be presented in the OE form, however the derivations of those equations in the first place from OE would not be feasible.

Finally, we would like to note the substantial conceptual difference between the Multi-term LSM approach and the OE approach. Indeed, the main idea of the multi-term LSM concept is the consideration that a priori estimates are equivalent to the measurements. In other words it is expected that a priori data can be characterized by their PDF and be treated equivalently to actual measurements. This is a different premise than the main idea of the Bayesian statistics approach, which assumes that the developer should always utilize a priori information about the entire solution before attempting the retrieval. Thus, this Bayesian assumption implicitly suggests having an a priori estimate for the entire unknown state vector as a starting point of developing the retrieval algorithm. With that in mind, the developer may feel a necessity to use at least some values as a priori estimates; for example, using his/her best guess about the possible solution even without having full justification for using
these estimates. As a result, if no objective link of the assumed a priori information has been established, the Bayesian approach becomes vulnerable to subjective assumptions of the developer (which somewhat contradicts to the principle of scientific objectivity). Thus, in spite the fundamental mathematical equivalence of Multi-term LSM and OE, the differences in the base concepts may have essential influence on forming specific methodological guidelines that may lead significantly different practical recommendations.

It is interesting that from a historical prospective, Ronald Fisher was among the opponents of the Bayesian statistics and promoted development of MML as an alternative strategy (see Fisher, 1956; Agresti and Hitchcock, 2005). Fatefully, the highly popular OE approach promoted by the textbook Rodgers (2000) is often considered as an equivalent of MML by remote sensing community (especially by inexperienced readers taking the text book as introductory coarse in retrieval methodology), and it somewhat promotes the Bayesian ideas.

## 2.2 Further Aspects of Inversion Optimization

The optimality of MML and LSM concept is a key methodological aspect used for realizing numerical inversion in GRASP. In this respect, it is necessary to remember that this optimality can be achieved under certain conditions that need to be respected; otherwise applying the MML and LSM concept may lead to unsatisfactory results.

### 2.2.1 Logarithmic Transformation as Rigorous Approach to Account for Non-Negativity of Measured and Retrieved Positively Defined Parameters

One of the situations where this aspect requires specific attention relates to securing the positive solutions in the retrieval of non-negative characteristics. Historically, the straightforward use of LSM-based linear constrained inversion could generate negative values for physically non-negative characteristics in some applications, and this created contradicting situations wherein empirical, largely intuitive algorithms (Chahine 1968; Twomey 1975) performed better than algorithms based on rigorous principles of statistical optimization. Such results reduce the trust and interest of the inverse algorithm developer community in studying the rather complex statistical estimation formalism and utilizing it in practical applications. At the same time, as was discussed by Dubovik and King (2000) and Dubovik (2004), after some adjustment of hypothesis about error statistic and correct application of the concept, the MML and LSM provide very fruitful solution recipes for such situations.

Specifically, as suggested by Dubovik et al. (1995), Dubovik and King (2000) and Dubovik (2004) the non-negativity constraints can be included in the statistical estimation concept by applying lognormal noise assumption in the retrieval optimization. This assumption of lognormal noise (instead of conventional normal noise) leads to implementing inversions in logarithmic space, i.e., to employing a logarithmic transformation of the forward model. Retrieval of logarithms of a physical characteristic (instead of absolute values) is an obvious way to avoid negative values for positively defined characteristics.

However, the literature devoted to inversion techniques tends to consider this apparently useful tactic as an artificial trick rather than a scientific approach to optimize solutions.

Such above misconception is probably caused by the fact that the pioneering efforts on inversion optimization by Phillips (1962), Tikhonov (1963), Twomey (1963) and many of later theoretical considerations (e.g., Hansen, 1992) were devoted to solving the Fredholm integral equation of the first kind, i.e., a system of linear equations produced by a quadrature. Examples include the retrieval of aerosol size distribution (King et al., 1978) or temperature profile of the atmosphere (Rodgers, 1976) by inverting spectral dependence of optical thickness. Considering optical thickness as a function of the logarithm of the aerosol concentrations or temperature profile requires replacing initially linear equation $f = K a$ by nonlinear one $f_j = f_j(\ln a_i)$. On the face of it, such a transformation of linear problems to non-linear ones is not enthusiastically accepted by scientific community as an optimization. On the other hand, in cases when a forward model is a nonlinear function of the retrieved parameters (e.g., atmospheric remote sensing in cases when multiple scattering effects are significant), the retrieval of logarithms is more easily accepted as a logical approach. Besides, non-linear Chahine-like iterative methods where proven to be efficient for inverting linear systems while Dubovik et al. (1995) and Dubovik and King (2000) have shown that under some assumptions these methods are equivalent to LSM formulated for logarithms.

Rigorous statistical considerations also reveal some limitations of applying Gaussian function for modeling errors in the measurement of positively defined characteristics. Indeed, the curve of the normal distribution is symmetrical and therefore the assumption of a normal PDF is equivalent to the assumption of the principal possibility of obtaining negative results even in the case of physically non-negative values. For such non-negative characteristics as intensities, fluxes, etc., the choice of the log-normal distribution for describing the measurement noise seems to be more correct due to the following considerations:

1) log-normally distributed values are positively-defined;

2) there is a number of theoretical and experimental reasons showing that for positively defined characteristics the log-normal curve (multiplicative errors) is the best for modeling random deviations in non-negative values [see the discussion of statistical experiments in the textbook of Tarantola (1987)].

In addition, from the practical viewpoint, using the lognormal PDF for noise optimization does not require any revision of normal concepts and can be implemented by simple transformation of the problem to the space of normally distributed logarithms. Therefore, there is a clear basis for considering log-normal PDF for measurements of fundamentally non-negative characteristics $f$ (such as intensities, for example). Correspondingly, logarithmical transformation can be applied to the initial system of equations, i.e.;

$$f^* = f(a) + \Delta f^* \rightarrow \ln f^* = \ln f(a) + \Delta \ln f^* \quad (2.2.1)$$
where $\Delta \ln f^*$ is a vector of measurement errors $\Delta \ln f_j^* = \ln f_j^* - \ln f_j^{\text{real}}$ that have Normal PFD $P(\ln f(\hat{a}) | \ln f^*)$. In reality the actual measurements as those of intensity are often really obtained using receivers with logarithmic responses.

In principle, MML or LSM do not directly assume distribution of errors in the final solution and formally in the unknowns $a$. At the same time, based on known properties of MML solution (see Edie et al., 1971), the MML or LSM estimates are also asymptotically normally distributed. It is obvious then, that MML or LSM estimates $a_i$ have asymptotically normally distributed errors; thus estimate based on normal PDF $P(\hat{a} | \ln f^*)$ cannot provide zero probability for $a_i < 0$ even if $a_i$ are positively defined by nature. On the other hand, the retrieval of logarithms $\ln a_i$ instead of absolute values $a_i$ eliminates the above contradiction, because the LSM estimates of $\ln a_i$ would have the normal distribution of $P(\ln \hat{a} | \ln f^*)$ (i.e., lognormal distribution of $a_i$ that assures positivity of non-negative $a_i$). In fact, the necessary conditions for optimality of MML and LSM is that the first derivatives of the measurements with the respect to the unknowns should be limited in the whole considered range of their variability (see Edie et al., 1971). In this respect, if measured $f$ and retrieved $a$ characteristics are positively defined $df_i/da_i$ are not limited when $a_i \to 0$ or $f_j \to 0$ and do not exist for negative $a_i$ and $f_j$. Correspondingly, application of MML and LSM in such conditions would not lead to asymptotically optimum solutions.

Thus, if the unknowns are positively defined parameters, there is a clear rationale in retrieving logarithms of unknowns instead of their absolute values, since the log-normal statistics (multiplicative errors) are more natural for them. Moreover, if both measured $f$ and retrieved $a$ characteristics are positively defined, the function $\ln f(\ln a)$ should likely be “well-behaved”, i.e., both the function and the derivatives $d\ln f_j/d\ln a_i$ should exist and be limited in whole range of the variability of both $\ln f$ and $\ln a$. Correspondingly the application of MML shouldn't be questioned. It can be noted here, once again, that Dubovik et al. (1995) and Dubovik and King (2000) demonstrated that the efficient (Chahine 1968; Twomey 1975) iterative procedures can be derived using LSM in logarithmic space.

Thus, accounting for non-negativity of solutions and/or non-negativity of measurements in GRASP is implemented by using of logarithms of unknowns ($a_i \to \ln a_i$) and/or measurements ($f_j \to \ln f_j$). Some empirical or semi-empirical parameters that can have negative values but bounded by a value, e.g. "-a", are made positively defined by adding a constant value "shift" - "b", so that the a+b is always >0 and logarithmic transformation for these parameters is possible.

## 2.2.2 Non-Linearity of Forward Problem and Levenberg-Marquardt Optimization of Iterative Solution

Most of theoretical derivations for inverse problems in general, and especially for statistically optimized inversion, are discussed for linear functional dependencies such as those proposed by the Fredholm integral equation of the first kind. At the same time the majority of practical remote sensing applications deal with the interpretation of observed characteristics that have non-linear dependence with sought parameters. Therefore, for algorithm

users it is often difficult to comprehend the potentials of different optimization approaches and identify the most appropriate retrieval development strategy. In this regard, most (and nearly all) considerations of the statistical optimization are based on the fact that in a small vicinity of the actual solution, variations of the observations are related linearly with the variation of the atmospheric parameters to be retrieved. Indeed, since the uncertainties in the observations are expected to be small, the optimization of statistical properties of the solution is done in the linear approximation even in non-linear algorithms. Moreover, the general numerical solution of a non-linear problem is formulated as an iterative procedure wherein the solution corrections are obtained by solving the problems in the linear approximation. Therefore, there is a very significant similarity in the equations used for inverse algorithms in linear and non-linear cases. The main difference is that non-linear inversion algorithms are almost always iterative and need to converge to a solution from a chosen initial guess. The achievement of convergence is an inherent necessity for the development of successful iterative solutions. Although the iterative solutions can be used for solving both linear and non-linear systems, the successful convergence of the algorithms is a more profound and a broader issue for non-linear inverse problems.

The solution of system Eq. 2.0.1 provided by Eq 2.0.6, 2.0.7 is written via iterations that can be used in both situations when $f_i^a(a)$ are linear or non-linear. In fact the square system $f(a) = f^*$ for the case of a non-linear $f_i^a(a)$ can be solved using Newtonian iterations: $a^{p+1} = a^p - K_p^{-1}(f(a^p) - f^*)$, where $K_p$ is the Jacobian of $f(a^p)$. The minimum of Eq. 2.0.4 corresponds a zero gradient $\nabla\Psi(a) = 0$ and can be found as $a^{p+1} = a^p - (K_{\nabla\Psi,p})^{-1}\nabla\Psi(a^p)$, where $K_{\nabla\Psi,p}$ is the Jacobian of $\nabla\Psi(a^p)$. For Eq. 2.1.5 one can write:

$$ \nabla\Psi(a^p) = K_p^T C^{-1} (f(a^p) - f^*), \quad (2.2.2) $$

and

$$ K_{\nabla\Psi,p} = K_p^T C^{-1} K_p + D_p \approx K_p^T C^{-1} K_p. \quad (2.2.2) $$

The elements of matrix $D_p$ include second derivatives of $f(a^p)$ and therefore can be neglected. Thus, under the above assumption Eq 2.0.6, 2.0.7, provide a so-called Quasi-Newtonian solution minimizing the quadratic function Eq. 2.0.4 for the case with non-linear functions $f(a)$ (a more detailed discussion can be found in Dubovik (2004) and various text books on numerical methods, (e.g., Tarantola, 1987):

$$ a^{p+1} = a^p - (K_p^T C^{-1} K_p)^{-1} K_p^T C^{-1} (f(a^p) - f^*). \quad (2.2.3) $$

Implementing a non-linear inversion by Newton-like methods requires assurance of iteration convergence. Iteration by Eq 2.0.6, 2.0.7 may not converge at all or converge to a wrong solution. The convergence difficulties may be caused by inadequate choice of the initial guess and/or limitations of the linear approximation used for the iterative guess correction. Indeed, for strongly non-linear functions $f_j(a)$, the minimized form $\Psi(a)$ may have a complex structure with several minima. The analysis of this
structure is desirable prior to inversion. However, when three or more unknowns are to be retrieved, such analysis is practically infeasible. Usually, researchers repeat retrievals with a set of initializations and select the best solution. The initializations and the criteria for selecting the best solution are commonly established based on the physical constraints of the application, experience, and intuition of the developers. Also, a convergence of non-linear solutions can be improved by modifying Eq. 2.0.7. The most established modification of Gauss-Newton iterations is widely known as the Levenberg-Marquardt method (Ortega and Reinboldt 1970; Press et al., 1992):

$$a^{p+1} = a^p - t_p (K_p^T C^{-1} K_p + \gamma D_p)^{-1} K_p^T C^{-1} (f(a^p) - f^*), \quad (2.2.4)$$

where matrix $D_p$ and the coefficients $0 < t_p \le 1$ and $\gamma \ge 0$ are selected empirically to provide convergence. The matrix $D_p$ is predominantly diagonal (unity matrix is often chosen as $D_p$) and addition of the term $\gamma D_p$ to $K_p^T C^{-1} K_p$ in Eq. 2.2.1 is analogous to using a priori estimates in linear inversions. Specifically, the matrix $K_p^T C^{-1} K_p$ can be singular on some of p-th iterations even if it is non-singular in the vicinity of the solution. Adding the term $\gamma D_p$ to $K_p^T C^{-1} K_p$ helps to pass the iteration process through the areas of $K_p^T C^{-1} K_p$ singularities. As pointed out by Press et al. (1992), the Levenberg-Marquardt formula generalizes the steepest descent method. Namely, Eq. 2.2.4 can be reduced to the steepest descent method by defining matrix $D_p$ in Eq. 2.2.4 as the unity matrix I and prescribing a large value to the parameter $\gamma$, i.e., $\gamma D_p = \gamma I \gg K_p^T C^{-1} K_p$. Thus, Eq. 2.2.4 always converges with appropriate selection of $\gamma$.

The multiplier $0 < t_p \le 1$ in Eq. 2.2.4 is invoked mainly to decrease the length of $\Delta a^p$, because the linear approximation may overestimate the correction $\Delta a^p$. Usually, $t_p$ is decreased by a factor (e.g., by 2, as done in GRASP) until a condition $\Psi(a^{p+1}) < \Psi(a^p)$ is satisfied. Underestimation of $\Delta a^p$ does not lead to a convergence failure and may only slow down the arrival to a solution. The addition of the term $\gamma D_p$ also reduces $\Delta a^p$. Correspondingly, using both $\gamma D_p$ ($\gamma > 0$) and $0 < t_p \le 1$ in the same iteration may seem redundant because both operations reduce $\Delta a^p$. However, using the multiplier $t_p$ is straightforward (compared to adding $\gamma D_p$) and sufficient in the application with moderately non-linear forward model with non-singular $K_p^T C^{-1} K_p$. On the other hand, if the matrix $K_p^T C^{-1} K_p$ is singular in some points $a^p$, using $t_p \le 1$ does not help and the inclusion of constraining term $\gamma D_p$ is necessary. Thus, the use of both $t_p$ and $\gamma D_p$ modifications in Levenberg-Marquardt (Eq. 2.2.4) complement each other in practice.

The above-mentioned similarity of predominantly diagonal matrix $D_p$ in Eq. 2.2.4 to a priori estimates term on the left side of Eq. 2.1.18 is used in GRASP following earlier methodological consideration by Dubovik and King (2000) and Dubovik (2004). These considerations are based on the fact that employing linear approximations for non-linear functions $f(a)$ in Eq. 2.0.6, 2.0.7 may result to a converge failure. Therefore, in case of inversion of Multi-term equations (as in Eq. 2.1.8, 2.1.14), if some of $f_k(a)$ are linear, it seems logical to expect less problems with convergence. This idea is elaborated by considering linear constraints applied to non-linear iterations (Dubovik (2004)). In principle, although, Eq. 2.1.8, 2.1.14 are constrained by a priori (presumably linear)

terms, solution $\hat{a}^p$ may fail to converge similarly to basic Newton-Gauss iterations by Eq. 2.2.23, 2.2.3. This is because at initial iterations ($p = 1, 2, \dots$) the influence of a priori terms in the constrained non-linear iteration are minor and the constrained iterations do not significantly differ from a non-constrained iteration. This is especially clear from Eq. 2.1.8, 2.1.9 which are written via weighting matrices and shown explicitly for single-pixel inversion in Eq. 2.1.17. Indeed, if the initial guess is far from the solution, the values ($f_j(a^p) - f_j^*$) are large and the measurement term dominates in the minimized quadratic form $\Psi(a^p)$ over a priori terms because the values of Lagrange multipliers $\gamma_2$ and $\gamma_3$ (see Eq. 2.1.9) are typically small. A priori terms start to matter only when fitting differences ($f_k(a^p) - f_k^*$) reach the level of measurements accuracy characterized by the standard deviation $\epsilon_1$. Therefore, some increase of a priori terms at initial iterations may improve the performance of Eq. 2.0.6, 2.0.7 and from Eq. 2.1.8, 2.1.9. This idea is exploited by the following considerations. Each p-th iteration in Eq. 2.1.8, 2.1.9 assumes the solution of the following overdetermined linear system:

$$\begin{align}
f_k (\Delta a^p) &= f_k^*(a^p) - f_k^* + \Delta f_k^* \Rightarrow \nonumber \\
f_k (\Delta a^p) &\approx K_{k,p} (\Delta a^p) = f_k^*(a^p) - f_k^* + \Delta f_k^{mes} + \Delta f_k^{lin}, \nonumber \\
(k = 1, \dots, K) \tag{2.2.5}
\end{align}$$

where $K_{k,p}$ is the Jacobi matrix of the first derivatives from $f_k(a)$ in the vicinity of $a^p$ and $\Delta f_k^{lin}$ denotes the errors of using the linear approximation of ($f_k(a^p) - f_k^*$) in the vicinity of $a^p$.

The difference of Eq. 2.2.5 with Eq. 2.1.2 is that it includes linearization errors $\Delta f_k^{lin}$ and therefore accounting for $\Delta f_k^{lin}$ can be introduced into Multi-term LSM. If $\Psi'(a^p)$ is defined using weighting matrices (Eq. 2.1.8, 2.1.9) then using $\epsilon_1^2 + (\epsilon_{1,lin})^2$ instead of $\epsilon_1^2$ in Lagrange multipliers definition can be used for optimizing the contributions of measurements and a priori terms. The value of the $\Delta f_k^{lin}$ variance is not known in each point $a^p$, but can be estimated from the value of the residual, i.e., analogously to Eq. 2.1.11 one can write:

$$\tilde{\epsilon}_1^2 + (\tilde{\epsilon}_{1,lin}(a^p))^2 \approx \frac{2\Psi'(a^p)}{\sum_{k=1,...,K} N_k - N_a} = \frac{2\Psi'(a^p)}{N_f^* + N_a - n}. \quad (2.2.6)$$

Using this equation, Eq. 2.1.9 can be re-written for non-linear iterations:

$$\gamma_k(a^p) = \frac{\tilde{\epsilon}_1^2 + (\tilde{\epsilon}_{1,lin}(a^p))^2}{\epsilon_k^2} \approx \frac{2\Psi'(a^p)}{\epsilon_k^2(\sum_{k=1,...,K} N_k - N_a)}. \quad (2.2.7)$$

This definition of Lagrange multiplier accounts for higher linearization errors at earlier iterations. In the close vicinity of the solution $a'$, where $\Delta f_k^{lin}$ is close to zero, Eq. 2.2.7 is reduced to Eq. 2.1.9. Hence, utilizing “adjustable” $\gamma_k$ ($k \ge 2$) in Eq. 2.1.9 improves the convergence while the final solution retains the same statistical properties.

Practice of using GRASP has shown that the convergence optimization employing adjustable Lagrange multipliers by
**Eq. 2.2.7** is often quite efficient. At the same time, it is not always sufficient; for example, there are situations when there is no linear or not at all a priori constraints in Eq. 2.1.13 and correspondingly in Eq 2.1.14, 2.1.17, etc. Nonetheless, it is clear that $\Delta\mathbf{a}^P$ should be limited, especially at the initial iterations when the linearization error is the largest. In such case, the determination of $\Delta\mathbf{a}^P$ in the iterative procedure can be considered as the solution of the following modified (compared to Eq. 2.2.5) system as:

$$ \left\{ \begin{array}{l} \mathbf{K}_{k,p} \Delta \mathbf{a}^p = \mathbf{f}_k^* (\mathbf{a}^p) - \mathbf{f}^* + \Delta \mathbf{f}^* + \Delta \mathbf{f}_p^{\text{lin}} \\ \Delta \mathbf{a}^{p,*} = \mathbf{0}^* + \Delta \mathbf{a} \end{array} \right. \qquad (2.2.8) $$

Correspondingly, using this additional requirement on limiting $\Delta\mathbf{a}^P$, an additional term will be introduced in Eq. 2.1.8:

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

where matrix $\mathbf{D}_{p,\Delta a}$ is diagonal matrix with the elements:

$$ \{\mathbf{D}_{p,\Delta a}\}_{ii} = \gamma_{p,\Delta a_i} = \frac{\varepsilon_{\Delta a_i}^2}{\varepsilon_{\Delta a_i}^2}. \qquad (2.2.10) $$

The variance $\varepsilon_{\Delta a_i}^2$ 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}^2$, i.e., $a_i, \max - a_i, \min \approx 3\varepsilon_{\Delta a_i}$. Also, in practical applications, the impact of the correction $\Delta\mathbf{a}^P$ is always scaled by a factor $t_p$: $\mathbf{a}^{p+1} = \mathbf{a}^{p+1} - t_p \Delta\mathbf{a}^P$, similarly as it is applied in Eq. 2.2.4.

Thus, the inversion in case of non-linear $f_k(\mathbf{a})$ and/or $f_i^a(\mathbf{a})$ is implemented within the frame of Eq. 2.0.6, 2.0.7 and include Levenberg-Marquardt like optimization of convergence using coefficient $t_p$ shown in Eq. 2.2.4. Also, for optimizing the convergence of solution weights of the a priori terms can be enhanced at earlier iterations (as shown in Eq. 2.2.7 and additional term $\mathbf{D}_{p,\Delta a}$ can be used as shown in Eq 2.2.9, 2.2.10.

## 2.2.3 Consistency of Multiple a Priori Constraints

The simultaneous utilization of multiple a priori constraints certainly relies on the assumption that all a priori constraints are consistent and do not contradict each other. For example, if a group of parameters $a_i$ ($i = 1, ..., N_i$) represent a physical function $a_i = y(x_i)$, one cannot use an irregular (highly unsmooth) a priori estimate function $a_i^* = y^*(x_i)$ as a priori estimate and simultaneously apply smoothness constraints on the retrieved solution. Obviously, if all a priori constraints are defined based on readily available data or information, no contradictory constraints should appear. However, if some of the constraints are set based on some ill-defined considerations, some erroneous a priori assumptions could occur in practical situations.

In order to detect the possible anomalies in the a priori assumptions or appearance of strong biases in the forward modeling, the value of the remaining miss-fit function $\Psi(\mathbf{a})$ given by Eq. 2.1.5 should be monitored. As a result of the inversion its resulting value should agree with the expected

noise assumptions and converge as provided by Eq 2.1.7, 2.1.11. In case of using Multi-term LSM formulation, the value $\Psi_k(\mathbf{a})$ is a sum of several terms:

$$ \begin{aligned} 2\Psi(\mathbf{a}) &= 2\sum_{k=1}^{K} \Psi_k(\mathbf{a}) = \sum_{k=1}^{K} (\mathbf{f}_k(\mathbf{a}) - \mathbf{f}_k^*)^\mathrm{T} \mathbf{C}_k^{-1} (\mathbf{f}_k(\mathbf{a}) - \mathbf{f}_k^*) \\ &= \min \rightarrow \sum_{k=1}^{K} N_k - n. \end{aligned} \qquad (2.2.11) $$

The value of each term is expected to reach its asymptotic limit:

$$ \Psi_k(\mathbf{a}) \approx N_k \frac{\sum_{k=1}^{K} N_k - n}{\sum_{k=1}^{K} N_k} \xrightarrow{\left(\sum_{k=1}^{K} N_k \gg n\right)} N_k. \qquad (2.2.12) $$

Obtaining significantly different values of $\Psi_k(\mathbf{a})$ from the expected value based on the assumptions made should considered as an indicator of anomaly that is likely related with in incorrect assumption made about known covariance matrix $\mathbf{C}_k$ or with other unidentified problems. Indeed, both overestimation and underestimation of accuracy of a priori assumption may have negative effects on the retrieval results. For example, if accuracy of a priori information is overestimated the desired accuracy of observation fit may not be achievable. On the other hand, if a priori information underestimated, the corresponding a priori terms would become very small and overall miss-fit function may have reasonably small values while some other terms can be under-fitted. In such situation the overall retrieval results and/or estimations of the accuracy of these obtained results may be inadequate.

## 2.2.4 Accounting for Measurements Redundancy

Accounting for data redundancy is another complex issue in implementing optimized inversion that is somewhat related to the consistency of the assumptions. This issue is of high practical importance, although it is not often addressed in inversion methodologies. For example, infinite enhancement of spectral or/and angular resolutions in remote sensing observations will not lead to accuracy improvements in the retrievals above a certain level. Based on common sense a simple increase in the number of observations $N_k$ leads to an increase in number of redundant measurements that do not help to improve the retrievals. However, theoretical considerations do not assume any “redundant” or “useless” observations. Performing N straightforward repetitions of the same observation (with established unchanged accuracy), simply means that the variance of this particular observation $f_i$ should decrease by a factor of N. Accordingly, the j elements of the corresponding covariance matrix C should decrease and the errors of the retrieved parameters should decrease appropriately, which contradicts to sensibleness. For the Multi-term LSM approach accounting for data redundancy is of particular importance. Indeed, individual data points from observations of the same type are usually comparable in accuracy. Therefore, it is unlikely that inverting single source of data should lead to a discrimination of some individual observations. In a multi-source inversion, the situation is different because an increase of the number of the observations in one of several inverted data
sets would lead to an increase in the weight of this data set even if the added observations were redundant from a practical point of view. Indeed, in the minimized quadratic form of Eq. 2.2.11 the higher the value of the k-th term $\Psi_k$, the stronger the contribution of k-th data set into the solution. In order to eliminate this obvious dependence of $\Psi_k$ on $N_k$, Dubovik and King (2000) suggested that the accuracy of a single measurement degrades as $1/N_k$ for redundant observations if several measurements are taken simultaneously, i.e.;

$$ \varepsilon_k^2(\text{multiple}) = N_k \varepsilon_k^2(\text{single}), \qquad (2.2.13) $$

where the term “multiple” indicates that several analogous measurements are taken simultaneously. Correspondingly, Eq. 2.1.9 can be written via accuracy of “single” measurement as follows:

$$ \gamma_k = \frac{N_1 \varepsilon_1^2(\text{single})}{N_k \varepsilon_k^2(\text{single})}. \qquad (2.2.14) $$

This definition of $\gamma_k$ makes the relative contributions of the terms $\gamma_k \Psi_k$ in Eq. 2.2.7 independent of $N_k$, and therefore equalizes the data sets with different numbers of observation. The relationship (2.2.13) assumes that $\varepsilon_k$ increases as $\sqrt{N_k}$ for data sets with “redundant” observations. Such an increase can be caused by the fact that the number of various sources of random errors may increase proportionally to a number of simultaneous measurements. Indeed, an excessive increase of spectral or/and angular resolution in remote sensing measurements may likely to results into decrease of the quality of each single measurement due to the increased complexity of the instrumentation and its calibration, in a sense that each single measurement taken as a part of the such complex observation likely has overall lower accuracy that if it would be taken alone or as a part of smaller set of observation. However, the assumption given by Eq. 2.2.13 is of intuitive character and may need to be verified in extensive practical experience.

Thus, identification of the measurements redundancy in practice is a difficult effort that strongly relies on the experience of the developer. Nevertheless, it can be advisable to consider data redundancy as a practical factor that may affect the retrieval. Namely, if Eq. 2.2.13 gives a value much higher than the level of expected measurement errors, then it is likely that noise assumptions need to be verified. In such case the ratios $N_1/N_k$ can be good indications of magnitude and directions of required adjustments in $\varepsilon_k^2$ in order to address the possible excessive domination of the large inverted data sets over smaller ones. For example, the assumption (2.2.13) was successfully employed in aerosol remote sensing retrievals (Dubovik and King, 2000), where it helped to harmonize the contribution of large sets of angular sky radiance measurements with much fewer observations of spectral optical thickness. A similar principle was used in earlier studies (Oshchepkov and Dubovik, 1993).

### 2.2.5 Asymptotical Character of MML and LSM Optimum Properties

The Least Squares Method (LSM) is probably the most frequently used numerical procedure for implementing the statistically

optimized solution of a system of equations. The LSM is used in diverse applications and is the basis of a vast family of methods interpreting indirect observations. As a result, most of the practical equations for solutions and error estimates used in the algorithms are directly or indirectly related to the LSM formalism. At the same time, LSM is a specific case of MML when the random noise in the data has a Gaussian distribution.

According to MML, the best estimates $\hat{a}$ of unknowns correspond to the maximum of likelihood function (PDF) $P(f(a)|f^*)$. The obtained solution $\hat{a}$ is statistically optimal in many senses as discussed in the textbooks devoted to fundamental statistics (Fourgeaut and Fuchs 1967; Edie et al., 1971). For the MML solution:

- $\hat{a}$ is asymptotically non-biased $\langle \hat{a} \rangle \rightarrow a^{\text{real}}$;

- $\hat{a}$ is asymptotically consistent $\rightarrow a^{\text{real}}$;

- $\hat{a}$ is asymptotically efficient, i.e., variance of $\hat{a}$ converges to the smallest possible value;

- $\hat{a}$ has asymptotically normal distribution $\hat{a}_N \xrightarrow[N \to \infty]{N(a^{\text{real}}, I_{\hat{a}}^{-1})}$, where $I_{\hat{a}}$ is Fisher information matrix.

It should be noted however, that ML method has above optimal characteristics if certain conditions are satisfied, e.g., first and second derivatives of the functions should exist and be finite (e.g., see Fourgeaut and Fuchs, 1967; Edie et al., 1971). For example, importance of satisfactions of these conditions in practice was outlined in Section 2.2.1.

In addition, the MML solution keeps many optimal characteristics even when there are a limited number of observations. The optimum properties of MML are closely connected with the Fisher information determination. For instance, in the case of several estimated parameters ($a$ is vector-column) Fisher information matrix is used. This matrix is formulated as a covariance matrix of logarithmic partial derivatives of PDF; i.e., the Fisher matrix has the following elements:

$$ \{I_{\hat{a}}\}_{ii'} = \int_{R_a} \dots \int_{R_a} \frac{\partial \ln P}{\partial a_i} \frac{\partial \ln P}{\partial a_{i'}} P da_1 \dots da_{N_a} = \left\langle \frac{\partial \ln P}{\partial a_i} \frac{\partial \ln P}{\partial a_{i'}} \right\rangle_{R_a} \qquad (2.2.15) $$

where $P = P(f(a)|f^*)$ and the integration is done over all space $R_a$ of possible $a$. Here only the case of unbiased estimation is considered (discussion including biased cases can be found elsewhere, e.g., Edie et al., 1971). The Fisher matrix defines the accuracy limits for estimates $\hat{a}$. Namely, if we use vector $\hat{a}$ for estimating any other value $b$ linearly dependent on $a$ (i.e., $b = ga$, $g$ is a vector-line of the coefficients - first derivatives), we have the estimates $\hat{b} = b^{\text{real}} + \Delta b$ with the variance limited by the Cramer-Rao inequality:

$$ \langle (\Delta b)^2 \rangle = g C_{\hat{a}} g^T \geq g I_{\hat{a}}^{-1} g^T, \qquad (2.2.16) $$

where $C_{\hat{a}}$ is the covariance matrix of $\hat{a}$. The MML produces the vector $\hat{a}$ of estimates $\hat{a}_i$ which are jointly effective in the sense that their covariance matrix $C_{\hat{a}}$ corresponds to the equality in Eq. 2.2.16. Thus, Fisher information has a clear practical meaning
and it can be easily evaluated using MML estimates covariance matrix (as it is given by Eq. 2.2.15), i.e.,

$$I_{\hat{a}} = C_{a, \text{MML}}^{-1}. \qquad (2.2.17)$$

There are some exceptions when Eq. 2.2.17 is not correct, but those cases are rather of interest for theoretical investigations than for practice [for the details and discussion see Edie et al. (1971)]. It should be noted, that the Fisher definition of information is more suitable for using in the optical and remote sensing fields than the popular Shanon's definition. For instance, the entropy (information in the random value about itself) of each particular estimate $\hat{a}_i$ is often considered for optimization of inversion algorithms; namely, the solution with the maximum of entropy is accepted as optimum. Let us consider the conceptual definitions of entropy in Shanon's and Fisher's methodologies. The entropy according to Shanon [details see in Shannon and Weaver (1949) or in review by Kolmagorov (1987)] is defined as

$$\hbar(P(\hat{a}_i)) = \int (-\log_2 P(\hat{a}_i)) P(\hat{a}_i) d\hat{a}_i \rightarrow N_{\text{bits}}. \qquad (2.2.18)$$

where $P(\hat{a}_i)$ is PDF of $\hat{a}_i$, $N_{\text{bits}}$ is the number of bits (binary digits) needed to represent the number of distinct estimates that could have been obtained. This formalism of information quantity uses logarithm with the base 2 and has a very suitable meaning for digitalization of information. This Shanon information formalism is widely used in the modern remote sensing and applied optics literature for evaluating information content of the measurements (e.g. see Peckham 1974; Rodgers 2000; Purser and Huang 1993, etc.). However, the experimental physics, the meaning given by Fisher information formalism sounds more useful. According to this formalism the entropy of $\hbar(P(\hat{a}_i))$ shows the accuracy of $\hat{a}_i$ determination and instead of Eq. 2.2.18 one can write:

$$\begin{aligned} \hbar(P(\hat{a}_i)) &= \int \left( \frac{\partial \ln P(\hat{a}_i)}{\partial \hat{a}_i} \right)^2 P(\hat{a}_i) d\hat{a}_i \\ &= \int \left( -\frac{\partial^2 \ln P(\hat{a}_i)}{\partial \hat{a}_i^2} \right) P(\hat{a}_i) d\hat{a}_i = \frac{1}{\langle (\Delta \hat{a}_i)^2 \rangle}. \end{aligned} \qquad (2.2.19)$$

Correspondingly, if we consider all possible estimates $\hat{a}_i$ obtained from initial data $f^*$, the estimate given by MML will have the maximum entropy. (NOTE that the second derivative of PDF is used in Eq. 2.2.19 for the emphasis of the similarity to Eq. 2.2.18, although in mathematical sense this definition is equal to Eq. 2.2.15 (see Linnik 1962; Edie et al., 1971).

It should be noted that all above properties of MML and LSM have asymptotical character. Strictly speaking this means that the measurement $f^*$ is repeated N time for observing exactly the same situation characterized by parameters $a_{\text{real}}$, and with $N \rightarrow \infty$ the asymptotical properties can be observed. Correspondingly, the PDF $P(\hat{a} | f^*)$ is a composition of $(N \rightarrow \infty)$ PDF of single observations $f_i^*$:

$$P(\hat{a} | f^*) = P(\hat{a} | f_1^*) P(\hat{a} | f_2^*) \dots P(\hat{a} | f_N^*) = \prod_{i=1, \dots, N} P(\hat{a} | f_i^*). \qquad (2.2.20)$$

In practical remote sensing, it is very difficult to meet the situation when exactly the same measurements can be repeated large number of times for unchanging atmosphere. For

example, in satellite remote sensing, a single observation would correspond only to the one state of the atmosphere, and any next satellite overpass over same site would strictly speaking observe a different atmosphere. In such situation, one can probably consider that each satellite measurement $f^*(x_j, y_k, t_i, \dots)$ is already a results of large number of records of electromagnetic interaction acts between the sensor and the atmosphere, to the extend that the measurement errors were significantly averaged in $f^*(x_j, y_k, t_i, \dots)$ to be inverted and the situation meet the requirements necessary for approaching asymptotical limits.

On the another hand if one studies the climatology of the atmospheric state $a(x,y,t, \dots)$ from extensive observations, it is possible to consider that the asymptotical properties can be expected in the resulting climatology that is based on the well-established extensive record of the observations $f^*(x_i, y_j, t_k, \dots)$ obtained in generally the same conditions of measurement noise formation:

$$P(\hat{a} | f^*(x, y, t, \dots)) = \prod_{N_{(i,j,k)} \rightarrow \infty} P(\hat{a} | f^*(x_i, y_j, t_k, \dots)). \qquad (2.2.21)$$

## 2.2.6 Optimization of Time Performance (Jacobian Calculations, etc.)

Some practical features were implemented in the GRASP algorithm for optimizing speed of the retrieval. A particular attention was devoted to the optimization of the calculations of the Jacobians that is the most time consuming component of Newton's retrieval algorithms. For example, for those parameters that didn't change at the previous iteration the Jacobians may not need to recalculated at each iteration (as this realized in GRASP for some applications). Also, as was discussed by Dubovik and King (2000) and Dubovik et al. (2011) the successful retrieval can be achieved by using the approximate and quick calculations of the Jacobians. In this respect GRASP can implement the calculations of atmospheric radiances with different accuracy, for example, by using analytical single-scattering approximation or multiple scattering regime of different accuracy [e.g., by adjusting the number of the terms in the expansion of the phase matrix and in the quadrature of directional integration, as discussed by Dubovik et al. (2011)]. Correspondingly, the simulation of Jacobians can be executed in GRASP with lower accuracy using faster calculations.

Additionally, as mentioned above the implementation of the PC analysis in the multi-pixel formalism (Xu et al., 2019) will be integrated into GRASP. This approach can also reduce significantly the calculation time in some situations, e.g., in processing of high spatial resolution satellite observations.

## 2.3 Error Estimates

Estimations of the retrieval errors in GRASP are based on LSM equations expressed for the case of Multi-term solutions written via weighting matrixes. Both the contribution of random and systematic error components are estimated as follows (see Dubovik, 2004):
$$C_{\Delta\hat{a}} = C_{\Delta\hat{a}_{\text{ran}}} + \hat{a}_{\text{bias}}\hat{a}_{\text{bias}}^T \quad (2.3.1)$$

$$C_{\Delta\hat{a}_{\text{ran}}} = \langle \Delta\hat{a}_{\text{ran}} (\Delta\hat{a}_{\text{ran}})^T \rangle = \left( \sum_{k=1}^{K} \mathbf{K}_k^T \mathbf{C}_k^{-1} \mathbf{K}_k \right)^{-1} \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 (2.3.2)$$

$$\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 (2.3.3)$$

where Jacobians $\mathbf{K}_k$ are calculated in the small vicinity of the point $\mathbf{a}^{\text{solution}}$, $\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 miss-fit of the data using Eq. 2.1.11.

Estimation of not only random retrieval error but also error retrieval bias $\hat{a}_{\text{bias}}$ is very important for the adequate evaluation of retrieval uncertainty, especially when multiple a priori constraints are used. For example, for the case of single-pixel retrieval given by Eq. 2.1.17 $C_{\Delta\hat{a}_{\text{ran}}}$ is expressed as:

$$C_{\Delta\hat{a}_{\text{ran}}} \approx \left( \sum_{k=1}^{K_1} \gamma_k \mathbf{K}_k^T \mathbf{W}_k^{-1} \mathbf{K}_k + \gamma_a \mathbf{W}_a^{-1} + \sum_{n=1}^{K_2} \gamma_n \mathbf{\Omega}_n \right)^{-1} \hat{\varepsilon}_1^2. \quad (2.3.4)$$

Analyzing this equation one can see rather obvious tendency: the higher the contributions of the second and the third terms the smaller the random errors are. Correspondingly, the more a priori constraints are used the lower the random errors of the retrieval. However, in practice a priori constraints can be unintentionally inadequate and introduce some systematic uncertainties, i.e., biases. In principle, there is no guaranteed approach for detecting those biases unless 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(\mathbf{a}^{\text{solution}}) - f_k^*$ or misfit of a priori constraints. For example, for Eq. 2.1.17 the bias can be introduced by a priori estimate $\hat{a}_{\text{bias}}^* = (\mathbf{a}^{\text{solution}} - \mathbf{a}^*)$ or unsmooth features in the retrieved solution: $\mathbf{a}_{n,\text{bias}}^{\text{smooth}} = \mathbf{\Omega}_n \mathbf{a}^{\text{solution}} \neq 0$. Correspondingly, the bias for single-pixel retrieval is estimated as:

$$\hat{a}_{\text{bias}} \approx \left( \sum_{k=1}^{K_1} \gamma_k \mathbf{K}_k^T \mathbf{W}_a^{-1} \mathbf{K}_k + \gamma_a \mathbf{W}_a^{-1} + \sum_{n=1}^{K_2} \gamma_n \mathbf{\Omega}_n \right)^{-1} \left( \sum_{k=1}^{K_1} \gamma_k \mathbf{K}_k^T \mathbf{W}_a^{-1} \Delta f_k^{\text{bias}} + \gamma_a \mathbf{W}_a^{-1} \hat{a}_{\text{bias}}^* + \sum_{n=1}^{K_2} \gamma_n \mathbf{a}_{n,\text{bias}}^{\text{smooth}} \right). \quad (2.3.5)$$

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 OE approaches) for constraining the retrieval. In a similar way, the a priori biases are estimated in the case when multi-pixel a priori constraints are used.

The Levenberg-Marquardt optimization of the convergence, discussed in Section 1.3.2 may also introduce a bias. Indeed, this optimization makes the iterations converge from given initial guess to fit the data even if the basic linear system is singular. Therefore, once 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 Eq 2.3.2-Eq 2.3.3 are modified as follows when Levenberg-Marquardt optimization is used in the retrieval, where is defined by Eq 2.2.9, 2.2.10:

$$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}_{p,\Delta a} \right)^{-1} \hat{\varepsilon}_1^2, \quad (2.3.6)$$

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}_{p,\Delta a} \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 (\mathbf{a}^{\text{solution}} - \mathbf{a}^{p=0}) \right). \quad (2.3.7)$$

Finally, in practice the users may not need directly the retrieved parameters but their functions $m(a)$ that can be calculated from the retrieved parameters. For example, GRASP retrieves parameters of aerosol microphysics (particle sizes, refractive indices, etc.) but users need Aerosol Optical Depth AOD. For such situation, GRASP is designed to provide a set of such diverse indirect characteristics with the possibilities of providing the uncertainties calculated as:

$$C_{\Delta\hat{m}} \approx \mathbf{M}(C_{\Delta\hat{a}_{\text{ran}}} + \hat{a}_{\text{bias}}\hat{a}_{\text{bias}}^T)\mathbf{M}^T = \mathbf{M}C_{\Delta\hat{a}_{\text{ran}}}\mathbf{M}^T + \mathbf{M}\hat{a}_{\text{bias}}(\mathbf{M}\hat{a}_{\text{bias}})^T \\ = C_{\Delta\hat{m}_{\text{ran}}} + \hat{m}_{\text{bias}}\hat{m}_{\text{bias}}^T, \quad (2.3.8)$$

where $\mathbf{M}$ - is the matrix of first derivatives $\{\mathbf{M}\}_{ji} = \frac{\partial m_j}{\partial a_i}|_{a_i^{\text{solution}}}$, $\mathbf{m}_j$ is j-th element of vector $\mathbf{m}$.

Finally, the effect of biases in the measurements on the solution bias $\hat{a}_{\text{bias}}$ is accounted for in Eq. 2.3.5 based on the assumption that the presence of biases is manifested in the non-zero miss-fits $\Delta f_k^{\text{bias}}$. Indeed, in many cases when systematic errors are present in the inverted measurements or in the physical forward model used in the inversion, the accurate matching of inverted data can't be achieved [e.g., see illustrations provided by numerical sensitivity tests for AERONET retrievals by Dubovik et al. (2000)]. At the same time, there are many situations when biases in the measurements may not significantly affect the residual (Eq. 2.1.11) and the miss-fits $\Delta f_k^{\text{bias}}$. For example, it is rather evident that the retrieval of aerosol SSA from AERONET ground-based measurements is highly sensitive to the calibration biases in the direct Sun measurements, while the fitting of these direct measurements is always quite accurate [see 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 of AERONET retrieval errors for the operational products is estimated using a series of ~27 numerical proxy inversion tests with the sets of perturbations in both input measurements and auxiliary input parameters (Sinyuk et al., 2020). A similar strategy can be used for evaluation of potential effects of undetected biases. Specifically, the bias term $(\hat{a}_{\text{bias}}) (\hat{a}_{\text{bias}})^T$ in Eq. 2.3.1 can be estimated as:
<figure><img src="image_7.png" /><figcaption>FIGURE 4 | General organization of forward model and its connection with the numerical inversion in the GRASP algorithm.</figcaption></figure>

$$a_{bias} a_{bias}^{T} \rightarrow \langle a_{bias} a_{bias}^{T} \rangle_{bias \text{ proxy set}} \quad (2.3.9)$$

where the values of the retrieval biases are estimated as an average effect from a preselected set of possible biases in measurements and auxiliary inputs.

## 3 GRASP FORWARD MODEL

The “forward model” module in the code implements simulations of the inverted remote sensing observations. The GRASP “forward model” is rather universal, i.e., can simulate a large variety of remote sensing observations (passive and active observations obtained from ground and space). The technical realization of different parts of forward model is already discussed in details in precedent publications (e.g., Dubovik, et al., 2011; Lopatin et al., 2013; Lopatin et al., 2021; Torres et al., 2017; Derimian et al., 2016; Li et al., 2019, etc.), therefore in this Section the discussion will be focused on different organizational aspects of the overall structure of GRASP forward model with the objective of showing users the main retrieval possibilities that have been already realized or can be relatively easy added to the software if needed. GRASP forward model consists from several distinct modules (Figure 4): Multiple Scattering, Aerosol single scattering columnar/volume properties, Aerosol vertical profile, Surface reflectance and Gas absorption calculations. These blocks are semi-independent in the sense that each block can be changed or entirely replaced with no or minimal effect on the other parts of the “forward model” routine. For example, GRASP “forward model” allows for the choice of phenomenological and physical approaches/models used for simulating surface reflectance (see Dubovik et al., 2011). As will be discussed below there is a very high flexibility in setting aerosol models. Moreover, even the radiative transfer block accounting for multiple scattering has a

rather standardized inputs/outputs and can be replaced if required by other radiative transfer routines. At the same time, each module has some very convenient customized features that may not be available in other analogous software packages and scientific codes.

Radiative transfer calculation accounting for multiple scattering effects in GRASP is implemented by on-line radiative transfer calculations using 1-D Successive Order of Scattering method. The scientific basis of this approach is documented in the paper by Lenoble et al. (2007). At the same time, some custom features are implemented in the radiative transfer module of GRASP as described by Dubovik et al. (2011). For example, the truncation procedure described by Waquet and Herman (2019) has been realized. The truncation procedure consists in removing the forward scattering peak observed in the phase function of large particles (few microns) or cloud droplets, due to diffraction. This approximation allows faster calculations for both total and polarized radiances and, as shown by Waquet and Herman (2019) maximal errors due to the used approximations do not exceed 0.001 for the degree of linear polarization for optical thickness smaller than 2.0, which is a sufficient accuracy for the most applications based on polarimetric measurements. In addition, several possibilities of using a number of trade-offs between the accuracy and the speed of the successive order of scattering method has been realized in GRASP (see Dubovik et al., 2011): including the possibilities of changing the number of terms M used in the expansion of the phase matrix into Legendre polynomials, the number of terms N used in Gaussian quadrature for zenith integration, number of numerical layers in vertical atmosphere properties integrations, etc. Moreover, very recently the atmospheric correction approach was realized in the GRASP forward model. This approach decouples the radiative transfer effects of the atmosphere and surface in a rather original way
and provides a simplified but still an accurate way of solving the radiative transfer equation (Litvinov et al., 2018). It provides sufficient accuracy for the atmospheric correction and high flexibility for further developments such as including adjacency effect, surface topology effect as well as the effect of cloud shadowing.

In addition, recently the radiative transfer module of GRASP has been extended by taking into account for the radiation coming from the thermal emission of the different atmospheric components and the surface (Herreras et al., 2020). In this realized implementation of the original 1-D Successive Order of Scattering numerical scheme by Lenoble et al. (2007) the source functions were extended to include Planck emission function. This modification allows for accounting the radiation coming from both scattering and emission origin in single and multiple scattering regimes. This new RT scheme accounting for emitted radiation has been validated against the reference radiative transfer code based on the discrete-ordinated method by Stamnes et al. (1988). The difference in the radiance, expressed in terms of brightness temperature, obtained from both codes under the same conditions is on average below 0.006 K. This difference is substantially below of the radiance accuracy of the most common instruments operating in thermal IR spectral range, like IASI (Blumstein et al., 2004) or CALIPSO IIR (Garnier et al., 2012) with an accuracy around 0.1 K. Thus, this realized RT development opened opportunity to apply GRASP approach for realizing atmospheric retrieval using both photometric and hyperspectral observations in visible and thermal IR spectral range.

The total surface reflectance matrix BRDM in GRASP code is modelled using Surface reflectance BRDF and BPDF. Specifically, as described in detail by Dubovik et al. (2011), 4 by 4 BRDM matrix is modelled using two terms $M = M_{\text{diff}} + M_{\text{spec}}$, where diffuse unpolarized term $M_{\text{diff}}$ has only one non-zero element $\{M_{\text{diff}}\}_{11}$ modelled using semi-empirical BRDF function and specular reflection term is represented by Fresnel reflection function scaled by selected BPDF coefficient. Nevertheless, GRASP can operate also with more physically based models of the BPDF (Litvinov et al., 2012) when physical constraints can be imposed by the surface structure and composition. In general, both BRDF and BPDF can be calculated in GRASP using a variety of subroutines representing different models (Dubovik et al., 2011; Litvinov et al., 2011a; Litvinov et al., 2011b; Litvinov et al., 2012). The retrieval relies on water-land mask and uses either models of land or water surface reflectance. In the mixed observed pixels, both the properties of land and water surface reflectance models are retrieved and summed up in the 1-D radiative transfer using land/water fraction.

For the ocean surface the reflection is mainly governed by the wind speed at sea level as suggested by the Cox-Munk model (Cox and Munk, 1954). This model is employed in most applications including polarimetric observation (e.g., Deuzé et al., 1999; Deuzé et al., 2001; Herman et al., 2005; Xu et al., 2016; Xu et al., 2017, etc.) and also used in GRASP (see Dubovik et al., 2011 and some details in; Fougnie et al., 2019).

In a contrast, the reflection matrix from land surfaces may differ very strongly from location to location. Therefore, in the many algorithms interpreting satellite observation over land, the key aspect is correct determination of appropriate surface reflectance model and appropriate parameters. A number of BRDF models developed for surface reflectance description from remote sensing measurements are included into GRASP algorithm: Rahman-Pinty-Vestarte (RPV) model (Rahman et al., 1993), kernel-driven semi-empirical models (Ross-Li sparse, Ross-Li dense, Ross-Roujen models (Ross, 1981; Li and Strahler, 1992; Roujean et al., 1992; Wanner et al., 1995), physically-based models for bare soil and vegetated surfaces (Litvinov et al., 2012) as well as physically-based models for snow and ice (Kokhanovsky and Zege, 2004; Kokhanovsky and Breon, 2012). For polarimetric remote sensing these BRDF models are combined in GRASP algorithm with models for surface polarized reflectance (BPDF models). Most of the existent BPDF models are based on the Fresnel equations of light reflection from the surface. For example, Nadal and Bréon (1999) have proposed simple two-parameter non-linear function of the Fresnel reflection for characterization of atmospheric aerosol over land surfaces based on POLDER observations of land surface reflectance. Maignan et al. (2009) have introduced a new linear BPDF model with only one free parameter and demonstrated that this simple model allows a similar fit to the POLDER measurements as a more complex non-linear model by Nadal and Bréon (1999). For accurate description of polarimetric measurements like RSP (Research Scanning Polarimeter) airborne instrument, a three-parameter semi-empirical model was proposed by Litvinov et al. (2011a) and Litvinov et al. (2011b) (e.g., Rondeaux and Herman 1991, Nadal and Bréon, 1999; Maignan et al., 2004; Maignan et al., 2009; Waquet et al., 2009a; Waquet et al., 2009b; Litvinov et al., 2010; Litvinov et al., 2011a; Litvinov et al., 2011b). Analysis of number of airborne polarimetric measurements show minor spectral dependence of polarized reflectance in the visible and infrared spectral (e.g., Rondeaux and Herman 1991, Nadal and Bréon, 1999; Maignan et al., 2004; Maignan et al., 2009; Waquet et al., 2009a; Waquet et al., 2009b; Litvinov et al., 2010; Litvinov et al., 2011a; Litvinov et al., 2011b). At present time the parametrization of surface polarized reflectance with Fresnel-based BPDF models looks acceptable for existent polarimetric space-borne VIS and NIR observations.

BRDF and BPDF models included in GRASP are capable to reproduce reasonably the surface total and polarized reflectance (Maignan et al., 2004; Maignan et al., 2009; Litvinov et al., 2011a; Litvinov et al., 2011b) and have already been used for interpreting observations by MISR, MODIS, POLDER and other instruments (Justice et al., 1998; Martonchik et al., 1998; Govaerts et al., 2010; Wagner et al., 2010).

The actual choice of BRDF or BPDF models in practical application and retrieval approaches is always a subject of the discussion and, very often, depends on user preference. For global processing of different remote sensing measurements (PARASOL, MERIS, OLCI, S5p/TROPOMI and other) with GRASP algorithm the optimal balance between speed, accuracy linearity and number of parameters was provided
<figure><img src="image_8.png" /><figcaption>**FIGURE 5** | Illustration of modeling properties of each aerosol component in GRASP algorithm that is represented as a mixture of homogeneous particles that may have different sizes, spheroidal shapes, complex refractive indices and vertical profiles of concentrations.</figcaption></figure>

by the combination of Ross-Li sparse BRDF model (Ross, 1981; Li and Strahler, 1992; Wanner et al., 1995) and one parametric Maignon-Breon model (Maignan et al., 2009). Nevertheless, other possible combinations of different BRDF and BPDF models are possible in GRASP algorithm and are always the subject of the studies on increasing retrieval performance.

Aerosol single scattering module is probably the most elaborated module of the GRASP algorithm that proposes a really large spectrum of different possibilities of approaching modeling optical properties of atmospheric aerosols. Indeed, in spite of the fact that GRASP is generally pursuing the retrieval of both aerosol and surface properties, it deeply relies on the heritage of aerosol retrieval advances (Dubovik et al., 1995; Dubovik et al., 2000; Dubovik and King, 2000; Dubovik et al., 2002; Dubovik, 2004; Dubovik et al., 2006) implemented for AERONET (see Holben et al., 1998) a worldwide network of currently over 600 radiometer sites that generate the data used to validate nearly all satellite observations of atmospheric aerosols. Currently, for all remote sensing applications aerosol can be modeled as a mixture of small polydisperse particles of various shapes and composition (Dubovik et al., 2006). Specifically, the optical properties of homogeneous layer of aerosol are defined by layer scattering and extinction optical thickness and by the elements of the scattering matrix $P_{ii}(\lambda; \Theta; h)$ that can be modeled with one or several aerosol components using microphysical properties of each component. Namely, as illustrated by **Figure 5** each component aerosol component represents a mixture of particles with:

a) different size by defining the volume size distribution $\frac{dV_k(r)}{d\ln r}$;

b) the different shapes using the mixture of randomly oriented spheroids with the distribution of axis or ratios $\frac{dn_k(\epsilon)}{d\ln\epsilon}$;

c) spectral complex index of refraction $n_k(\lambda)$, $K_k(\lambda)$ and

d) vertical profile of aerosol component volume concentration $\frac{dV_k(h)}{d\ln h}$.

Correspondingly, $\tau_{scat/ext}$ and $P_{ii'}(\lambda; \Theta; h)$ of each homogeneous layer are modeled as the following:

$$ \tau_{scat/ext} = \sum_{k=1}^{K} \left( \int_{\ln h_{\min}}^{\ln h_{\max}} \int_{\ln \epsilon_{\min}}^{\ln \epsilon_{\max}} \int_{\ln r_{\min}}^{\ln r_{\max}} \frac{C_{scat/ext}(n_k(\lambda); k_k(\lambda); h; \epsilon; r)}{V(r)} \frac{dV_k(h)}{d\ln h} \frac{dn_k(\epsilon)}{d\ln \epsilon} \frac{dV_k(r)}{d\ln r} d\ln h \, d\ln \epsilon \, d\ln r \right) \quad (3.1) $$

and

$$ \tau_{scat} P_{ii'}(\lambda; \Theta; h) = \sum_{k=1}^{K} \left( \int_{\ln h_{\min}}^{\ln h_{\max}} \int_{\ln \epsilon_{\min}}^{\ln \epsilon_{\max}} \int_{\ln r_{\min}}^{\ln r_{\max}} \frac{C_{ii'}(\Theta; n_k(\lambda); k_k(\lambda); h; \epsilon; r)}{V(r)} \frac{dV_k(h)}{d\ln h} \frac{dn_k(\epsilon)}{d\ln \epsilon} \frac{dV_k(r)}{d\ln r} d\ln h \, d\ln \epsilon \, d\ln r \right) \quad (3.2) $$

Here, as illustrated by **Figure 5**, each component is described by the size distribution $\frac{dV_k(r)}{d\ln r}$, shape distribution $\frac{dn_k(\epsilon)}{d\ln\epsilon}$, spectral complex index of refraction $n_k(\lambda)$, $K_k(\lambda)$ and vertical profile $\frac{dV_k(h)}{d\ln h}$, where $\lambda$ denotes wavelength, $\Theta$ - scattering angle, $h$ - height of the layer, $\epsilon$ - axis ratios of spheroid, and $r$ denotes radius of volume equivalent sphere and $C(n_k(\lambda); K_k(\lambda); h; \epsilon; r)$ is a cross-section of extinction, total or angular scattering of aerosol particle.

All these characteristics of aerosol component can be modeled using approaches of different complexity.

Aerosol size distribution $\frac{dV_k(r)}{d\ln r}$ - can be modeled:

* using a representation of the size distribution as a superposition of triangular $N_i$ bins by retrieving concentration of each bin $\frac{dV_k(r_i)}{d\ln r}$ for each radii $r_i$. This approach is used in AERONET retrieval by Dubovik and King (2000);

* using a representation of the size distribution as a superposition of log-normal $N_i$ bins by retrieving concentration of each log-normal bin $\frac{dV_k(r_i)}{d\ln r}$. This approach is described in Dubovik et al. (2011). Each log-normal bin may have different position and width. This approach is more appropriate when only very limited number of bins is used;

* using a bi-modal log-normal approximation of the size distribution. In this case the parameters of bi-modal log-normal function are retrieved (see Dubovik et al., 2011). This approximation is extensively used by Torres et al. (2017) for inverting measured AOD.

* using the size distribution with the fixed shape of size distribution (using any of above modeling approaches). In this case only total volume of the aerosol component is retrieved.

Aerosol shape distribution $\frac{dn_k(\epsilon)}{d\ln\epsilon}$ - is simulated assuming aerosol as mixture of randomly oriented spheroids using of
pre-computed spheroid kernels software described by Dubovik et al. (2006). The following approaches for modeling shape distribution are realized using detailed shape distribution defined as $\frac{dn_k(\epsilon_i)}{d\ln\epsilon} = q_{k,i}$, where.

* $q_{k,i}$ is a fraction of k-th aerosol component particles with axis $\epsilon_i$.

* $q_{k,i}$ is a fraction of particles of k-th aerosol component with j-th predetermined shape distribution. The definition of two or more fractions can be defined in nearly arbitrary way within the pre-calculated spheroid kernels by Dubovik et al. (2006).

For example, for AERONET (Dubovik et al., 2006) and POLDER (Dubovik et al., 2011) retrievals only two fractions used $q_{k,\epsilon=1}$ - fraction of spherical particles and $q_{k,\epsilon \neq 1} = 1 - q_{k,\epsilon=1}$ - fraction of non-spherical particles (with fixed shape distribution). In that case only one parameter $q_{k,\epsilon=1}$ is retrieved. The non-spherical component of aerosol in that case is simulated using axis ratio obtained by Dubovik et al. (2006) from fitting the desert dust phase matrix measurements by ensemble of spheroids [see detailed description in the paper by Dubovik et al. (2006) and Dubovik et al. (2011)].

Aerosol vertical distribution $\frac{dV_k(h)}{d\ln h}$ can be defined in the following ways:

* the profile is defined as an exponential profile $\frac{dV_k(h)}{dh} \sim e^{-h/\alpha}$. The profile is normalized to unity and scale height $\alpha$ is retrieved parameter.

* the profile is defined as Gaussian profile $\frac{dV_k(h)}{dh} \sim e^{-\frac{1}{2}(h-h_0)^2}$, The profile is normalized to unity, the mean height $h_0$ and standard deviation $\sigma_h$ can be retrieved.

* the profile is defined as a superposition of N$_i$ triangular bins by retrieving concentration of each bin $\frac{dV_k(h_i)}{d\ln h} = C_{h,i}$ for each height $h_i$. This is similar approach as used above for modeling size distribution, more details and illustration of the approach are provided by Lopatin et al. (2013) and Lopatin et al. (2021).

It should be noted that as illustrated in Figure 5 the aerosol shape $\frac{dn_k(\epsilon_i)}{d\ln\epsilon} = q_{\epsilon,i}$ and vertical $\frac{dV_k(h_i)}{d\ln h} = C_{h,i}$ distributions are always retrieved in form of fractions. For example, for vertical profile $\frac{dV_k(h_i)}{d\ln h} = C_{h,i}$ it is always assumed:

$$ \sum_{i=1,...,N_h} C_{h,i} = 1 \qquad (3.3) $$

Using the above constraint, the following $N_h-1$ values are retrieved in GRASP:

$$ a_i = \frac{C_{h_i+1}}{C_{h_i}} (i = 1, ..., N_h - 1) \qquad (3.4) $$

Correspondingly,

$$ C_{h,1} = \frac{1}{\left( 1 + \sum_{i=1,...,N_h-1} a_i \right)} \quad \text{and} \quad C_{h,i+1} = C_{h,1} a_i (i = 1, ..., N_h - 1) \qquad (3.5) $$

The above is written for the case when the profile grid points are constant, i.e.

$$ \Delta_i = \ln h_{i+1} - \ln h_i = \text{const} \qquad (3.6) $$

If $\Delta_i$ is not constant then Eq. 3.5 should be rewritten as

$$ C_{h,1} = \frac{2}{1 + \sum_{i=1,...,N_h-1} a_i (\Delta_i + \Delta_{i+1})} \qquad (3.7) $$

Size distribution can also be retrieved using fractions of particle different sizes:

$$ a_i = \frac{\frac{dV_k(r_{i+1})}{d\ln r}}{\frac{dV_k(r_i)}{d\ln r}} (i = 1, ..., N_r - 1) \qquad (3.8) $$

However, in that case, together with $N_r-1$ $a_i$ parameters the total concentration $C_r$ is also retrieved:

$$ C_r = \sum_{i=1,...,N_r} \frac{dV_k(r_i)}{d\ln r} \qquad (3.9) $$

Retrieving $C_r$ and $N_r-1$ $a_i$ parameters allows for separating the parameter defining the total amount of particles with the set of parameters $a_i$ defining shape of size distribution.

Complex index of refraction: $n(\lambda)$-real and $k(\lambda)$-imaginary parts of the refractive index can be defined using several different strategies:

* the spectral values $n(\lambda_i)$ and $k(\lambda_i)$ are retrieved for each wavelength in similar way as it is realized for AERONET (Dubovik and King, 2000) and POLDER (Dubovik et al., 2011) retrievals;

* the spectral values $n(\lambda_i)$ and $k(\lambda_i)$ are modeled as a homogeneous mixture of different aerosol components mixed internally. At present there are two possibilities in GRASP [see details in Li et al. (2019)]:

i)- several aerosol components with known spectral dependencies $n(\lambda)$ and $k(\lambda)$ are mixed internally using Maxwell-Garnett effective medium approximation: the particles are composed of black carbon, brown carbon, absorbing insoluble, non-absorbing insoluble embedded in a soluble host.

ii)- several aerosol components with known spectral dependencies $n(\lambda)$ and $k(\lambda)$ are mixed internally proportionally to the volume fraction of each component including black carbon, brown carbon, iron oxide, non-absorbing mineral dust, etc.

* the spectral values $n(\lambda_i)$ and $k(\lambda_i)$ can be assumed and fixed for each aerosol component.

The aerosol modeling options described above are suitable for retrievals when such aerosol parameters as size, shape, spectral refractive index and vertical distribution can be explicitly retrieved. However, in the situations with limited information content retrieval of those detailed characteristics could be challenging or impossible. In such a situation, the aerosol single scattering properties can be modelled as an external mixture of several
aerosol components and the columnar properties of each component:

$$ \tau_{scat/ext} = \sum_{k=1}^{K} C_{v}^{k} \rho_{scat/ext}^{k}(\lambda) \qquad (3.10) $$

and

$$ \tau_{scat} P_{ii'}(\lambda; \Theta; h) = \sum_{k=1}^{K} C_{v}^{k} \rho_{scat/ext}^{k}(\lambda) P_{ii'}^{k}(\lambda; \Theta; h) \qquad (3.11) $$

where $\rho_{scat/ext}^k(\lambda)$ and $P_{ii'}^k(\lambda; \Theta; h)$ denote the scattering/ extinction per unit of volume and phase matrix of each aerosol component that are pre-calculated using complex refractive index, size and shape distributions assumed for each aerosol component. Correspondingly, only K concentrations drive the modeling of columnar properties of aerosol. The number of the retrieved parameters is significantly reduced in this approach, which is helpful for the observations with limited sensitivity to the size, shape and refractive index of the aerosol particles. For example, as will be discussed below this approach was successfully employed in processing MERIS and POLDER satellite data by GRASP. Also, as discussed by Lopatin et al. (2021) this approach is rather suitable for processing vertically resolved observations by lidar or airborne radiosonde data. In applications to MERIS and POLDER satellite data, the vertical aerosol distribution was assumed the same for all the components and modeled as Gaussian or exponential function, while in processing of vertically resolved observations a detailed separate vertical distribution can be retrieved for each aerosol component. In general, the aerosol components are associated with optically distinct types of aerosol based on particle sizes, scattering and absorption capabilities, etc. For example, in GRASP applications to MERIS, POLDER, lidar and radiosonde data, the aerosol components were defined based on AERONET climatology by Dubovik et al. (2002) and Smirnov et al. (2002), the exact definitions are discussed by Lopatin et al. (2021). At the same time, these components can be easily redefined and modified in frame of GRASP software.

Thus, aerosol in GRASP retrieval can be represented as rather sophisticated mixture of one or several aerosol components that can differ by particle size distribution, shape distribution, vertical profile and complex index of refraction. All these characteristics can be either assumed and fixed or retrieved. The complexity of modeling each of these characteristics depends on the information content of observations used in each specific GRASP application.

Gaseous absorption can also be fully accounted in GRASP forward modeling of atmospheric radiances. The possibility of rigorous accounting for atmospheric gases absorption was recently incorporated with the help of the team from Free University of Berlin with an objective to explore the synergy between photometric and spectrometric observations and combined retrieval of properties of both atmospheric

aerosol and gases. The rationale of such an approach is in improving the accuracy of atmosphere components characterization. Indeed, the gases and aerosol affect radiation very differently, and, therefore, very different measurements are used for their retrieval. For example, Figure 6 illustrates a typical spectral dependence of extinction by atmospheric gases, and fine (radius < ~ 0.5 micron) and coarse (radius > ~ 0.5 micron) aerosol particles using complex refractive index of quartz. It is clear from this illustration that absorption of atmospheric gases has much finer spectral features than aerosol, especially in visible spectral range.

Therefore, in practice most commonly the retrievals of aerosol and gases are separated. In general, radiometers are used to make measurements in only a few selected spectral “window” channels with minimal gaseous absorption with the purpose to characterize aerosol properties. In a contrast, the gaseous retrievals tend to rely on the hyper-spectral features of gaseous absorption observed by spectrometers, with aerosol contribution subtracted as smooth “background”. In these regards, the current version of GRASP allows joint retrieval of both aerosol and gases from joint spectrometric and photometric observations.

**Equation 3.12** illustrates the modeling of the transmitted radiances on the properties of atmospheric aerosol, gases and molecular scattering for cloud-free atmosphere.

$$ \begin{aligned} I(\lambda, h) &\sim \int_{\lambda_{min}}^{\lambda_{max}} \Omega(\lambda) I_0(\lambda) \left( \begin{array}{l} -m \left( \int_0^h \left( \sum_{i=1}^{N_{abs}} N_i(h) C_{i,abs}^{gas}(\lambda, T(h), P(h)) + C_{ext}^{gas}(\lambda, h) + C_{ext}^{mol}(\lambda, h) \right) dh \right) + \text{mult scat} \end{array} \right) d\lambda = \\ &= \int_{\lambda_{min}}^{h_{max}} \left( \Omega(\lambda) I_0(\lambda) \left( e^{-m \left( C_{i,abs}^{gas}(\lambda) + C_{ext}^{gas}(\lambda) + C_{ext}^{mol}(\lambda) \right)} + \text{mult scat} \right) \right) d\lambda \end{aligned} \qquad (3.12) $$

where $N_i(h)$ is profile of i-th gas concentration, the cross-section of gaseous absorption $C_{i,abs}^{gas}(\lambda, T(h), P(h))$ depends on temperature ($T$) and pressure ($P$) and $\Omega(\lambda)$ denotes spectral filter function. At the moment, GRASP forward model can account in monochromatic calculations for any gas in any spectral range for a maximum of ten different gaseous species simultaneously using a line-by-line approach. The calculations of the radiative flux calculations realized in GRASP [as described by Derimian et al. (2016)] also account for such details as dependence of water vapor absorption cross section on its concentration due to continuum (e.g., see Kato et al., 1999). The vertical profiles of pressure and temperature are to be provided as an additional and necessary input for line-by-line calculations. In order to speed up the calculations, the look-up-tables with the absorption spectral information for each considered gas species are used. These look-up-tables of absorption values are usually pre-calculated for a wide range of pressures, temperatures and wavelengths and stored into look-up-tables that are later called in the retrieval. As a result, only a renormalization to the required concentration and a simple interpolation to the defined pressure and temperature profile are
<figure><img src="image_9.png" /><figcaption>**FIGURE 6** | The displayed extinction for fine (radius < ~ 0.5 micron) and coarse (radius > ~ 0.5 micron) aerosol particles was calculated using complex refractive index of quartz.</figcaption></figure>

needed to obtain the gas absorption profile. In order to obtain these pre-calculated values the CGASA line-by-line model (developed by the Institute for Space Sciences, Free University Berlin) is used. Notwithstanding, the look-up-table format is readable, easy to understand and it is open to the users to include their own look-up-tables which fulfill their requirements.

However, direct spectral integration of line-by-line fine structure requires consideration of a very large number of channels with full accounting for multiple scattering for each channel. Such direct implementation of the spectral integration has high accuracy but is highly time consuming. Therefore, an alternative methodology called a “k-distribution approach” (see Doppler et al., 2014) is also been integrated into GRASP forward model to speed up the integration time. The accuracy of this methodology is somewhat reduced compared to line-by-line procedure but sufficient for most remote sensing applications. The k-distribution methodology relies on the same basic assumption of smooth aerosol spectral behavior and reduces the number of multiple scattering runs for the spectral integration to about 10 or even less instead of thousand. Specifically, the approach sorts the gaseous absorption coefficients “k” within the spectral function to a limited set of bins (~10) by the magnitude of the coefficients and then implements only a single multiple scattering run for each “k” bin.

It should be noted that the calculation of the k-distribution kernels themselves are generated externally and need to be prepared as input prior using k-distribution modeling in GRASP. The dependence of GRASP on external sources to use this technique can be considered as a drawback. At the same time, the fully standardized format of this input information in GRASP allows the utilization of a vast variety of different principles (Uncorrelated, TOA-correlated, Layer-Correlated, etc.) and algorithms that can be used to calculate k-distribution parameters. The user can select the specific methodology to match the specific requirements in each situation and is not forced to use one specific approach.

Thus, depending on available measurements and information both the profiles of gaseous concentrations as well as temperature and pressure can be determined.

# 4 MULTI-TERM LSM: RETRIEVAL PRACTICE USING GRASP

GRASP algorithm was designed as generalized approach that could be used in diverse remote sensing applications (see introduction by Dubovik et al. (2014)]. Indeed, the Multi-term LSM approach employed for numerical inversion and forward model allowing simulations of transmitted and reflected radiation in the atmosphere allows an efficient exploitation of the GRASP algorithm for diverse applications for atmospheric remote sensing including retrieval of detailed properties of aerosol particles from laboratory and *in situ* observations, as well as from ground-based, airborne and satellite passive and active remote sensing observations. GRASP also retrieves detailed properties of underlying surface reflectance from down-looking satellite and airborne observations. In addition, due to very general principles realized for inversion and forward model calculations, GRASP application scope is expected to be significantly expanded, so that it can be used for applying both radiometric and hyper spectral observations in a wide spectral range from ultra violet (UV) to thermal infrared (TIR) and retrieving not only aerosol and surface properties but also properties of atmospheric gases and clouds. Some applications for retrieving properties of hydrosols can also be envisioned.

## 4.1 Laboratory and *In Situ* Observations of Aerosol Single Scattering

Laboratory and *in situ* instrumentation generally can provide more detailed and elaborated measurements of radiances transmitted and diffused by a sample of air volume or surface area compared to the remote sensing using ground-based, airborne or satellite measurements. For example, using laboratory and *in situ* instruments most of the single scattering characteristics, $τ_{scat/ext}$ and the elements of scattering matrix $P_{ii'}(λ; Θ; h)$ can be directly measured. Using the approach described by Eq 3.1, 3.2 these characteristics can be simulated by
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**FIGURE 7** | The results of a study using GRASP to characterize the size and refractive index of monodisperse polystyrene spheres from multiwavelength polar nephelometer measurements. (A): scattering matrix element measurements (circles) at 473, 532 and 671 nm obtained from a PSL sample along with the corresponding GRASP fits (solid lines). (B): retrieved real part of the refractive index, alongside three prior measurements of PSL refractive index (Ma et al., 2003; Sultanova et al., 2003; Jones et al., 2013). The embedded subplot in the upper-right corner shows the GRASP retrieved size distribution (blue) alongside the manufacturer's specified central radius (red dashes) and ±0.5σ distribution width (red dots).</figcaption></figure>

GRASP based on known microphysical aerosol properties (particle size distribution, shape parameter, etc.). Inversely, if measurements of $τ_{scat/ext}$ and/or $P_{ii'}(λ; Θ; h)$ are available, the microphysical properties of the aerosol can be retrieved by GRASP.

The retrieval of aerosol microphysical properties from the measurements of the phase matrix $P_{ii'}(λ; Θ; h)$ elements at several wavelengths (together with or without $τ_{scat/ext}$) is the base for more challenging remote sensing satellite and ground-based observations. Indeed, such angular remote sensing observations from satellite and ground-based instruments are less detailed and affected by complex multiple scattering effects.

The angular measurements of $P_{ii'}(λ; Θ; h)$ in a wide range of scattering angles (~from 0° to 180°) at several wavelengths contain practically maximal information about aerosol microphysical properties, especially if these angular measurements are complimented with spectral measurements of total scattering or extinction $σ_{scat/ext}$ of the aerosol sample. Therefore, the retrieval of practically all aerosol characteristics can be expected from such detailed measurements. For example, in studies by Dubovik et al. (2006) such properties as detailed size and shape (axis or aspects ratios) distributions were retrieved together with spectral complex index of refraction from laboratory measurements of aerosol samples phase matrices. In the framework of GRASP software, the inversion of full scattering matrix (6 elements) at any wavelengths and for any angular range selection together with the spectral measurements of $σ_{scat/ext}(λ, h)$ is realized. In addition, the measurements of any simple functions of $P_{ii'}(λ; Θ; h)$ and $σ_{scat/ext}(λ, h)$ can be easily added if required. For example, the measurements of total absorption $σ_{abs}(λ, h)$, backscattering $β(λ, 180°, h)$ or total scattering $σ_{scat}(λ, Θ_{min}, Θ_{max}, h)$ in the selected range of scattering angles: $Θ_{min} ≤ Θ ≤ Θ_{max}$ can be included in the set of inverted single scattering properties with

GRASP. From the practical viewpoint the inversion of the phase function $P_{11}(λ; Θ; h)$ (or angular scattering $σ_{scat}(λ, h)P_{11}(λ; Θ; h)$) and linear polarization $P_{12}(λ; Θ; h)$ at several wavelengths (together with or without $τ_{scat/ext}$) is one of the most popular applications. Indeed, most of passive remote sensing techniques using angular information based on the measurement of atmospheric radiation characterized by $I(λ; Θ; h)$ element of Stocks vector (e.g., see Mishchenko et al., 2006) and linear polarization characterized by $Q(λ; Θ; h)$ and $U(λ; Θ; h)$ elements of Stocks vector. In single scattering approximation, $I(λ; Θ; h)$, $Q(λ; Θ; h)$ and $U(λ; Θ; h)$ measured at the top or bottom of the atmosphere depends only on $P_{11}(λ; Θ; h)$ and $P_{12}(λ; Θ; h)$. Some other elements of phase matrix contribute to $I(λ; Θ; h)$, $Q(λ; Θ; h)$ and $U(λ; Θ; h)$ only via rather complex multiple scattering interactions of solar radiation with the atmosphere. Therefore, the single scattering contributions of $P_{11}(λ; Θ; h)$ and $P_{12}(λ; Θ; h)$ often dominate passive radiometric and polarimetric remote sensing measurements. Correspondingly, the retrieval of aerosol microphysical properties from the measured phase function $P_{11}(λ; Θ; h)$ together with (or without) linear polarization $P_{12}(λ; Θ; h)$ at several wavelengths can be considered as a base for a more challenging remote sensing satellite and ground-based observations.

Therefore, the development of aerosol property retrievals from angular measurements of the phase function and linear polarization helps to obtain understanding of a retrieval's potential and to design relevant remote sensing approaches. In these regards, in the framework of GRASP the potential for a wide diversity of retrievals can be evaluated on the basis of single scattering sensitivity tests. For example, there are a number of cases where the phase function and linear polarization was measured directly and used for the retrieval of detailed aerosol properties with the help of the GRASP approach, including studies by (Dolgos and Martins, 2014; Espinosa et al., 2017;
<figure><img src="image_12.png" /></figure>

Espinosa et al., 2019 W. R.; Puthukkudy et al., 2017; Schuster et al., 2019, etc.). A demonstration of one such retrieval study using measurements of polystyrene (PSL) spheres from Espinosa et al. (2017) is shown in **Figure 7**.

In addition, such remote sensing observations such as aerosol optical depth or lidar backscattering measured from the ground are sensitive mainly to single scattering aerosol properties and generally modeled using the single scattering approximation. The application of GRASP for these observations will be considered below.

## 4.2 Ground-Based Passive Observations

GRASP approach is developed based on the heritage of AERONET retrieval developments by Dubovik and King (2000), Dubovik et al. (2000) and Dubovik et al. (2006). Therefore, GRASP can be used for inversion of AERONET-like radiometric ground-based observation data in a similar way as it is done in AERONET operational aerosol retrieval. Therefore, applying GRASP to such data should provide essentially similar retrieval results. At the same time, even though GRASP has adapted the main conceptual elements

from the AERONET retrieval, the actual algorithm is different, therefore some minor differences can be observed even in the identical application.

The potential of applying GRASP to radiometric ground-based observations lies in the possibilities of trying and evaluating new retrieval concepts. Indeed, as described above in *Theory of Multi-Term LSM Inversion*, GRASP allows one to use a variety of the retrieval setups. For example, in difference with AERONET operational retrieval, it is possible to attempt deriving more detailed shape distribution, to employ multi-component aerosol mixture with different refractive indices for each component, to approximate size distribution by log-normal functions, etc. The examples of such diverse AERONET inversions can be found in a number of published studies. Li et al. (2019) demonstrated the retrieval of aerosol properties from AERONET observations assuming aerosol as a mixture of different aerosol components. As it was demonstrated this approach allows adequate modeling of aerosol AERONET observation and the retrieved results provide some inside on aerosol composition. Torres et al. (2014) made extensive analysis of aerosol retrieval sensitivity to geometrical configuration of
<figure><img src="image_13.png" /><figcaption>**FIGURE 9** | The illustration of the inversion dynamics evolution through at the different iterations. (A): fitting of direct Sun hyperspectral observations of Pandora. (B): evolution of the retrieved AOD values.</figcaption></figure>

ground-based sun/sky-radiometer observations. This study was conducted using GRASP algorithm, and employed different assumption of aerosol vertical variability in the series of synthetic tests. In more recent studies the GRASP based methodology was established for retrieving aerosol properties from only direct Sun AERONET observation of optical thickness. Torres et al. (2017) describe the “GRASP-AOD” approach in details and Torres and Fuertes (2020) demonstrate extensive application of the method to AERONET AOD observations. GRASP-AOD uses the assumptions of bi-modal lognormal size distribution and provides AOD of fine and coarse aerosol modes (AODF and AODC). The results of comparisons by Torres and Fuertes (2020) have shown that GRASP-AOD provides AODF and AODC comparably good and sometimes even better accuracy than the established retrieval approaches by O'Neill et al. (2003) and Perez-Ramirez et al. (2015). In addition, GRASP-AOD provides information about aerosol size distribution including concentrations, standard deviations, median and effective radii of total, fine and coarse modes of aerosol size distributions. These retrieved aerosol properties agree well with the results of full AERONET retrieval performed for the same AOD observations.

GRASP also can be used for inverting diverse non-AERONET observations from ground or AERONET observations that generally included in AERONET processing. For example, Torres and Fuertes (2020) demonstrated aerosol retrieval using direct Sun AOD observation combined with the sky-radiances measured in solar aureole. Such retrieval showed to provide superior about aerosol size distribution compare to GRASP-AOD only retrieval and can be applied to the AERONET data acquired in partially cloudy environment. Roman et al. (2017, 2021) applied GRASP for aerosol retrieval from the measurements by lunar aureole with a sky camera. Torres et al. (2017) and Popovici et al. (2018) have applied GRASP-AOD approach for interpretation of AOD observations by lunar-photometer. The upper panel of **Figure 8** illustrates the size distribution by GRASP-AOD from AOD observations by lunar-photometers. The lower panel of **Figure 9** illustrates

the improvements in the aerosol retrieval by inverting radiances measured in solar aureole together with AOD data compare to inversion of AOD data only (Torres and Fuertes, 2020).

It is important to underline that GRASP allow rather straightforward use of ground-based radiometric observations in the retrievals using synergy with other ground-based, satellite or air-borne observations. For example, the combined processing of ground-based radiometric and lidar observations as will be discussed in next sections.

Finally, as described in *Theory of Multi-Term LSM Inversion*, several new modules were integrated into GRASP forward model for making it possible to simulate both photometric and hyperspectral observations in visible and thermal IR part of the spectrum. Therefore, at present GRASP is ready for a processing not only photometric but also spectrometric ground-based observations in wide spectral range. Specifically, very recently first tests have been done with application of GRASP to combined inversion AERONET type radiometric observation together with CLIMAT thermal TIR radiometer (Legrand et al., 2000; Brogniez et al., 2003) or PANDORA spectrometer (Herman et al., 2015). These developments are yet at early stages, nonetheless the first test indicate high potential for these new applications. For example, AERONET/CLIMAT retrieval will allow the retrieval of a consistent aerosol optical characteristic in wide visible-TR spectral range and synergy inversion of joint AERONET and PANDORA observations is promising for improved simultaneous retrieval of aerosol and such atmospheric gases as NO2 and Ozone (Herreras-Giralda et al. in preparation). These applications were overall realized and tested with synthetic data, where NO2 are retrieved simultaneously with aerosol. **Figure 9** illustrates the synthetic retrieval applied to a combined data set including almost 150 Pandora channels ranging from 400 to 440 nm, and the four standard AERONET channels (440, 675, 870 and 1,020 nm). In this application GRASP algorithm successfully operates with the hyperspectral features of gas absorption and the smooth aerosol optical characteristics. The real data inversion results of this combined retrieval of aerosol and NO2 have been done using DIVA (see *Synergetic Retrievals and Other Diverse GRASP*
<figure><img src="image_14.png" /><figcaption>**FIGURE 10** | An example of GaRRLiC/GRASP retrieval columnar and vertical aerosol properties from a combination of ground-based radiometer and multi-wavelength lidar, LR denotes lidar ratio, RRI and IRI denote real and imaginary refractive indices respectively.</figcaption></figure>

Applications) platform and initial validation has been done against independent algorithms as AERONET, for aerosol, and PGN (Pandonia Global Network) (Herman et al., 2009; Tzortziou et al., 2012) for gas related products.

Thus, GRASP is being expanded for retrieval of also atmospheric gases in addition to aerosol and surface. Moreover, some very recent developments are initiated to include cloud retrieval in GRASP framework. Theretofore, to reflect this development, it is appropriate to rename GRASP as Generalized Retrieval of Atmosphere and Surface Properties instead of Generalized Retrieval of Aerosol and Surface Properties.

## 4.3 Ground-Based Active Observation and Synergy With Passive Observations

Initially the possibilities of processing active lidar observation were integrated into GRASP by Lopatin et al. (2013) as part synergetic GARRLiC (Generalised Aerosol Retrieval form Radiometer and Lidar Combination) retrieval that inverts simultaneously the photometric and lidar observations. Since then GARRLiC/GRASP has been employed and discussed in numerous studies (Benavent-Oltra et al., 2017; Benavent-Oltra et al., 2019; Benavent-Oltra et al., 2021; Hu et al., 2019; Tsekeri et al., 2013; Tsekeri et al., 2017; Roman
et al., 2018; Titos et al., 2019; Molero et al., 2020; Parajuli et al., 2020; Konsta et al., 2021 etc.). Moreover, it is currently being employed for operational processing of such combined data in the frame of the European ACTRIS infrastructure (https://www.actris.eu/). The standard GaRRLiC approach by Lopatin et al. (2013) inverts spectral lidar data together with sun/sky scanning AERONET like radiometer, and the detailed columnar size distributions, and spectral complex index of refraction can be retrieved for two fine and coarse aerosol components together with two vertical profiles of their concentration can be derived as illustrated by Figure 10.

While the original GARRLiC scheme was limited to the specific observational set of multi-wavelength elastic scattering lidar together with AERONET-like sun/sky-radiometer observations, the area of GRASP application was significantly extended during last years and now covers a variety of active and vertically resolved observations that can be processed. For example, at present GRASP can be applied for interpretation of such vertical observations of atmosphere as profiles of extinction, backscatter, normalized elastic and inelastic lidar signals, profiles of volume and particle depolarization. In addition, a possibility of lidar only retrieval was introduced in GRASP together with several technical retrieval details improving processing of vertically resolved observations (Lopatin et al., 2021).

Thus, at present GRASP allows simulation and inversion of point or vertical observations of the profiles of backscattering, extinction, 6 elements of scattering matrix. In addition, several characteristics that are simple functions of the above scattering characteristics were introduced as a potential input data for GRASP. Specifically, both the elastic and inelastic lidar measurements can be inverted by GRASP. The lidar equation for ground-based observations is realized as the following:

$$L(\lambda, h) = A(\lambda)\beta(\lambda, h) \exp\left(-2 \int_{Z_{\min}}^{Z(h)} \sigma(\lambda, z) dz\right), \quad (4.3.1)$$

where $\sigma(\lambda, h)$ denote extinction and $\beta(\lambda, h)$ backscattering profiles. The extinction profile includes aerosol, molecular scattering and gaseous absorption components $\sigma(\lambda, h) = \sigma_a + \sigma_m + \sigma_g$ and backscattering includes aerosol and molecular components $\beta(\lambda, h) = \beta_a + \beta_m$, $A(\lambda)$ is the instrumental constant, $z(h)$ is lidar path related with the atmospheric altitudes of target $h$, ground level $h_{BOA}$ and zenith angle of lidar inclination $\Theta_L$ as $z(h) = \frac{h-h_{BOA}}{\cos(\Theta_L)}$. Here both backscattering and extinction profiles depend of aerosol properties. The ambiguity in separation of the backscattering by the layer and extinction of the lidar signal by underlying layers is considered as the key challenge in the interpretation of elastic lidar signals. The lidar systems using inelastic scattering address that ambiguity by measuring the following signal:

$$L_{inel}(\lambda, h) = A(\lambda, \lambda')\beta_{inel}(\lambda, h) \exp\left(-\int_{Z_{\min}}^{Z(h)} (\sigma(\lambda', z) + \sigma(\lambda, z)) dz\right) \quad (4.3.2)$$

where scattering, $\lambda'$ - wavelength of exciting impulse that triggers inelastic backscatter at wavelength $\lambda$, $\beta_{inel}(\lambda, h)$ inelastic backscattering of atmosphere, $A(\lambda, \lambda')$ is known calibration constant. The shift $\lambda' \to \lambda$ in inelastic backscattering $\beta_{inel}(\lambda, h)$ could be a result of gaseous molecules emission frequency shifts due rotations and vibrations of molecules and can be rather accurately estimated based on known characteristics of emitted lidar impulse and atmospheric gases. Therefore, neglecting the differences in aerosol extinction at $\lambda'$ and $\lambda, \sigma_a(\lambda, h)$ can be considered as the only fully unknown characteristic in Eq. 4.3.2 and can be obtained from $L_{inel}(\lambda, h)$ by rather straightforward transformations. The measurements of such advanced lidar systems as HSRL are usually converted to the backscatter $\beta... (\lambda, h)$ and extinction profiles $\sigma... (\lambda, h)$ (Hair et al., 2008; Rogers et al., 2009, etc.) that can be used as an input to GRASP algorithm for conducting a full aerosol retrieval.

Another characteristic measured by advanced lidars with polarimetric capabilities is the profile of volume and aerosol particle depolarization. The profile of particle depolarization can be estimated from the ratio of lidar returns $L_\perp$ and $L_\parallel$ obtained with from emitted polarized light beams as following (Freudenthaler et al., 2009):

$$\delta = \frac{L_{\perp}(\lambda, h)}{L_{\parallel}(\lambda, h)} = \frac{\beta_{\perp}(\lambda, h)}{\beta_{\parallel}(\lambda, h)} = \frac{P_{11}(\lambda, 180^{\circ}, h) - P_{22}(\lambda, 180^{\circ}, h)}{P_{11}(\lambda, 180^{\circ}, h) + P_{22}(\lambda, 180^{\circ}, h)}, \quad (4.3.3)$$

where the subscripts “$\perp$” and “$\parallel$” indicate cross- and co-polarized components correspondingly. Both profiles of aerosol and atmospheric volume depolarization ratio can be used for GRASP.

Thus, diverse vertically resolved observations discussed above from a single instrument or in a combination with other observations (e.g., a passive instrument) can be inverted by GRASP. The number and type of aerosol parameters retrieved depend on a chosen configuration of employed aerosol forward model (see Theory of Multi-Term LSM Inversion) and a priori constraints that should be chosen in accordance with the information content of inverted data. For example, if only observations of a single wavelength lidar are inverted, only profile of one aerosol component can be retrieved using a priori assumptions about most of other aerosol characteristics. In case of inversion of observation from multi-wavelength advanced lidar system, a number of parameters and scope of retrieved aerosol information can be significantly extended. For example, Lopatin et al. (2021) demonstrated simultaneous retrieval of four profiles of different aerosol components from such data (aerosol was modeled using Eq 3.10, 3.11). Lopatin et al. (2021) provided more detailed discussion of GRASP applications to inversion of vertically resolved observations and numerous illustrations.

## 4.4 Satellite and Airborne Passive Observations

One of the original motivations for initializing GRASP developments was an idea of creating enhanced retrieval of
aerosol from satellite observations by multi-angle polarimeter such as POLDER. In fact, numerous theoretical and practical studies have concluded that polarimetry is an approach that can provide an accurate characterization of aerosols with the detail and precision sufficient for many important applications (Mishchenko et al., 1997; Mishchenko and Travis, 1997; Hasekamp and Landgraf, 2005; Hasekamp and Landgraf, 2007; Kokhanovsky et al., 2010; Knobelspiesse et al., 2012). On the other hand, as discussed by Dubovik et al. (2019), the interpretation of multi-angular multi-spectral polarimetric (MAP) data is quite challenging from the fundamental point of view. Indeed, polarimetry is highly sensitive to a large number of atmospheric parameters, and accounting adequately for all these sensitivities in the retrieval algorithm is very demanding, especially in the satellite applications where large volumes of data need to be processed with a minimal delay, and for quite a long time there were no available operational MAP product demonstrating the practical advantages of MAP retrieval. In these regards, GRASP was adapting AERONET retrieval methodology and, in difference with conventional look-up-table (LUT approaches) attempted highly advanced statistically optimized fitting of satellite data implementing rigorous search for the solution in a continuous space of solutions. Such approaches are considered as the most promising for satellite retrieval and with some technical differences are deployed in other algorithms of new generations developed for interpretation of MAP observations (Hasekamp and Landgraf, 2007; Hasekamp et al., 2019; Fu et al., 2020; Dubovik et al., 2011; Xu et al., 2016; Xu et al., 2017, etc.). The challenges of realizing such state-of-the-art algorithms lie in the fundamental difficulties to adapt them optimally to the high sensitivities of MAP instruments and in specific technical issues, such as, a need of significantly more computational time than LUT algorithms to simulate the satellite signal and the Jacobian derivatives matrices online. In these regards, the developed structure of GRASP algorithm allows realizing rigorous inversion of MAP satellite observations. In addition, the GRASP software was extensively optimized to reduce time of calculation and adapted to practical application to real satellite observations.

At present, 18 months of POLDER-1 and -2 and 9 years of POLDER-3 observations have been processed and several versions of the retrieval product have been archived at the AERIS/ICARE Data and Services Center (http://www.icare.univ-lille1.fr) and GRASP-OPEN site (https://www.grasp-open.com). For POLDER, GRASP utilizes radiance and polarization observations from all available spectral channels with minor gaseous absorption: 5 and 6 channels for POLDER-1 and -2, and POLDER-3, correspondingly and 3 polarized radiances spectral channels for all instruments. The retrieval uses a unique global set of constraints (no location-specific assumptions) and a single initial guess globally. The radiative transfer computations accounting for multiple interactions of the scattered solar light in the atmosphere are performed on-line. All GRASP retrievals were performed at POLDER native resolution POLDER-1 and -2

at ~7 km and POLDER-3 at ~6 km. In the saved data archives, the original POLDER/GRASP retrievals are stored at Level-1, Level-2 and -3 products and are publicly available in the form of daily, monthly, seasonal, yearly and climatological datasets. The Level-2 data contain full resolution data filtered following the established quality criteria. Level 3 data is aggregated into a 0.1° and 1° grid box using the sinusoidal projection from Level-2 data.

GRASP allows a variety of different possibilities on modeling aerosol scattering, surface reflectance and implementing atmospheric RT calculations and due to various reasons several different configurations of the atmospheric forward model were used to generate aerosol products. For example, the full POLDER-3/PARASOL data archive was processed by GRASP using the three following retrieval configurations: POLDER-3/GRASP «optimized», «high-precision» and «models». The observations of POLDER-1 and -2, at present, were processed using only a single GRASP/Models approach. The «optimized» and «high-precision» differ only by the precision of the RT calculations, while they use the same aerosol model driven by aerosol size distribution, spectral values of complex index of refraction, fraction of spherical particles and the aerosol scale height (Dubovik et al., 2011). The «models» approach uses the assumption of an external mixture of several aerosol components (see Eq 3.10-3.11) and directly retrieved parameters include aerosol concentrations and a scale height. In addition, recently POLDER-3/GRASP aerosol product has been generated using “component” approach (Li et al., 2019; Li et al., 2020a; Li et al., 2020b). This approach retrieves the size resolved fractions of aerosol components representing the different composition species, such as black carbon, brown carbon, fine/coarse mode non-absorbing soluble and insoluble, coarse mode absorbing and aerosol water. The retrieved fractions drive the aerosol spectral index of refraction in modeling of atmospheric radiances. **Figure 11** illustrates the climatology of aerosol component columnar mass concentration derived from POLDER-3 over East Asia region by the GRASP/Component algorithm.

Independently of which aerosol model approach all retrieval data products contain the aerosol main aerosol characteristics including spectral aerosol optical depth (AOD), aerosol absorption optical depth (AAOD) and single scattering albedo (SSA) as well as Ångström exponent (AE), spectral fine mode AOD (AODF) and coarse mode AOD (AODC). **Figure 12** illustrate the climatology of these parameters.

The main aerosol POLDER-3/GRASP products including AOD, AE, AODF, AODC, SSA and AAOD were extensively evaluated using validations against AERONET and comparisons with the original POLDER algorithm (PARASOL/Operational), and MODIS Collection 6 aerosol products. For example, Schutgens et al. (2021) have evaluated GRASP/Models Level3 1-degree SSA against AERONET and compared with other satellite SSA products. The studies recognized GRASP/ Models as one of reliable and most extensive data SSA set. Chen et al. (2020) have conducted the detailed validation of Level 3 0.1 degree products. The studies have shown that the POLDER-3/ GRASP retrieval provided reliable aerosol products. Specifically, POLDER-3/GRASP spectral products including AOD for six
<figure><img src="image_15.png" /></figure>

**FIGURE 11** | Climatology of aerosol component columnar mass concentration derived from POLDER-3 over East-Asia by the GRASP/Component approach: (A) fine mode black carbon, (B) fine mode brown carbon, (C) coarse mode mineral dust.

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

**FIGURE 12** | Climatological values of AOD (565), Angstrom exponent, SSA(670) and scale height for winter 2009 provided by POLDER-3/GRASP optimized product.
<figure><img src="image_17.png" /><figcaption>**FIGURE 13** | The illustrations of the POLDER/GRASP product comparisons with AERONET data. Upper panel: the correlations of POLDER-3/Models AOD with AERONET for several selected wavelengths (440, 550, 870 nm). Middle panel: the correlations of POLDER-1,2/GRASP AOD with AERONET for several selected wavelengths (440, 550, 870 nm). Lower panel: the correlations of the detailed aerosol parameters retrieved by POLDER-3/GRASP/HP product, from left to right: AE, AODF(550), AODC(550), SSA(865).</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>**FIGURE 14** | The correlations of the detailed aerosol parameters retrieved by POLDER-3/GRASP/Component product, from left to right: AOD(550), AODF(550), AODC(550), AE.</figcaption></figure>

wavelengths in the range 443–1,020 nm agree well with the AERONET AOD measurements, e.g., for POLDER-3/Models AOD correlation coefficients R are ≥0.86 over land and ≥0.94 over ocean with BIAS not exceeding 0.01 over land and 0.02 over ocean for all wavelengths. The upper panel of **Figure 13** demonstrates the correlations of satellite AOD with AERONET for several selected wavelengths.

The data from POLDER-1 and -2/GRASP retrievals were not included in the analysis by Chen et al. (2020), while the limited comparisons are shown in middle panel of **Figure 13** demonstrate that the quality of the POLDER-1 and -2/GRASP retrievals are expected to be close to those of POLDER-3/GRASP retrievals.

The comparisons with MODIS aerosol products showed that the POLDER-3/GRASP AOD retrievals are very coherent with popular MODIS data while also exhibit some important advancements. For example, POLDER-3/GRASP retrievals provide more reliable detailed aerosol parameters such as AE, AODF and AODC especially over land and such parameters as SSA and AAOD that are generally not available from MODIS-like instruments. The validation of POLDER-3/GRASP products by Chen et al. (2020) showed a robust correlation of the retrieved SSA and AAOD spectral values with AERONET (440–1,020 nm), correlations increase for the retrievals corresponding to the events with higher AOD. For AAOD retrievals overall the BIAS did not exceed 0.01, suggesting that POLDER-3/GRASP products can be used for making global estimations of AAOD at such level of uncertainty. The lower panel of **Figure 13** demonstrates the correlations of the detailed POLDER-3/GRASP/HP products. The detailed discussion can be found in the study by Chen et al. (2020).

One key finding of Chen et al. (2020) analysis is that the best retrieval of total AOD is provided by the simplest approach (GRASP/ Models) in which the retrieval is restrained to a superposition of predefined aerosol components, significantly reducing the number of free parameters for retrieval. The AOD retrieved from more complex GRASP/HP and GRASP/Optimized approaches over land has

notable bias (~0.06–0.07 at 500). In these regards, the GRASP/ Component [that was not considered by Chen et al. (2020)] provided apparently overall the most coherent total and detailed aerosol properties. Indeed, the validation by Zhang et al. (2021) of GRASP/Component POLDER-3 products against AERONET showed that total AOD from GRASP/Component is close in accuracy to GRASP/Models and higher in accuracy than AOD form GRASP/HP and GRASP/Optimized. Also, the accuracy of most of detailed aerosol products (AOF, AOC, AE) from GRASP/ Component approach is overall higher or close to that of best results of GRASP/HP and GRASP/Optimized. For example, **Figure 17** illustrates results for AOD, AODF, AODC and AE over land that can be compared in **Figures 13, 14**.

Evidently, that future efforts on improving the GRASP retrieval will be aimed at achieving accurate retrievals within one approach, however the situation also reveals the challenge of developing a unique approach that can provide a retrieval of all parameters with highest accuracy from MAP observations. Indeed, multi-angular polarimetric observations have sensitivity to different aerosol properties, and therefore the MAP algorithms tend to be designed for the retrieval of large number of parameters, while in the situations with low aerosol presence the information from observations may not be sufficient to retrieve all parameters reliably. Moreover, POLDER like MAP observations have high potential for essentially helping to improve extensive monitoring air quality parameter that are vital for many evaluating the dynamic of environment. For example, Wei et al. (2020) demonstrated higher capacity of POLDER product than single-view MODIS data for characterization of PM2.5 from space and Wei et al. (2021) presented a methodology of POLDER/GRASP products for deriving PM10 the characteristics that is generally even more difficult to obtain from remote sensing than PM2.5. Moreover, MAP and specifically POLDER/GRASP data provide were used as additional constraints for improving global emissions of the atmospheric components in chemical transport models (e.g., Chen et al., 2018; Chen et al., 2019; Elguindi et al., 2020). This
<figure><img src="image_19.png" /><figcaption>**FIGURE 15** | Climatological values (9 years averages) of surface reflectance properties retrieved from POLDER-3. (A): the land BRDF parameters. (B): the ocean surface reflectance parameters.</figcaption></figure>

supports the highly promising concept of synergizing satellite with available modeled information for advancing satellite remote sensing.

In all above processings, the underlying surface reflectance was retrieved simultaneously with the aerosol properties and the detailed BRDF (Bi-directional Reflectance Distribution Function) and BPDF (Bi-directional Polarization Distribution Function) of ocean and land surfaces were derived. Specifically, the parameters of “Ross-Li model” BRDF and Maignan et al. (2009) BPDF models were retrieved over land. **Figure 15** demonstrates the climatological surface reflectance properties retrieved from POLDER-3.

The preliminary validations of POLDER surface reflectance show a generally good agreement with other surface products as those from MODIS, especially for such general parameters as surface albedo. At the same time, for detailed BRDF properties some differences are obvious for the Ross-Li model parameter related with angular anisotropy of surface BRDF. For example, **Figure 16** illustrates the correlations of three Ross-Li model BRDF parameters provided by POLDER-3/GRASP and retrieved by MODIS. The correlations for 2-nd and 3-rd parameters are significantly lower than for the first one. This can be explained by high sensitivity of both POLDER and MODIS observations to the total reflectance of the land surface, while the sensitivity to the BRDF anisotropy, driven by 2nd and 3rd BRDF parameters, is evidently higher for the multi-viewing POLDER than for single-view MODIS measurements.

The reflective properties of ocean surface are modeled analogously to earlier developments (Deuzé et al., 2001;

Herman et al., 2005; Tanré et al., 2011). The Fresnel's reflection on the agitated sea surface is considered by using the Cox and Munk model (Cox and Munk, 1954). The water leaving radiance term and the white cap reflection are taken into account for by Lambertian unpolarized reflectances (Voss et al., 2007). The whitecaps reflectance are dependent on the wind speed at sea surface (Koepke 1984). The seawater reflectance at short wavelengths is not negligible and depends on the properties of the oceanic waters. Thus, in POLDER-3/GRASP, the magnitude of seawater reflectance at each wavelength and the wind speed could be retrieved together with aerosol. The lower panel of **Figure 15** illustrates climatology of ocean surface parameters retrieved from POLDER-3. More illustrations of the chlorophyll, water leaving radiances and wind speed retrievals and their comparison with MODIS results and ECMWF reanalysis can be found in the paper by Frouin et al. (2019). The detailed analysis of both land and ocean surface reflectance properties provided by POLDER/GRASP are ongoing and expected to be discussed in the separate paper.

It is also important to emphasize that POLDER-3/GRASP retrievals are based on rigorous optimized inversion that searches for statistically optimized fitting in a continuous space of solution without using widely used Look-up-Tables, unlike as most of the conventional satellite retrievals. As a result, POLDER/GRASP aerosol and surface reflectance products are globally-consistent based on exactly the same aerosol modeling approach over land and ocean, unique set of a priori constraints and initial guess, while retrieving surface reflectance properties simultaneously with aerosol. As
<figure><img src="image_20.png" /><figcaption>FIGURE 16 | The correlations of the land BRDF monthly average parameters retrieved for September 2008 from POLDER-3 and MODIS.</figcaption></figure>

discussed by more in Dubovik et al. (2019) the similar type of approaches is expected to become common and evolve further in the coming era of multiple MAP instruments. In these regards, there are several GRASP based developments for the future polarimetric mission. For example, an operational aerosol retrieval algorithm for future 3MI/EPS-SG polarimeter (Fougnie et al., 2019) has been developed using GRASP (Marbach et al., 2020). There are also some efforts to adapt GRASP for processing the observations by the future DPC, Aerosol-UA, HARP2 missions. For example, **Figure 17** illustrates the application of the GRASP algorithm for processing AirHARP (airborne Hyper Angular Rainbow Polarimeter) observations made during the NASA ACEPOL 2017 campaign (Puthukkudy et al., 2020a). GRASP retrieved AOD using the AirHARP data is validated using the collocated AERONET and HSRL2 measurements. **Figure 17C** shows a good agreement of GRASP AirHARP retrieved AOD with the collocated AERONET AOD observation with a maximum bias of -0.009 and mean absolute error of 0.015 at 870 nm band. Comparison of GRASP retrieved AOD with measurement of HSRL2 lidar for a forest fire smoke plume shows good correlation and agreement with a correlation coefficient of 0.940 at 532 nm (see **Figure 17D**). Furthermore, preliminary aerosol retrievals based on the HARP CubeSat data using the GRASP was presented at the AGU Fall meeting 2020 (Puthukkudy et al., 2020b). There are also successful GRASP developments for SGLI/GICOM-C and MISR observations (Fuertes et al., 2020).

GRASP approach has also been successfully applied to several single-view satellite radiometers. For example, GRASP has been used to process 10 years (2002–2012) of MERIS (MEDium Resolution Imaging Spectrometer). MERIS is an imaging spectrometer operated from Envisat satellite platform that scanned the Earth’s surface by the so-called push-broom method. Linear CCD arrays provided spatial sampling in the across-track direction, with the satellite’s

motion provided scanning in the along-track direction (Rast et al., 1999). MERIS/GRASP retrieval provided aerosol and surface reflectance retrieval product at 8 used channels 413, 443, 490, 510, 560, 665, 760 and 870 nm. The data were inverted at resolution of 10 km in cloud-free conditions as determined by original MERIS cloud-mask algorithm at latitudes below 60 degrees north and south. The retrieval of aerosol and surface properties is conducted simultaneously. An external mixture (see Eq. 4.3.3) of four aerosol components was used for modeling aerosol. The surface BRDF as well as general retrieval approach was similar to the POLDER/GRASP processing e.g., multi-pixel concept was applied. The full list of the retrieved parameters, as well as, Level 2 and Level 3 products of all derived parameters are available at GRASP-OPEN website (https://www.grasp-open.com).

The MERIS/GRASP products were validated against AERONET and compared to other satellite products. The analysis has shown that MERIS/GRASP retrievals provided robust and reasonably accurate results. For example, as illustrated in **Figure 18** MERIS/GRASP spatial distribution of aerosol hot spots agrees qualitatively well with POLDER/ GRASP products. **Figure 19** shows the global correlations of MERIS/GRASP AOD(560) with AERONET data. The correlations coefficients are ~0.76 over land and ~0.84 (over ocean) with RMSE of 0.155 over land and ~0.06 over ocean. There is a rather small systematic bias of ~0.03 over ocean. However, over land the bias significantly higher: ~0.09. The origin of such high bias over land is being investigated and expected to be addressed in ongoing efforts on updating MERIS/GRASP product. The land surface products are generally in good agreement with the products from other satellite instruments, such as MODIS and POLDER. The lower part of **Figure 19** illustrates the global comparisons of
<figure><img src="image_21.png" /><figcaption>**FIGURE 17** | The application of GRASP algorithm for processing AirHARP observations during the ACEPOL 2017 campaign (Puthukkudy et al., 2020a): (A) RGB image of a forest fire smoke scene captured on 27th October 2017 18:16 UTC; (B) AOD map retrieved for different spectral bands data for the pixels shown in the red rectangle; (C) Comparison of AirHARP-GRASP retrieved AOD with collocated AERONET; (D) Scatterplot of the AOD measured by HSRL2 with GRASP-AirHARP retrieved AOD for a smoke plume on 7th November 2017 19:31 UTC. [Adapted version of the plots originally published in Puthukkudy et al. (2020b), CC BY 4.0].</figcaption></figure>

retrieved DHR (Direct Hemispheric Reflectance) with values provided by MODIS.

The GRASP is been also applied to other single or dual-view radiometers. For example, at present there are several ingoing GRASP based developments of algorithms and retrieval products for Sentinel-3, Sentinel-4 and Sentinel-S5P, VENUS, PRISMA, HIMAWARI, SGLI-2, etc.

It should be noted that while the estimations of the retrieval errors were not generated as a part of POLDER/GRASP and MERIS/GRASP products, but were included in all new and recent developments. For example, Figure 20 illustrates the result of the numerical tests, were the AOD was retrieved from synthetic observation of S5P/GRASP AOD simulated over Banizoumbou AERONET site (the random noise of ~3% was added to the simulations). It can be seen that the error estimates of satellite data obtained based on Eq 2.3.1-Eq 2.3.7 and denoted by the error bars agree well with assumed data. The detailed discussion

on performance of GRASP error estimation concept is provided by Herrera et al. (2020).

## 4.5 Synergetic Retrievals and Other Diverse GRASP Applications

As described above the GRASP is a versatile algorithm that have been already used in a wide variety of applications including diverse passive and active ground-based and satellite observations as well as *in situ* and laboratory observations. Figure 21 shows the overview diagram of GRASP applications. The several key and mature applications have been described above in this section. At the same time, there are some other different but relevant applications of GRASP that have been realized by different studies. For example, GRASP has been used for processing observations by moon-light sky-camera by Roman et al.