On the retrieval of aerosol optical depth over cryosphere using passive remote sensing

The lack of aerosol information over the cryosphere introduces large uncertainties to our understanding of phenomenon, known as the Arctic Amplification (AA) and its feedback mechanisms. The aerosol optical depth (AOD) describes the optical characteristics of aerosol loading. This paper describes a novel algorithm, which retrieves AOD above snow-covered regions from the measurements of the up-welling radiation at the top of atmosphere, observed by the Advanced Along-Track Scanning Radiometer (AATSR) and the Sea and Land Surface Temperature Radiometer (SLSTR) instruments. The algorithm optimizes the generic eXtensible Bremen Aerosol/cloud and surfacE parameters Retrieval (XBAER) approach for longer wavelengths over the cryosphere. The algorithm utilizes the characteristics of solar bidirectional distribution properties of snow and aerosol at wavelength 3.7 μm to derive above-snow-AOD. Since the impact of fine-mode aerosol on 3.7 μm is ignorable, the retrieved AOD in this manuscript represents mainly coarse-mode dominated part. A novel method to extract the solar reflection part at 3.7 μm is presented and used in the surface parameterization. Two aerosol types (sea saltdominated and dust-dominated) are used and the best-fit type is derived by an iterative procedure, using a LookUp-Table (LUT) approach. Sensitivity studies of the impact on the retrieved AOD using XBAER algorithm, which investigate the impacts of aerosol type, snow surface emissivity and potential cloud contamination under typical AATSR observation conditions, are presented. The sensitivity studies show that the surface parametrization and aerosol typing are suitable for the retrieval of above-snow-AOD over the Arctic snow-covered region. AOD observations retrieved in this study from AATSR (2002–2012) observation collocated with those from the Aerosol Robotic Network (AERONET) sites over Greenland show good agreement. 72.1% of the match-ups fall into the expected error envelope of (± 0.15AOD ± 0.025). The AATSR derived above-snow-AOD at 0.55 μm research product has also been compared with Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) aerosol product, the Mineral Aerosols Profiling from Infrared Radiances (MAPIR) derived Infrared Atmospheric Sounding Interferometer (IASI) AOD research product, and the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2) AOD simulations over Greenland on April 2011. The comparison reveals that all datasets show similar patterns for the AOD above Greenland. The AOD is smaller in central Greenland and larger over the coastline regions. The XBAER derived above-snow-AOD has improved coverage, as compared to that of the existing AATSR aerosol product. The transition between above-snow-AOD and AOD derived over surrounding ocean surfaces does not indicate any systematic errors. Two aerosol transport events have been well-captured by the XBAER derived above-snow-AOD research product. The new algorithm is also applied to the SLSTR onboard Sentinel-3 demonstrating new SLSTR above-snow-AOD data products, and its value for research in the changing AOD during the period of Arctic Amplification.

and is a complex multiphase phenomenon. The positive feedback on temperature is determined by and affects a variety of atmospheric physical, chemical and biogeochemical processes (Serreze and Francis, 2006;Wendisch et al., 2019). Aerosol plays a significant role in the physico-chemical processes at high latitudes (Abbatt et al., 2019). Aerosol amounts and distributions are influenced by and have impact on the AA. This impact is direct through the absorptions and scattering of solar and thermal infrared radiation (Miller and Tegen, 1998;Samset et al., 2014) and indirect through the influence of aerosol on cloud properties (Haywood and Boucher, 2000;Sassen et al., 2003). In addition, deposition of aerosol, in particular that continuing black carbon, onto the Arctic snow/ice surfaces reduces the surface albedo (Qian et al., 2015;Yasunari et al., 2015). This increases the absorption of solar energy at the surface and thus is a positive feedback.
The chemical composition of Arctic Haze (Quinn et al., 2007) is largely determined by variety of process. Aerosol with Brown and Black Carbon is produced to a minor extend by the oxidation of natural sand anthropogenic missions of hydrocarbons and to a major extend by fires/ biomass burning (Heintzenberg et al., 1981;Quinn et al., 2007;Nielsen et al., 2019). Black carbon has been widely discussed in the literature (Hansen and Nazarenko, 2004;Ramanathan and Carmichael, 2008;Bond et al., 2013;Jacobi et al., 2019). However, other important coarse mode dominated aerosol types, such as mineral dust and sea salt, have not yet received sufficient attention (Zwaaftink et al., 2016;Frey et al., 2019). Hesaraki et al. (2017) highlighted the need to qualify better the impacts of Arctic coarse-mode dominated aerosols.
The dust mass column loading in the atmosphere is approximately 0.11 million metric tons between 60°and 90°N (Takemura et al., 2009). The deposition of dust in the Arctic region is estimated to be 6.80 million metric tons per year (Vincent, 2018). Atmospheric dust in the Arctic consists of Asian dust transported by the westerlies (~38%), African dust carried initially by the trade winds (~32%) as well as local dust (~27%) (Takemura et al., 2009;Zwaaftink et al., 2016;van der Does et al., 2018). The local contribution of dust aerosol to the total amount of aerosol in the Arctic is more important than previously thought (Zwaaftink et al., 2016). Local dust contribution, particularly in autumn, may be due to the exposure of fine sediments to the atmosphere, as a result of retreating of ice masses (Vincent, 2018) and other periglacial processes (Bullard, 2013).
The impact of dust on the Arctic cryosphere is poorly-understood, because of the lack of knowledge and measurements (Boy et al., 2019). The reduction of glacier albedo over Greenland and Iceland is partially ascribed to dust deposition (Dumont et al., 2014;Wittmann et al., 2017). The accuracy of satellite sea/ice surface temperature products is also affected by dust aerosols in the Arctic, because dust aerosol affects the satellite observed surface emissivity in thermal infrared channels (Vincent, 2018). Better knowledge of the radiative impacts of mineral dust is required to assess its role in the Arctic Amplification (Lambert et al., 2013;Kylling et al., 2018).
Sea salt is a dominant primary aerosol source in Arctic regions. These particles originate from bubble bursting over the Arctic ocean (Nilsson et al., 2001), mobilized saline snow (Huang and Jaeglé, 2017) and blooming of frost flowers by wind over sea-ice-covered areas (Xu et al., 2016). Frost flowers grow on imperfections in the surface ice during periods of sub-zero temperatures around −22°C. Spiky structures form, which have been found to house microorganisms. The reduction of sea-ice cover during the AA may increase the amount of sea salt aerosols in the Arctic atmosphere (Struthers et al., 2011;May et al., 2016). Sea salt particles are a major source for Arctic background aerosol, especially during summer (Deshpande and Kambra, 2014). Model simulations show that a positive response of sea salt aerosol emissions to the decrease of sea ice cover will lead to an increase of 23% of natural AOD due to the increase of sea salt aerosols by a factor of 2-3 (Struthers et al., 2011). However, Arctic sea salt aerosol from blown snow over sea ice is a missing natural source in model simulations, especially during winter (Frey et al., 2019). Wind mobilization of snow salinated by interacting with frost flowers also contributes to the Arctic sea salt aerosol (Obbard et al., 2009). May et al. (2016) reported that sea salt aerosols have the potential to change Arctic cloud formation and snowpack. Sea salt aerosol debromination may be a dominant source of tropospheric bromine (Zhu et al., 2019;Domine et al., 2004) and affects atmospheric bromine cycle in the polar regions (Hara et al., 2018).
In-situ measurements and model simulations are the two major ways to characterize aerosol properties, especially coarse model dominated aerosols, in the Arctic. For dust aerosols, the High Latitude and Cold Climate Dust (LCCD) community (http://www.hlccd.org/) through a worldwide collaboration plans to harmonize the spatiotemporal limited ground-based measurements to improve understanding of high latitude dust emissions. Individual measurements have been obtained to understand the characteristics and radiative effects of Arctic dust (Intrieri and Shupe, 2004;Schnell and Jefferson, 2015). Modellers highlight the importance of mineral dust for the Arctic aerosol budget and climate (Zwaaftink et al., 2016;Sand et al., 2017;Kylling et al., 2018). However, the current aerosol-transport models (in particular global models with coarse resolution, personal communication with Dr. Bernd Heinold) tend to underestimate local, high-latitude dust emissions and they underestimate/overestimate long-travelled mineral dust from the Sahara or Gobi Desert (Sand et al., 2017;Weinzierl et al., 2017). The underestimation is attributed most probably to a lack of dust sources in the models (Koven and Fung, 2008;Tegen et al., 2013). This is because the dust sources and other aerosol sources occur at sub-grid scales, which cover relatively small areas (Sato et al., 2016). Transport of Sahara or Gobi dust tends to be overestimated or underestimated. The former is because of insufficient wet removal process (Joussaume, 1990;Choobari et al., 2014), and the latter due to a lack of a source term in the models.
The ACLOUD/PASCAL ("Arctic CLoud Observations Using airborne measurements during polar Day/"Physical feedbacks of Arctic boundary layer, Sea ice, Cloud and AerosoL) campaign provides shipbased measurements of aerosol during 24 May to 20 July 2017 from Bremerhaven to Longyearbyen and Tromsö (Wendisch et al., 2019). The first tropospheric sea-salt aerosol dataset obtained during the Atmospheric Tomography (ATom) mission has been presented by Murphy et al. (2019) and high sea-salt aerosol concentration has been observed over the Arctic ice-covered regions.
Global knowledge about aerosol optical properties under cloud-free conditions (King et al., 1999;Kaufman et al., 2002) is retrieved from passive remote sensing instrumentation, making spectral (Levy et al., 2013;Hsu et al., 2013;Mei et al., 2017aMei et al., , 2017bLyapustin et al., 2018), angular (Popp et al., 2016), temporal (Govaerts and Luffarelli, 2018;Gupta et al., 2019) and/or polarization (Tanré et al., 2011) measurements. Recently, the retrieval of aerosol optical depth (AOD) has been achieved above cloud (Jethva et al., 2013;Meyer et al., 2015;Sayer et al., 2016;Mei et al., 2019b). Up to the present, the retrieval of aerosol properties using passive remote sensing over the polar regions using visible-NIR spectral range (Tomasi et al., 2015) is limited. This is because of the combination of bright surfaces, low aerosol loading under large sun zenith angle observation condition makes deconvolving the AOD from the surface reflectance challenging.
Research algorithms for AOD have been discussed using the dual viewing capabilities of AATSR (Advanced Along-Track Scanning Radiometer) (Istomina et al., 2011;Mei et al., 2013a), synergy of MODIS (Moderate Resolution Imaging Spectroradiometer) observations onboard Terra and Aqua (Mei et al., 2013b) and PARASOL (Mei et al., 2019b). However, these proposed algorithms have not yet been used operationally to produce AOD over the Arctic snow/ice covered regions due to remaining cloud identification, surface parametrization and aerosol typing issues.
The detection/retrieval of surface/atmospheric properties over the Arctic snow/ice covered regions frequently uses observations in thermal infrared channels, e.g. 3.7 μm. Allen et al. (1990) proposed the use of the reflected part of 3.7 μm for snow/cloud discrimination. A similar idea has been used for sea and lake ice detection (Dorofy et al., 2016). Roger and Vermote (1998) proposed a method to extract the solar reflection from the reflected and emitted component, observed by Advanced very-high-resolution radiometer (AVHRR). This method also takes into account atmospheric effects, by using empirical corrections. The method from Roger and Vermote (1998) can be used to extract the solar reflective part of 3.7 μm over both land and ocean conditions. The reflected part of 3.7 μm has also been used for the retrieval of surface, cloud and aerosol properties. There are several advantages of using the 3.7 μm channel to retrieve surface properties. This is because the channel is sensitive to surface properties, but the reflectance of most aerosol types is negligible (Kim et al., 2008). Relevant and more detailed explanations are found in Kaufman and Remer (1994). Platnick et al. (2001) uses the 3.7 μm reflectance to retrieve the cloud optical thickness and droplet size over snow and ice surfaces. Because the contribution of snow/ice to the TOA reflectance is relatively small, the use of 3.7 μm reflectance to retrieve atmospheric parameters can significantly reduce the impact of snow/ice surfaces (Platnick et al., 2001). This idea has been implemented in the MODIS cloud products (Platnick et al., 2017). Examples of the use of the 3.7 μm reflectance for the retrieval of AOD are found in previous publications (Kim et al., 2008;Istomina et al., 2011;Mei et al., 2014).
Two reasons make the retrieval of coarse-mode dominated AOD using 3.7 μm possible: (1) the contribution of snow/ice to the TOA reflectance being relatively small (dark surface); (2) the 3.7 μm channel is insensitive to aerosols, except for the coarse-mode dominated aerosol types such as dust (Kim et al., 2008). In this study we use the 3.7 μm channel to derive above-snow-AOD at 0.55 μm (hereafter referred to above-snow-AOD) over the snow/ice covered regions. For this study the eXtensible Bremen Atmospheric and surfacE parameter Retrieval (XBAER) algorithm has been extended to derive AOD over snow (see Mei et al., 2017a;Mei et al., 2018a). This followed our investigation of the maximum "information content" obtained by combining the observations of the solar and thermal infrared channels with an improved knowledge of cloud (Mei et al., 2017b). In addition, the surface parametrization (Mei et al., 2017a) and aerosol typing (Levy et al., 2013) were optimized for use in this algorithm. The standard XBAER algorithm has been used to derive cloud properties (Mei et al., 2018b) and both cloud and aerosol properties simultaneously (Mei et al., 2019b).
This manuscript is organized as follow. AATSR instrument characteristics and the pre-processing ("parallax effect" correction and cloud identification) of AATSR before AOD retrieval are both presented in Section 2. Other aerosol data products, used for the comparison with XBAER derived above-snow-AOD, and the collocation method are explained in Section 2. The AOD retrieval algorithm (surface parameterization and aerosol typing) is presented in Section 3. Section 4 focuses on the sensitivity study, undertaken to assess and understand the impact of aerosol typing, surface parameterization and potential cloud contamination on the AOD retrieval. A preliminary comparison of XBAER-derived AOD with the Mineral Aerosols Profiling from Infrared Radiances (MAPIR) derived InfraRed Atmospheric Sounding Interferometer (IASI), Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2), Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP) AOD products and the validation using Aerosol Robotic Network (AERONET) are presented and discussed in Section 5. A case study of Sahara dust transport to Greenland using the XBAER-derived AOD is presented in Section 6. The application of the XBAER algorithm to the observations of the Sea and Land Surface Temperature Radiometer (SLSTR) instrument is described in Section 7. The above-snow-AOD data product, retrieved using XBAER, demonstrates the capability and advantages of this new dataset. Section 8 provides a summary and the conclusion for this study.

AATSR instrument and pre-processing
The AATSR onboard European Space Agency (ESA) satellite ENVISAT (ENVIroment SATEllite) is the successor of ATSR-1 (Along Track Scanning Radiometer 1) and ATSR-2 (Along Track Scanning Radiometer 2) on board ERS-1 (European Remote Sensing 1) and ERS-2 (European Remote Sensing 2). The ATSR-2 measurements on the ERS-2 were extended as compared to those on ERS-1 to include the visible. AATSR was further optimized and provides observations over the period 2002-2012. The successor to AATSR is SLSTR (Sea and Land Surface Temperature Radiometer) on board Sentinel-3A and Sentinel-3B, launched during February 2016 and April 2018, respectively.
AATSR is a dual-view (nadir and 55°forward) radiometer, with a resolution of 1 km at nadir observations. The swath width of AATSR is 500 km. AATSR provides global coverage in 3 days for IR channels and 6 days for visible channels. There are seven wavelengths of AATSR instruments for both nadir and forward observations: 0.55, 0.66, 0.87, 1.6, 3.7, 11, and 12 μm. The ATSR family of instruments and SLSTR have the primary goal of observing sea surface temperature.
In addition, the dual-viewing observation strategy provides unique capabilities to derive AOD globally (Popp et al., 2016). The dualviewing observation capability enable the surface properties to be better characterized. This enables the Top Of the Atmosphere, TOA, radiance to be separated into that arising from scattering by aerosol and that from surface. The two viewing observations have scenes collocated at ground level. Consequently, there is a shift of observation scene (in the along-track direction). This shift depends on the surface elevation and, or, target (e.g. cloud) position. This behavior is called the AATSR "parallax effect". The "parallax effect" is used to derive cloud (Carbajal Henken et al., 2014) and aerosol (Virtanen et al., 2014) layer height. However, for the retrieval of aerosol over elevated regions such as Greenland, the removal of the "parallax effect" is needed before the retrieval.
In the XBAER algorithm, the "parallax" shift estimation method, described in Virtanen et al. (2014), is used. The scene shift is estimated by the term (tan(θ f ) − tan(θ n )) × h, where θ f and θ n are the viewing zenith angles for forward and nadir observations respectively, and h is the surface elevation. The National Oceanic and Atmospheric Administration (NOAA) provides the required surface elevation data (https:// www.ngdc.noaa.gov/mgg/topo/globe.html, last access: 4 June 2019). To minimize further the AATSR "parallax effect" for the retrieval of above-snow-AOD in this study, the retrieval is performed on a 9 × 9 scene or scene averaged box.
The XBAER algorithm (described later in Section 3) uses the anisotropic properties of surface and aerosols at 3.7 μm. Accurate knowledge of the AATSR instrument observation geometry is necessary to assess the "information content" (Tanre et al., 1996) of AATSR measurements for a given set of observation angles. We performed a statistical analysis of solar zenith angle (SZA), viewing zenith angle (VZA) and relative azimuth angle (RAA) over Greenland (60°N-85°N, 75°W-10°W) for April and September 2011, using the AATSR instrument. April and September are reported to be months having Arctic haze and background aerosol in the AERONET climatology AOD dataset (Mei et al., 2019a). The statistical analysis was performed as follow: (1) Excluding AATSR observations with SZA smaller than 0°; (2) The ranges of SZA, VZA and RAA are [0°,90°], [0°,90°], [0°,360°]; (3) Statistical analysis for SZA, VZA and RAA for every 5°interval. Large solar zenith angles (SZA > 50°) characterize the passive remote sensing observations in the Arctic, with the most frequently observed ranges of angles being [65][66][67][68][69][70][71][72][73][74][75] and [70][71][72][73][74][75][76][77][78][79][80], for April and September respectively. Consequently, we selected 70°as a typical SZA in the Arctic region. This will be used later to investigate the surface and atmospheric anisotropic properties. The SZA distributions for forward and nadir are the same ( Fig. 1(a) and (b), (c) and (d)). This is because the measurements have the same time of observation (ignoring the~130 second delay of forward observations) over the same geographic location. Fig. 1(e)-(h) shows the available pairs of (VZA, RAA), which provide valuable information to assess the AATSR observed surface and atmospheric anisotropic properties. Theoretically, the ranges for VZA and RAA are [0°-90°] and [0°-360°], respectively. The BRDF (bidirectional reflectance distribution function) is used to describe the surface and atmospheric anisotropic properties (Lucht et al., 2000). The BRDF is typically shown on a polar plot, as shown in Fig. 1(e)-(h). The large values (red) in Fig. 1(e)-(h) show the satellite has a higher probability of observing the surface and atmosphere at the selected (VZA, RAA) observations. Fig. 1(e) and (g) (AATSR nadir observation for April and September), (f) and (h) (AATSR forward observation for April and September) shows similar patterns, indicating that AATSR has stable observation geometries. The AATSR nadir and forward observations show different patterns, indicating that the dual-viewing observations is used to separate the aerosol signal from the surface contribution observed by AATSR at TOA.

SLSTR onboard Sentinel-3
The European Commission's Copernicus Sentinel-3A/3B Sea and Land Surface Temperature Radiometer (SLSTR) is the successor of AATSR onboard ENVISAT. Sentinel-3A and 3B were launched on the 16th of February 2016 and 25th of April 2018, respectively. SLSTR is built and designed using the heritage of AATSR. In particular, the dualviewing technique, facilitates the separation and retrieval of surface and atmospheric properties. The SLSTR provides a backward (oblique) observation rather than the forward observation in AATSR. The SLSTR instrument has some new and advanced features, such as more spectral bands and a better spatial resolution of 0.5 km for visible and SWIR channels. Moreover, SLSTR provides an OLCI overlapped swath coverage (1400 km) for nadir observation and an increased dual-view swath width of 740 km. In Table 1 AATSR and SLSTR characteristics are compared.

IASI MAPIR dust AOD research product
The Infrared Atmospheric Sounding Interferometer (IASI) instruments fly on board three satellite platforms launched successively in October 2006 (Metop-A), September 2012 (Metop-B) and November 2018 (Metop-C). From their polar sun-synchronous orbit they each offer an almost global Earth coverage twice a day, at about 9:30 and 21:30 local solar time. IASI is a nadir-viewing high-resolution thermal infrared Fourier transform spectrometer. It has a wide swath of 2200 km (satellite viewing angles up to 48.3°on both sides) containing 30 elementary fields of view each composed of 4 pixels. The pixel size and shape depend on the viewing angle: from a circle of 12 km ground diameter at nadir to an ellipse of 39 km by 20 km for the biggest viewing angles. IASI records the radiance from 645 cm −1 to 2760 cm −1 at a resolution of 0.5 cm −1 after apodization. The radiometric noise in the thermal infrared atmospheric window (800-1200 cm −1 ) is about 0.2 K (Clerbaux et al., 2009).
The Mineral Aerosols Profiling from Infrared Radiances (MAPIR) algorithm retrieves vertical profiles of dust aerosol concentrations from IASI cloud-free measurements in the thermal infrared atmospheric window. Callewaert et al. (2019) describes its version 4.1, which was used to process of the IASI/Metop-A data, for latitudes between 60°S and 60°N. MAPIR uses the Optimal Estimation Method (Rodgers, 2000) and the Radiative Transfer model for TOVS version 12 (RTTOV v12, Saunders et al., 2017). The MAPIR retrieval state vector is composed of the surface temperature and the dust aerosol concentration in 7 layers centered every 1 km from 0.5 to 6.5 km altitude. The dust aerosols are modelled as spherical particles, with a log-normal distribution with 2 μm effective radius and the "dust-like" refractive index from GEISA-HITRAN. The 10 μm AOD is calculated by the integration of the concentration profiles and multiplication by the dust extinction coefficient at 10 μm. The 550 nm AOD is obtained by using a conversion factor, based on the knowledge of the dust extinction at both wavelengths.
The standard version 4.1 of MAPIR processes only the cloud-free IASI spectra, for cloud free scenes determined by the EUMETSAT IASI level 2 cloud fraction. That product often flags dust plumes as clouds, especially in areas where they rarely occur, therefore preventing the standard retrieval from being performed. Another limitation is that, the standard retrieval only runs when a signature of the dust presence is observed. This approach misses low AOD plumes but was necessary for the processing, in terms of computing resources. Finally, the standard retrieval does not cover polar regions. Consequently, we processed all scenes (i.e. including cloudy scenes and not using the dust presence prefilter) in the target area (20°N-90°N, 60°W-10°E) for the period April 2011, in preparation for the comparison with XBAER derived abovesnow-AOD.
Two additional changes to MAPIR were made for this analysis. Firstly, the minimal a priori value for the dust concentration (used when the a priori climatology shows no dust aerosols, which is the case for polar areas) was increased from 0.1 to 2 particles/cm 3 , in each retrieval layer (corresponding to a 10 μm column AOD of about 0.06). Secondly, the standard deviation on the a priori Ts was decreased to 1 K (instead of 5 K over ocean and 15 K over land). Both those changes are specific to this analysis and answer to the attempt to detect a small dust amount under difficult conditions (very low surface temperature in polar areas means very low signal in the thermal infrared). The standard deviation on the a priori Ts for the standard MAPIR retrievals is much higher because the a priori used is the Ts retrieved with the EUMETSAT IASI algorithm. The latter shows significant bias over deserts. Consequently, the MAPIR dust retrievals need a low constraint on Ts. In the specific case of polar areas, the sensitivity of the thermal infrared seems to be significantly lower, making it difficult (or impossible) to retrieve both aerosols and Ts with the weak constraints, used in the standard MAPIR. Consequently, we increased the constraint on the Ts retrieval, and we increased the a priori aerosols concentration, to boost the sensitivity to the aerosols. This results in the EUMETSAT IASI Ts being wrong. The adapted MAPIR retrieval is not able to correct the bias but instead modifies the aerosol content.

CALIOP
The Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO), the joint satellite mission between United States (NASA) and France (Centre National d'Etudes Spatiales/CNES), is a sun-synchronous polar-orbit satellite, part of the international Afternoon-Train (A-Train) constellation of Earth-observation satellites (Stephens et al., 2018), since June 2006 (Winker et al., 2010). In order to fulfill its main science objective, to provide the three-dimensional observations of aerosol and clouds on a global scale, CALIPSO comprises three co-aligned near-nadir viewing instruments, namely, an Imaging Infrared Radiometer (IIR), a Wide Field Camera (WFC) and the Cloud-Aerosol Lidar with Orthogonal Polarization (CALIOP), an elastic backscatter Nd:YAG lidar (Winker et al., 2007). CALIOP operates a twowavelength transmitter that emits simultaneously linearly polarized light at 0.532 μm and 1.064 μm, and a three-channel receiver, that separately measures the backscatter by aerosol and clouds signals, and specifically, the 1.064 μm and the 0.532 μm parallel and perpendicular backscattered components, relative to the polarization plane of the transmitted 0.532 μm beam .
Through utilizing CALIOP, CALIPSO provides high-resolved aerosol and cloud profiles of backscatter coefficient (at 0.532 μm and 1.064 μm) and particulate depolarization ration (at 0.532 μm), along the CALIPSO orbit track and in addition the feature type and aerosol subtype classification of the detected atmospheric layers ). To be more specific, CALIOP Level 2 (L2) Version 4 (V4) algorithm classifies the detected atmospheric features into distinct categories, namely, "clear air", "cloud", "tropospheric aerosol", "surface", "subsurface", "totally attenuated" and "low/no confidence" . In case of detected features classified as "tropospheric aerosols", the algorithm attempts further classification based on the optical properties of the detected atmospheric layers (e.g. integrated attenuated backscatter at 0.532 μm and estimated particulate depolarization ratio), the altitude of the detected layer (layer Top/Base) and the sublayer type of surface (land/ocean). Accordingly, the aerosol classification algorithm classifies the detected aerosol features into distinct aerosol subtypes, including among others, the "marine", "dust", "dusty marine" and "polluted dust" aerosol classes Kim et al., 2018).
To focus only the pure-dust observations of CALIPSO, we utilize the "Pure-Dust product". The methodology of the pure-dust product was initially established in the framework of the European Aerosol Research Lidar Network (EARLINET; Pappalardo et al., 2014; https://www. earlinet.org; last access: 21/05/2019), for the extraction of the puredust component of the total aerosol load ). Accordingly, the developed methodology was implemented to CALIPSO, in the framework of the European Space Agency (ESA) "LIdar climatology of Vertical Aerosol Structure for space-based lidar simulation studies (LIVAS; http://lidar.space.noa.gr:8080/livas/; last access: 21/ 05/2019)" project (Amiridis et al., 2015). The suggested methodology assumes the CALIOP "dust", "dusty marine" and "polluted dust" aerosol subtypes as external mixtures of two aerosol types, and specifically, of a pure-dust component and a non-dust component. Accordingly, CALIOP particle linear depolarization ratio measurements are used in conjunction with the strongly depolarizing optical property of mineral dust Hofer et al., 2017;Mamouri and Ansmann, 2017) to retrieve the CALIPSO-based backscatter coefficient profiles of pure-dust at 0.532 μm. In addition, through the implementation of suitable dust Lidar-Ratio (LR) values to backscatter coefficient profiles of pure-dust (Baars et al., 2016), the extinction coefficient profiles of pure-dust at 0.532 μm along the CALIPSO orbit-track are obtained .
For the needs of the study, the CALIPSO L2 pure-dust profiles are further processed following the CALIPSO Level 3 (L3) quality-assurance methodology Tackett et al., 2018), in order to produce the CALIPSO Cloud-Free L3-equivalent pure-dust Optical Depth product in 1°× 1°grid spatial resolution, through the vertical integration of the pure-dust extinction coefficient profile at 0.532 μm.
The CALIPSO-based pure-dust product has been validated against AERONET over North Africa and Europe in Amiridis et al. (2013), while the methodology-related uncertainties are extensively discussed in Marinou et al., 2017. The pure-dust product is used for the assessment of dust events (e.g. Kosmopoulos et al., 2017;Solomos et al., 2018), to describe the three-dimensional dust climatology and dust transport pathways (e.g. Proestakis et al., 2018), and for the validation and evaluation of dust models (e.g. Konsta et al., 2018). In this paper, the pure-dust product is used as a reference dataset of Dust Optical Depth (DOD) and, for the evaluation of the under-development algorithm for aerosol retrievals over high-reflectivity surfaces.

AERONET and collocation method
AERONET (AErosol RObotic NETwork) is a ground-based aerosol observation network (Holben et al., 1998). AERONET includes permanent sites, which provide long-term, continuous and quality-controlled aerosol optical/microphysical/radiative dataset, and observations provided by campaign measurements (Holben et al., 1998). The quality of AERONET derived aerosol properties are grouped into three levels: Level 1.0 (unscreened), Level 1.5 (cloud-screened), and Level 2.0 (cloud-screened and quality-assured). The Level 2 data is recommended for the validation of satellite aerosol products. The AERONET aerosol product has updated to version 3, in which a better-quality control and cloud screening have been included compared to version 2.0 (Giles et al., 2019), although potential problem may occur in presence of cirrus clouds (Smirnov et al., 2018). This study uses the AERONET version 3 Level 2 product.
The AERONET-Satellite spatial-temporal collocation method has been widely used for the aerosol community. The method is described in details in Ichoku et al. (2002). The same collocation method is used in this study. XBAER-derived above-snow-AOD are averaged over an area of ± 25 km around the AERONET site, and the AERONET observed AOD are averaged over a temporal window of ± 30 min around the satellite overpass time (Levy et al., 2013). At least five XBAER-derived AOD over the ± 25 km square area (~20%) and two sun photometer measurements within ± 30 minute time window are required (Mei et al., 2011;Levy et al., 2013). The expected error (EE) envelope, defined as ( ± 0.15AOD ± 0.05) (Remer et al., 2005), is also adapted for this study.
There are four permanent sites over Greenland: Thule (76.516°N, 68.769°W), Ittoqqortoormiit (70.485°N, 21.951°W), Kangerlussuaq (66.996°N, 50.621°W), and Narsarsuaq (61.156°N, 45.419°W). All these four sites are used for the validation of XBAER derived above-snow-AOD. However, these four sites are located near the coastline, therefore typically snow-free. Taking the aerosol geographical homogeneity (no strong source) into account, the following steps are performed before the collocation. (1) Taking AERONET location as the central point, if all pixels situated in the ± 25 km square area are 100% snow covered, then the standard collocation between AERONET and XBAER derived above-snow-AOD will be performed; (2) if (1) is not the case, the nearest 100% snow covered pixel (compared to the AERONET site) is determined based on the MODIS monthly snow cover product (Hall et al., 2016); (3) taking the nearest full snow covered pixel as the lower left corner of the 50 km × 50 km square area, if all pixels in the 50 km × 50 km square area are 100% snow covered region, the nearest pixel will be used as the "new" AERONET position and the standard collocation between "new" AERONET position and XBAER derived above-snow-AOD will be performed. Fig. 2 shows the CALIOP orbits for April 2011 and corresponding AERONET positions for all months based on the collocation method above.

Algorithm
The above-snow-AOD in the XBAER algorithm is retrieved using our knowledge of the different bidirectional properties of the solar reflection of the measured radiation in the channel 6 (3.7 μm) of AATSR/ SLSTR instruments (hereafter referred to as AATSR) for snow and aerosols. The absolute solar reflection for spectral channels in the shortinfrared (e.g. 2.1 μm) is assumed to have no strong bidirectional properties (absolute values are small) for selected satellite observation geometries over a snow surface (Platnick et al., 2001). The information (see Tanre et al., 1996 and references therein), used to decouple the contribution of snow reflectance from the TOA reflectance from satellite observations, is the difference between surface and atmosphere bidirectional properties at the AATSR observation geometries and wavelengths (see Fig. 1). Fig. 3 shows the bidirectional properties of surface (snow) and three typical aerosol types (soot, dust and sea-salt) (Quinn et al., 2007). These three aerosol types are chosen, because the typical transport aerosol from middle-low latitude are soot and dust while the background aerosol in the Arctic is frequently sea-salt dominated. The simulations of TOA reflectance at 3.7 μm have been performed using the radiative transfer package SCIATRAN (Rozanov et al., 2014). The forward and backward scattering in this study are defined by RAA = 0°and 180°, respectively. The snow surface was considered as a layer with an optical thickness of 1000 and a geometrical thickness of 1 m composed of ice particles and positioned above a black surface. The snow layer is assumed to be vertically and horizontally homogeneous without any surface roughness and composed of monodisperse snow particles. Their single-scattering properties are described by sparsely packed aggregates of eight hexagonal solid columns as described in Yang et al. (2013). The impact of snow impurities is neglected. The radiative transfer calculations were performed ignoring the contribution of thermal emission, assuming the black surface, and the atmosphere consisting of dust-dominated (Dubovik et al., 2006), soot-dominated and sea-salt-dominated aerosol types (Hess et al., 1998). The details of the settings are listed below: (1) Snow surface has been simulated using Yang et al. (2013) database with snow particle shape of 8403 elements according to recent investigation (Järvinen et al., 2018); (2) Snow grain size is set to be 100 μm according to the statistical investigation of MODIS snow product derived from MODIS snow covered-area and grain size retrieval algorithm (Painter et al., 2009) over Greenland; (3) Aerosol type: dust-dominated type (Dubovik et al., 2006) and soot-dominated type sea-salt-dominated type (Hess et al., 1998); (4) Solar zenith angle is set to be 70°; (5) AOD at 0.55 μm is set to be 0.5, the maximal value used later in Look-Up- Table; (6) Meteorological conditions are set to April at 75°N latitude from B2D Chemical Transport Model (CTM) ).
The analysis of snow BRDF properties for visible channels is described in previous publications (Clémence et al., 2015;Gatebe and Poudyal, 2018;Jiao et al., 2019). In the visible channels, the snow single scattering albedo is close to 1, thus the snow BRDF is not significantly influenced by the single scattering parameter, which is used to describe the surface anisotropic properties (Aoki et al., 2000). This study focuses on the snow properties at 3.7 μm. According to Fig. 3(a), we see that the anisotropic reflection of a snow layer is significant at 3.7 μm, which is consistent with previous publications (Leroux et al., 1997;Aoki et al., 2000) because the snow particles are more absorptive at 3.7 μm compared to visible channels. Despite the existing reflectance asymmetry of snow layer, for the AATSR observation geometries (see Fig. 1), the snow BRDF properties (normalized to nadir observation) are not so strong, especially compared to certain aerosol particles (e.g. dust). Since the new algorithm uses a specific wavelength (3.7 μm), which is not the same as classical aerosol retrieval algorithm using visible channels (e.g. Levy et al., 2013;Mei et al., 2017aMei et al., , 2017b, the retrieval is not sensitive to certain aerosol types, such as fine-mode dominated aerosol types (e.g. biomass burning). Fig. 3(c) shows an example of the BRDF feature of soot-dominated aerosol type (polluted continental aerosol as defined in Hess et al., 1998). We can see that the contribution of aerosols at 3.7 μm to almost all geometries is ignorable, thus our retrieval algorithm is only sensitive to coarse-mode dominated aerosol types, such as sea salt dominated and dust dominated aerosol types in the Arctic. Fig. 3(b) and (d) shows the BRDF patterns of those two aerosol types at 3.7 μm. Both aerosol types show stronger anisotropic reflection for the AATSR observation angles compared to snow surface, indicating the aerosol reflectance dominates the angular-depended feature of the TOA reflectance. Thus, the change of satellite L. Mei, et al. Remote Sensing of Environment 241 (2020) 111731 observed TOA reflectance at two different observation directions is explained in large part by the aerosol loading. Similar conclusions can be drawn for other solar zenith angles (e.g. 60°and 80°).

Cloud identification
The cloud identification is performed before the AOD retrieval. A three steps method is used in this study to minimize cloud contamination in the above-snow-AOD research product. A cloud identification algorithm to detect cloud over snow/ice regions has been descried in Istomina et al. (2010). This algorithm uses the observations of the AATSR visible (VIS) and thermal infrared (TIR) channels. The difference in the spectral response of snow and cloud yields a set of relative thresholds. The algorithm has been used for the retrieval over snow/ice covered regions as described in Istomina et al. (2011) and Mei et al. (2013a). This cloud identification algorithm (Istomina et al., 2010) is readily implemented and suitable for removing scenarios with cloud for subsequent retrieval of AOD in cloud free scenes for large geographic scale (e.g. global scale). Consequently, this cloud identification method is used as the first-step of the cloud identification method in the XBAER algorithm. However, potential cloud contamination of scenes remains (Mei et al., 2013a). The second step for the cloud identification in this study is the use the corresponding MERIS (MEdium Resolution Imaging Spectrometer) observations (AATSR and MERIS are both onboard ENVISAT, thus they observe the same position at the same time, ignoring the AATSR parallax effect and~130 second time delay of AATSR forward observation). The cloud identification algorithm for MERIS to separate cloud from snow is described in Mei et al. (2017b) (Fig. 16 in Mei et al. (2017b)). A similar idea has been used as the pre-processing step for the retrieval of snow albedo (Zege et al., 2015). In order to avoid the cloud adjacency effect (Koren et al., 2007), if a pixel is marked as cloud, the surrounding 5 × 5 pixels will be automatically marked as cloud as well. In addition, a similar method involving manual-checking (Dr. Vincent, personal communication) as descried in Vincent (2018) is used as the third step. Fig. 4 shows an example of the comparison between different cloud identification methods. Fig. 4(a) shows the RGB composition figure of an AATSR orbit (33187) over Greenland on the 5th of July 2008. Fig. 4(a) indicates that almost all central part of this orbit is covered by cloud. Fig. 4(b) shows the cloud identified using the approach from Istomina et al. (2010), which is the first step of the cloud identification used in XBAER. Unfortunately,~30% of the cloud-covered scenes are identified as cloud free scenes. This introduces bias in the retrieved AOD, if ignored. Fig. 4(c) shows the result of the three-steps cloud identification method used in XBAER algorithm. Cloud contamination has been reduced to being of negligible significance. Fig. 4(d) shows the final cloud free and clouded ground scenes taking the cloud adjacency effect (Koren et al., 2007) into account.

Surface reflectance estimation
Measurements of the upwelling radiation at 3.7 μm comprise both solar radiation scattered at the surface and the thermal emission from the surface, which is strongly temperature dependent. The results of sensitivity tests support the use of the solar reflection part of 3.7 μm in cloud free scenes to derive above-snow-AOD. This section describes the approach used to extract the solar reflection part of 3.7 μm.
The clear-sky radiance at channel 3.7 μm at the TOA is given by: where L(θ, θ 0 , φ) is the measured radiance from satellite observation, L s (θ, θ 0 , φ) and L t (θ, φ) are the solar reflected and thermal radiation of the surface-atmosphere system. θ, θ 0 and φ are the viewing zenith angle, solar zenith angle, and relative azimuth angle, respectively. Following the method suggested by Roger and Vermote (1998), we can extract the solar reflective part. In particular, we assume that the thermal emission of the atmosphere can be neglected and second term in Eq.
(2) is given by: where ε is the surface emissivity, B(T s ) is the Plank function with surface temperature T s . adapted XBAER cloud mask + adjacency pixels (if a pixel is classified as cloud, the surrounding 5 × 5 pixels will be considered to have cloud present).
The atmospheric thermal emission caused by gaseous/aerosol absorption for cloud-free condition is assumed to be negligible at 3.7 μm. The impact will be estimated in Section 3.4.
Following Spangenberg et al. (2001), brightness temperature measured at 11 μm (T 11 ) instead of T s is used in the Planck function, and the reflected solar component can be calculated as following: The solar reflected radiation at TOA is given in the path radiance representation (Chandrasekhar, 1950;Kaufman et al., 1997) where R s 0 (θ, θ 0 , φ) is the TOA reflectance assuming black surface, ξ(θ, θ 0 ) is the total (diffuse and direct) transmittance from the sun to the surface and from surface to the satellite, s is the atmospheric spherical albedo, A is the Lambertian surface albedo.
Accounting for the following relationship between radiance and reflection function: and substituting Eq. (6) into Eqs. (4) and (5), where = E E π 0 is the normalized solar irradiance. Based on the Kirchoff relation for a Lambertian reflector, we have: and substituting Eq. (8) into Eq. (7), we obtain the quadratic equation with respect to A: Here, the coefficients are given as: 0 and arguments θ, θ 0 , φ are omitted for simplification. The positive solution of Eq. (9) is the required estimation of surface reflection. T 11 , μ 0 , E, L are the direct measurements provide by satellite observations. s, ξ, R s 0 are parameters pre-calculated and kept in the LUT.

Aerosol type and Look-Up-Table
In a similar manner to the standard XBAER algorithm, this version of XBAER algorithm uses the Look-Up- Table (LUT) method to derive above-snow-AOD. The LUT is used to calculate the parameters in Eq. (9). The LUT includes the atmospheric reflectance (black underlying surface with surface albedo equal to 0), upward and downward transmittance and spherical albedo, simulated using the radiative transfer package SCIATRAN (Rozanov et al., 2014) for given wavelengths and observation geometries. In order to calculate LUTs, the information of aerosol particle properties, aerosol vertical profiles and corresponding atmosphere status are needed.
One important issue for AOD retrieval in the Arctic is the change of atmospheric reflectance as a function of the elevation, especially over Greenland. Thus, we set the typical geometry combinations of (SZA, VZA, RAA) to (70°, 0°, 30°) (70°, 55°, 30°) for nadir and forward viewing, respectively. Other atmospheric parameters (O 3 , NO 2 , H 2 O, CO 2 , CO, Ch 4 , O 2 , NO) and atmosphere status (pressure, temperate profiles) are obtained from Bremen 2D model for April and 75°N latitude (Aschmann et al., 2009).
We compare the calculation of atmospheric reflectance, transmittance and spherical albedo for sea salt aerosol and dust aerosol for elevation height of 0 km and 3.694 km (highest elevation of Greenland) with AOD at 0.55 μm of 0.5 (largest AOD in LUT), respectively. For sea salt aerosol, the differences of (atmospheric reflectance, transmittance downward, transmittance upward and spherical albedo) are (0.042%, 0.076%, 0.015%, 0.021%) for nadir viewing and (0.006%, 0.076%, 0.0014%, 0.021%) for forward direction viewing. For dust aerosol, the differences of (atmospheric reflectance, transmittance downward, transmittance upward and spherical albedo) are (0.009%, 0.0018%, 0.0009%, 0.0258%,) for nadir viewing (0.024%, 0.0018%, 0.0013%, 0.0258%) for forward direction viewing. This investigation indicates that we can ignore the elevation effect for the retrieval of aerosol using the 3.7 μm. A previous sensitivity study demonstrated in Mei et al. (2017a) shows negligible AOD error from the assumed aerosol height (vertical shape). Consequently, in this work, an "exponential vertical distribution" of aerosol number density with aerosol geometric height equal to 3.0 km is used.

Inverse problem
For the AATSR nadir (n) and forward (f) observations, the surface reflectance at 3.7 μm, A n and A f , are obtained after solving Eq. (9). They depend on a selected AOD for a given aerosol type. Thus, an optimal AOD is found during an iterative process. The latter is stopped when the surface reflectance differences at 3.7 μm between two observation angles reaches predefined criteria.
We then minimize the following equations for dust-dominated and sea salt dominated aerosol types.
where τ ss and τ du are the AOD at 0.55 μm for sea salt dominated type and dust dominated type, respectively. The final output is AOD at 0.55 μm for the selected ε. Eq. (12) is solved using the pre-calculated LUT iteratively employing the Brent's method (Brent, 1973). Fig. 5 shows the logic flowchart for the XBAER algorithm, used to derive AOD over snow surfaces. The algorithm including the pre-processing steps: parallax effect correction, cloud identification. The extraction of solar reflection at 3.7 μm for both nadir and forward observations. The minimalization steps to find the best-fit AOD.

Sensitivity study
Cloud identification, surface parameterization are key issues for any L. Mei, et al. Remote Sensing of Environment 241 (2020) 111731 AOD retrieval algorithm (Mei et al., 2017a). We now investigate the impact of the aerosol typing and surface parametrization in Sections 4.1 and 4.2. In Section 4.3 the impact of cloud contamination is presented.

Impact of aerosol typing
To investigate the impact of aerosol typing in XBAER algorithm, we simulated the AATSR observations (brightness temperature) at 3.7 μm and 11 μm for the aerosol types: dust-dominated and sea salt dominated. The radiative transfer calculations were performed in this case, L. Mei, et al. Remote Sensing of Environment 241 (2020) 111731 including the contribution of thermal emission and reflection of solar light with the following settings: (1) Snow layer was considered as Lambertian reflector and snow emissivity was set to be 0.964 at 3.7 μm according to Spangenberg et al. (2001); (2) SZA, VZA and RAA were set to 70°, (0°, 55°) and 30°, respectively, according to AATSR observation geometry presented in Fig. 1; (3) AOTs were set to be [0.01, 0.05, 0.1, 0.15, 0.2, 0.3, 0.5], which are consistent with the LUT settings; (4) Other settings are the same as given in Section 3. Fig. 6 shows the impact of aerosol typing in XBAER algorithm. In the case of adequate aerosol type (same aerosol type is used in both forward simulation and retrieval process), XBAER algorithm retrieves the abovesnow-AOD successfully for both aerosol types, with maximum error of < 5%. A slight overestimation can be seen for both aerosol types, especially for large AODs and dust-dominated aerosol type. The observed retrieval errors characterize the potential accuracy of the XBAER algorithm and may be explained by: (1) The assumption that the atmospheric thermal emission is negligible in the derivation of Eq. (3); (2) The use of brightness temperature measured at 11 μm instead of the surface temperature in Eq. (4); (3) The assumption of Lambertian reflection of snow layer used in Eq. (5).
We have performed different SCIATRAN simulations to investigate the three issues raised. To assess the impact of thermal emission on AOD (point 1), we switched the aerosol thermal emission off in the forward simulations and repeated the retrieval process. This resulted in the AOD retrieval results having an error of about 1% for both aerosol types. To investigate the impact of using brightness temperate at 11 μm (point 2), we have used the surface temperature based on the temperature vertical profile defined in B2D CTM in the retrieval process, rather than the SCIATRAN simulated TOA brightness temperature. We found that an~5% error in the estimated surface reflectance, is caused by using brightness temperate at 11 μm, depending on AOD and observation geometries. However, this 5% in surface reflectance can be largely mitigated by using the dual-viewing observations to find the "effective surface reflectance" and optimal AOD in iterations. The AOD uncertainty caused by using brightness temperate at 11 μm is < 2%. The impact of potential uncertainty using Eq. (5) (point 3) can be evaluated using a given snow property database (Yang et al., 2013). We have simulated the brightness temperature at both 3.7 and 11 μm for dual-viewing observations with the settings of a snow layer as Fig. 3(a). Positive biases for AOD at 0.55 μm in a range of [0.01, 0.03] can be found depending on AOD. The error introduced by Eq. (5) is caused by the requirement of "effective" Lambertian assumption of surface reflectance, which is not exactly the same as surface BRDF reflectance as we use in the retrieval. Potential correction due to this "inconsistency" has been proposed in Tanré et al. (1979), but is out of this scope for this study.
The green lines in both Fig. 6(a) and (b) indicate the large impacts of using a wrong aerosol type in the retrieval. An overestimation of~2 times occurs when sea salt type is used in retrieval and dust type in forward simulations. And an underestimation of~2 times occurs when dust type is used in retrieval and sea salt type in forward simulations. This is explained by the fact that the same amount of aerosol (same AOD), dust aerosol provides~2 times larger reflectance at the typical AATSR observation geometries, thus for the same atmospheric reflectance, 2 times less AOD is needed to produce the same TOA reflectance (see Fig. 3). Even in a "pristine environment" as the Arctic, the aerosol type is always a mixture of different aerosol components, thus an under/overestimation of AOD can be found in the retrieval of satellite measurements, unless a mixture for type of aerosol if fitted.

Impact of snow surface emissivity
The snow surface emissivity is an important issue in the surface parameterization described in Section 3. To investigate the impact of the snow surface emissivity on the AOD retrieval, we have performed the radiative transfer calculations setting the snow emissivity at 3.7 μm to 0.978, 0.964 and 0.962 and all other settings were used as given in Section 4.1. We note that the minimal and maximal values of snow emissivity are selected according to estimations of Kondo and Yamazawa (1986). Fig. 7 shows the impact of the snow emissivity on the AOD retrieval for both dust and sea salt aerosol types. We can see that the impact of snow emissivity for the typical emissivity values/ranges is ignorable. The three lines in both Fig. 7(a) and (b) are overlapped. This can be explained by the use of Eq. (3) when extracting the solar reflection at 3.7 μm. The dynamical emissivity during the iterative steps enable XBAER to find the best-fit surface and aerosol properties.

Impact of cloud contamination
Although we used a three-steps cloud identification method, potentially sub scene scale cloud and thus cloud contamination may still be present in a scene. On average of 60% of Arctic is covered by cloud (Platnick et al., 2017) and higher probably of cloud contamination is expected, compared to middle-low latitude regions.
We have performed a statistical analysis of the cloud properties over the full Arctic circle (> 60°N) using the MODIS cloud product (version 6.1) for the period 2000 to 2018 (Platnick et al., 2017) and obtained Fig. 6. Impact of aerosol typing on the AOD retrieval accuracy. Reference AOD and XBAER AOD are the input AOD and retrieved AOD. (a) Dust-dominated aerosol type is used in the retrieval; (b) sea salt-dominated aerosol type is used in the retrieval. In (a), Dust-Dust (red line) means Dust-dominated aerosol type is used in both forward and retrieval; Dust-Seasalt (green line) means Sea salt dominated aerosol type is used in forward simulation and Dust-dominated aerosol type is used in retrieval; in (b), Seasalt-Seasalt (red line) means Seasalt-dominated aerosol type is used in both forward and retrieval; Seasalt-Dust (green line) means dust dominated aerosol type is used in forward simulation and Sea saltdominated aerosol type is used in retrieval. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) reasonable cloud properties (COT, CTH, CBH, and effective radius) for the investigation of cloud contamination on above-snow-AOD retrieval in the Arctic. In additional to the cloud settings described in the two section above, we use additional settings on cloud as follow: (1) a Lambertian surface with constant surface albedo is used and snow emissivity is set to be 0.964; (2) the cloud phase is set to be ice and ice cloud particle shape was selected to be the aggregated solid columns (Baum et al., 2014); (3) the cloud optical thicknesses are set to be [0.01, 0.1, 0.25, 0.5]; (4) Effective radius of ice crystals is set to be 35 μm according to the statistical analysis; (5) Cloud top height and bottom height are 6 km and 5 km, respectively. Fig. 8 shows the impacts of cloud contamination on the above-snow-AOD retrieval. Unlike the typical AOD retrieval over middle-low latitude regions where the underlying surfaces are relative dark and the cloud contamination always introduces positive biases in the retrieved AOD (Popp et al., 2016), the impact of cloud on above-snow-AOD retrieval is more complicated, depending mainly on the surface reflectance and the aerosol/cloud properties. Both panels (a) and (b) of Fig. 8 show similar patterns of cloud contamination on above-snow-AOD retrieval for both aerosol types. Potential thin cloud (e.g. cirrus cloud in the Arctic regions), which is a great challenge for any cloud identification method (Mei et al., 2017b), may introduce positive biases for low AODs and negative biases on high AODs. The reason for this is that cloud contamination under relative clear atmosphere (e.g. AOD < 0.1) enhances the brightness temperature at 3.7 μm due to dominated enhancement of the scattering part. For relative hazy case, the absorption of cloud is enhanced because the photon path length increases due to increase of multiple-scattering caused by aerosol-cloud particle interaction. And the turning point from positive bias to negative bias (hereafter as turning point) depends on the aerosol absorption. A lower absorption aerosol layer will enhance the cloud absorption faster, thus a smaller turning point can be expected. Single scattering albedo at 0.55 μm are 0.932, 0.904 for dust and sea salt aerosol types used in the retrieval, respectively. The turning point are 0.07 and 0.15 for dust and sea salt respectively.

Results and evaluation
This section presents the XBAER derived above-snow-AOD and comparisons with in-situ measurements and other existing aerosol (research) products. The structure of this section includes: (1) Compare XBAER-derived above-snow-AOD with AERONET to quantitatively qualify the accuracy of the research product; (2) based on the evaluation of (1), compare the AOD spatial distribution patterns from existing products over a selected regional scale (Greenland) to understand the performance of XBAER; (3) extend the comparison of (2) to the full Arctic region; (4) with a comprehensive understanding of the performance/accuracy of XBAER derived above-snow-AOD research product based on the output of (1), (2), and (3), a detailed cased study is presented to illustrate potential application in next section. A direct comparison with AERONET observation shows the overall quality of the above-snow-AOD research product, derived from XBAER. This is shown Fig. 7. Impact of surface emissivity on the AOD retrieval accuracy. (a) Dust-dominated aerosol type is used in the retrieval; (b) sea salt-dominated aerosol type is used in the retrieval. Three difference emissivity (0.978 (red line), 0.964 (green line) and 0.962 (blue line)) are used for forward simulations. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 8. Sensitivity study for cloud contamination. (a) Dust-dominated aerosol type is used in the retrieval; (b) sea salt-dominated aerosol type is used in the retrieval. The typical COT values (0.01, 0.1, 0.25, 0.5) are used in the forward simulations and cloud-free condition in the retrieval. L. Mei, et al. Remote Sensing of Environment 241 (2020) 111731 in Fig. 9. Fig. 9 shows the comparison for AOD at 0.55 μm. Only coarsemode dominated conditions base on AERONET fine mode fraction are used for the comparison. The comparison collocates the full AATSR mission (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) with four existing AERONET sites over Greenland (Thule, Ittoqqortoormiit, Kangerlussuaq, Narsarsuaq). The observation periods for those sites are April 2008-present, October 2009-present, March 2009-present and March 2013. Due to (1) polar night and (2) large cloud cover, the AERONET observations over the Arctic regions (e.g. Greenland) have fewer measurements compared to middle-low latitude regions. 369 (N as shown in Fig. 9) collocations between AATSR and AERONET were found for the whole AATSR mission. The color of each grid (0.01 × 0.01 increment) represents the number of match-ups. The AOD values are required to be between 0.01 and 0.5. AOD values outside of this region are excluded in the comparison. Taking the small AOD values and AOD variability into account, an expected error (EE) envelope (Levy et al., 2013) is defined as ( ± 0.15AOD ± 0.025). Although there are intensive discussions about the use of linear regression for the validation of satellite aerosol products in the aerosol satellite community (e.g. in the AeroCom/AeroSAT meetings), no other better way is currently widely accepted by the community. This work still follows the classical validation strategy, even though this validation method may not really show the capability of the retrieval algorithm, especially over the Arctic, where AODs have very small values inside a small variability range. The regression equation is Y (x) = (0.764 ± 0.055) x + (0.029 ± 0.004) with a correlation coefficient of 0.615. About 72% of the match-ups fall into the predefined EE and a clear overestimation of XBAER derived abovesnow-AOD can be found (22% match-ups are above EE). Unlike the classical aerosol retrieval over dark-moderately bright surface (e.g. standard XBAER algorithm as presented in Mei et al. (2017a)), the overestimation of AOD over snow can't be explained simply by cloud contamination as clouds may introduce underestimation as presented in Section 3. Even though clouds still introduce positive biases for the relatively low AOD cases, the absolute biases are typically < 0.05 for AOD smaller than 0.1. A double overestimation may also be due to an improper aerosol typing or possible change of surface conditions. The overestimations frequently occur during snow-melting season or heavy snow-fall periods, during which the surface properties changed rapidly. However, currently we cannot find a physical-robust way to evaluate the error introduced by the surface parameterization (e.g. the anisotropic properties at the dual-viewing observations). The majority of the match-ups (68.3%) are fall into low aerosol loadings (AOD < 0.1). The comparison between XBAER above-snow-AOD and AERONET observations shows a promising quality of the research product.
The validation using AERONET sites shows the quality of XBAER for the limited regions around the AERONET sites. In order to have an overall comparison of XBAER derived above-snow-AOD research product for a larger spatial coverage, Fig. 10 shows the monthly AOD of L. Mei, et al. Remote Sensing of Environment 241 (2020) 111731 April 2011 for CALIOP , MERRA (Rienecker et al., 2011), IASI (MAPIR) (Callewaert et al., 2019) and AATSR (XBAER). In order to highlight the spatial distribution patterns, the MERRA simulations have been increased by a factor of 1.2. Due to the different temporal and spatial resolutions of these four datasets, we can expect similar but different patterns. And the differences between AOD products can sometimes be the dominant features. Fig. 10 shows that the four datasets show similarities in general. Large differences in AOD absolute values can be found especially between MERRA and other three AOD products. Greenland shows a majority AOD value of < 0.05, especially for the elevated central Greenland. The coastline regions show higher AOD values. CALIOP, IASI (MAPIR) and AATSR (XBAER) have similar AOD values around coastline regions, with values of about 0.15-0.2, MERRA shows lower AODs. Hardenberg et al. (2012) showed that the AOD from model simulations are significantly underestimated due to lack of aerosol sources in the Arctic regions. IASI (MAPIR) and AATSR (XBAER) show a good agreement for both the spatial patterns and AOD magnitude on a monthly scale over Greenland. Above analyses demonstrate the considerable potential of the XBAER algorithm to retrieve above-snow-AOD for the full Arctic circle (60°-90°N). Tomasi et al. (2015) presented two representative days (29 March and 3 May 2006) for the Arctic clear day and haze day. Fig. 11 shows the AATSR derived above-snow-AOD. Fig. 11(a) and (b) is the Swansea aerosol product produced in the European Space Agency's Climate Change Initiative (CCI) project, which has been evaluated to be the best AATSR AOD product (Popp et al., 2016). According to Fig. 11(a) and (b), only very limited "open water" regions can be used to retrieve AOD and almost no retrieval can be found for the whole Arctic land region due to large cloud/snow/sea ice coverage. Fig. 11(c) and (d) shows the above-snow-AOD derived by XBAER using AATSR. The AATSR AOD coverage has been largely extended by XBAER compared to the Swansea retrieval (Popp et al., 2016) and the AODs for 3 May are in general larger compared to 29 March, which shows good spatial consistency with the AOD figures presented in Tomasi et al. (2015) using the method proposed in Mei et al. (2013a). Large AODs for 3 May can be seen especially over north Russia and Canada. The pollution over Europe and Russia has been observed to transport to Spitsbergen (Treffeisen et al., 2007;Stone et al., 2014;Shi et al., 2019). Although very limited regions are covered by the AATSR AODs provided by the Swansea retrieval, a smooth transition between Swansea retrieved AOD over open water and XBAER above-snow-AOD can be seen, indicating the promising quality of above-snow-AOD derived by XBAER.

Case study
In this section, we applied the XBAER derived above-snow-AOD for an aerosol transport scenario to check if the AOD research product can capture the transport event. A recent publication (Francis et al., 2018) proposed a new mechanism of long-distance transport of Saharan dust reaching the Arctic. A dust transport event on April 2011 has been analyzed in detail and the dust transport was associated with a Saharan cyclone formed over southern Morocco due to the intrusion into subtropics in early period of April 2011 (Francis et al., 2018). The intense cyclone created a strong near surface wind with a speed of 20 m·s −1 . This blew the Sahara dust poleward (Francis et al., 2018). Fig. 12 shows a time series of the dust event between the 6th and 9th of April 2011, retrieved from observations of satellite instrumentation and model simulations. Fig. 12  because the MERRA simulations have assimilated the MODIS aerosol product (Rienecker et al., 2011). This agreement between MODIS and MERRA simulations shows that MERRA can be used to capture and quantify this dust event. Since the MODIS aerosol product does not cover cloud/snow covered regions, the plume was not observed on the 9th of April, but it has been modelled by MERRA. On the 9th of April, the MERRA simulations enable us to track the trajectory of the dust plume. The dust plume reached East Greenland, where its AOD is~0.3 ( Fig. 12(l)). Fig. 13 shows the dust event observed over Greenland. Limited swath with of AATSR limits the coverage over Greenland. Days when AATSR has coverage over East Greenland and is thus able to capture the transported dust plume described in Fig. 12 were selected. Four days (3,9,12,20 April) were chosen and the MERAA simulations with the shortest time interval with respect to AATSR overpass time are used. We also show the observations from the MAPIR derived IASI dust aerosol research product for those four days as support information to understand this dust event. All three data sources show good agreement for the background aerosol conditions over Greenland, with AOD of < 0.05. Both XBAER and MAPIR provide AOD under cloud free conditions, thus "data gaps" are presented in Fig. 13. For the 3rd of April, most regions over Greenland show low AOD, except the west coast. Both MERRA and MAPIR observed the "high" aerosol pattern over the west coast. However, AATSR did not make measurements in these regions on the 3rd of April. On the 9th of April, the Saharan plume reached east Greenland and its AOD is approximately 0.3. MERRA and XBAER AODs are in good agreement, while MAPIR AOD only observed that part of the plume which covered the southern part of Greenland. The local time differences of the three (research) products does not explain the differences of the spatial AOD patterns. Both XBAER and MAPIR AOD have "noisy" features, in particular MAPIR. We attribute this to potential (broken) cloud and related contamination for both research products. The "red points", which are not smooth high AOD, in the MAPIR retrievals are dubious although they do pass the quality filters, potential cloud contamination may still remain. The sensitivity for MAPIR retrieval is low with low AODs on cold surfaces, which is responsible for the noisy results in addition to the plausible presence of remaining clouds. And in XBAER, the effect of Sastrugi on snow may also contribute to the non-smooth AOD features (Leroux and Fily, 1998). The two "bluish" patterns over east part of Greenland observed by XBAER are likely due to the cloud contamination, which is consistent with the sensitive studies presented above. For 12 and 20 April, both MERRA and XBAER show limited aerosol information over Greenland while MAPIR still catches some aerosol information over the southern part of Greenland. The different patterns between MAPIR and MERRA/ XBAER may due to the low sensitivity of MAPIR, the a-priori (EUME-TSAT IASI level 2 retrieval) of surface temperature (Ts) shows unreasonably high Ts in the areas where MAPIR retrieves dust. Good agreement between MAPIR and MERRA over Northeast Greenland is observed on 20 April, where unfortunately no XBAER data is available due to cloud cover. The over-lapped regions between XBAER and MAPIR show acceptable agreement, especially over the coast line, where MERRA shows underestimation as above. The above analysis shows that XBAER derived above-snow-AOD research product over snow is capable to catch aerosol in the Arctic snow-covered regions.

Applying XBAER to SLSTR measurements
The XBAER algorithm has been applied to the data from the SLSTR instrument on board Sentinel-3AFor this study, we chose observations over the Kamchatka regions in eastern Russia, where a large group of volcanoes are situated and most regions are covered by snow. Volcanic eruptions in this region produce large amounts of ash typically with height of 3-4 km and the ash plumes travel southeast, south, and southwest (Girina et al., 2018). A volcanic eruption (56°39′12″N 161°21′42″E) occurred on 19th of February 2019. Fig. 14(left) shows  (a, b, c, d) shows the MODIS RGB composition for 6-9 April 2011 (red lines with arrow shows the wind speed (at level 60 in ECMWF dataset) and wind direction: the length of the lines and the arrows show wind speed and direction, respectively); Middle plane (e, f, g, h) shows the corresponding MODIS collection 6.1 AOD (no information for cloud and snow condition) for 6-9 April 2011; Lower plane (i, j, k, l) shows the MERRA simulated AOD for 6-9 April 2011. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) the RGB composition figure from MODIS. The reason to use MODIS RGB rather than SLSTR RGB figures is due to limited channels of SLSTR compared to MODIS, MODIS provides a better and direct figure of "true color" of the world. Fig. 14(left) shows that this region is mostly covered by snow. Fig. 14(right) shows the XBAER derived AOD at 0.55 μm on the 19th of February 2019. Relatively large AOD values are observed around the volcano, in particular over its southern slope. The AOD at 0.55 μm of is about 0.3. A fresh ash plume is seen travelling in the northeast direction, which is well-captured in SLSTR. This is the yellowish pixel northeast of the volcano shown in Fig. 14(right).
The above demonstrates the successful use of XBAER to derive above-snow-AOD using SLSTR measurements. The features of SLSTR AOD appear "noisier" than that retrieved from AATSR data. This may be explained by the following two reasons: (1) the SLSTR instrument provides dual-viewing observations for nadir and backward (oblique) directions (not forward as AATSR), which provides weaker information to distinguish aerosols from underlying snow surfaces; (2) the calibration for SLSTR associated with a new instrument is not final. Potential calibration issues may also contribute to the "noisy" features observed here.

Conclusion
A new retrieval algorithm to retrieve AOD over snow covered regions using passive satellite instruments is presented. The new AOD retrieval algorithm is an optimized version of the generic XBAER algorithm. It exploits the unique information, resulting from the dualviewing observation capability of ESA's AATSR and SLSTR instruments. The XBAER uses the differences between the anisotropic properties of solar reflection of aerosol and snow surface at 3.7 μm, in which the Fig. 13. AOD i) top panel (a-d) simulated using MERRA, ii) middle panel (e-h) retrieved using XBAER from AATSR observations and iii) bottom panel (i-l) derived using MAPIR from MODIS for 3, 9, 12 and 20 April for the dust event shown in Fig. 12.
Fig. 14. XBAER derived AOD over snow on 2019-02-18 when Kamchatka volcano (red star) erupted. using SLSTR onboard Sentinel 3A. For a better illustration, we use MODIS RGB (a) and show corresponding XBAER derived AOD over snow. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) solar reflection contribution from snow is weak. The XBAER algorithm yields above-snow-AOD for coarse-mode dominated aerosol types.
The "parallax effect" of AATSR/SLSTR instruments is accounted for, before retrieval. XBAER uses a cloud identification algorithm from Istomina et al. (2010) and Mei et al. (2017b) in a pre-processing step. A new surface parameterization has been used in the XBAER algorithm to decouple the solar reflection from thermal emission at 3.7 μm. A sensitivity study taking aerosol typing, surface impacts and potential cloud contamination, into account has been performed. The main findings are as follows: (1) The XBAER algorithm when applied to AATSR/SLSTR observations yields above-snow-AOD for the appropriate aerosol type and minimized cloud contamination; (2) an inappropriate aerosol type results in significant errors in retrieved above-snow-AOD in XBAER algorithm; (3) inaccuracy in the knowledge of surface properties (e.g. surface emissivity) play a minor role in the retrieval; (4) cloud contamination introduces both positive and negative biases depending on the aerosol absorbing and cloud properties.
The validation using AERONET shows good agreement between XBAER derived above-snow-AOD and AERONET observations, with about 72% match-ups that fall into the EE of ( ± 0.15AOD ± 0.025). The monthly mean AODs over Greenland show good agreement of the spatial distribution patterns with other satellite-derived and model-simulated products such as CALIOP, MERRA and IASI. The XBAER-derived above-snow-AOD research product has significantly better coverage than the current existing AATSR AOD product and no obvious snow-sea contrast has been observed. In two case studies, the new research product identifies the dust event, which occurred on the 9th of April 2011 and the volcanic eruption on the 19th of February 2019.
The current XBAER algorithm uses combined cloud identification algorithms. However, potential cloud contamination may still exist, which may introduce positive or negative biases depending on the aerosol and cloud properties. The criteria used in the current cloud identification method will be further linked to the cloud particle size (cloud effective radius) based on existing databases (Yang et al., 2013) and other theory (e. g. asymptotic radiative transfer) (Kokhanovsky et al., 2018). Thin clouds, such as cirrus, are an important source of cloud contamination in the XBAER approach. The ice crystals in a cirrus cloud are much smaller than the snow grain size of the underlying snow surface, some simplified analytical theories to derive an acceptable ice cloud effective radius (Mei et al., 2018b) and snow grain size will be helpful to distinguish the thin clouds from snow-covered surface. The updated cloud mask may reduce the potential cloud contamination in the XBAER algorithm.
The surface parameterization coupled with the path radiance representation (Eq. (5)) may introduce some uncertainties in the abovesnow-AOD from inaccurate BRDF treatments. The path radiance representation requires an "effective Lambertian albedo" and the atmospheric correction (Eq. (9)) gives a direction reflectance for given geometries. Our investigation shows a small positive difference between the "effective Lambertian albedo" and the direction reflectance for nadir observation, and a small negative difference is found for forward observation using AATSR instrument. Although a reasonable AOD can still be retrieved during the iteration steps by minimizing the positive and negative differences to approach a zero difference, a better way to taking these differences into account is still needed. In the next version of XBAER algorithm, we plan to link the surface bidirectional reflectance to an angular independent parameter: snow grain size. The above-snow-AOD and snow grain size will be retrieved simultaneously.
Our investigation of the impacts of aerosol typing on above-snow-AOD retrieval shows the current "weakly absorbing" aerosol type (Levy et al., 2013;Mei et al., 2017a;Mei et al., 2019b) is an optimal aerosol parametrization for the Arctic regions. This aerosol parametrization enables a "single" type to represent the whole Arctic due to the dynamical relationship between aerosol optical properties to AOD itself. The XBAER algorithm may also be applicable to the above-sea-ice-AOD, using a similar idea. This requires an investigation of the angular and spectral BRDF features of snow and sea-ice (Clémence et al., 2015). The XBAER has been demonstrated to deliver AOD for coarse dominated aerosol types and data from instruments such as AATSR/SLSTR. In the future, the AOD retrieved from AATSR/SLSTR instruments for all types of particles will be investigated.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.