The Determination of Snow Albedo from Satellite Measurements Using Fast Atmospheric Correction Technique

: We present a simpliﬁed atmospheric correction algorithm for snow / ice albedo retrievals using single view satellite measurements. The validation of the technique is performed using Ocean and Land Colour Instrument (OLCI) on board Copernicus Sentinel-3 satellite and ground spectral or broadband albedo measurements from locations on the Greenland ice sheet and in the French Alps. Through comparison with independent ground observations, the technique is shown to perform accurately in a range of conditions from a 2100 m elevation mid-latitude location in the French Alps to a network of 15 locations across a 2390 m elevation range in seven regions across the Greenland ice sheet. Retrieved broadband albedo is accurate within 5% over a wide (0.5) broadband albedo range of the (N = 4155) Greenland observations and with no apparent bias.


Introduction
There is a decreasing trend in both the extent and the reflective power of the terrestrial cryosphere with important climate change feedbacks [1][2][3][4]. The solar light reflectance from snow and ice has a bi-directional character, depending on the direction of illumination and on the observation direction, which can be measured using ground, airborne, and satellite optical instruments. Climate models utilize snow spectral plane albedo, which provides total reflected solar light power for a given wavelength and a given solar incidence angle that depends on location and time. Satellite measurements, of particular importance for studies of polar environment [4], are usually performed with a fixed observation geometry. Therefore, special procedures are needed to convert satellite-measured reflectance to a plane albedo [5,6]. The broadband plane albedo (BBA) can be derived using various parameterizations or by integration of the spectral plane albedo with account for the spectral snow irradiance at the snow surface [7]. The optical signals measured by satellite are influenced not just by light reflected from the surface but by atmospheric extinction, scattering and absorption. Therefore, atmospheric effects need to be removed to obtain accurate surface retrievals. Similarly, remote sensing of the atmosphere requires removal of the surface contribution to the observed signal. Atmospheric remote sensing is more readily made in cases of dark underlying surfaces, such as the ocean. In polar regions, the retrievals of atmospheric aerosol load over bright snow and ice surfaces are challenging and often hardly possible because the signal is dominated by the bright surface and not by atmospheric aerosol.
A key task of this study is to provide an accurate determination of spectral and broadband plane albedo of snow and ice using satellite observations amid the challenge of atmospheric absorption by ozone, molecular light scattering and light scattering, and absorption by atmospheric aerosols. It is assumed that aerosol optical properties are known a priori, i.e., from aerosol climatology, forecasts, or ground measurements for the case of polluted snow/atmosphere. In the case of clean snow and atmosphere, we do not rely on any a priori information on atmospheric aerosol loading and properties in our snow and ice albedo retrieval technique. In any case, the generally low polar aerosol loading [4] reduces the influence of the aerosol contribution to the retrieved surface albedo. In the case of polluted snow, retrievals are performed outside strong atmospheric absorption bands (e.g., O 2 and H 2 O). The ozone absorption effects are fully accounted for in the retrieval framework. While the algorithm is easily portable to other multi-spectral instruments observing the cryosphere from space, we present an application to data from the Ocean and Land Colour Instrument (OLCI) on board the European Union Copernicus Sentinel-3A satellite. The theoretical modeling of spectral snow reflectance is performed as in [6]. The earlier atmospheric correction used in [6], which appears in OLCI Snow Properties module incorporated in the European Space Agency (ESA) SeNtinel Application Platform (SNAP), can be biased in case of strong atmospheric pollution episodes (arctic haze, etc.) because it neglects scattering and absorption by liquid and solid particles suspended in atmosphere. This shortcoming of the previous algorithm as presented in [6] is eliminated in this study.

Theory
We perform the retrievals using separate retrieval chains for clean snow (Case 1 snow) and for polluted snow (Case 2 snow). Here, we use the analogy with the classification of Case 1 and Case 2 water as proposed in [8] (see also [9]). Case 1 water corresponds to relatively clean water where most of the absorption is due to phytoplankton, and Case 2 water contains other impurities including mineral particles. In our application, we define Case 1 as the situation where snow properties are determined just by snow grains without significant interference from impurities or living matter (cells, algae, etc.). The snow Case 1 is often met in Antarctica-far from any significant aerosol sources and limited algal populations. The areal extent of the clean dry snow areas on Greenland and Antarctica ice sheets makes the Case 1 snow dominant on a global scale. Additionally, a simplified atmospheric correction is possible in this case [6]. The selection of clean snow pixels is performed as follows. First, we check the reflectance in OLCI band 1. If it is larger than the dynamic threshold value (THV), it is assumed that the ground scene is covered by unpolluted snow (the majority of pixels in the terrestrial cryosphere). The THV is derived from the synthetic radiative transfer calculations for the assumed (default: 0.1) aerosol optical thickness at 550 nm (see Appendix A).

Case 1 snow
The simplified atmospheric correction for Case 1 snow is described in [6] and summarized below. It is based on the fact that the pure snow spherical albedo can be accurately parameterized using the following equation: where α(λ) = 4πχ/λ is the bulk ice absorption coefficient for a given wavelength λ, χ (see e.g., https://atmos.washington.edu/ice_optical_constants/, last access: 07/01/2020) is the imaginary part of ice refractive index, and l is the effective absorption length. The snow spectral reflectance R s is related to the snow spherical albedo, which is a three dimensional integral of R s with respect to solar Remote Sens. 2020, 12, 234 3 of 18 and viewing zenith angles and relative azimuthal angle (RAA) [6] via the following approximate equation [6]: where x is a geometrical correction coefficient depending on R 0 and on the angular function u [6] evaluated at the cosine of the solar zenith angle (SZA) µ 0 or at the cosine of the viewing zenith angle (VZA) µ: and we use the following approximation for the angular function [6]: The value of R 0 gives the non-absorbing underlying surface reflectance (r s = 1). One can use OLCI measurements at 865 and 1020 nm to determine both effective absorption length and R 0 from Equation (2) under the assumption that the atmosphere does not affect the satellite signal at these channels [6].
The determined value of effective absorption length makes it possible to derive the snow spherical albedo at any wavelength using Equation (1). The plane albedo r p is defined via the integral of the azimuthally averaged reflection function with respect to the viewing zenith angle [6]. As a matter of fact, r p can be also derived from the spherical albedo using the following simple approximation [6]: or with account for Equation (1): Also, one can derive the underlying snow spectral reflectance function using Equations (1) and (2). Therefore, the procedure for the determination of Case 1 spectral albedo from space is straightforward. It was validated in [6]. Generally, the errors in the retrieved albedo are below 1-3% depending on the wavelength λ.

Case 2 snow
The retrievals for the Case 2 snow are more complicated. In this case, the satellite measurements of snow spectral reflectance in the visible are influenced by various pollutants or living matter (cells, algae, etc.). Therefore, there is no way to estimate snow spectral reflectance/albedo in the visible using measurements in the near infrared as it is done for the Case 1 snow (see above). Then, we use yet another approach described below.
The top-of-atmosphere reflectance for the atmosphere-underlying snow system can be presented in the following way [10,11]: where R ag is the atmospheric contribution to the measured signal, r ag is the spherical albedo of the atmosphere, r s is the bottom-of-atmosphere snow spherical albedo, and T ag is atmospheric transmittance from the top-of-atmosphere to the underlying surface and back to the satellite position. In the case of Lambertian underlying surfaces, the underlying surface reflectance does not depend on solar and viewing observation directions, and Equation (7) is valid with r s = R s , where R s is the underlying Lambertian surface reflectance. The snow is not exactly the Lambertian reflector; therefore, we replace r s in the numerator of Equation (7) by the snow reflectance [see Equation (2)]. Such an approximation makes it possible to have the correct limit for the top-of-atmosphere reflectance [see Equation (2)] in the case of absence of atmosphere. The term in the dominator of Equation (7)  reflections between snow and atmosphere, and the account for the snow reflectance directional nature in the dominator of this equation is of secondary importance. Then, it follows: The reflectance of non-absorbing snow R 0 in Equation (8) is calculated using simple analytical approximation, as discussed in [6]. We do not derive the value of R 0 from OLCI measurements themselves because such a derivation for the polluted snow can be influenced by the type and the load of pollutants.
We use channels that are not influenced by water vapor and oxygen absorption effects, although we account for the ozone absorption effects. Equation (8) is very general and valid outside and inside molecular absorption bands. We account for the ozone absorption in a simplified way. Namely, we derive free of ozone absorption top-of-atmosphere reflectance R c using the following equation: , where T O3 is the atmospheric transmittance with account for the ozone absorption (see Appendix A). Then, Equation (8) is transformed to a simplified approximation: where the functions R a , r a , T a (see Appendix A) have the same meaning as R ag r ag , T ag , respectively, except for atmosphere not influenced by gaseous absorption processes (e.g., ozone absorption). The spherical albedo of underlying snow surface can be found from Equation (9) provided that the aerosol model is known. In this case, the snow spherical albedo r s is the only unknown parameter in Equation (9) and can be readily calculated, solving the transcendent Equation (9) with respect to r s . For the wavelengths where the aerosol contribution is low and can be neglected, R a ∼ 0, r a ∼ 0, T a ∼ 1, and an analytical solution of Equation (9) is possible: where the analytical expression for R 0 is given in [6]. The functions R a , T, and r a depend on aerosol and molecular scattering parameters and can be stored in look-up-tables for various aerosol models. Because aerosol load is weak in the Arctic and Antarctica, various approximations for the functions mentioned above can be used. In particular, we calculate these functions in the framework of approximations described in the Appendix A. We solve the transcendent Equation (9) with respect to r s for all OLCI wavelengths free of water vapor and oxygen absorption in the Case 2 snow. The broadband albedo (BBA), either plane or spherical, is calculated from the spectral plane or the spherical albedo using the integration between the wavelengths λ a and λ b as shown below [7]: where F(λ) is the incident solar flux at the snow surface, and r p,s (λ) is plane (p) or spherical (s) albedo depending on whether plane or spherical BBA r p,s (λ a , λ b ) is to be calculated. The indices a and b signify the wavelengths λ used. We assume that the incident solar flux can be approximated by the following analytical function: where we ignore rapid oscillations of F(λ), which are due to gaseous absorbers. This is possible because r p,s (λ) is a continuous function, which acts as a filter of high frequencies. The coefficients in Equation (12) are derived from the fit of F(λ) calculated using the Santa Barbara DISORT Radiave Remote Sens. 2020, 12, 234

of 18
Transfer (SBDART) code [12] to Equation (12) in the spectral range 0.3-2.4 µm. They are given in Table 1. The calculations of F(λ) are performed at the parameters listed in Table 2 for the rural aerosol model [13]. Clearly, F(λ) depends on the location and the time. We find that the choice of aerosol model in the calculation of F(λ) only weakly influences the calculations of BBA (see Equation (11)). The spectral snow albedo needed as input for SBDART is calculated assuming clean snow with the effective diameter of spherical ice grains equal to 0.25 mm. Generally, the results are only weakly sensitive to the variation of the function F(λ) [7]. We therefore assume solar flux independent from the location of the retrieval and from solar zenith angles. Table 1. The coefficients of approximation given by Equation (12).  For the Case 1 snow, the broadband albedo is calculated numerically using Equations (1), (5), (11), and (12) in the spectral range 0.3-2.4 micrometers. Also, other limits of integration can be used (say, to derive visible or near-infrared BBA).
For the Case 2 snow, the spherical albedo is known only for selected OLCI channels as derived from Equation (9). Therefore, we use interpolation to get the spherical albedo between the measurement points needed for the evaluation of integral (11). For the spectral range below 865 nm, we use: While, for wavelengths larger than 865 nm, we use: We use the dependencies as shown in Equations (13) and (14) because we find that the measurements can be approximated by the second order polynomial for the spectral range below 865 nm and the exponential function for the wavelengths above 865 nm. The coefficients (a, b, c) are found separately for the intervals 400-709 nm and 709-865 nm using the following wavelength triplets: (400, 560, 709 nm) and (709, 753, and 865 nm), respectively.
The coefficients ( , σ) are derived from OLCI measurements at 865 and 1020 nm at the values of R meas (1020 nm) equal to or smaller than 0.5. Otherwise, Equation (1) [and not Equation (14)] is used at λ > 865 nm with the effective absorption length derived from the value of spherical albedo at 1020 nm. We use different approaches for the pixels with small and large values of R meas (1020 nm) because the case of comparatively large values of R meas (1020 nm) corresponds to snow. Otherwise, ice or extremely dirty snow is present. Then, Equation (1) is not valid.
Integral (11) for the spherical broadband albedo with account for Equations (12)-(14) can be evaluated analytically. The answer is: Remote Sens. 2020, 12, 234 Here, the coefficients a j , b j , c j are the same as presented in Equation (13) with j = 1 for the first spectral interval (0.3-0.709 microns) and j = 2 for the second spectral interval (0.709-0.865 microns).
is the integral given in the dominator in Equation (11) [evaluated analytically with account for Equation (12)] and: At the R (1020 nm) equal to or above 0.5, the analytical expression for the BBA cannot be derived (because one accounts for Equation (1) and not Equation (13) in Equation (11)). Then, the numerical integration procedure is followed.
The broadband plane albedo is calculated in a similar way as a broadband spherical albedo using Equation (5) for the transformation of spherical to plane albedo.
This concludes the description of this new fast radiative transfer Snow and ICE surface albedo retrieval (SICE) that accounts for atmospheric scattering and absorption effects. The SICE algorithm can be considered as an update of the previous version of the algorithm (called S3Snow [6]) that appeared in the Snow Properties module of SNAP.

Snow Spectral Albedo
The validation of spectral albedo for the case of clean snow is reported in [6], where a detailed description of ground and satellite measurements, not repeated here, may be found. In the case of polluted snow, we follow a different procedure than that suggested in [6]. Here, we use the improved atmospheric correction, which explicitly accounts for molecular/aerosol light scattering and absorption effects. The results for the French Alps, Col du Lautaret validation site (45.041288N, 6.410557E, 2100 m a.s.l.) on 17 April, 2018 (Figure 1a) confirm that the current SICE planar albedo retrieval has a higher accuracy as compared to the earlier S3Snow algorithm [6] for the cases studied. Also, unlike the S3Snow retrievals, there is a possibility to vary aerosol load in the framework of the updated retrieval, which is currently not the case for S3Snow plane albedo retrieval results. As it follows from Figure 1b, the variation of aerosol optical thickness (AOT) (500 nm) in the range 0.07-0.35 does not change the plane albedo retrieval accuracy considerably (above 3% for the case studied). Note that the AOT (500 nm) obtained from the Copernicus Atmospheric Monitoring Service near-real-time forecast product (https://www.ecmwf.int/en/about/what-we-do/environmental-services/copernicusatmosphere-monitoring-service, last access: 07/01/2020) is 0.125 for the case studied. We conclude that Remote Sens. 2020, 12, 234 7 of 18 the precise information on the spectral aerosol optical thickness is not needed for the accurate snow spectral albedo retrievals for the low aerosol load characteristic for the Col du Lautaret validation site located at 2,100 m a.s.l. [4]. It should be pointed out that the discrepancy of satellite and ground plane albedo measurements is not solely due to the retrieval but also partially due to surface spatial inhomogeneity (local scale effects vs. much broader scale satellite pixel effects), time difference between satellite and ground measurements, and influence of 3-D effects from surrounding mountains.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 18 characteristic for the Col du Lautaret validation site located at 2,100 m a.s.l. [4]. It should be pointed out that the discrepancy of satellite and ground plane albedo measurements is not solely due to the retrieval but also partially due to surface spatial inhomogeneity (local scale effects vs. much broader scale satellite pixel effects), time difference between satellite and ground measurements, and influence of 3-D effects from surrounding mountains.

Snow Broadband Albedo
We validate broadband albedo measured in the 0.3-2.4 micron wavelength range using ground measurements from fifteen Programme for the Monitoring of the Greenland Ice Sheet (PROMICE) automatic weather stations, with albedo data already described with more detail in [6]. In comparisons presented here, the closest hourly observations are considered. Occasionally for northern sites (KPC_L, KPC_U), there are multiple comparisons each day. Clear sky conditions are estimated from downward longwave irradiance data after [6]. A total of 4,146 individual comparisons are made. The PROMICE BBA data include a correction for measurement platform obstruction of the radiometer field of view [14] that increases average PROMICE albedo by 0.034.
At the PROMICE SCO_U location (Figures 2 and 3, Table 3) on the eastern ice sheet, where clear sky conditions are common, we observe in three May-September years (2017, 2018, 2019) over a wide (0.52) BBA range (from 0.33 indicative of impurity rich bare ice to 0.85 indicative of dry clean snow cover) a very similar temporal pattern in both the ground observations and from the OLCI retrievals.

Snow Broadband Albedo
We validate broadband albedo measured in the 0.3-2.4 micron wavelength range using ground measurements from fifteen Programme for the Monitoring of the Greenland Ice Sheet (PROMICE) automatic weather stations, with albedo data already described with more detail in [6]. In comparisons presented here, the closest hourly observations are considered. Occasionally for northern sites (KPC_l, KPC_U), there are multiple comparisons each day. Clear sky conditions are estimated from downward longwave irradiance data after [6]. A total of 4,146 individual comparisons are made. The PROMICE BBA data include a correction for measurement platform obstruction of the radiometer field of view [14] that increases average PROMICE albedo by 0.034.
At the PROMICE SCO_U location (Figures 2 and 3, Table 3) on the eastern ice sheet, where clear sky conditions are common, we observe in three May-September years (2017, 2018, 2019) over a wide (0.52) BBA range (from 0.33 indicative of impurity rich bare ice to 0.85 indicative of dry clean snow cover) a very similar temporal pattern is found in both the ground observations and from the OLCI retrievals.    (Table 3) in comparison to the BBA retrievals from the current retrieval (SICE), the earlier Ocean and Land Colour Instrument (OLCI) retrieval (S3SNOW) after [6] and for the National Aeronautics and Space Administration (NASA) Moderate Resolution Imaging Spectroradiometer (MODIS) MOD10A1 product [15]. The symbols have the following meanings: PROMICE (crosses), SICE (circles), S3SNOW (rhombus), MOD10A1 (triangles). The several outliers for the SICE product (too high BBA) are related to the problem with our automatic cloud detection procedure.
Among the fifteen PROMICE locations spanning a wide spatial scale [2076 km north-south (18.9° latitude) and 2390 m in elevation], SICE BBA agreement is as high as is realistic to expect with unattended automatic weather station (AWS) observations even though we filter the sample to select cases when AWS tilt recordings are under 1 degree and downward longwave-derived cloud index is under 0.3 [6] (Table 3). It is very encouraging to find: regression slopes averaging insignificantly from unity; an average multi-site correlation coefficient of 0.869 and an average rootmean-square difference (RMSD) of 0.056. The relatively low correlation at the EGP site is more due to the relatively small (~0.1) seasonal fluctuation than a higher error at that site. In idealized circumstances, i.e., with tilt well under 1 degree and evidence of absolutely clear sky conditions from smooth downward shortwave diurnal curves, agreement on a case by case basis can be better than 0.02 (not shown), but such high agreement is not to be expected from a large sample of automatically-selected cases.  (Table 3) in comparison to the BBA retrievals from the current retrieval (SICE), the earlier Ocean and Land Colour Instrument (OLCI) retrieval (S3SNOW) after [6] and for the National Aeronautics and Space Administration (NASA) Moderate Resolution Imaging Spectroradiometer (MODIS) MOD10A1 product [15]. The symbols have the following meanings: PROMICE (crosses), SICE (circles), S3SNOW (rhombus), MOD10A1 (triangles). The several outliers for the SICE product (too high BBA) are related to the problem with our automatic cloud detection procedure.
Among the fifteen PROMICE locations spanning a wide spatial scale [2076 km north-south (18.9 • latitude) and 2390 m in elevation], SICE BBA agreement is as high as is realistic to expect with unattended automatic weather station (AWS) observations even though we filter the sample to select cases when AWS tilt recordings are under 1 degree and downward longwave-derived cloud index is under 0.3 [6] (Table 3). It is very encouraging to find: regression slopes averaging insignificantly from unity; an average multi-site correlation coefficient of 0.869 and an average root-mean-square difference (RMSD) of 0.056. The relatively low correlation at the EGP site is more due to the relatively small (~0.1) seasonal fluctuation than a higher error at that site. In idealized circumstances, i.e., with tilt well under 1 degree and evidence of absolutely clear sky conditions from smooth downward shortwave diurnal curves, agreement on a case by case basis can be better than 0.02 (not shown), but such high agreement is not to be expected from a large sample of automatically-selected cases. The current approach, advanced from that reported in [6] by the improved atmospheric correction, has increased agreement with the PROMICE data and an apparent accuracy that also exceeds that in the comparison with the National Aeronautics and Space Administration (NASA) Moderate Resolution Imaging Spectroradiometer (MODIS) MOD10A1 [15] albedo product (Figures 4 and 5). Examples are made for the southern Greenland ice sheet QAS_L PROMICE location and the northwestern ice sheet THU_L PROMICE location (Figure 4, right). The current approach, advanced from that reported in [6] by the improved atmospheric correction, has increased agreement with the PROMICE data and an apparent accuracy that also exceeds that in the comparison with the National Aeronautics and Space Administration (NASA) Moderate Resolution Imaging Spectroradiometer (MODIS) MOD10A1 [15] albedo product ( Figures  4 and 5). Examples are made for the southern Greenland ice sheet QAS_L PROMICE location and the northwestern ice sheet THU_L PROMICE location (Figure 4, right).

Discussion
Standard underlying surface albedo retrieval algorithms based on single view observations can be corrected for surface anisotropy effects using multiple day observations of reflected solar light for a given site to cover necessary illumination/observation geometries needed for the respective integration procedures with respect to the corresponding zenith, viewing, and relative azimuth angles. In our approach, we use the analytical relationship between the bottom-of-atmosphere reflectance and the spherical albedo for clean snow underlying surface (Equations (1) and (2)) in the near infrared (865 and 1020 nm), where atmospheric contribution to the signal as registered on a satellite is small, to derive the snow spherical albedo from measurements at a fixed illumination/observation geometry. We underline that our major assumption is that atmospheric influences on OLCI measurements at 865 and 1020 nm in polar regions are weak and can be neglected. Then, top of atmosphere (TOA) reflectance almost coincides with bottom of atmosphere (BOA) reflectance (see Equation (2)). In the framework of our technique, the multiple day observations for the same site are not required, and the snow albedo for a given place can be derived in approximately one hour after the satellite acquisition time. In the case of polluted snow, the spherical albedo is found from Equation (9) for an assumed aerosol model. The technique also incorporates the calculation of plane spectral and broadband albedo. We find that the errors of the Case 1 snow are usually in the range 1-2% in the visible as compared to the ground measurements [6]. They can increase to 3-5% for the spectral albedo in the near IR and for polluted snow. The retrievals for the dark pixels are less accurate because the underlying theory is more accurate for the case of bright pixels [6]. The retrieval method presented here is the extension of the technique described in [6] for the cases where the aerosol load cannot be neglected.
Water exists in three thermodynamic phases (liquid, solid, gas) both in the atmosphere and in the underlying surface. The separation of clean (Case 1) and polluted (Case 2) waters has been useful in oceanic remote sensing using spaceborne observations. We show that a similar separation of satellite retrievals for clean and polluted snow areas (Case 1 and Case 2 snow) is useful in remote sensing of snow from space. Actually, a similar separation of cases is of importance in cloud remote sensing, where modern cloud remote sensing algorithms are based on the assumption of clean (Case 1) clouds. The polluted (Case 2) clouds exist, but up to now, their study is much less advanced.
In this paper, we propose fast snow albedo retrieval techniques for both Case 1 and Case 2 snow. The results for the clean snow are more accurate and robust. The retrievals for the Case 2 snow are less accurate and are based on the simplified atmospheric correction procedure specified in Equation (9) and the general relationship between reflectance and albedo given by Equation (10).

Discussion
Standard underlying surface albedo retrieval algorithms based on single view observations can be corrected for surface anisotropy effects using multiple day observations of reflected solar light for a given site to cover necessary illumination/observation geometries needed for the respective integration procedures with respect to the corresponding zenith, viewing, and relative azimuth angles. In our approach, we use the analytical relationship between the bottom-of-atmosphere reflectance and the spherical albedo for clean snow underlying surface (Equations (1) and (2)) in the near infrared (865 and 1020 nm), where atmospheric contribution to the signal as registered on a satellite is small, to derive the snow spherical albedo from measurements at a fixed illumination/observation geometry. We underline that our major assumption is that atmospheric influences on OLCI measurements at 865 and 1020 nm in polar regions are weak and can be neglected. Then, top of atmosphere (TOA) reflectance almost coincides with bottom of atmosphere (BOA) reflectance (see Equation (2)). In the framework of our technique, the multiple day observations for the same site are not required, and the snow albedo for a given place can be derived in approximately one hour after the satellite acquisition time. In the case of polluted snow, the spherical albedo is found from Equation (9) for an assumed aerosol model. The technique also incorporates the calculation of plane spectral and broadband albedo. We find that the errors for the Case 1 snow are usually in the range 1-2% in the visible as compared to the ground measurements [6]. They can increase to 3-5% for the spectral albedo in the near IR and for polluted snow. The retrievals for the dark pixels are less accurate because the underlying theory is more accurate for the case of bright pixels [6]. The retrieval method presented here is the extension of the technique described in [6] for the cases where the aerosol load cannot be neglected.
Water exists in three thermodynamic phases (liquid, solid, gas) both in the atmosphere and in the underlying surface. The separation of clean (Case 1) and polluted (Case 2) waters has been useful in oceanic remote sensing using spaceborne observations. We show that a similar separation of satellite retrievals for clean and polluted snow areas (Case 1 and Case 2 snow) is useful in remote sensing of snow from space. Actually, a similar separation of cases is of importance in cloud remote sensing, where modern cloud remote sensing algorithms are based on the assumption of clean (Case 1) clouds. The polluted (Case 2) clouds exist, but up to now, their study is much less advanced.
In this paper, we propose fast snow albedo retrieval techniques for both Case 1 and Case 2 snow. The results for the clean snow are more accurate and robust. The retrievals for the Case 2 snow are less accurate and are based on the simplified atmospheric correction procedure specified in Equation (9) and the general relationship between reflectance and albedo given by Equation (10). We find that the influence of the aerosol load on the retrieval of the snow surface albedo is weak in the case of small atmospheric aerosol optical thickness characteristic for Arctic and alpine areas. As a matter of fact, Equation (10) performs well not only for snow but also for other types of weakly absorbing and strongly scattering media, such as clean and polluted bare ice. Therefore, the technique proposed here can be used to study the albedo of terrestrial bare ice surfaces, as demonstrated in Figures 3-5, where low albedo values correspond not to snow but to ice underlying surfaces.

Conclusions
Through comparison with independent ground observations, the proposed fast atmospheric correction technique is shown to perform accurately in wide a range of conditions from a 2100 m elevation mid-latitude location in the French Alps to a Greenland ice sheet network of 15 locations spanning a 2076 km north-south, 18.9 degrees latitude, and 2390 m in elevation. It should be pointed out that snow albedo satellite retrievals are often biased due to the assumed shapes of ice grains (spheres, columns, fractal particles, etc.) used in the retrieval process. We use the notion of the effective absorption length in this work. It makes it possible to include all shape-dependent constants in the value of effective absorption length determined from the satellite measurements themselves. This reduces the snow grain shape effect on the retrievals (at least in the OLCI spectral range). The atmospheric correction is performed assuming the aerosol model and the aerosol optical thickness ahead of retrievals. The associated errors do not lead to considerable errors in the retrieved snow albedo in the case of low aerosol load, as demonstrated in Figure 1.
The current approach, advanced from that reported in [6] by the improved atmospheric correction, has increased agreement with ground observations and an apparent accuracy that also exceeds that of the NASA MODIS MOD10A1 [15] broadband albedo product. A next step is to process the full OLCI catalogue from Sentinel-3A and B satellites over 100% snow or ice covered areas of our planet using this new algorithm. The product would offer the climate research community a new and enhanced quality snow and ice albedo product, which will lead to the advancement of our knowledge of snow albedo effects on the terrestrial climate change [1,2].   [22,23]) data, Level 1.5 retrievals with the residual error < 5%, and AOT (440 nm) < 0.20). Total number of retrievals is 5316 divided in 12 groups with 443 asymmetry parameters in each group. Table A1. The spectral dependence of ozone vertical optical thickness τ ozone in terrestrial atmosphere at the ozone load equal to 405 Dobson units (DU). The results are derived assuming particular shapes of temperature, pressure, and ozone concentration vertical distribution as discussed in [24]. We find that the variation of profiles does not change the value of τ ozone significantly.