Aerosol retrieval over snow using the RemoTAP algorithm

In order to conduct accurate aerosol retrieval over snow, the Remote Sensing of Trace Gases and Aerosol Products (RemoTAP) algorithm developed by SRON Netherlands Institute for Space Research is extended with a bi-directional reflection distribution function (BRDF) for snow surfaces. The capability of the extended algorithm is validated with both synthetic measurements and real satellite measurements from the Polarization & Anisotropy of Reflectances for Atmospheric Sciences coupled with Observations from a Lidar (PARASOL), and a comparison has been made to retrievals with the baseline RemoTAP (without a snow kernel). For retrievals with real PARASOL observations, we use pixels over Aerosol Robotic Network (AERONET) stations for validation and we use the MODIS snow cover products to identify pixels over snow. We evaluate the retrieved aerosol optical thickness (AOT) at 550 nm (τ₅₅₀), single-scattering albedo (SSA) at 550 nm (ω₅₅₀) and Ångström exponent (AE) for 440–870 nm (AE_440−870). The experiments with both synthetic and real data show that the extended RemoTAP maintains capability for snow-free pixels and has obvious advantages in accuracy and the fraction of successful retrievals for retrieval over snow, especially over surfaces with snow cover >75 %. According to the real-data experiment, we find that the retrieval algorithm has difficulty in fitting the PARASOL 1020 nm band, where snow reflectance is significantly lower than that for the visible bands. When we perform a four-band retrieval (490, 565, 670, 865 nm) with the extended RemoTAP, we obtain a good retrieval result for τ₅₅₀, ω₅₅₀ and AE_440−870. Therefore, the four-band retrieval with the extended RemoTAP is recommended for aerosol retrieval over snow.

1 Introduction

Global climate change is greatly influenced by aerosol through aerosol–cloud interaction (e.g., Twomey, 1974; Li et al., 2011; Rosenfeld et al., 2014; Hasekamp et al., 2019a; Gryspeerdt et al., 2020; Quaas et al., 2020) and aerosol–radiation interaction (e.g., Koren et al., 2004; Yu et al., 2006; Myhre, 2009; Guo et al., 2017; Lacagnina et al., 2017; Witthuhn et al., 2021). These effects cause a significant radiative forcing of climate (Bellouin et al., 2020; Haywood et al., 2021). According to the Sixth Assessment Report (AR6) of the Intergovernmental Panel on Climate Change (IPCC), aerosols still represent the largest uncertainty in our quantification of global climate change (IPCC, 2023). Therefore, in order to better study climate change, it is essential to give an accurate estimate of global aerosol properties. Multi-angle polarimetric (MAP) satellite measurements provide the richest set of information on aerosol properties from a passive remote sensing point of view. So far, the only MAP instrument that has provided a multi-year data set has been POLDER-3 (Polarization and Directionality of the Earth's Reflectances 3) on board PARASOL (Polarization & Anisotropy of Reflectances for Atmospheric Sciences coupled with Observations from a Lidar; the whole thing hereafter referred to as PARASOL), which was active from 2005–2013. In the near future a number of new MAP instruments will be launched, including SPEXone (Hasekamp et al., 2019b) and HARP-2 (McBride et al., 2020) on the NASA PACE mission (Werdell et al., 2019) and the 3MI instrument (Fougnie et al., 2018) on MetOp-SG.

After years of development, there exist a number of aerosol retrieval algorithms which are available for multi-angle and multi-spectral polarization sensors. The algorithms can be classified into two groups: algorithms based on lookup tables (LUTs) and full inversion algorithms. LUT-based algorithms (Herman et al., 1997; Deuzé et al., 2000; Waquet et al., 2016) provide faster calculation but are less accurate than full inversion algorithms. Among the full inversion algorithms are the microphysical aerosol properties from polarimetry (MAPP) algorithm (Stamnes et al., 2018), the algorithm developed at the Jet Propulsion Laboratory (JPL) (Xu et al., 2016, 2017, 2018), the Generalized Retrieval of Aerosol and Surface Properties (GRASP) algorithm (Dubovik et al., 2011, 2014, 2021; Chen et al., 2020), and the Remote Sensing of Trace Gases and Aerosol Products (RemoTAP) algorithm developed by SRON Netherlands Institute for Space Research (Hasekamp et al., 2011; Fu and Hasekamp, 2018; Fu et al., 2020; Lu et al., 2022). Till now, only the RemoTAP algorithm and the GRASP algorithm have demonstrated capability at the global scale.

In order to generate global aerosol products, aerosol retrieval over snow remains an important challenge for the abovementioned algorithms. Until now, studies on aerosol retrieval over snow have focused mainly on intensity sensors (Istomina et al., 2011; Mei et al., 2013) and on aerosol optical thickness (AOT) retrievals in the Arctic region.

In our paper, we extend the RemoTAP algorithm to carry out aerosol retrieval over snow from MAP measurements of PARASOL. We evaluate the capability of the extended RemoTAP algorithm using synthetic observations as well as real PARASOL retrievals which are validated by the Aerosol Robotic Network (AERONET). The paper is organized as follows: Sect. 2 introduces the methodology of RemoTAP for aerosol retrieval over snow; Sect. 3 describes the satellite data and ancillary data used in our real-data retrieval and the validation data; Sect. 4 shows results for synthetic retrievals; Sect. 5 shows the results for real PARASOL retrievals and provides a recommended routine for aerosol retrieval over snow from PARASOL; Sect. 6 discusses the results of the paper and future research and concludes the paper.

2 Methodology

2.1 Forward model

For aerosol retrieval over snow, detailed information for the extended RemoTAP algorithm is described below.

The aim of the RemoTAP algorithm is to retrieve a state vector x from measurement vector y by inverting the forward model F:

$$\begin{array}{}\text{(1)}& \mathbf{y}=\mathbf{F}\left(\mathbf{x}\right)+{\mathbf{e}}_{\mathbf{y}},\end{array}$$

where e_y refers to the error vector including measurement error and modeling error.

The measurement by polarization sensors is described by the intensity vector $\left[I,Q,U,V\right]$ (Stokes, 1851). For Earth observation, the V parameter can be ignored in most cases. In this study, the Stokes parameters are calculated by the SRON radiative transfer model LINTRAN v2.0 (Schepers et al., 2014). The measurement vector y contains parameters including the top-of-atmosphere (ToA) apparent reflectance ρ and degree of linear polarization (DoLP) at different spectral bands and different observation geometries:

$$\begin{array}{}\text{(2)}& \mathit{\rho }={\displaystyle \frac{I}{E_{\mathrm{0}}}},\end{array}$$

$$\begin{array}{}\text{(3)}& \mathrm{DoLP}={\displaystyle \frac{\sqrt{Q^{\mathrm{2}}+U^{\mathrm{2}}}}{I}},\end{array}$$

where E₀ is the ToA solar irradiance perpendicular to the solar beam.

The state vector x contains parameters related to aerosol and surface characteristics. To describe aerosol properties in the state vector, following Lu et al. (2022), three log-normal modes are applied: (1) mode 1 is a fine mode for which the state vector includes aerosol column number (N_aer), effective radius (r_eff), effective variance (v_eff), spherical fraction (f_sph), aerosol layer height (z_aer here refers to the altitude of the aerosol layer center) and the refractive index coefficients corresponding to the standard spectra (D'Almeida et al., 1991; Kirchstetter et al., 2004) of inorganic aerosol, black carbon (imaginary part) and organic carbon (imaginary part). (2) Mode 2 is an insoluble coarse mode, which consists of non-spherical dust. The state vector includes N_aer, r_eff, v_eff, z_aer and refractive index coefficients corresponding to the standard spectra of dust (imaginary part). The fixed parameter is f_sph=0 and the coefficient for the real part of the dust refractive index is fixed to 1. z_aer is assumed to be the same for mode 1 and mode 2. (3) Mode 3 is a soluble coarse mode. The state vector includes N_aer, r_eff and v_eff, and the refractive index coefficient corresponds to the standard spectra of inorganic aerosol. The fixed parameters are f_sph=1 and z_aer=500 m.

Table 1 Parameters in the state vector x utilized in the retrieval. c₁, c₂ and c₃ correspond to coefficients for standard refractive index of inorganic aerosol (INORG), black carbon (BC), organic carbon (OC) and dust (DU). BPDF denotes the bi-directional polarization distribution function.

Table 2 Spectral band information for PARASOL.

To describe the surface, the surface reflection matrix ${\mathbf{R}}_{\mathrm{surf}}\left(\mathit{\lambda },{\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)$ is given by

$$\begin{array}{}\text{(4)}& {\mathbf{R}}_{\mathrm{surf}}\left(\mathit{\lambda },{\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)=r_{\mathrm{11}}\left(\mathit{\lambda },{\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)\mathbf{D}+{\mathbf{R}}_{\mathrm{pol}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right),\end{array}$$

where λ is the wavelength; θ_s and θ_v are the sun zenith angle and view zenith angle, respectively; φ is the relative azimuth angle; D is a null matrix except D₁₁=1, $r_{\mathrm{11}}\left(\mathit{\lambda },{\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)$ is described by the Ross–Li bi-directional reflection distribution function (BRDF) model, extended by a snow kernel:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{r}_{\mathrm{11}}& \left(\mathit{\lambda },{\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)=A\left(\mathit{\lambda }\right)\left[\mathrm{1}+k_{\mathrm{geo}}{f}_{\mathrm{geo}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)\right.\\ & \left.+k_{\mathrm{vol}}{f}_{\mathrm{vol}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)+k_{\mathrm{snow}}{f}_{\mathrm{snow}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)\right],\end{array}\end{array}$$

where A(λ) is the isotropic reflectance; $f_{\mathrm{geo}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)$ and $f_{\mathrm{vol}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)$ are the geometric (Li-Sparse) kernel and volumetric (Ross-Thick) kernel function, respectively (Wanner et al., 1995); $f_{\mathrm{snow}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)$ is the snow kernel function (Jiao et al., 2019); and k_geo, k_vol and k_snow are the coefficients for the Li-Sparse, Ross-Thick and snow kernel, respectively. It is important to note that the inclusion of $k_{\mathrm{snow}}{f}_{\mathrm{snow}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)$ is the major difference between the baseline RemoTAP and the extended RemoTAP of the present work. The Li-Sparse kernel function $f_{\mathrm{geo}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)$, Ross-Thick kernel function $f_{\mathrm{vol}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)$ and snow kernel function $f_{\mathrm{snow}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)$ are given in Eqs. (6) and (7):

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{f}_{\mathrm{geo}}& \left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)=O\left({\mathit{\theta }}_{\mathrm{s}}{}^{\prime },{\mathit{\theta }}_{\mathrm{v}}{}^{\prime },\mathit{\phi }\right)-\left[\sec {\mathit{\theta }}_{\mathrm{s}}{}^{\prime }+\sec {\mathit{\theta }}_{\mathrm{v}}{}^{\prime }\right.\\ & \left.-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(\mathrm{1}-\cos \mathrm{\Theta }{}^{\prime })\sec {\mathit{\theta }}_{\mathrm{s}}{}^{\prime }\sec {\mathit{\theta }}_{\mathrm{v}}{}^{\prime }\right],\end{array}\end{array}$$

$$\begin{array}{}\text{(7)}& f_{\mathrm{vol}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)={\displaystyle \frac{\left(\mathrm{\Theta }-\mathit{\pi }/\phantom{\mathit{\pi }\mathrm{2}}\mathrm{2}\right)\cos \left(\mathit{\pi }-\mathrm{\Theta }\right)+\sin \left(\mathit{\pi }-\mathrm{\Theta }\right)}{\cos {\mathit{\theta }}_{\mathrm{s}}+\cos {\mathit{\theta }}_{\mathrm{v}}}}-{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}},\end{array}$$

where Θ is the scattering angle ($\cos \mathrm{\Theta }=-\cos {\mathit{\theta }}_{\mathrm{s}}\cos {\mathit{\theta }}_{\mathrm{v}}-\sin {\mathit{\theta }}_{\mathrm{s}}\sin {\mathit{\theta }}_{\mathrm{v}}\cos \mathit{\phi })$, θ_s^′ and θ_v^′ are the equivalent zenith angles (the transformation is $\mathit{\theta }{}^{\prime }={\tan }^{-\mathrm{1}}\left(b/r\cdot \tan \mathit{\theta }\right)$ for the sun zenith angle and view zenith angle, respectively, where b and r are vertical and horizontal crown radius, respectively), and $O\left({\mathit{\theta }}_{\mathrm{s}}{}^{\prime },{\mathit{\theta }}_{\mathrm{v}}{}^{\prime },\mathit{\phi }\right)$ is the overlap function given by Li and Strahler (1992).

The snow kernel function $f_{\mathrm{snow}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)$ is given by Jiao et al. (2019):

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& f_{\mathrm{snow}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)=R_{\mathrm{0}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)\left[\mathrm{1}-\mathit{\alpha }\right.\\ & \left.\cdot \cos \left(\mathit{\pi }-\mathrm{\Theta }\right)\cdot {\mathrm{e}}^{-\cos \left(\mathit{\pi }-\mathrm{\Theta }\right)}\right]+\mathrm{0.4076}\cdot \mathit{\alpha }-\mathrm{1.1081},\end{array}\end{array}$$

where $R_{\mathrm{0}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)$ is the reflectance for a semi-infinite, non-absorbing snow layer at zero absorption and α is an empirical parameter used to correct $R_{\mathrm{0}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)$ for the underestimation of the reflectance in the forward-scattering direction. $R_{\mathrm{0}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)$ is given by Kokhanovsky et al. (2005):

$$\begin{array}{}\text{(9)}& R_{\mathrm{0}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)={\displaystyle \frac{\begin{array}{c}{K}_{\mathrm{1}}+K_{\mathrm{2}}\cdot \left(\cos {\mathit{\theta }}_{\mathrm{s}}+\cos {\mathit{\theta }}_{\mathrm{v}}\right)\\ +K_{\mathrm{3}}\cdot \cos {\mathit{\theta }}_{\mathrm{s}}\cdot \cos {\mathit{\theta }}_{\mathrm{v}}+P\left(\mathrm{\Theta }\right)\end{array}}{\mathrm{4}\left(\cos {\mathit{\theta }}_{\mathrm{s}}+\cos {\mathit{\theta }}_{\mathrm{v}}\right)}},\end{array}$$

$$\begin{array}{}\text{(10)}& P\left(\mathrm{\Theta }\right)=\mathrm{11.1}\cdot {\mathrm{e}}^{-\mathrm{0.087}\cdot \mathrm{\Theta }}+\mathrm{1.1}\cdot {\mathrm{e}}^{-\mathrm{0.014}\cdot \mathrm{\Theta }},\end{array}$$

where K₁, K₂ and K₃ are three constants. In our algorithm, α=0.3, K₁=1.247, K₂=1.186 and K₃=5.157 are fixed when calculating the snow kernel, as suggested by Jiao et al. (2019) and Kokhanovsky et al. (2005).

${\mathbf{R}}_{\mathrm{pol}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)$ in Eq. (4) is given by

$$\begin{array}{}\text{(11)}& {\mathbf{R}}_{\mathrm{pol}}\left({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi }\right)=B_{\mathrm{pol}}{\displaystyle \frac{\exp \left[-\tan \left(\frac{\mathit{\pi }-\mathrm{\Theta }}{\mathrm{2}}\right)\right]\exp (-v){\mathbf{F}}_p(m,\mathrm{\Theta })}{\mathrm{4}(\cos {\mathit{\theta }}_{\mathrm{s}}+\cos {\mathit{\theta }}_{\mathrm{v}})}},\end{array}$$

where B_pol is the free linear parameter, F_p(mΘ) is the Fresnel scattering matrix and m is the refractive index (Maignan et al., 2009). In our experiment, m=1.5 and v=0.1 are fixed when calculating R_pol.

To characterize the surface properties, we include A(λ), k_geo, k_vol, k_snow and B_pol as the surface parameters in the state vector. Theoretically, the surface model still maintains the ability to depict snow-free surfaces because k_snow is fit in the retrieval (for snow-free surfaces, it should retrieve k_snow=0). In our algorithm, A(λ) is fit separately for each wavelength, which provides full flexibility to represent any spectral shape of the snow albedo. An alternative method to deal with the spectral dependence is discussed by Kokhanovsky et al. (2023), where A(λ) is parameterized with the effective absorption length L, which is valid for snow with a different microstructure and pollution level.

The aerosol and surface parameters in the state vector x are shown in Table 1.

2.2 Inversion algorithm

To retrieve the state vector x from the measurement y, a damped Gauss–Newton iteration method with Phillips–Tikhonov regularization is employed (Hasekamp et al., 2011; Fu and Hasekamp, 2018). We shortly summarize the method here. The aim of the inversion algorithm is to find the solution $\widehat{\mathbf{x}}$ of the minimization–optimization problem:

$$\begin{array}{}\text{(12)}& \widehat{\mathbf{x}}=\left[{∥{\mathbf{S}}_y^{-\frac{\mathrm{1}}{\mathrm{2}}}\left(\mathbf{F}\left(\mathbf{x}\right)-\mathbf{y}\right)∥}^{\mathrm{2}}+{\mathit{\gamma }}^{\mathrm{2}}{∥{\mathbf{W}}^{-\frac{\mathrm{1}}{\mathrm{2}}}\left(\mathbf{x}-{\mathbf{x}}_a\right)∥}^{\mathrm{2}}\right],\end{array}$$

where x_a is the a priori state vector, W is the diagonal weight matrix in order to remove the order-of-magnitude difference of each state parameter ($W_{ii}=\mathrm{1}/x_{a,i})$, S_y is the diagonal measurement error covariance matrix which is related to the sensor and γ is the regularization parameter (Hasekamp et al., 2011).

Figure 1 Reference reflectance for vegetation, soil and snow. The snow data were downloaded from the National Snow and Ice Data Center (NSIDC) (https://nsidc.org/data/hma_sbrf/versions/1, last access: 8 June 2023), and the soil and vegetation data were downloaded from ASTER spectral library (https://speclib.jpl.nasa.gov/, last access: 8 June 2023).

Since the forward model F is nonlinear, an iterative strategy is utilized to conduct the inversion. For each iteration n, the linear approximation of the forward model is given in Eq. (13):

$$\begin{array}{}\text{(13)}& \mathbf{F}\left(\mathbf{x}\right)\approx \mathbf{F}\left({\mathbf{x}}_n\right)+\mathbf{K}(\mathbf{x}-{\mathbf{x}}_n),\end{array}$$

where K is the Jacobian matrix containing the derivatives of the forward model F with respect to the parameters in state vector x and $K_{ij}=\frac{\partial F_i}{\partial x_j}\left({\mathbf{x}}_n\right)$.

With linear approximation, Eq. (12) can be simplified:

$$\begin{array}{}\text{(14)}& {\stackrel{\mathrm{̃}}{\mathbf{x}}}_{n+\mathrm{1}}=\left[{∥\stackrel{\mathrm{̃}}{\mathbf{K}}\left(\stackrel{\mathrm{̃}}{\mathbf{x}}-{\stackrel{\mathrm{̃}}{\mathbf{x}}}_n\right)-\stackrel{\mathrm{̃}}{\mathbf{y}}∥}^{\mathrm{2}}+{\mathit{\gamma }}^{\mathrm{2}}{∥\stackrel{\mathrm{̃}}{\mathbf{x}}-{\stackrel{\mathrm{̃}}{\mathbf{x}}}_a∥}^{\mathrm{2}}\right],\end{array}$$

where $\stackrel{\mathrm{̃}}{\mathbf{x}}={\mathbf{W}}^{-\frac{\mathrm{1}}{\mathrm{2}}}\mathbf{x}$, $\stackrel{\mathrm{̃}}{\mathbf{y}}={\mathbf{S}}_y^{-\frac{\mathrm{1}}{\mathrm{2}}}\left(\mathbf{y}-\mathbf{F}\left({\mathbf{x}}_n\right)\right)$ and $\stackrel{\mathrm{̃}}{\mathbf{K}}={\mathbf{S}}_y^{-\frac{\mathrm{1}}{\mathrm{2}}}{\mathbf{KW}}^{\frac{\mathrm{1}}{\mathrm{2}}}$.

The solution to Eq. (14) for iteration step n is given by

$$\begin{array}{}\text{(15)}& {\stackrel{\mathrm{̃}}{\mathbf{x}}}_{n+\mathrm{1}}={\stackrel{\mathrm{̃}}{\mathbf{x}}}_n+\mathrm{\Lambda }{\left({\stackrel{\mathrm{̃}}{\mathbf{K}}}^T\stackrel{\mathrm{̃}}{\mathbf{K}}+{\mathit{\gamma }}^{\mathrm{2}}\mathbf{I}\right)}^{-\mathrm{1}}\left[{\stackrel{\mathrm{̃}}{\mathbf{K}}}^T\stackrel{\mathrm{̃}}{\mathbf{y}}-{\mathit{\gamma }}^{\mathrm{2}}\left({\stackrel{\mathrm{̃}}{\mathbf{x}}}_n-{\stackrel{\mathrm{̃}}{\mathbf{x}}}_a\right)\right],\end{array}$$

where I is the identity matrix and Λ is the filter factor ($\mathrm{0}\le \mathrm{\Lambda }\le \mathrm{1})$, which is utilized to control the step size per iteration. For each iteration, the optimal Λ and γ are chosen via goodness-of-fit assessment by comparing χ² given in Eq. (16):

$$\begin{array}{}\text{(16)}& {\mathit{\chi }}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{N_{\mathrm{meas}}}}{\sum }_{i=\mathrm{1}}^{N_{\mathrm{meas}}}{\left({\displaystyle \frac{y_i-F_i}{e_i}}\right)}^{\mathrm{2}},\end{array}$$

where N_meas is the length of measurement vector y; y_i and F_i refer to the measurements and the results of forward model, respectively; and e_i is the measurement uncertainties determined by the instrument. After the final iteration step, we use χ²<5 as the threshold to determine whether the retrieval was successful or not. For a set of retrievals, the fraction of successful retrievals (FoSR) is determined based on this filter.

Table 3 Observation geometry, aerosol properties and surface properties used to create synthetic PARASOL observations. c_veg, c_soil and c_snow are the fraction of vegetation, soil and snow, respectively. Distribution “linear” refers to $X\sim \mathrm{U}(X_{\min },X_{\max })$ and distribution “logarithmic” refers to $\ln X\sim \mathrm{U}(\ln X_{\min },\ln X_{\max })$, where X is the property value and X_min and X_max are the minimum and maximum, respectively.

Table 4 Descriptions for different synthetic measurements.

3 Data

3.1 PARASOL data

The microsatellite Polarization & Anisotropy of Reflectances for Atmospheric Sciences coupled with Observations from a Lidar (PARASOL), equipped with the Polarization and Directionality of the Earth's Reflectances 3 (POLDER-3) instrument, was launched on 18 December 2004 (Fougnie et al., 2007; Lier and Bach, 2008). PARASOL includes nine spectral bands, and three bands contain polarization information (the details of the bands are shown in Table 2).

PARASOL provides multi-directional and multi-spectral data from December 2004 to December 2013, which have been used for aerosol retrieval (e.g., Tanré et al., 2011; Dubovik et al., 2011; Fu and Hasekamp, 2018; Chen et al., 2020), and the resolution is approximately 6 km in the nadir. In our experiments, we consider five bands (490, 565, 670, 865, 1020 nm) because the 443 nm band suffers from stray light which may be in particular problematic over bright snow surfaces (Fougnie et al., 2007), 763 and 765 nm are mainly used to retrieve cloud oxygen pressure, and 910 nm is usually used to retrieve water vapor (Leroy et al., 1997).

3.2 Ancillary data

As input to the retrieval, ancillary data are needed, including surface pressure and profiles of air temperature and relative humidity after a moist ozone mass mixing ratio. These data are obtained from the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) (Gelaro et al., 2017). Additionally, the cloud fraction is from the Moderate Resolution Imaging Spectroradiometer (MODIS). In order to conduct a reliable cloud screening, only pixels with a cloud fraction under 0.2 are used in our experiments.

Figure 2 Synthetic-data retrievals of τ₅₅₀, ω₅₅₀ and AE_440−870 among the baseline RemoTAP (red bar) and extended RemoTAP (blue bar) over different surfaces (see Table 4). Panels (a1, a2), (b1, b2), (c1, c2) and (d) show the bar plots of τ₅₅₀, ω₅₅₀, AE_440−870 and the fraction of successful retrievals (FoSR), respectively. Panels (a1, b1, c1) and (a2, b2, c2) show RMSE and bias, respectively.

3.3 Snow cover data

During the process of retrieval and validation, pixels with snow cover are selected based on the MODIS/Aqua Snow Cover Daily L3 Global 500 m SIN Grid (MYD10A1) product (Hall et al., 2019). Each tile is generated from the MODIS/Aqua Snow Cover 5-Min L2 Swath 500 m (MYD10_L2) product. Given the row and column in a MODIS tile, the 500 m MODIS sinusoidal grid is converted to a PARASOL sinusoidal grid. It is important to note that the MODIS snow cover product has been pre-processed with a cloud filtering procedure.

3.4 AERONET data

In this paper, the main concern is aerosol retrieval over snow. The RemoTAP-retrieved aerosol properties are validated with AERONET data (Holben et al., 1998). AERONET provides three levels of data quality: level 1.0, level 1.5 and level 2.0. Level 1.0 provides unscreened data which have rarely been used for validation. Level 1.5 provides near-real-time, cloud-screened data with instrument quality control. Level 2.0 provides quality-assured data on the basis of level 1.5 by applying pre-field and post-field calibrations. The retrieved AOT and Ångström exponent (AE) are validated with AERONET direct-sun level-2.0 AOT data (Giles et al., 2019). The retrieved single-scattering albedo (SSA) is validated with AERONET inversion level-2.0 SSA data (Sinyuk et al., 2020).

Figure 3 Synthetic-data retrievals of τ₅₅₀, ω₅₅₀ and AE_440−870 among retrievals of the extended RemoTAP versus synthetic truth over pure snow surfaces (c_snow=100 %). Panels (a), (b) and (c) show the scatterplot of τ₅₅₀, ω₅₅₀ and AE_440−870, respectively. The color indicates the density of data points, where yellow indicates high density and blue–purple low density (viridis color map).

3.5 Data pre-processing

In the real-data experiments, the PARASOL data are pre-processed in steps. The first step is to match global PARASOL L1 measurement data with global AERONET validation data (AOT, SSA and AE). For each successfully matched group of pixels (hereafter referred to as a colocation), the difference between the measurement time of AERONET data and PARASOL pixels is within 1 h and the distance between AERONET data and PARASOL pixels is within 20 km. The second step is to match MODIS snow cover data with AERONET-colocated PARASOL data and thus divide the colocated PARASOL data into snow pixels and snow-free pixels.

4 Synthetic-data experiments

The forward model of RemoTAP is used to generate the synthetic PARASOL measurement, and noise is subsequently added according to a Gaussian distribution. For ToA reflectance, the simulated noise (1 standard deviation) is 1 %, and for DoLP it is 0.007 (absolute). The set of synthetic measurements contains 1000 pixels with randomly generated input land properties, aerosol properties, auxiliary properties and geometry properties. The settings for the properties are shown in Table 3.

Figure 4 Real-data retrievals of τ₅₅₀ among five-band (panel a1, a2), four-band (panel b1, b2) and three-band (panel c1, c2) RemoTAP versus AERONET over snow-dominant surfaces (c_snow>75 %). Panels (a1), (b1) and (c1) show the scatterplot of the baseline RemoTAP, and panels (a2), (b2) and (c2) show the scatterplot of the extended RemoTAP. The statistics with threshold required by the Global Climate Observing System (GCOS) are also shown in the bottom right of each panel.

The surface properties in the synthetic data set are created by mixing the contribution of the surface reflection by vegetation, soil and snow. By controlling the fraction of vegetation, soil and snow, four sets of synthetic measurements are created, and the detailed information for these four synthetic measurements is listed in Table 4. The isotropic reflectance A(λ) is calculated with equation $A\left(\mathit{\lambda }\right)=c_{\mathrm{veg}}{A}_{\mathrm{veg}}\left(\mathit{\lambda }\right)+c_{\mathrm{soil}}{A}_{\mathrm{soil}}\left(\mathit{\lambda }\right)+c_{\mathrm{snow}}{A}_{\mathrm{snow}}\left(\mathit{\lambda }\right)$. A_veg(λ), A_soil(λ) and A_snow(λ) refer to the reference reflectance spectra for vegetation, soil and snow, respectively (shown in Fig. 1). For the kernel coefficients of Li-Sparse ($k_{\mathrm{geo}}=\mathrm{0.087}{c}_{\mathrm{veg}}+\mathrm{0.158}{c}_{\mathrm{soil}})$ and Ross-Thick ($k_{\mathrm{vol}}=\mathrm{0.688}{c}_{\mathrm{veg}}+\mathrm{0.547}{c}_{\mathrm{soil}})$, the constant values we use are found by Litvinov et al. (2011).

Figure 5 Real-data retrievals of τ₅₅₀, ω₅₅₀ and AE_440−870 for the baseline four-band RemoTAP (red bar), extended four-band RemoTAP (blue bar), baseline three-band RemoTAP (green bar) and extended three-band RemoTAP (brown bar) over snow surfaces ($\mathrm{0}\mathit{\%}<c_{\mathrm{snow}}\le \mathrm{100}\mathit{\%})$ and snow-free surfaces (c_snow=0 %). Panels (a1, a2), (b1, b2) and (c1, c2) show the bar plots of τ₅₅₀, ω₅₅₀ and AE_440−870, respectively. Panels (a1, b1, c1) and (a2, b2, c2) show RMSE and bias, respectively.

These four sets of synthetic measurements are taken as the input for RemoTAP to conduct the retrieval with the baseline RemoTAP setup (without snow BRDF) and the extended RemoTAP setup (with snow BRDF). During validation, τ₅₅₀, ω₅₅₀ and AE_440−870 are chosen as the main performance indicators. In our experiment, a five-band (490, 565, 670, 865, 1020 nm) retrieval is conducted and the validation results for different synthetic situations are shown in Fig. 2.

Comparing the RMSE, bias and fraction of successful retrievals between the baseline and extended RemoTAP in Fig. 2, we can conclude that the baseline RemoTAP, which utilizes the Ross–Li model to characterize the ground surfaces, has poor capability to retrieve aerosol properties over snow. For synthetic measurements of pure snow (snow_pure) and dominated by snow (snow_domi), the fraction of successful retrieval is low (34.8 % and 51.6 %, respectively) and the retrieval accuracy for the successfully retrieved pixels is low, with an RMSE of 0.268 and 0.197, respectively, for τ₅₅₀. The extended version of RemoTAP with the Ross–Li snow model on the other hand has much better performance with a high fraction of successful retrievals (89.1 % to 97.4 %) and an RMSE in τ₅₅₀ of 0.084 and 0.097, respectively. As expected, the performance over snow-covered pixels is still worse than for snow-free pixels because the high signal from the bright snow surface overwhelms the aerosol signal, leading to reduced aerosol information content. For the synthetic measurements in snow_rand, in which snow is not dominant, the baseline RemoTAP achieves acceptable performance, but the performance is still worse than the extended RemoTAP in both accuracy and the fraction of successful retrievals. For ω₅₅₀ retrieval, the performance of the baseline is slightly better than the extended RemoTAP, but it should be noted that only a few pixels are left after the χ² filtering. Overall, the results demonstrate the importance of extending RemoTAP with the snow kernel in order to achieve the capability for aerosol retrieval over snow.

Looking at the validation of the synthetic measurement snow_free in which there is no snow present, the extended RemoTAP maintains a similar capability to the baseline RemoTAP in terms of RMSE, bias and fraction of successful retrievals. For synthetic measurement with snow (snow_pure, snow_domi and snow_rand), especially when snow is dominant among the land cover (snow_pure and snow_domi), the extended RemoTAP has good performance for aerosol retrieval over snow. The scatterplot for the synthetic measurement snow_pure retrieved with the extended RemoTAP is shown in Fig. 3. The accuracy for τ₅₅₀, ω₅₅₀ and AE_440−870 retrievals is good, with an RMSE of 0.084 for τ₅₅₀ retrieval, an RMSE of 0.021 for ω₅₅₀ retrieval and an RMSE of 0.293 for AE_440−870 retrieval.

Figure 6 Fraction of successful retrievals (left y axis) for real-data retrievals among baseline four-band RemoTAP (red bar) and extended four-band RemoTAP (blue bar) over different snow cover intervals. The number of the matched snow pixels (right y axis) over different snow cover intervals is plotted with green bars.

The synthetic results suggest that PARASOL measurements have enough information to allow the inclusion of the snow kernel in the retrieval state vector. Also, they show that by inclusion of the snow kernel, the extended RemoTAP is capable of performing aerosol retrievals over snow in a consistent setup. The next section will show the performance of real measurements.

5 Real-data experiments

In our experiments, five-band (490, 565, 670, 865, 1020 nm), four-band (490, 565, 670, 865 nm) and three-band (490, 565, 670 nm) retrievals are conducted with the extended RemoTAP, as well as with the baseline RemoTAP. The validation of τ₅₅₀ for different retrieval setups over snow-dominant surfaces (c_snow>75 %) is shown in Fig. 4. Figure 4 uses the same color map as Fig. 3.

Figure 7 Example of PARASOL measurements and the RemoTAP results for ToA reflectance (panels a1, b1) and ToA DoLP (panels a2, b2). Solid lines in panels (a1, a2, b1, b2) refer to PARASOL measurements, dotted lines in panel (a1, a2) refer to the baseline RemoTAP results, and dotted lines in panel (b1, b2) refer to extended RemoTAP results. The pink lines refer to the 490 nm band, the blue lines refer to the 670 nm band and the yellow lines refer to the 865 nm band. The pixel is located at (43.75^∘ N, 96.65^∘ W), and the snow cover of this pixel is 68.41 %.

According to the validation with AERONET for τ₅₅₀ retrieval over snow-dominant surfaces shown in Fig. 4, we find the following:

1. For the five-band retrieval, the extended RemoTAP offers a better performance than the baseline RemoTAP, but the accuracy is still far from acceptable. This is caused by the fact that the retrieval algorithm has difficulty in fitting the 1020 nm band, where snow reflectance decreases significantly, making the reflectance much lower than that of the visible bands (see Fig. 1). In order to obtain a better result than the five-band retrieval, we investigate retrievals with a reduced number of spectral bands. First of all, the 1020 nm band is excluded because the algorithm has difficulties in fitting the strong spectral change between the visible bands and the 1020 nm band (synthetic experiments indicate that this only works with a very accurate first-guess value for the BRDF in the 1020 nm band). We also investigate an even more reduced set of wavelength bands including only the 490, 565 and 670 nm bands (three-band retrieval). The performance of the three-band and four-band retrieval is much better than that of the five-band retrieval (e.g., for the four-band retrieval the RMSE is reduced from 0.307 to 0.133). Interestingly, for the three-band retrieval, the RMSE is further reduced to 0.101. We do not show the scatterplot validation figures for SSA and AE retrievals over snow-dominant surfaces (c_snow>75 %) because not enough AERONET data are available for validation for snow-dominated pixels.

2. The extended RemoTAP has a significantly better agreement with AERONET than the baseline RemoTAP for the three-, four- and five-band retrievals. We would also like to emphasize that the fraction of pixels which meet the GCOS accuracy requirement also has a huge increase compared to the baseline RemoTAP of up to 45 percentage points for three- and four-band retrievals. This demonstrates the importance of adding the snow kernel to our BRDF model for retrievals over snow.

Figure 6 summarizes the RMSE and bias in τ₅₅₀, ω₅₅₀ and AE_440−870 for different retrieval setups for both snow-free (c_snow=0 %) and snow-covered pixels (0 % $<c_{\mathrm{snow}}\le \mathrm{100}\mathit{\%})$. For SSA, we would like to note again that the statistics over snow have limited value and no strong conclusions should be drawn from these numbers. For AE, although there are not enough statistics over snow-dominant surfaces (c_snow>75 %), there are enough pixels for validation when $\mathrm{0}\mathit{\%}<c_{\mathrm{snow}}\le \mathrm{100}\mathit{\%}$.

According to the validation results for retrievals over snow-free surfaces (c_snow=0 %) and snow surfaces ($\mathrm{0}\mathit{\%}<c_{\mathrm{snow}}\le \mathrm{100}\mathit{\%})$ shown in Fig. 5, we can conclude the following:

1. For snow-free surfaces (c_snow=0 %), the extended RemoTAP offers similar performance to the baseline RemoTAP for three-band and four-band retrievals, respectively. The only exception is the bias in ω₅₅₀ retrieval, where extended RemoTAP retrievals have larger bias than the baseline RemoTAP, but the difference is small (<0.01), especially also given the AERONET uncertainty in SSA and the limited number of validation points for SSA.

2. For snow-covered surfaces ($\mathrm{0}\mathit{\%}<c_{\mathrm{snow}}\le \mathrm{100}\mathit{\%})$, the baseline RemoTAP four-band and three-band retrievals fail to provide useful results with good accuracy, while both three-band and four-band extended RemoTAP retrievals are able to have good performance.

3. Comparing the extended RemoTAP for three- and four-band retrievals, the results of τ₅₅₀ retrievals are similar, but the four-band extended RemoTAP retrieval has slightly better performance for ω₅₅₀ and AE_440−870 retrievals (comparing the blue and brown bars). Therefore, the four-band extended RemoTAP retrieval is recommended for aerosol retrieval over snow from PARASOL.

In addition to the advantage in accuracy, taking four-band RemoTAP as an example, there is also an obvious increase in the fraction of successful retrievals (χ²<5), as shown in Fig. 6 for different snow fractions bins. On average, inclusion of the snow kernel leads to an increase in successful retrievals by 5.8 percentage points, and when the snow cover is large this increases to 18.9 percentage points, compared to the baseline, resulting in retrievals that are more successful by a factor ∼ 2. This indicates that the extended RemoTAP provides better performance in terms of both accuracy and goodness of fit for aerosol retrieval over snow.

Figure 7 shows an example of the difference between PARASOL measurements and the RemoTAP forward model simulations (after convergence) for ToA reflectance and DoLP for aerosol retrieval over a snow-dominated pixel. Comparing the ToA reflectance for the baseline and extended RemoTAP in Fig. 7a1 and b1, the difference between PARASOL measurements and the results of the extended RemoTAP forward model is much smaller than that of the baseline RemoTAP. Comparing the ToA DoLP in Fig. 7a2 and b2, the performance is quite comparable, and there is a slight advantage at the 865 nm band for the extended RemoTAP. It is important to note that there are some unphysical oscillations in the PARASOL measurements which are probably caused by interpolation error in the PARASOL level-1C processing for inhomogeneous scenes.

6 Conclusions

We extended the RemoTAP algorithm with a snow kernel in the BRDF model to carry out aerosol retrieval over snow from PARASOL MAP measurements. We conducted synthetic retrievals to investigate the necessity and advantage of extending RemoTAP and conducted retrievals on real measurements. For PARASOL retrievals with four bands (490, 565, 670, 865 nm), the extended RemoTAP retrieval results agree well with AERONET of optical properties, for retrievals both over snow and for snow-free areas. The RMSE, bias and fraction of retrievals within the GCOS requirements are 0.118 %, 0.023 % and 57.2 %, respectively, for τ₅₅₀ retrieval over partly snow-covered surfaces. This is much improved compared to the baseline RemoTAP (without the snow kernel) that yields an RMSE, bias and fraction of retrievals within the GCOS requirements of 0.221 %, 0.071 % and 48.2 %, respectively. Furthermore, the fraction of successful retrievals also improves by up to 18.9 percentage points compared with the baseline. The improvement is most striking for surfaces that have snow cover > 75 %, where the number of successful retrievals more than doubles and the fraction of retrievals that are within the GCOS requirement increases from 3.4 % to 48.1 %. The performance of the extended RemoTAP in terms of accuracy and goodness of fit is in good agreement with the expectation from the synthetic-data experiments. A limitation of the extended RemoTAP is that it is not able to fit PARASOL measurements at 1020 nm, where the snow albedo is substantially lower than at lower wavelengths. Therefore, the four-band extended RemoTAP is recommended as the best choice for aerosol retrieval over snow surfaces considering the performance for AOT, SSA and AE retrieval, from both synthetic and real retrievals.