Journal cover Journal topic
Atmospheric Measurement Techniques An interactive open-access journal of the European Geosciences Union
Journal topic
Atmos. Meas. Tech., 11, 6589-6603, 2018
https://doi.org/10.5194/amt-11-6589-2018
Atmos. Meas. Tech., 11, 6589-6603, 2018
https://doi.org/10.5194/amt-11-6589-2018

Research article 14 Dec 2018

Research article | 14 Dec 2018

# Joint retrieval of surface reflectance and aerosol properties with continuous variation of the state variables in the solution space – Part 1: theoretical concept

Combined Inversion of Surface and AeRosol
Yves Govaerts and Marta Luffarelli Yves Govaerts and Marta Luffarelli
• Rayference, 1030 Brussels, Belgium
Abstract

This paper presents a new algorithm for the joint retrieval of surface reflectance and aerosol properties with continuous variations of the state variables in the solution space. This algorithm, named CISAR (Combined Inversion of Surface and AeRosol), relies on a simple atmospheric vertical structure composed of two layers and an underlying surface. Surface anisotropic reflectance effects are taken into account and radiatively coupled with atmospheric scattering. For this purpose, a fast radiative transfer model has been explicitly developed, which includes acceleration techniques to solve the radiative transfer equation and to calculate the Jacobians. The inversion is performed within an optimal estimation framework including prior information on the state variable magnitude and regularisation constraints on their spectral and temporal variability. In each processed wavelength, the algorithm retrieves the parameters of the surface reflectance model, the aerosol total column optical thickness and single-scattering properties. The CISAR algorithm functioning is illustrated with a series of simple experiments.

1 Introduction

Radiative coupling between atmospheric scattering and surface reflectance processes prevents the use of linear relationships for the retrieval of aerosol properties over land surfaces. The discrimination between the contribution of the signal reflected by the surface and that scattered by aerosols represents one of the major issues when retrieving aerosol properties using space-borne passive optical observations over land surfaces. Conceptually, this problem can be modelled to solve a radiative system composed of at least two sets of layers, where the upper layers include aerosols and the bottom ones represent the soil–vegetation strata. The problem can be further complicated by the intrinsic anisotropic radiative behaviour of natural surfaces due to the mutual shadowing of the scattering elements, which is also affected by the amount of incident radiation . In most cases, an increase in aerosol concentration is responsible for an increase in the fraction of diffuse sky radiation which, in turn, smooths the effects of surface reflectance anisotropy. Though multispectral information is critical for the retrieval of aerosol properties, the spectral dimension alone does not allow full characterisation of the underlying surface reflectance, which often offers a significant contribution to the total signal observed at the satellite level. In this regard, the additional information contained in multispectral and multi-angular observations has proven essential to characterising aerosol properties over land surfaces.

pioneered the development of a retrieval method dedicated to the joint retrieval of surface reflectance and aerosol properties based on the inversion of a physically based radiative model. This method has been subsequently improved to permit the processing of any geostationary satellites accounting for their actual radiometric performance . This new versatile version of Pinty's algorithm has permitted the generation of a global surface albedo product from archived data acquired by operational geostationary satellites around the globe . These data included observations acquired by an old generation of radiometers with only one broad solar channel on board the European Meteosat First Generation satellite, the US Geosynchronous Operational Environmental Satellite (GOES) and the Japanese Geostationary Meteorological Satellite (GMS). It is now routinely applied in the framework of the Sustained and COordinated Processing of Environmental satellite data for Climate Monitoring (SCOPE-CM) initiative for the generation of essential climate variables . An improved version of this algorithm has been proposed by to take advantage of the multispectral capabilities of Meteosat Second Generation Spinning Enhanced Visible and Infrared Imager (MSG SEVIRI) operated by EUMETSAT and includes an optimal estimation (OE) inversion scheme using a minimisation approach based on the Marquardt–Levenberg method .

The strengths and weaknesses of the algorithm proposed by are discussed in Sect. 2. In their approach, the solutions of the radiative transfer equation (RTE) are pre-calculated and stored in look-up tables (LUTs) for a limited number of state variable values. Aerosol properties are limited to six different models dominated either by fine or coarse particles. Two major drawbacks result from the use of predefined aerosol models stored in precomputed LUTs. Firstly, only a limited region of the solution space is sampled as a result of the reduced range of variability for state variables stored in the LUTs. For instance, in order to reduce the size of the LUTs, limit the maximum aerosol optical thickness to 1. Secondly, the use of predefined aerosol models constitutes a major drawback since the solution space is not continuously sampled.

and , among others, demonstrated the advantages of a retrieval approach based on continuous variations of the aerosol properties as opposed to a LUT-based approach relying on a set of predefined aerosol models. Even considering a large number of aerosol models, LUT-based approaches underperform compared to methods with multivariate continuity in the solution space .

A new joint surface reflectance–aerosol properties retrieval approach is presented here that overcomes the limitations resulting from precomputed RTE solutions stored in LUTs.

The advantages of a continuous variation of the aerosol properties in the solution space against a LUT-based approach is discussed in Sect. 3. The proposed method expresses the single-scattering albedo and phase function values as a linear mixture of basic aerosol models. The forward radiative transfer model that includes the Jacobians, i.e. the partial derivative, is described in Sect. 4. With the exception of gaseous transmittance, this model no longer relies on LUTs, and the RTE is explicitly solved. The inversion method is described in Sect. 5. Finally, the ability to express aerosol single-scattering properties as a linear combination is in illustrated Sect. 6 with simulated data representing various scenarios, including small and large particles. Practical aspects of the application of the CISAR algorithm for the retrieval of both surface and aerosol properties from actual satellite data are addressed in (hereafter referred to as Part 2).

2 Lessons learned from previous approaches

proposed an algorithm for the joint retrieval of surface reflectance and aerosol properties to demonstrate the possibility of generating essential climate variables (ECVs) from data acquired by operational weather geostationary satellites. Due to limited operational computational resources available at that time in the EUMETSAT ground segment, where the data were processed, the development of this algorithm was subject to strong constraints. The RTE solutions were precomputed and stored in LUTs with a very coarse resolution, limiting the maximum aerosol optical thickness (AOT) to 1, which represented a severe limitation over the Sahara region where AOT values can easily exceed such a limit. Furthermore, the radiative coupling between aerosol scattering and gaseous absorption was not taken into account. This algorithm, referred to as geostationary surface albedo (GSA), has been subsequently modified by to include an estimation of the retrieval uncertainty. This updated version has permitted the generation of a global aerosol product derived from observations acquired by operational weather geostationary satellites . Since then, it has been routinely applied in the framework of the SCOPE-CM initiative to generate a climate data record (CDR) of surface albedo .

The GSA algorithm has been further improved for the processing of SEVIRI data on board MSG for the retrieval of the total column AOT from observations acquired in three solar bands centred at 0.6, 0.8 and 1.6 µm . The method developed by these authors relies on an OE approach wherein surface reflectance and daily aerosol load are simultaneously retrieved. The inversion is performed independently for each aerosol model and the one with the smallest cost function is selected. A physically based radiative transfer model accounting for non-Lambertian surface reflectance and its radiative coupling with atmospheric scattering is inverted against daily accumulated SEVIRI observations. However, this land daily aerosol (LDA) algorithm suffers from two major limitations: (i) the use of predefined aerosol models, and (ii) the algorithm delivers only one mean aerosol value per day when applied to MSG SEVIRI data . This latter issue has been addressed by , who retrieve an aerosol optical thickness value for each SEVIRI observation. The former issue prevents a continuous variation of the state variables characterising the aerosol single-scattering properties as required to find the minimum of the cost function. A consistent implementation of such an approach is not straightforward since aerosol models are defined as prior knowledge of the observed medium but uncertainties cannot be easily assigned to this decision. Consequently, the estimated retrieval uncertainty is inconsistent as it does not account for the use of prior information and its associated uncertainties.

Figure 1Aerosol dual-mode models based on in the {g,ω0} space derived from the aggregation of aerosol single-scattering properties retrieved from AERONET observations . Classes 1 to 3 are dominated by the fine mode and 4 to 6 by the coarse one.

demonstrated the advantages of a retrieval method based on continuous variations of aerosol single-scattering properties in the solution space as opposed to a LUT-based approach derived for a limited number of predefined aerosol models. proposed an original method for the retrieval of aerosol microphysical properties, which also does not necessitate the use of predefined aerosol models. This method retrieved more than 100 state variables, therefore requiring a considerable number of observations, such as those provided by multi-angular and -polarisation radiometers like Polarisation et Anisotropie des Réflectances Au SOmmet de l'Atmosphère (PARASOL) or the future Multi-viewing Multi-channel Multi-polarization Imaging (3MI) instrument on board EUMETSAT's Polar System Second Generation . Instruments delivering such a large number of observations are rather scarce, as most of the current passive optical sensors do not offer instantaneous multi-angular observation capabilities nor information on polarization. The primary objective of this paper is to address the limitations resulting from conventional approaches based on LUTs and/or a limited number of predefined aerosol models by proposing a method that can be applied to observations acquired by single or multi-view instruments.

3 Continuous variation of aerosol properties in the solution space

Aerosol single-scattering properties include the single-scattering albedo ω0 and the phase function Φ in RTE. explained the benefits of representing predefined aerosol models in a two-dimensional solution space composed of these aerosol single-scattering properties. For the sake of clarity, they limited the phase function in that 2-D space to the first term of the Legendre coefficients, i.e. the asymmetry parameter g. However, one should keep in mind that the reasoning applied in this section should be applied to the entire phase function Φ. These aerosol single-scattering properties are themselves determined by aerosol microphysical properties, such as the particle size distribution, shape and their complex index of refraction. The objective of retrievals that assume aerosol models is to provide a reasonable sampling of the {g,ω0} space. Omitting areas of that space may produce biased retrievals, as discussed in . The inversion process proposed by in this paper relies on a set of six models which have been defined from AErosol RObotic NETwork (AERONET) data aggregation . That choice of models was intended to provide a sampling of solution space representative of real-world conditions. The inversion is repeated for each aerosol class and the result with the best fit is reported, rather than having continuously varying aerosol properties, as would be preferable.

Figure 2Example of sensitivity of aerosol single-scattering properties to particle median radius (green arrows) and imaginary part of the refractive index (red arrows) at 0.44 and 0.87 µm for fine-mode, F (rmf=0.1µm), and coarse mode, C (rmc=2.0µm). The length of the arrows reflects the magnitude of the change.

A visual inspection of Fig. 1 based on reveals that aerosol models occupy different regions in the {g,ω0} space according to the dominant particle size distribution, i.e. fine or coarse. Within that space, an aerosol model is defined by the spectral behaviour of the {g(λ),ω0(λ)} pairs, where λ indicates the wavelength. The proposed fine-mode models vary mostly as a function of ω0, which is largely determined by the imaginary part of the refractive index ni. Conversely, aerosol models dominated by coarse particles show little dependency on g and are therefore organised parallel to the ordinate axis. The main parameter discriminating these latter models is the median radius rm, which essentially determines the asymmetry parameter value at a given wavelength.

To illustrate the dependence of g and ω0 on the median radius rm (http://eodg.atm.ox.ac.uk/user/grainger/research/aerosols.pdf, last access: (9 December 2018) and imaginary part of the refractive index ni, fine- and coarse-mono-mode aerosol models were generated with rm=0.15 and 2.0 µm respectively. The other microphysical values have been fixed to σr=0.5µm, nr=1.42 and ni=0.008, where σr is the radius standard deviation and nr the real part of the refractive index. These values were selected to ease the explanation of the organisation of the aerosol models in Fig. 1. Black dots in Fig. 2 show the corresponding location of {g,ω0} at 0.44 and 0.87 µm. The magnitude of the red arrows illustrate the sensitivity to a ni change of ±0.0025 and the green ones to a rm change of ±25 %. For the fine mono-mode (F), changes in ni essentially translate in displacement along the ω0 axis, while changes in rm result in changes almost parallel to the g axis. There is also a clear relationship between the particle size and g for that mode. A change in the particle size results in a change in g, while ω0 remains relatively unchanged. The situation is quite different for the coarse mono-mode, where changes in both ni and rm induce displacement parallel to the ω0 axis with limited impact on g values.

Figure 3Example of a region (light-blue area) in the {g,ω0} solution space at 0.44 µm defined by four aerosol vertices: single fine-mode non-absorbing (FN), single fine-mode absorbing (FA), coarse mode with small radius (CS) and coarse mode with large radius (CL). The isolines show the probability that the aerosol single-scattering properties derived from AERONET observations with the method of fall within the delineated spaces.

The actual extent of solutions in the {g,ω0} space for a given spectral band can be outlined by a series of vertices defined by aerosol single-scattering properties (Fig. 3). Following Fig. 2, these vertices are defined by absorbing and non-absorbing fine mono-mode models with small radii of about 0.1 µm, labelled respectively FA and FN, and by two coarse mono-modes with different radii, i.e. large (1 µm) and small (0.3 µm), labelled respectively CL and CS. In Sect. 4, we will see how any pair of single-scattering albedo and phase function values can be expressed as a linear combination of the vertex properties.

The position of these vertices is critical as they should encompass the most likely aerosol single-scattering properties that could be observed at a given time and location. Different approaches could be used to define the position of these vertices. The positions can be derived from the analysis of typical aerosol single-scattering properties available in databases such as the Optical Properties of Aerosols and Clouds (OPAC) . Alternatively, it is also possible to follow a similar approach to the one proposed in , who analysed the single-scattering albedo and phase function values derived from AERONET observations acquired in a specific region of interest for a given period . The red isoline in Fig. 3 delineates the area [g,ω0] of the solution space where 99.7 % of the aerosol single-scattering properties derived by from AERONET observations are located. The green and blue lines respectively show the 95 % and 68 % probability regions. These values have been derived using all available level 2 AERONET observations since 1993. Finally, the model proposed by can be used to determine the spectral variations of the single-scattering properties outside the spectral bands measured by AERONET. This work relies on simulated data and the aerosol vertices have been positioned to sample the solution space in a realistic way. When processing actual satellite data over a specific region or period, it is advised to calculate the isolines corresponding to that region of interest from AERONET observations and to adjust the position of the aerosol vertices accordingly as performed in Part 2.

## 4.1 Overview

The forward model, named FASTRE, simulates the TOA bidirectional reflectance factor (BRF) ${y}_{m}\left(\mathbf{x},\mathbf{b};\mathbf{m}\right)$ as a function of the independent parameters m defining the observation conditions and a series of state variables x describing the state of the atmosphere and underlying surface. Model parameters b represent variables such as total column water vapour that influence the value of ${y}_{m}\left(\mathbf{x},\mathbf{b};\mathbf{m}\right)$ but cannot be retrieved from the processed space-based observations due to the lack of information. The independent parameters m include the illumination and viewing geometries 0v) and the wavelength dependence. The RTE is solved with the matrix operator method optimised by for a limited number of quadrature points.

The model simulates observations acquired within spectral bands $\stackrel{\mathrm{̃}}{\mathit{\lambda }}$ characterised by their spectral response. Gaseous transmittances in these bands are precomputed and stored in LUTs. The model computes the contributions from single and multiple scattering separately, the latter being solved in Fourier space. In order to reduce the computation time, the forward model relies on the same atmospheric vertical structure as in , i.e. a three-level system containing two layers (Fig. 4). The lowest level, Z0, represents the surface. The lower layer, La, ranging from levels Z0 to Za, contains the aerosol particles. Molecular scattering and absorption also take place in that layer, which is radiatively coupled with the surface for both single and the multiple scattering. The upper layer, Lg, ranging from Za to Zs, is only subject to molecular absorption.

Figure 4Atmospheric vertical structure of the FASTRE model. The surface is at level Z0 and radiatively coupled with the lower layer La extending from level Z0 to Za. This layer includes scattering and absorption processes. The upper layer, Lg, runs from level Za to Zs and only accounts for gas absorption processes.

The surface reflectance ${r}_{\mathrm{s}}\left({\mathbf{x}}_{\mathbf{s}},\mathbf{b};\mathbf{m}\right)$ over land is represented by the so-called RPV (Rahman–Pinty–Verstraete) model characterised by four parameters, ${\mathbf{x}}_{\mathbf{s}}=\mathit{\left\{}{\mathit{\rho }}_{\mathrm{0}},k,\mathrm{\Theta },{\mathit{\rho }}_{c}\mathit{\right\}}$, that are all wavelength dependent . The ρ0 parameter, included in the [0,1] interval, controls the mean amplitude of the BRF and strongly varies with wavelengths. The k parameter is the modified Minnaert's contribution that determines the bowl or bell shape of the BRF and typically varies between 0 and 2. The asymmetry parameter Θ of the Henyey–Greenstein phase function varies between −1 and 1. The ρc parameter controls the amplitude of the hotspot due to the “porosity” of the medium. This parameter varies between −1 and 1. For the simulations over the ocean, the Cox–Munk model is used as implemented in .

Aerosol single-scattering properties in the layer La are represented by an external mixture of a series of predefined aerosol vertices as explained in Sect. 4.2. The Lg layer only contains absorbing gas not included in the scattering layer, such as high-altitude ozone, the part of the total column water vapour not included in layer La and a few well-mixed gases.

The FASTRE model expresses the TOA BRF in a given spectral band $\stackrel{\mathrm{̃}}{\mathit{\lambda }}$ as a sum of the single ${I}_{\mathrm{s}}^{↑}$ and multiple ${I}_{\mathrm{m}}^{↑}$ scattering contributions as in

$\begin{array}{}\text{(1)}& {y}_{m}\left(\mathbf{x},\mathbf{b};\mathbf{m}\right)={T}_{{L}_{g}}\left(\mathbf{b};\mathbf{m}\right)\frac{{I}_{\mathrm{s}}^{↑}\left(\mathbf{x},\mathbf{b};\mathbf{m}\right)+{I}_{\mathrm{m}}^{↑}\left(\mathbf{x},\mathbf{b};\mathbf{m}\right)}{{E}_{\mathrm{0}}^{↓}\left(\mathbf{m}\right){\mathit{\mu }}_{\mathrm{0}}},\end{array}$

where

• ${I}_{\mathrm{s}}^{↑}\left(\mathbf{x},\mathbf{b};\mathbf{m}\right)$ is the upward radiance field at level Za due to the single scattering,

• ${I}_{\mathrm{m}}^{↑}\left(\mathbf{x},\mathbf{b};\mathbf{m}\right)$ is the upward radiance field at level Za due to the multiple scattering,

• ${T}_{{L}_{g}}\left(\mathbf{b};\mathbf{m}\right)$ denotes the total transmission factor in the Lg layer,

• ${E}_{\mathrm{0}}^{↓}\left(\mathbf{m}\right)$ denotes the solar irradiance at level Zs corrected for the Sun–Earth distance variations.

The single-scattering contribution is written

$\begin{array}{ll}& {I}_{\mathrm{s}}^{↑}\left(\mathbf{x},\mathbf{b};\mathbf{m}\right)=\\ \text{(2)}& & \frac{{E}_{\mathrm{0}}^{↓}\left(\mathbf{m}\right){\mathit{\mu }}_{\mathrm{0}}}{\mathit{\pi }}\mathrm{exp}\left(\frac{-\mathit{\tau }{}_{{L}_{a}}}{{\mathit{\mu }}_{\mathrm{0}}}\right){r}_{\mathrm{s}}\left({\mathbf{x}}_{\mathbf{s}},\mathbf{b};\mathbf{m}\right)\mathrm{exp}\left(\frac{-\mathit{\tau }{}_{{L}_{a}}}{{\mathit{\mu }}_{\mathrm{v}}}\right),\end{array}$

where $\mathit{\tau }{}_{{L}_{a}}$ is the total optical thickness of layer La. μ0 and μv are the cosine of the illumination and viewing zenith angles respectively.

The multiple-scattering contribution, ${I}_{\mathrm{m}}^{↑}\left(\mathbf{x},\mathbf{b};\mathbf{m}\right)$, is solved in the Fourier space for all illumination and viewing directions of the quadrature directions Nθ for 2Nθ−1 azimuthal directions. The contribution ${I}_{\mathrm{m}}^{↑}\left(\mathbf{x},\mathbf{b};\mathbf{m}\right)$ in the direction 0v) is interpolated from the surrounding quadrature directions. Finally, the Jacobian ${\mathbf{k}}_{{x}_{i}}=\frac{\partial {y}_{m}\left({\mathbf{x}}_{\mathbf{i}},\mathbf{b};\mathbf{m}\right)}{\partial {x}_{i}}$ of ${y}_{m}\left(\mathbf{x},\mathbf{b};\mathbf{m}\right)$ for parameter xi are calculated as finite differences.

## 4.2 Scattering layer La properties

The layer La contains a set of mono-mode aerosol models v characterised by their single-scattering properties, i.e. the single-scattering albedo ${\mathit{\omega }}_{\mathrm{0},v}\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)$ and phase function ${\mathrm{\Phi }}_{v}\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }},{\mathrm{\Omega }}_{g}\right)$ where Ωg represents the scattering angle. The different vertices are combined into this layer according to their respective optical thickness ${\mathit{\tau }}_{v}\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)$ with the total aerosol optical thickness ${\mathit{\tau }}_{a}\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)$ of the layer being equal to

$\begin{array}{}\text{(3)}& {\mathit{\tau }}_{a}\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)=\sum _{v}{\mathit{\tau }}_{v}\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right).\end{array}$

The phase function ${\mathrm{\Phi }}_{v}\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }},{\mathrm{\Omega }}_{g}\right)$ is characterised by a limited number Nκ of Legendre coefficients equal to 2Nθ−1. The decision to use this number results from a trade-off between accuracy and computational time. When Nκ is too small, the last Legendre moment is often not equal to zero and the delta-M approximation is applied (Wiscombe1977). In this case, the αd coefficient of the delta-M approximation is equal to Φv(Nκ). The Legendre coefficients κj become

$\begin{array}{}\text{(4)}& {c}_{j}=\frac{{\mathit{\kappa }}_{j}-{\mathit{\alpha }}_{d}}{\mathrm{1}-{\mathit{\alpha }}_{d}},\end{array}$

and the truncated phase function is denoted by ${\mathrm{\Phi }}_{v}^{\prime }$. The corrected optical thickness ${\mathit{\tau }}_{v}^{\prime }\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)$ and single-scattering albedo ${\mathit{\omega }}_{\mathrm{0},v}^{\prime }\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)$ of the corresponding aerosol model become

$\begin{array}{}\text{(5)}& {\mathit{\tau }}_{v}^{\prime }\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)=\left(\mathrm{1}-{\mathit{\omega }}_{\mathrm{0},v}{\mathit{\alpha }}_{d}\right){\mathit{\tau }}_{v}\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)\end{array}$

and

$\begin{array}{}\text{(6)}& {\mathit{\omega }}_{\mathrm{0},v}^{\prime }\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)=\frac{\mathrm{1}-{\mathit{\alpha }}_{d}}{\mathrm{1}-{\mathit{\omega }}_{\mathrm{0},v}{\mathit{\alpha }}_{d}}\phantom{\rule{0.125em}{0ex}}{\mathit{\omega }}_{\mathrm{0},v}\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right).\end{array}$

The layer total optical thickness, ${\mathit{\tau }}_{{L}_{a}}$, is the sum of the gaseous, τg, the aerosol, ${\mathit{\tau }}_{\mathrm{a}}^{\prime }$ and the Rayleigh, τr, optical depth

$\begin{array}{}\text{(7)}& {\mathit{\tau }}_{{L}_{a}}\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)={\mathit{\tau }}_{\mathrm{g}}\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)+{\mathit{\tau }}_{\mathrm{a}}^{\prime }\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)+{\mathit{\tau }}_{\mathrm{r}}\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right),\end{array}$

with ${\mathit{\tau }}_{\mathrm{a}}^{\prime }\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)={\sum }_{v}{\mathit{\tau }}_{v}^{\prime }\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)$. The single-scattering albedo of the scattering layer is equal to

$\begin{array}{}\text{(8)}& {\mathit{\omega }}_{\mathrm{0}}^{\prime }\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)=\frac{{\sum }_{c}{\mathit{\omega }}_{\mathrm{0},v}^{\prime }\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right){\mathit{\tau }}_{v}^{\prime }\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)}{{\mathit{\tau }}_{\mathrm{a}}^{\prime }\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)}\end{array}$

and the layer average phase function

$\begin{array}{}\text{(9)}& {\mathrm{\Phi }}^{\prime }\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }},{\mathrm{\Omega }}_{g}\right)=\frac{{\sum }_{c}{\mathrm{\Phi }}_{v}^{\prime }\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }},{\mathrm{\Omega }}_{g}\right){\mathit{\tau }}_{v}^{\prime }\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)}{{\mathit{\tau }}_{\mathrm{a}}^{\prime }\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)}.\end{array}$

## 4.3 Gaseous layer properties

It is assumed that only molecular absorption takes place in layer Lg. The height of level Za is used to partition the total column water vapour and ozone concentration in each layer assuming a US76 standard atmosphere vertical profile. This height is not retrieved and is therefore a model parameter of FASTRE, which should be derived from some climatological values. ${T}_{{L}_{g}}$ denotes the total transmission of that layer.

Table 1Relative bias and root mean square error in percentage between FASTRE and the reference RTM in various spectral bands.

## 4.4 FASTRE model accuracy

The simple atmospheric vertical structure composed of two layers is the principal assumption of the FASTRE model. In order to evaluate the accuracy of FASTRE, a similar procedure to that in has been applied. The outcome of FASTRE has been evaluated against a more elaborated 1-D radiative transfer model (RTM) (Govaerts2006) for sun and viewing angles varying from 0 to 70, for various types of aerosols, surface reflectance and total column water vapour values. This reference RTM represents the vertical structure of the atmosphere with 50 layers. The mean relative bias and relative root mean square error (RMSE) between the reference model and FASTRE have been estimated in the main spectral bands used for aerosol retrievals. The relative RMSE, Rr, is estimated as

$\begin{array}{}\text{(10)}& {R}_{\mathrm{r}}=\sqrt{\frac{\mathrm{1}}{N}\sum _{N}{\left(\frac{{y}_{m}\left(\mathbf{x},\mathbf{b};\mathbf{m}\right)-{y}_{\mathrm{r}}\left(\mathbf{x},\mathbf{b};\mathbf{r}\right)}{{y}_{\mathrm{r}}\left(\mathbf{x},\mathbf{b};\mathbf{r}\right)}\right)}^{\mathrm{2}}},\end{array}$

where ${y}_{\mathrm{r}}\left(\mathbf{x},\mathbf{b};\mathbf{m}\right)$ is the TOA BRF calculated with the reference model. In this paper, the FASTRE model solves the RTE using 16 quadrature points Nθ, which provides a good compromise between speed and accuracy. Results are shown in Table 1. The relative RMSE between FASTRE and the reference model is typically in the range of 1 %–3 %. Another comparison of FASTRE has been made against actual Project for On-Board Autonomy-Vegetation (PROBA-V) observations . These comparisons show an RMSE in the range [0.024–0.038].

5 Inversion process

## 5.1 Overview

Surface reflectance characterisation requires multi-angular observations ${\mathbf{y}}_{\mathrm{\Omega }\stackrel{\mathrm{̃}}{\mathrm{\Lambda }}}$, the acquisition of which can take between several minutes, as is the case for the Multi-angle Imaging SpectroRadiometer (MISR) instrument, and several days, as is the case for the Ocean and Land Colour Instrument (OLCI) on board Sentinel-3 or the Moderate Resolution Imaging Spectroradiometer (MODIS). In the former case, data are often assumed to be acquired almost instantaneously, i.e. with the atmospheric properties remaining unchanged during the acquisition time. Such a situation considerably reduces the calculation time required to solve the RTE, as the multiple scattering term ${I}_{\mathrm{m}}^{↑}\left(\mathbf{x},\mathbf{b};\mathbf{m}\right)$ needs to be estimated only once per spectral band. In the latter case, atmospheric properties cannot be assumed to be invariant and the multiple scattering contribution needs to be solved for each observation. When geostationary observations are processed, the accumulation period is often reduced to 1 day, and the assumption that the atmosphere does not change can be converted into an equivalent radiometric uncertainty . Strictly speaking, it should be assumed that atmospheric properties have changed when the accumulation time exceeds several minutes .

The retrieved state variables in each spectral band $\stackrel{\mathrm{̃}}{\mathit{\lambda }}$ are composed of the xs parameters characterising the state of the surface and the set of aerosol optical thicknesses τv for the aerosol vertices that are mixed in layer La. Prior information consists of the expected values xb of the state variables x characterising the surface and the atmosphere on one side, and regularisation of the spectral and/or temporal variability of τv on the other side. Uncertainty matrices Sx are assigned to this prior information. Finally, uncertainties in the measurements Sy are assumed to be normally distributed with zero mean. The inversion process of the FASTRE model will be herein referred to as the Combined Inversion of Surface and AeRosol (CISAR) algorithm.

## 5.2 Cost function

The fundamental principle of OE is to maximise the probability $P=P\left(\mathbf{x}|{\mathbf{y}}_{\mathrm{\Omega }\stackrel{\mathrm{̃}}{\mathrm{\Lambda }}},{\mathbf{x}}_{\mathbf{b}},\mathbf{b}\right)$ with respect to the values of the state vector x, conditional to the value of the measurements and any prior information. The conditional probability takes on the quadratic form (Rodgers2000):

$\begin{array}{ll}P\left(\mathbf{x}\right)\propto & \phantom{\rule{0.25em}{0ex}}\mathrm{exp}\left[-\left({y}_{m}\left(\mathbf{x},\mathbf{b};\mathbf{m}\right)-{\mathbf{y}}_{\mathrm{\Omega }\stackrel{\mathrm{̃}}{\mathrm{\Lambda }}}{\right)}^{T}{\mathbf{S}}_{y}^{-\mathrm{1}}\left({y}_{m}\left(\mathbf{x},\mathbf{b};\mathbf{m}\right)-{\mathbf{y}}_{\mathrm{\Omega }\stackrel{\mathrm{̃}}{\mathrm{\Lambda }}}\right)\right]+\\ & \mathrm{exp}\left[-\left(\mathbf{x}-{\mathbf{x}}_{b}{\right)}^{T}{\mathbf{S}}_{x}^{-\mathrm{1}}\left(\mathbf{x}-{\mathbf{x}}_{b}\right)\right]+\\ & \mathrm{exp}\left[-{\mathbf{x}}^{T}{\mathbf{H}}_{a}^{T}{\mathbf{S}}_{a}^{-\mathrm{1}}{\mathbf{H}}_{a}\mathbf{x}\right]+\\ \text{(11)}& & \mathrm{exp}\left[-{\mathbf{x}}^{T}{\mathbf{H}}_{l}^{T}{\mathbf{S}}_{l}^{-\mathrm{1}}{\mathbf{H}}_{l}\mathbf{x}\right],\end{array}$

where the first two terms represent weighted deviations from measurements and the prior state parameters respectively, the third represents the AOT temporal smoothness constraints and the fourth represents the AOT spectral constraint, with respective uncertainty matrices Sa and Sl. The algorithm proposed by implements similar temporal and spectral smoothness constraints. The two matrices Ha and Hl, representing respectively the temporal and spectral constraints, can be written as block diagonal matrices:

$\begin{array}{}\text{(12)}& \mathbf{H}=\left(\begin{array}{ccccc}{\mathbf{H}}^{{\mathit{\rho }}_{\mathrm{0}}}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbf{H}}^{k}& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathbf{H}}^{\mathit{\theta }}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& {\mathbf{H}}^{{\mathit{\rho }}_{c}}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& {\mathbf{H}}^{\mathit{\tau }}\end{array}\right),\end{array}$

where the four blocks ${\mathbf{H}}^{{\mathit{\rho }}_{\mathrm{0}}}$, Hk, Hθ and ${\mathbf{H}}^{{\mathit{\rho }}_{c}}$ express the spectral constraints between the surface parameters. Their values are set to zero when these constraints are not active. The submatrix ${\mathbf{H}}_{a}^{\mathit{\tau }}$ can also be written using blocks ${\mathbf{H}}_{a;\stackrel{\mathrm{̃}}{\mathit{\lambda }},v}^{\mathit{\tau }}$ along the diagonal. For a given spectral band $\stackrel{\mathrm{̃}}{\mathit{\lambda }}$ and aerosol vertex v, the block ${\mathbf{H}}_{a;\stackrel{\mathrm{̃}}{\mathit{\lambda }},v}^{\mathit{\tau }}$ is defined as follows:

$\begin{array}{}\text{(13)}& {\mathbf{H}}_{a;\stackrel{\mathrm{̃}}{\mathit{\lambda }},v}^{\mathit{\tau }}\phantom{\rule{0.25em}{0ex}}{\mathbit{\tau }}_{\stackrel{\mathrm{̃}}{\mathit{\lambda }},v}=\left(\begin{array}{ccccc}\mathrm{1}& -\mathrm{1}& \mathrm{0}& \mathrm{\dots }& \mathrm{\dots }\\ \mathrm{0}& \mathrm{1}& -\mathrm{1}& \mathrm{0}& \mathrm{\dots }\\ \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }\\ \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{1}& -\mathrm{1}\\ \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{0}\end{array}\right)\phantom{\rule{0.25em}{0ex}}\left(\begin{array}{c}{\mathit{\tau }}_{\stackrel{\mathrm{̃}}{\mathit{\lambda }},v,\mathrm{1}}\\ {\mathit{\tau }}_{\stackrel{\mathrm{̃}}{\mathit{\lambda }},v,\mathrm{2}}\\ \mathrm{⋮}\\ {\mathit{\tau }}_{\stackrel{\mathrm{̃}}{\mathit{\lambda }},v,{N}_{t}-\mathrm{1}}\\ {\mathit{\tau }}_{\stackrel{\mathrm{̃}}{\mathit{\lambda }},v,\mathrm{1},{N}_{t}}\end{array}\right).\end{array}$

In the same way, the submatrix ${\mathbf{H}}_{l}^{\mathit{\tau }}$ can be written using blocks ${\mathbf{H}}_{l;v,t}^{\mathit{\tau }}$. For a given aerosol vertex v and time t, the block ${\mathbf{H}}_{l;v,t}^{\mathit{\tau }}$ is defined as follows:

$\begin{array}{}\text{(14)}& {\mathbf{H}}_{l;v,t}^{\mathit{\tau }}\phantom{\rule{0.25em}{0ex}}{\mathbit{\tau }}_{v,t}=\left(\begin{array}{ccccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}\\ -\frac{{\mathit{ϵ}}_{\mathrm{2}}}{{\mathit{ϵ}}_{\mathrm{1}}}& \mathrm{1}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{0}& -\frac{{\mathit{ϵ}}_{\mathrm{3}}}{{\mathit{ϵ}}_{\mathrm{2}}}& \mathrm{1}& \mathrm{\dots }& \mathrm{0}\\ \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& \mathrm{\ddots }& \mathrm{0}\\ \mathrm{\dots }& \mathrm{\dots }& \mathrm{\dots }& -\frac{{\mathit{ϵ}}_{{N}_{\mathit{\lambda }}}}{{\mathit{ϵ}}_{{N}_{\mathit{\lambda }}-\mathrm{1}}}& \mathrm{1}\end{array}\right)\phantom{\rule{0.25em}{0ex}}\left(\begin{array}{c}{\mathit{\tau }}_{\mathrm{1},v,t}\\ {\mathit{\tau }}_{\mathrm{2},v,t}\\ {\mathit{\tau }}_{\mathrm{3},v,t}\\ \mathrm{⋮}\\ {\mathit{\tau }}_{{N}_{\stackrel{\mathrm{̃}}{\mathit{\lambda }}},v,t}\end{array}\right),\end{array}$

where the ϵl represents the uncertainties associated with the AOT spectral constraints of the individual vertex v, bounding the solution space. The spectral variations of τv between band ${\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{l}$ and ${\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{l+\mathrm{1}}$ are assumed to vary, as

$\begin{array}{}\text{(15)}& \frac{{\mathit{\tau }}_{{\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{l},v}}{{\mathit{\tau }}_{{\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{l+\mathrm{1}},v}}=\frac{{e}_{{\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{l}}}{{e}_{{\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{l+\mathrm{1}}}},\end{array}$

where ${e}_{{\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{l}}$ is the extinction coefficient in band ${\stackrel{\mathrm{̃}}{\mathit{\lambda }}}_{l}$.

Maximising the probability function in Eq. (11) is equivalent to minimising the negative logarithm

$\begin{array}{}\text{(16)}& J\left(\mathbf{x}\right)={J}_{y}\left(\mathbf{x}\right)+{J}_{x}\left(\mathbf{x}\right)+{J}_{a}\left(\mathbf{x}\right)+{J}_{l}\left(\mathbf{x}\right),\end{array}$

with

$\begin{array}{}\text{(17)}& {J}_{y}\left(\mathbf{x}\right)& =\left({y}_{m}\left(\mathbf{x},\mathbf{b},\mathrm{\Omega }\right)-{\mathbf{y}}_{\mathrm{\Omega }\stackrel{\mathrm{̃}}{\mathrm{\Lambda }}}\right){\mathbf{S}}_{y}^{-\mathrm{1}}\left({y}_{m}\left(\mathbf{x},\mathbf{b},\mathrm{\Omega }\right)-{\mathbf{y}}_{\mathrm{\Omega }\stackrel{\mathrm{̃}}{\mathrm{\Lambda }}}{\right)}^{T}\text{(18)}& {J}_{x}\left(\mathbf{x}\right)& =\left(\mathbf{x}-{\mathbf{x}}_{b}\right){\mathbf{S}}_{x}^{-\mathrm{1}}\left(\mathbf{x}-{\mathbf{x}}_{b}{\right)}^{T}\text{(19)}& {J}_{a}\left(\mathbf{x}\right)& ={\mathbf{x}}^{T}{\mathbf{H}}_{a}^{T}{\mathbf{S}}_{a}^{-\mathrm{1}}{\mathbf{H}}_{a}\mathbf{x}\text{(20)}& {J}_{l}\left(\mathbf{x}\right)& ={\mathbf{x}}^{T}{\mathbf{H}}_{l}^{T}{\mathbf{S}}_{l}^{-\mathrm{1}}{\mathbf{H}}_{l}\mathbf{x}.\end{array}$

Notice that the cost function J is minimised with respect to the state variable x, so that the derivative of J is independent of the model parameters b. The need for angular sampling to document the surface anisotropy leads to an unbalanced dimension of nx and ny with ny>nx, where ny and nx represent the number of observations and state variables respectively. According to , these additional observations should improve the retrieval as, from a statistical point of view, repeating similar observations implies that the variance should decrease. Accordingly, the magnitude of the elements of the covariance matrix should decrease as $\mathrm{1}/\sqrt{{n}_{y}}$. Thus, repeating similar observations results in some enhancements of retrieval accuracy, which is proportional to the ratio nynx. Hence, the cost function which is actually minimised is ${J}_{s}\left(\mathbf{x}\right)={J}_{y}\left(\mathbf{x}\right)+{n}_{y}/{n}_{x}\left({J}_{x}\left(\mathbf{x}\right)+{J}_{a}\left(\mathbf{x}\right)+{J}_{l}\left(\mathbf{x}\right)\right)$.

## 5.3 Retrieval uncertainty estimation

The retrieval uncertainty is based on the OE theory, assuming linear behaviour of ${y}_{m}\left(\mathbf{x},\mathbf{b};\mathbf{m}\right)$ in the vicinity of the solution $\stackrel{\mathrm{^}}{\mathbf{x}}$. Under this condition, the retrieval uncertainty ${\mathit{\sigma }}_{\stackrel{\mathrm{^}}{\mathbf{x}}}$ is determined by the shape of J(x) at $\stackrel{\mathrm{^}}{\mathbf{x}}$

$\begin{array}{ll}& {\mathit{\sigma }}_{\stackrel{\mathrm{^}}{\mathbf{x}}}^{\mathrm{2}}={\left(\frac{{\partial }^{\mathrm{2}}{J}_{s}\left(\mathbf{x}\right)}{\partial {\mathbf{x}}^{\mathrm{2}}}\right)}^{-\mathrm{1}}=\\ \text{(21)}& & {\left({\mathbf{K}}_{x}^{T}\phantom{\rule{0.125em}{0ex}}{\mathbf{S}}_{y}^{-\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathbf{K}}_{x}+{\mathbf{S}}_{x}^{-\mathrm{1}}+{\mathbf{H}}_{a}^{T}{\mathbf{S}}_{a}^{-\mathrm{1}}{\mathbf{H}}_{a}+{\mathbf{H}}_{l}^{T}{\mathbf{S}}_{l}^{-\mathrm{1}}{\mathbf{H}}_{l}\right)}^{-\mathrm{1}},\end{array}$

where Kx is the Jacobian matrix of ${y}_{m}\left(\mathbf{x},\mathbf{b};\mathbf{m}\right)$ calculated in $\stackrel{\mathrm{^}}{\mathbf{x}}$. Combining Eqs. (21) and (8), the uncertainty in the retrieval of ω0 in band $\stackrel{\mathrm{̃}}{\mathit{\lambda }}$ is written

$\begin{array}{}\text{(22)}& {\mathit{\sigma }}_{\stackrel{\mathrm{^}}{{\mathit{\omega }}_{\mathrm{0}}}}^{\mathrm{2}}\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)=\sum _{v}{\left(\frac{{\mathit{\omega }}_{\mathrm{0},v}\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)-{\mathit{\omega }}_{\mathrm{0}}\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)}{{\mathit{\tau }}_{a}\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)}\right)}^{\mathrm{2}}{\mathit{\sigma }}_{\stackrel{\mathrm{^}}{{\mathit{\tau }}_{v}}}^{\mathrm{2}}\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right).\end{array}$

A similar equation can be derived for the estimation of ${\mathit{\sigma }}_{g}^{\mathrm{2}}$.

## 5.4 Acceleration methods

The minimisation of Eq. (16) relies on an iterative approach with ${y}_{m}\left(\mathbf{x},\mathbf{b};\mathbf{m}\right)$ and the associated Jacobians Kx being estimated at each iteration. In order to reduce the calculation time dedicated to the estimation of ${y}_{m}\left(\mathbf{x},\mathbf{b};\mathbf{m}\right)$ and Kx, a series of methods have been implemented. All quantities that do not explicitly depend on the state variables, such as the observation conditions m, model parameters b, quadrature point weights, etc. are computed only once prior to the optimisation.

When solving the RTE, the estimation of the multiple scattering term is by far the most time-consuming step. Hence, during the iterative optimisation process, when the change $\mathrm{\Delta }{\mathit{\tau }}_{a}\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right)$ between iteration j and j+1 is small, the multiple scattering contribution at iteration j+1 is estimated with

$\begin{array}{ll}& {I}_{\mathrm{m}}^{↑}\left({\mathit{\tau }}_{a}\left(j+\mathrm{1},\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right),\mathbf{b};\mathbf{m}\right)={I}_{\mathrm{m}}^{↑}\left({\mathit{\tau }}_{a}\left(j,\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right),\mathbf{b};\mathbf{m}\right)\\ \text{(23)}& & +\frac{\partial {I}_{\mathrm{m}}^{↑}\left({\mathit{\tau }}_{a}\left(j,\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right),\mathbf{b};\mathbf{m}\right)}{\partial {\mathit{\tau }}_{a}}\mathrm{\Delta }{\mathit{\tau }}_{a}\left(\stackrel{\mathrm{̃}}{\mathit{\lambda }}\right).\end{array}$

This approximation is not used in two successive iterations to avoid inaccurate results, and the single-scattering contribution is always explicitly estimated.

Table 2List of aerosol properties used for the simulations. The parameters rmf and rmc are the median fine- and coarse-mode radii expressed in µm. Their respective standard deviations are ${\mathit{\sigma }}_{{r}_{\mathrm{mf}}}$ and ${\mathit{\sigma }}_{{r}_{\mathrm{mc}}}$. The parameters nr and ni are the real and imaginary part of the refractive index in the indicated bands. Nf and Nc are the fine- and coarse-mode particle concentration in number of particles per cm3.

6 Algorithm performance evaluation

## 6.1 Experimental set-up

A simple experimental set-up based on simulated data has been defined to illustrate the performance of the CISAR algorithm as a function of the chosen solution space. More specifically, the capability of CISAR to continuously sample the {g,ω0} solution space is examined in detail. For the sake of simplicity, a noise-free multi-angular observation vector, ${\mathbf{y}}_{\mathrm{\Omega }\stackrel{\mathrm{̃}}{\mathrm{\Lambda }}}$, where Ω expresses the illumination and viewing geometries, is assumed to be acquired instantaneously in the principal plane and in the spectral bands listed in Table 1. A radiometric uncertainty of 3 % is assumed to compose Sy. In this ideal configuration, the solar zenith angle (SZA) is set to 30. In these experiments, to concentrate on the retrieval of aerosol properties, the surface parameters prior values are set to the true values used in the simulation (Table 5) with an ascribed uncertainty of 0.03. The first guess values are randomly chosen within this uncertainty interval. Part 2 explains how prior information on the surface parameters can be derived. No prior information is assumed for the aerosol optical thickness; i.e. the prior uncertainty is set to very large values. Only regularisation of the spectral variations of τa is applied.

Table 3Microphysical parameter values for the four vertices (FA, FN, CS and CL) in the selected spectral bands. Radius are given in µm.

The CISAR algorithm performance evaluation is based on a series of experiments corresponding to different selections of aerosol properties, both for the forward simulation of the observations and their inversion. Three different aerosol models are used in the forward simulations: F0, which only contains fine-mode; F1, which contains a dual-mode particle size distribution dominated by small particles; and F2, composed of a dual-mode distribution dominated by the coarse particles. Table 2 contains the values of the size distribution and refractive indices of these aerosol models. Values for the four vertices (FA, FN, CL and CS) enclosing the solution space as illustrated in Fig. 3 are given in Table 3. When the observations simulated with aerosol types F0, F1 or F2 are inverted, the list of vertices actually used depends on the type of experiment as indicated in Table 4. For all these scenarios, an AOT of 0.4 at 0.55 µm is assumed.

Table 4List of experiments: the names are provided in the first column. The active vertices in each experiment are indicated with the × symbol. The last column indicates the name of the aerosol model used to simulate the observations.

Table 5Values of the true and retrieved surface RPV parameters for experiment F00. Wavelengths are given in µm.

Table 6Retrieved AOT error and uncertainties for the six experiments. The error ϵτ is calculated as the difference between the retrieved and the true values, δτ is the relative error in percent and στ is the retrieval uncertainty estimated with Eq. (21). Wavelengths are given in µm.

Values used for the RPV parameters in the four selected bands are indicated in Table 5. They correspond to typical BRF values that would be observed over a vegetated surface with a leaf area index value of 3 and a bright underlying soil.

Figure 5(a) Results of experiment F00 in the {g,ω0} space. The aerosol vertices used for the inversion are FN (blue) and FA (green). The forward aerosol properties are shown in black and the retrieved ones in red. Vertical and horizontal red bars indicate the uncertainty of the retrieved values. (b) Retrieved AOT (red circles). The retrieval uncertainty is shown with the vertical red lines. True values are indicated with black crosses. True and retrieved values are slightly staggered to ease the reading.

Figure 6Same as Fig. 5 but for experiment F10.

## 6.2 Results

### 6.2.1 Experiment F00

The purpose of the first experiment (F00) is to demonstrate that the CISAR algorithm can accurately retrieve aerosol properties in a simple situation, showing therefore that the inversion process works correctly. The F0 aerosol model used to simulate the observations is only composed of fine particles with a median radius of 0.08 µm, i.e. the same value as for the FN and FA vertices used for the inversion. Hence, the retrieval is limited to the imaginary part of the index of refraction, the real part being set to 1.4. With a retrieval configuration restricted to the use of only two vertices, the solution space for each wavelength is limited to a straight line between the two vertices.

Results are shown in Fig. 5 for the atmosphere and Table 5 for the surface. The asymmetry factor g and single-scattering albedo ω0 are almost exactly retrieved. There is practically no uncertainty in the retrieval of g because of the constraints imposed by the fact that the particle radius is the same as for the F0 aerosol model. The retrieved AOT is also in very good agreement with the true values as can be seen on the right panel. The retrieval error ϵτ is defined as the difference between the retrieved and the true AOT values. Results are summarised in Table 6. This first experiment demonstrates that it is possible to retrieve the properties of the aerosol model F0 as a linear combination of the vertices FA and FN when only the absorption varies, the particle median radius being constant.

A comparison between the true and retrieved values in Table 5 shows that the surface parameters are very accurately retrieved. As stated in Sect. 6.1, prior information on the magnitude of the RPV parameter is assumed to be unbiased with an uncertainty of 0.03. The corresponding posterior uncertainties exhibit a significant decrease for the ρ0 parameter at all wavelengths. Similar behaviour is not observed for the other parameters. As explained in , the k and Θ parameters, controlling the surface reflectance anisotropy, are strongly correlated with the amount of atmospheric scattering. Consequently, the retrieved uncertainties decrease with increasing wavelengths, i.e. as a function of the actual AOT. Despite the observations taking place in the principal plane, the posterior uncertainty on the hot spot parameter remains equal to the prior one as a result of atmospheric scattering. This fact is attributed to the relatively high value of the true AOT, and the consequent amount of scattering attenuating the hot spot. Results for the surface parameter retrieval exhibits very similar behaviour for the other experiments and will not be shown.

### 6.2.2 Experiment F10

Let us now examine the case in which both rm and ni differ from those of the vertices used for the inversion. The aerosol type F1 is used for the forward simulation with rmf=0.1µm for the predominant fine mode and rmc=0.93µm for the coarse mode. The same aerosol vertices as in experiments F00 are used for the inversion.

The results in Fig. 6 show that ω0 is reasonably well retrieved unlike the g parameter, which is systematically underestimated. At any given wavelengths, it is not possible to retrieve g values outside the bounds defined by the FA and FN vertices. Consequently, the retrieved AOT values are underestimated by about 10 % (Table 6). This example illustrates the retrieval failure when the actual solution lays outside the {g,ω0} space defined by the active vertices.

Figure 7Same as Fig. 5 but for experiment F11.

Figure 8Same as Fig. 5 but for experiment F12.

Figure 9Same as Fig. 5 but for experiment F13.

### 6.2.3 Experiments F11–F13

In order to improve the retrieval of the F1 aerosol model properties, the additional aerosol CS vertex with rm=0.3µm has been added for the inversion process. Results of experiment F11 are displayed in Fig. 7. Retrieved g values are no longer underestimated. The single-scattering albedo is slightly underestimated. It should be noted that the estimated uncertainty associated with g increases with wavelength and is particularly large at 0.87 µm, but rather underestimated at 0.44 µm. The improvement in the AOT retrieval accuracy is noticeable in the 0.44 and 0.55 µm bands where the magnitude of ϵr is reduced from 0.062 to 0.005 and from 0.042 to −0.021 respectively (Table 6). At larger wavelengths, the benefit of adding the CS vertex is less noticeable though the magnitude of ϵr remains below 0.05. Finally, the retrieval uncertainty slightly increases from 0.199 up to 0.239 for the 0.44 µm band because of the use of additional state variables τv associated with the inclusion of an additional vertex. Similar behaviour is observed in the other bands.

Figure 10Same as Fig. 5 but for experiment F21.

Figure 11Same as Fig. 5 but for experiment F22.

Figure 12Same as Fig. 5 but for experiment F23.

For experiment F12, the CS vertex is substituted by vertex CL which has a median radius of 1.0 µm. The use of this vertex instead of CS considerably improves the retrieval of g and of ω0 at large wavelengths (Fig. 8). As can be seen in Fig. 2, the sensitivity of aerosol single-scattering properties to particle median radius and imaginary part of the refractive index depends on the wavelength. Hence, a similar performance of the algorithm at all wavelengths should not be expected. The errors ϵτ in this experiment F12 are further reduced compared to experiment F11 with the exception of the 0.44 µm band. The CISAR algorithm retrieves total AODs consistent with the truth.

Finally, in experiment F13 the inversion was performed using all four vertices (Fig. 9). This additional “degree of freedom” translates into an increase of the estimated uncertainty ${\mathit{\sigma }}_{\stackrel{\mathrm{^}}{\mathit{\tau }}}$ as a result of the large number of possible way to combine these four vertices to retrieve the properties of the aerosol model F1. In other words, using these two coarse mode vertices does not improve the characterisation of F1. The actual benefit of adding this fourth vertex, i.e. expanding the solution space, is therefore not straightforward, and it should be noted that increasing the number of vertices impacts the computational time. This series of experiments has shown that, of the four considered, the use of the FN, FA and CL vertices provides the best combination for the retrieval of aerosol model F1. With this combination, the FN and FA vertices allow for control of the amount of radiation absorbed by the aerosols, while the effects of the particle size are controlled by the CL vertex.

### 6.2.4 Experiments F21–F23

The retrieval of aerosol model F2, a dual-mode particle size distribution dominated by coarse particles, is now examined. This model is composed of a fine mode with radius rmf=0.08µm and a coarse mode with radius rmc= 0.77 µm. As for the retrieval of the F1 aerosol model, three combinations of vertices have been explored, i.e. FN, FA and CS for experiment F21 (Fig. 10), FN, FA and CL for experiment F23 (Fig. 11) and FN, FA, CS and CL for experiment F22 (Fig. 11). Essentially the same conclusions hold as for the retrieval of aerosol model F1. The retrieval of the F2 model properties expressed as a linear combination of the FN, FA and CL vertices provides the best solution, with both g and ω0 being retrieved at all wavelengths.

7 Discussion and conclusion

This paper describes the CISAR algorithm designed for the joint retrieval of surface reflectance and aerosol properties. Previous attempts to perform a joint retrieval have been reviewed, discussing their advantages and weaknesses. The limitations due to a combined used of discrete and continuous state variables in retrieval methods based on OE as in are discussed in Sect. 2. The new method presented in this paper specifically addresses these limitations, allowing continuous variations of the aerosol single-scattering properties in the solution space without the aerosol microphysical properties explicitly appearing as state variables.

A fast-forward radiative transfer model has been designed which solves the radiative transfer equation without relying on precomputed look-up tables. This model considers two atmospheric layers. The upper layer only hosts molecular absorption. The lower layer accounts for both absorption and scattering processes due to aerosols and molecules and is radiatively coupled with the surface represented with the RPV BRF model. Single-scattering aerosol properties in this layer are expressed as a linear combination of the properties of vertices enclosing part of the solution space.

A series of different experiments has been devised to analyse the behaviour of the CISAR algorithm and its capability to retrieve aerosol single-scattering properties as well as the optical thickness. This discussion focuses on the retrieval of aerosol models dominated by the fine mode or coarse mode. These two models have pretty different spectral behaviour in the {g,ω0} space and yet the CISAR algorithm is capable of retrieving the corresponding single-scattering properties in both cases with estimated uncertainties of about 15 %.

These experiments illustrate the possibility of using Eqs. (8) and (9) for the continuous retrieval of the aerosol single-scattering albedo and phase function. These equations assume linear behaviour of ω0 and g in the solution space. Such assumptions have proven to be valid for the case addressed in experiment F00. This assumption is not exactly true for the retrieval aerosol models of a fine- and a coarse-particle size modes. However, the retrieved aerosol single-scattering properties are derived much more accurately than with a method based on a limited number of predefined aerosol models as in , where the single-scattering properties of only predefined models are retrieved. It thus represents a major improvement with respect to these types of retrieval approaches without requiring the use of a large number of state variables as in the method proposed by , where aerosol microphysical properties are explicitly included in the set of retrieved state variables.

The choice of vertices outlining the {g,ω0} solution space is critical. In these experiments, the best retrieval is obtained using three vertices, i.e. one vertex composed of small weakly absorbing particles (FN), one vertex composed of small absorbing particles (FA) and one vertex composed of large particles (CL). The use of a fourth vertex (CS) does not improve the retrieval and increases the estimated retrieval uncertainty.

This set of experiments represents ideal conditions, i.e. noise-free observations in the principal plane with no bias on the surface prior. This choice is motivated by the need to keep the result interpretation simple to illustrate how the new retrieval concept developed in this paper works. These experiments show the possibility of retrieving aerosol single-scattering properties within the solution space provided it is correctly bounded by the vertices. It is clear that adding noise to the observations will degrade the quality of the retrieval. Similar conclusions can hold if the observations are taking place far from the principal plane, where most of the angular variations occur. Part 2 addresses the performance of CISAR when applied to actual satellite data.

Data availability
Data availability.

Results presented in Sect. 6 are available from the authors upon request.

Supplement
Supplement.

It includes the plots of case F22, adding a 3 % Gaussian noise to the simulated TOA BRF for AOT = 0.05 (Fig. S1), AOT = 0.2 (Fig. S2), AOT = 0.4 (Fig. S3) and AOT = 0.8 (Fig. S4). Figure S5 shows the merged results in the {g,ω0} space of experiments F11, F12 and F13. Figure S6 shows the merged results in the {g,ω0} space of experiments F21, F22 and F23. The supplement related to this article is available online at: https://doi.org/10.5194/amt-11-6589-2018-supplement.

Author contributions
Author contributions.

YG developed the FASTRE model, conceived the experimental set-up and wrote the paper. ML contributed to the development of the inversion method.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

The authors would like to thank the reviewers for their fruitful suggestions.

Edited by: Andrew Sayer
Reviewed by: three anonymous referees

References

Cox, C. and Munk, W.: Measurement of the Roughness of the Sea Surface from Photographs of the Sun's Glitter, J. Opt. Soc. Am., 44, 838–850, https://doi.org/10.1364/JOSA.44.000838, 1954. a

Diner, D. J., Hodos, R. A., Davis, A. B., Garay, M. J., Martonchik, J. V., Sanghavi, S. V., von Allmen, P., Kokhanovsky, A. A., and Zhai, P.: An optimization approach for aerosol retrievals using simulated MISR radiances, Atmos. Res., 116, 1–14, https://doi.org/10.1016/j.atmosres.2011.05.020, 2012. a, b

Dubovik, O., Sinyuk, A., Lapyonok, T., Holben, B. N., Mishchenko, M., Yang, P., Eck, T. F., Volten, H., Munoz, O., Veihelmann, B., van der Zande, W. J., Leon, J. F., Sorokin, M., and Slutsker, I.: Application of spheroid models to account for aerosol particle nonsphericity in remote sensing of desert dust, J. Geophys. Res.-Atmos., 111, 11208–11208, 2006. a, b, c, d, e, f

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

Fischer, J. and Grassl, H.: Radiative transfer in an atmosphere-ocean system: an azimuthally dependent matrix-operator approach, Appl. Optics, 23, 1032–1039, 1984. a

Govaerts, Y. and Lattanzio, A.: Retrieval Error Estimation of Surface Albedo Derived from Geostationary Large Band Satellite Observations: Application to Meteosat-2 and -7 Data, J. Geophys. Res., 112, D05102, https://doi.org/10.1029/2006JD007313, 2007. a, b

Govaerts, Y., Luffarelli, M., and Damman, A.: Effects of Sky Radiation on Surface Reflectance: Implications on The Derivation of LER from BRF for the Processing of Sentinel-4 Observations, in: Living Planet Symposium 2016, 9–13 May, Prague, Czech Republic, 2016. a

Govaerts, Y. M.: RTMOM V0B.10 User's Manual, 2006. a

Govaerts, Y. M., Lattanzio, A., Taberner, M., and Pinty, B.: Generating global surface albedo products from multiple geostationary satellites, Remote Sens. Environ., 112, 2804–2816, https://doi.org/10.1016/j.rse.2008.01.012, 2008. a, b

Govaerts, Y. M., Wagner, S., and Dubovik, O.: Enhanced retrieval of loading and detailed micro-physics of atmospheric aerosol from MTG/FCI observations, in: EUMETSAT Meteorological User Conference, 20–24 September, Córdoba, Spain, 2010a. a

Govaerts, Y. M., Wagner, S., Lattanzio, A., and Watts, P.: Joint retrieval of surface reflectance and aerosol optical depth from MSG/SEVIRI observations with an optimal estimation approach: 1. Theory, J. Geophys. Res., 115, 10.1029/2009JD011779, https://doi.org/10.1029/2009JD011779, 2010b. a, b, c, d, e, f, g, h, i, j, k, l, m, n

Hess, M., Koepke, P., and Schult, I.: Optical properties of aerosols and clouds: The software package OPAC, B. Am. Meteorol. Soc., 79, 831–844, 1998. a

Kokhanovsky, A. A., Deuzé, J. L., Diner, D. J., Dubovik, O., Ducos, F., Emde, C., Garay, M. J., Grainger, R. G., Heckel, A., Herman, M., Katsev, I. L., Keller, J., Levy, R., North, P. R. J., Prikhach, A. S., Rozanov, V. V., Sayer, A. M., Ota, Y., Tanré, D., Thomas, G. E., and Zege, E. P.: The inter-comparison of major satellite aerosol retrieval algorithms using simulated intensity and polarization characteristics of reflected light, Atmos. Meas. Tech., 3, 909–932, https://doi.org/10.5194/amt-3-909-2010, 2010. a

Lattanzio, A., Schulz, J., Matthews, J., Okuyama, A., Theodore, B., Bates, J. J., Knapp, K. R., Kosaka, Y., and Schüller, L.: Land Surface Albedo from Geostationary Satelites: A Multiagency Collaboration within SCOPE-CM, B. Am. Meteorol. Soc., 94, 205–214, https://doi.org/10.1175/BAMS-D-11-00230.1, 2013. a, b

Liu, Q. and Ruprecht, E.: Radiative transfer model: matrix operator method, Appl. Optics, 35, 4229–4237, 1996. a

Luffarelli, M. and Govaerts, Y.: Joint retrieval of surface reflectance and aerosol properties with continuous variation of the state variables in the solution space: Part 2: Application to geostationary and polar orbiting satellite observations, Atmos. Meas. Tech. Discuss., https://doi.org/10.5194/amt-2018-265, in review, 2018.  a

Luffarelli, M., Govaerts, Y., and Damman, A.: Assessing hourly aerosol properties retrieval from MSG/SEVIRI observations in the framework of aeroosl-cci2, in: Living Planet Symposium 2016, 9–13 May, Prague, Czech Republic, 2016. a, b

Luffarelli, M., Govaerts, Y., Goossens, C., Wolters, E., and Swinnen, E.: Joint retrieval of surface reflectance and aerosol properties from PROBA-V observations, part I: algorithm performance evaluation, in: Proceedings of MultiTemp 2017, 27–29 June, Bruges, Belgium, 2017. a

Manolis, I., Grabarnik, S., Caron, J., Bézy, J.-L., Loiselet, M., Betto, M., Barré, H., Mason, G., and Meynart, R.: The MetOp second generation 3MI instrument, Proc. SPIE Remote Sens., 88890J, https://doi.org/10.1117/12.2028662, 2013. a

Marquardt, D.: An Algorithm for Least-Squares Estimation of Nonlinear Parameters, SIAM J. Appl. Math., 11, 431–441, 1963. a

Pinty, B., Roveda, F., Verstraete, M. M., Gobron, N., Govaerts, Y., Martonchik, J. V., Diner, D. J., and Kahn, R. A.: Surface albedo retrieval from Meteosat: Part 1: Theory, J. Geophys. Res., 105, 18099–18112, 2000a. a, b

Pinty, B., Roveda, F., Verstraete, M. M., Gobron, N., Govaerts, Y., Martonchik, J. V., Diner, D. J., and Kahn, R. A.: Surface albedo retrieval from Meteosat: Part 2: Applications, J. Geophys. Res., 105, 18113–18134, 2000b. a

Rahman, H., Pinty, B., and Verstraete, M. M.: Coupled surface-atmosphere reflectance (CSAR) model, 2. Semiempirical surface model usable with NOAA Advanced Very High Resolution Radiometer Data, J. Geophys. Res., 98, 20791–20801, 1993. a

Rodgers, C. D.: Inverse methods for atmospheric sounding, Series on Atmospheric Oceanic and Planetary Physics, World Scientific, Oxford, 2000. a

Schuster, G. L., Dubovik, O., Holben, B. N., and Clothiaux, E. E.: Inferring black carbon content and specific absorption from Aerosol Robotic Network (AERONET) aerosol retrievals, J. Geophys. Res., 110, 1017–1017, 2005. a

Serene, F. and Corcoral, N.: PARASOL and CALIPSO: Experience Feedback on Operations of Micro and Small Satellites, in: SpaceOps 2006 Conference, 21 June, American Institute of Aeronautics and Astronautics, https://doi.org/10.2514/6.2006-5919, 2006. a

Vermote, E. F., Tanré, D., Deuzé, J. L., Herman, M., and Morcrette, J. J.: Second simulation of the satellite signal in the solar spectrum, 6S: An overview, IEEE T. Geosci. Remote, 35, 675–686, 1997. a

Wagner, S. C., Govaerts, Y. M., and Lattanzio, A.: Joint retrieval of surface reflectance and aerosol optical depth from MSG/SEVIRI observations with an optimal estimation approach: 2. Implementation and evaluation, J. Geophys. Res., 115, D02204, https://doi.org/10.1029/2009JD011780, 2010. a, b

Wiscombe, W. J.: The Delta-M Method: Rapid Yet Accurate Radiative Flux Calculations for Strongly Asymmetric Phase Functions, J. Atmos. Sci., 34, 1408–1422, 1977. a