Arctic Sea Ice Albedo Estimation from Fengyun-3C/Visible and Infra-Red Radiometer

Abstract

The sea ice albedo can amplify global climate change and affect the surface energy in the Arctic. In this paper, the data from Visible and Infra-Red Radiometer (VIRR) onboard Fengyun-3C satellite are applied to derive the Arctic sea ice albedo. Two radiative transfer models, namely, 6S and FluxNet, are used to simulate the reflectance and albedo in the shortwave band. Clear sky sea ice albedo in the Arctic region (60°~90°N) from 2016 to 2019 is derived through the physical process, including data preprocessing, narrowband to broadband conversion, anisotropy correction, and atmospheric correction. The results are compared with aircraft measurements and AVHRR Polar Pathfinder-Extended (APP-x) albedo product and OLCI MPF product. The bias and standard deviation of the difference between VIRR albedo and aircraft measurements are −0.040 and 0.071, respectively. Compared with APP-x product and OLCI MPF product, a good consistency of albedo is shown. And analyzed together with melt pond fraction, an obvious negative relationship can be seen. After processing the 4-year data, an obvious annual trend can be observed. Due to the influence of snow on the ice surface, the average surface albedo of the Arctic in March and April can reach more than 0.8. Starting in May, with the ice and snow melting and melt ponds forming, the albedo drops rapidly to 0.5–0.6. Into August, the melt ponds begin to freeze and the surface albedo increases.

1. Introduction

Surface albedo is defined as the ratio between the upward reflected irradiance and the downward incident irradiance under the influence of solar radiation, representing the surface reflection ability. For the global average surface, only 13% of energy is reflected, and the rest is absorbed by the Earth system [1]. Compared with the land and ocean surface, sea ice is more capable of reflecting solar radiation. As a result of climate warming, the decrease in sea ice, glaciers, snow, and ice shelf reduces the albedo of the Earth’s surface, enhancing the surface absorption of solar radiation and causing global temperature to rise and sea ice to melt faster [2,3]. It is called the sea ice albedo positive feedback mechanism, which is the second main factor for the Arctic amplification phenomenon [4,5,6]. In recent years, sea ice extent in the Arctic has generally shown a decreasing trend, and the ice age and sea ice thickness have also been diminishing [7]. These changes can be seen from the trend of albedo change. Therefore, it is significantly meaningful to monitor sea ice albedo in the Arctic for the study of global climate change and surface energy budget.

The harsh environment in the Arctic makes it difficult to conduct continuous scientific research in the Arctic compared to most other regions. Fortunately, with the continuous development of satellite remote sensing technology, polar-orbiting meteorological satellites can ensure the acquisition of Arctic observation data with high spatial resolution and long time series. Based on a large number of satellite data, several retrieval methods have been developed to retrieve the sea ice albedo in the Arctic. The sea ice albedo products now available include the forecast albedo from ERA5, the AVHRR Polar Pathfinder-Extended (APP-x) albedo product, the CM SAF Cloud, Albedo And Surface Radiation dataset from AVHRR data–second edition (CLARA-A2) albedo product, and the MERIS/OLCI melt pond fraction and sea ice albedo (MPF) product.

The earlier albedo retrieval method is the parameter algorithm, exploring the relationship between sea ice albedo and parameters to estimate the sea ice albedo. However, there are spatial differences and seasonal variations of sea ice albedo, resulting in low accuracy [8]. Most studies are based on the Advanced Very High Resolution Radiometer (AVHRR) data to estimate the Arctic broadband sea ice albedo with a traditional retrieval algorithm [9,10,11,12,13,14]. The traditional algorithm is the most popular algorithm, simulating the transfer process of the solar signal, mostly including the narrowband to broadband conversion, anisotropy correction, and atmospheric correction. But it differs in different orders in various methods. AVHRR is carried on the NOAA series polar-orbiting meteorological satellites, including one visible channel, two near-infrared channels and three thermal infrared channels, making it possible to estimate long time series and the high temporal resolution of clear sky albedo.

Both APP-x and CLARA-A2 products are retrieved using traditional retrieval method based on AVHRR data, but there are differences in the specific retrieval process [12,13,14]. The APP-x product is compared with the in situ measurements from the Surface Heat Budget of the Arctic Ocean (SHEBA) project in 1998. The bias and the root mean square error of the difference are −0.028 and 0.08, respectively. The CLARA-A2 product is validated with the in situ measurements from the Tara polar ocean expedition in 2007, of which the bias and the root mean square error of the difference are 0.078 and 0.097, respectively [15]. However, the spatial resolution is 25 km in both products, which cannot embody the advantages of high resolution in the visible band and near-infrared band.

To realize a much faster calculation, a direct estimation method based on MODIS is developed [16], which is first applied for the land albedo with high resolution and relatively maturation [17,18]. The algorithm considers a pixel as a mixture of sea ice, ocean, and snow, building a linear empirical lookup table between the top of atmosphere reflectance and surface broadband albedo based on a bidirectional reflectance distribution function (BRDF) database. Compared with in situ measurements from Tara, the bias and the root mean square error of the difference are 0.008 and 0.068, respectively. The algorithm can estimate sea ice albedo quickly as long as it knows the solar and observation geometry; however, it heavily depends on the prior BRDF database, and is easily affected by the sensor noise and cloud detection error.

According to the fact that part of sea ice melts and melt ponds form on the surface in summer, the Melt Pond Detector (MPD) algorithm, which also follows the steps of the traditional algorithm, firstly implements a series of iterative procedure to recognize the melt ponds [19,20]. And for anisotropy correction, surface BRDF is calculated based on an optical model of the sea ice reflection of white ice and pond. Compared with the airborne measurements from Arctic CLoud Observations Using airborne measurements during polar Day project (ACLOUD) on 25 June 2017, the bias and the root mean square error of the difference are 0.006 and 0.02, respectively. However, this algorithm can only be applicable to areas with 100% sea ice concentration and the presence of melt ponds.

Table 1. Relevant information for main arctic sea ice albedo products and algorithms.

Product	Retrieval Algorithms	Spatial Resolution	Accuracy	Time Period	References
ERA5 reanalysis	Data assimilation	0.25°	\	1940–present	Hersbach et al. [21]
APP-x	Traditional algorithm	25 km	RMSE = 0.08	1982–present	Key et al. [9] Key et al. [16]
CLARA -A2	Traditional algorithm	25 km	RMSE = 0.069	1982–2015	Riihela et al. [13] Karlsson et al. [14]
MPF	Melt pond detection	12.5 km	RMSE = 0.02	2002–2011 2017–2023	Zege et al. [20] Pohl et al. [19]
\	Direct estimation	\	RMSE = 0.068	\	Qu et al. [16]

presents the relevant information for the main arctic sea ice albedo products and algorithms.

Fengyun-3 is the Chinese second-generation polar-orbiting meteorological satellite with the goal of achieving all-weather, multispectral, and three-dimensional observations of global atmospheric and geophysical elements. The third satellite of Fengyun-3 (FY-3C) was launched on 23 September 2013. The Visible and Infra-Red Radiometer (VIRR) on board has highly sensitive visible channels and three infrared atmospheric window channels, which can be used for sea ice albedo retrieval and cloud detection, respectively.

Since the surface of sea ice is often covered by a certain thickness of snow, the observed surface albedo is mostly a combined result of snow and sea ice. So, the retrieval algorithm and characteristic of surface albedo include the study of sea ice and snow. At present, there is no operational sea ice albedo product from FY-3C satellites. This study attempts to develop the algorithm appliable to the FY-3C/VIRR sensor for only one visible band and one near infrared band. And the traditional retrieval algorithm with the amount of water vapor and the aerosol optical depth suitable for the global scale is most widely used, but is not quite appropriate for the Arctic. This study obtains the retrieval coefficients from different radiative transfer models with new atmospheric composition range. Section 2 describes the data and method used. Section 3 presents the retrieval result and validates it with airborne observations, the APP-x albedo product, and the OLCI MPF product. Section 4 describes the discussion of spatial and temporal distribution. Section 5 presents the conclusions.

2. Materials and Methods

2.1. Materials

2.1.1. FY-3C/VIRR L1B Data

FY-3C/VIRR has five visible channels and five near-infrared and infrared channels. Limited by the dynamic range and noise required for retrieval, the reflectance of channels 1 and 2 and emissive radiance of channels 3, 4, and 5 were used in this study.

Table 2. Relevant information of arctic sea ice albedo algorithms.

Channel	Center Wavelength (μm)	Band Range (μm)	NER (%)/ NETD (300 K)	Dynamic Range ($\mathsf{\rho }$/K)	Application in Algorithm
1	0.630	0.58$–$0.68	0.1%	$0–$100%	Albedo calculation; Cloud detection
2	0.865	0.84$–$0.89	0.1%	$0–$100%	Albedo calculation; Cloud detection
3	3.740	3.55$–$3.93	0.3 K	180$–$350 K	Cloud detection
4	10.80	10.3$–$11.3	0.2 K	180$–$350 K	Cloud detection
5	12.00	11.5$–$12.5	0.2 K	180$–$350 K	Cloud detection

NER refers to Noise Equivalent Reflectance, and NETD refers to Equivalent Temperature Difference.

presents the information of these channels.

The VIRR Level 1B data contain not only the observed radiance, but also the longitude, latitude, satellite zenith angle, satellite azimuth angle, solar zenith angle, solar azimuth angle, digital elevation model, and land/sea mask of each pixel. The L1B data are generated every 5 min in an HDF5 format with nadir resolution of 1.1 km and a local descending time of 10:15. Its date range is 25 September 2013 to 5 February 2020 with June and July in 2015 missing, which will affect the continuous analysis. Therefore, only data from 2016 to 2019 are processed in this study.

2.1.2. ERA5 Hourly Data on Single Levels

ERA5 is a global reanalysis dataset of climate and weather. It provides multiple parameter products of the atmosphere from 1950 to present by data assimilation. It combines model data with observations and stays updated with new observations. In this study, the collections “total column water vapor” and “sea-ice cover” on single levels are used. The spatial resolution is 0.25° × 0.25°, and temporal resolution is 1 h [21].

2.1.3. Aircraft Measurements

The ACLOUD campaign and the Airborne measurements of radiative and turbulent FLUXes of energy and momentum in the Arctic boundary layer (AFLUX) campaign are both parts of Arctic Amplification: Climate Relevant Atmospheric and Surface Processes and Feedback Mechanisms (AC)3, which aims to observe clouds, aerosol particles, and some surface properties with two aircrafts, Polar 5 and Polar 6 in ACLOUD [22], and one aircraft, Polar5 in AFLUX [23]. The ACLOUD campaign was conducted from 23 May to 26 June 2017, which was exactly in the melting season. And the AFLUX campaign was conducted in March and April, 2019, which was in early spring. Kipp & Zonen CMP22 Pyranometers carried on the aircrafts were used to monitor the upward and downward broadband solar irradiance at flight level, between 0.2 and 3.6 μm in the north of Svalbard, Norway. The observed region was in the marginal sea ice zone, mainly formed by thin ice, with sea ice thickness of about 0.6–2 m [24]. The surface albedo can be calculated as the ratio of downward solar irradiance to upward solar irradiance.

2.1.4. APP-x Albedo Product

NOAA processes AVHRR observations and atmospheric reanalysis data and generates the APP-x product, which aims to obtain surface, cloud, and radiative transfer characteristics over the Arctic and Antarctic from 1982 to present [15]. The spatial resolution of albedo products is 25 km, and it includes twice daily products with local solar time 14:00 and 04:00 [12]. FY-3C is a morning satellite. Therefore, only APP-x product with the central local time of 1400 is considered in the comparison with results.

2.1.5. MPF V1.7 Data

As the temperature increases, the incomplete melting of snow and ice will form melt ponds on the ice surface. The areal fraction of sea ice covered by these melt ponds is called melt pond fraction. Based on the Ocean and Land Colour Instrument (OLCI) onboard Sentinel-3 and MEdium Resolution Imaging Spectrometer (MERIS) onboard Envisat, the daily clear sky melt pond fraction is estimated by iteratively changing the parameters in the area with 100% sea ice concentration and the daily clear sky spectral albedo is calculated based on the BRDF and fraction of sea ice and melt ponds [19,25]. The time range is from 15 May to 15 September, 2002 to 2011 (based on MERIS) and 2017 to 2020 (based on OLCI), with a spatial resolution of 12.5 km and temporal resolution of 1 day [26].

2.2. Methods

In this paper, considering the transfer process of solar radiation, the retrieval method is divided into the following steps. First, the necessary data processing and cloud detection is performed to recognize the clear sky sea ice surface. Second, the visible and near-infrared satellite observation data are converted into broadband reflectance using 6S radiative transfer model. The anisotropic factor correction is then implemented to correct the reflectance directional distribution. Lastly, the atmospheric correction coefficients are obtained from FluxNet radiative transfer model. The retrieval scheme of all steps is shown in Figure 1.

2.2.1. Cloud Detection

In the shortwave band, solar radiation will be reflected back by clouds in the atmosphere, which means that the value of albedo observed by satellite may be from clouds instead of the Earth’s surface. The cloud coverage is high over the Arctic region, accounting for more than 50% most of the year [27], contributing great differences to the results. Therefore, cloud detection is required.

The cloud detection method is the threshold detection. Since the cloud and snow/ice have similar high reflectance, the cloud detection tree is shown in Figure 2, which contains the following parameters: the brightness temperature difference of 11 μm and 12 μm (BTD45), the brightness temperature difference of 3.7 μm and 12 μm (BTD34), the reflectance in 0.63 μm (REF), and the uniformity detection. The BTD34 larger than BTD34_THRESH can detect the thick cloud. The emissivity of clouds increases as the cloud thickness increases. Depending on the wavelength and thickness of the cloud, the emissivity varies between 0 and 0.97. For thick clouds, an emissivity of 3.7 μm is about 0.8, and that of 11 μm is about 0.97 [28]. The emissivity of ice and snow does not change much in the three channels, and its emissivity is greater than 0.96 in the thermal infrared band, which is regarded as 0.99 in many applications [29,30]. Both ice and liquid clouds have higher transmittances at 11 µm than at 12 μm for their absorption and scattering, so the split window channels can detect the cirrus cloud and thin water cloud. And the BTD45_THRESH varies with brightness temperature from 220 K to 280 K, which is shown in

Table 3. BTD45_THRESH as a function of 11 μm brightness temperature.

11 μm BT (K)	220	230	240	250	260	270	280
BTD45_THRESH (K)	0.8	0.9	1.2	1.45	1.8	2.45	3.4

. The single channel of reflectance lower than REF_THRESH is used to detect thin clouds, since some thin clouds are too thin to appear apparent characteristic of brightness temperature difference. And then the uniformity detection is implemented to test the cloud edge. Figure 3 shows the clear sky reflectance map of channel 1 after cloud detection.

2.2.2. Narrowband to Broadband Conversion

Solar radiation can be divided into three parts, namely, visible, near-infrared, and shortwave infrared radiation. The wavelength range of visible band is about 0.35–0.75 μm, that of the near-infrared band is 0.75–2.8 μm, and that of the shortwave infrared is 2.8–4 μm. Accordingly, the broadband reflectance can be converted into three parts, too. When the wavelength is larger than 1.4 μm, the reflected radiation of ice and snow drops very fast [31], accounting for a small proportion of the total solar radiation; therefore, the spectral reflectance in the shortwave infrared band can be negligible.

The two channels of FY3C/VIRR are 0.58–0.68 μm and 0.84–0.89 μm, distributing properly in the wavelength range of visible band and near infrared band, respectively. Then, the conversion formula for broadband reflectance and two sensor reflectance values can be calculated approximatively as follows:

$${\rho }_{toa}=\mathrm{a}{\mathsf{\rho }}_{1,\mathrm{toa}}+\mathrm{b}{\mathsf{\rho }}_{2,\mathrm{toa}}+\mathrm{c},$$

(1)

where ${\mathsf{\rho }}_{1,\mathrm{toa}}$ and ${\mathsf{\rho }}_{2,\mathrm{toa}}$ denote the top of atmosphere reflectance of channels 1 and 2 from VIRR, respectively. ${\rho }_{toa}$ represents broadband reflectance at the top of atmosphere.

This conversion process is realized by 6S radiative transfer model [32]. The 6S model can simulate the radiation transmission in the wavelength range of 0.25–4 μm under different atmospheric and surface conditions. In terms of band range, the two visible spectral ranges of VIRR and their corresponding spectral response functions (SRF) are input into the model, respectively, and the spectral response of the broadband is 1.0. The in situ-measured reflectance spectra from MOSAIC expedition are chosen to represent different conditions of surface conditions. There are spectral reflectance curves under clear sky or nearly clear sky shown in Figure 4, whose surfaces mainly comprise ice, snow-covered ice, and pond-covered ice. These curves are input into the 6S model to simulate the top of the atmosphere reflectance. The range of measurement time is from June to September, which contains the new/old snow and new/old ice, etc. The missing part in 1300–1400 nm and 1800–2000 nm is smoothly interpolated because it is too noisy. And the wavelength larger than 2500 nm is assumed to be a constant which would not introduce deviation. For the observation geometry, the nadir angle for the view zenith angle and at increment of 5° intervals in the range of 0–85° for the solar zenith angle are set in the model. The subarctic winter and summer are chosen for the atmosphere model and four aerosol combinations are input with six aerosol optical depths at 550 nm (0.05, 0.1, 0.2, 0.3, 0.4, 0.5). The top of atmosphere reflectance in the visible and near-infrared bands and the broadband (short wave) under the above conditions are obtained by model simulation. The four aerosol combinations are 0% soot, 5% water-soluble, and 95% oceanic component; 5% soot, 5% water-soluble, and 90% oceanic component; 10% soot, 5% water-soluble, and 85% oceanic component; and 15% soot, 5% water-soluble, and 80% oceanic component, respectively.

Table 4. The narrowband to broadband conversion coefficients and relative bias with the variation of AOD.

AOD	a	b	c	Relative Bias
all	0.3110	0.5522	0.0124	\
0.05	0.2701	0.6028	0.0053	0.89%
0.10	0.2573	0.6153	0.0061	0.87%
0.20	0.2892	0.5772	0.0103	0.32%
0.30	0.3353	0.5259	0.0140	−0.36%
0.40	0.3956	0.4636	0.0155	−1.04%
0.50	0.3698	0.4770	0.0222	−0.93%

presents the narrowband to broadband conversion coefficients calculated by the least square method of six AODs from 0.05 to 0.5. The coefficients fitting from all six AODs are applied to the conversion, which could introduce errors. Table 4 also gives the relative bias that can be caused by the variation in AOD.

In the regression, two thirds simulated points are used for fitting and the remaining third is used for validation. For validation points, the mean bias and standard deviation between the simulated broadband albedo and calculated broadband albedo is 0.0 and 0.0119, respectively. Figure 5 shows the clear sky broadband reflectance at the top of atmosphere before anisotropy correction.

Here, only two narrowband channels are used to represent the shortwave broadband. There may be a question as to whether this is sufficient compared with other multi-channel sensors, such as MODIS, ASTER, and MISR. To answer the question, two more channels (0.45–0.55 μm and 0.74–0.79 μm) are applied to the process of narrowband to broadband conversion. Residual Standard Error (RSE) and Coefficient of Determination (R²) are introduced to evaluate result of the regression. Two atmosphere models, fifteen sensor zenith angles, and thirty surface types are considered in simulation. Two thirds are used for fitting, and one third is used for evaluation. For the one third, the RSE and R² between fitting results and simulation of four channels are 0.0107 and 0.9965, while the RSE and R² between fitting results and simulation of two channels are 0.0170 and 0.9959. It can be said that the results of four channels are better than those of two channels, but there is no significant deviation for two channels either.

2.2.3. Anisotropy Correction

Anisotropy correction is performed to correct the dependence of reflectance on the sun–satellite–surface geometry. For the Lambertian plane, the reflected radiance is the same for any direction. However, reflectance is directional under the earth surface. The forward scattering is stronger for the snow and ice surface. The anisotropic factors are used to correct the reflectance bidirectional properties of the surface and are obtained from a lookup table. It is derived from angular radiation models that are required for analysis of satellite measurements of Earth radiation. The models describe the angular variation in radiance through angular bins under 12 kinds of surface types. And the bins are varied with solar zenith angle, satellite zenith angle and relative azimuth angle [33]. The anisotropic factor is one of the shortwave model parameters, which can realize the bidirectional correction for radiation. It is calculated using the following equation [34]:

$$f=\frac{{\mathrm{L}}_{\mathrm{r}}\left({\mathsf{\theta }}_{\mathrm{i}},{\mathsf{\phi }}_{\mathrm{i}},{\mathsf{\theta }}_{\mathrm{r}},{\mathsf{\phi }}_{\mathrm{r}}\right)}{{\mathrm{L}}_{\mathrm{r}}^{\mathrm{i}\mathrm{d}}\left({\mathsf{\theta }}_{\mathrm{i}},{\mathsf{\phi }}_{\mathrm{i}}\right)}=\frac{{\mathrm{E}}_{\mathrm{i}}\left({\mathsf{\theta }}_{\mathrm{i}},{\mathsf{\phi }}_{\mathrm{i}}\right)}{{\mathrm{L}}_{\mathrm{r}}^{\mathrm{i}\mathrm{d}}\left({\mathsf{\theta }}_{\mathrm{i}},{\mathsf{\phi }}_{\mathrm{i}}\right)}·\frac{{\mathrm{L}}_{\mathrm{r}}\left({\mathsf{\theta }}_{\mathrm{i}},{\mathsf{\phi }}_{\mathrm{i}},{\mathsf{\theta }}_{\mathrm{r}},{\mathsf{\phi }}_{\mathrm{r}}\right)}{{\mathrm{E}}_{\mathrm{i}}\left({\mathsf{\theta }}_{\mathrm{i}},{\mathsf{\phi }}_{\mathrm{i}}\right)}=\frac{\mathsf{\pi }}{\alpha }\ast \frac{{\mathrm{L}}_{\mathrm{r}}\left({\mathsf{\theta }}_{\mathrm{i}},{\mathsf{\phi }}_{\mathrm{i}},{\mathsf{\theta }}_{\mathrm{r}},{\mathsf{\phi }}_{\mathrm{r}}\right)}{{\mathrm{E}}_{\mathrm{i}}\left({\mathsf{\theta }}_{\mathrm{i}},{\mathsf{\phi }}_{\mathrm{i}}\right)}=\frac{\rho }{\alpha }$$

(2)

where $f$ is the anisotropic factor, ${\mathrm{L}}_{\mathrm{r}}\left({\mathsf{\theta }}_{\mathrm{i}},{\mathsf{\phi }}_{\mathrm{i}},{\mathsf{\theta }}_{\mathrm{r}},{\mathsf{\phi }}_{\mathrm{r}}\right)$ is reflected radiance, ${\mathrm{L}}_{\mathrm{r}}^{\mathrm{i}\mathrm{d}}\left({\mathsf{\theta }}_{\mathrm{i}},{\mathsf{\phi }}_{\mathrm{i}}\right)$ is the radiance reflected by the ideal Lambertian surface, ${\mathrm{E}}_{\mathrm{i}}\left({\mathsf{\theta }}_{\mathrm{i}},{\mathsf{\phi }}_{\mathrm{i}}\right)$ is the incident irradiance, $\mathsf{\theta }$ is the zenith angle, and $\mathsf{\phi }$ is the azimuth angle. And the correction formula is as follows:

$${\mathsf{\alpha }}_{\mathrm{toa}}={\rho }_{toa}/f$$

(3)

where ${\mathsf{\alpha }}_{\mathrm{toa}}$ is the albedo at the top of atmosphere after anisotropy correction. Figure 6 shows the clear sky albedo at the top of atmosphere after anisotropy correction. The proportion of the correction to the origin value is from −20.38% to 13.80%.

2.2.4. Atmospheric Correction

On ideal conditions (no atmosphere), some of the incident radiation from the Sun is absorbed by the Earth’s surface while the rest is reflected off into space. Therefore, the surface properties directly determine the measured radiance. However, in an actual situation, some of the incident photons are absorbed or scattered by the atmosphere. Considering the process of atmospheric radiation, the incident radiation can be divided into three parts: part of the solar irradiance is absorbed by atmospheric components, part is absorbed while reaching the ground, and the remaining part is scattered by the atmosphere, and then the scattered signal is observed by the satellite sensors [35]. The relationship can be represented as the equation:

$${\mathrm{E}}_{\mathrm{toa}}={\mathrm{E}}_{\mathrm{abs}}+{\mathrm{E}}_{\mathrm{s}}\left(1-{\mathsf{\alpha }}_{\mathrm{s}}\right)+{\mathsf{\alpha }}_{\mathrm{toa}}{\mathrm{E}}_{\mathrm{toa}}$$

(4)

where ${\mathrm{E}}_{\mathrm{toa}}$ refers to the incident irradiance at the top of atmosphere, ${\mathrm{E}}_{\mathrm{abs}}$ refers to the irradiance absorbed by the atmosphere, ${\mathsf{\alpha }}_{\mathrm{s}}$ refers to the surface albedo, and ${\mathrm{E}}_{\mathrm{s}}$ refers to the incident irradiance at the surface. According to the above radiative transfer equation, the relationship between surface albedo and the top of atmosphere albedo can be expressed as follows:

$${\mathsf{\alpha }}_{\mathrm{toa}}=({\mathrm{E}}_{\mathrm{toa}}-{\mathrm{E}}_{\mathrm{abs}}-{\mathrm{E}}_{\mathrm{s}})/{\mathrm{E}}_{\mathrm{toa}}+{\mathsf{\alpha }}_{\mathrm{s}}{\mathrm{E}}_{\mathrm{s}}/{\mathrm{E}}_{\mathrm{toa}}$$

(5)

Nevertheless, the process of calculating the solar irradiance based on the above equation for each pixel is slow and complicated. The absorption of water vapor and ozone, as well as the scattering of aerosols are primary causes of the signal interference from the atmosphere in the visible and near-infrared bands. Therefore, these parameters can be modified within the atmospheric radiative transfer model in order to simulate the atmospheric transfer process. Equation (5) can be simplified to Equation (6):

$${\mathsf{\alpha }}_{\mathrm{toa}}=\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{c}\mathrm{e}\mathrm{p}\mathrm{t}+\mathrm{s}\mathrm{l}\mathrm{o}\mathrm{p}\mathrm{e}\times {\mathsf{\alpha }}_{\mathrm{s}}$$

(6)

A lookup table for atmospheric correction coefficients can be established to make the correction effective [36]. This process is achieved through the FluxNet radiative transfer model [37], which is the neural network version of the Streamer model. Compared with other radiative transfer models, this model can directly simulate the radiative irradiance both at the surface and the top of atmosphere, and is suitable for the calculation of a wide spectral range with the built-in snow and ice surface types. It has a fast calculation with the adequate accuracy for the study. For simulation, the variables total column ozone, aerosol optical depth, and total column water amount are input to FluxNet model to change the atmospheric conditions. The input ranges of above variables are 6.5–9.5 g/m², 0.05–0.45, and 0.05–2.55 g/cm², and the sampling intervals are 1.0 g/m², 0.1, and 0.5 g/cm², respectively.

Thus, the incident and reflected flux density at the top of atmosphere under different atmospheric conditions and solar zenith angles are simulated, as well as at the surface of snow/ice. The albedo at the top of atmosphere and surface are then calculated according to the following formula:

$${\mathsf{\alpha }}_{\mathrm{toa}}\left({\mathsf{\theta }}_0\right)={\mathrm{E}}_{\mathrm{t}}\left({\mathsf{\theta }}_0\right)/{\mathrm{M}}_{\mathrm{t}}$$

(7)

$${\mathsf{\alpha }}_{\mathrm{s}}\left({\mathsf{\theta }}_0\right)={\mathrm{E}}_{\mathrm{s}}\left({\mathsf{\theta }}_0\right)/{\mathrm{M}}_{\mathrm{s}}\left({\mathsf{\theta }}_0\right)$$

(8)

where${ \mathrm{M}}_{\mathrm{t}}$ represents incident solar irradiance, ${\mathrm{E}}_{\mathrm{t}}\left({\mathsf{\theta }}_0\right)$ represents reflected irradiance in the atmosphere, ${\mathrm{M}}_{\mathrm{s}}\left({\mathsf{\theta }}_0\right)$ represents the surface incident irradiance, and ${\mathrm{E}}_{\mathrm{s}}\left({\mathsf{\theta }}_0\right)$ represents the surface reflected irradiance.

Figure 7 shows the relationship between the surface albedo and the top of atmosphere albedo under a range of ozone, aerosol optical depth, and water vapor. From Figure 7a–c, it can be concluded that the amount of water vapor is the main cause in the atmospheric correction, while the total column ozone has a little influence on the corrected surface reflectance; therefore, the fixed value of ozone is set as 6.96 g/m². In addition, it is found that the simulation result of the radiative transfer model is not reliable when the solar zenith angle is large; therefore, the correction coefficients when the solar zenith angle is greater than 60° is taken as the coefficients of angle equalized to 60°. Since the satellite goes through the same place twice a day, there are actually two values into one pixel. The point with smaller solar zenith angles is selected for the projection and then the inversion, to avoid errors by large angles. In this way, there will be almost no points with larger solar zenith angles after projection, which can also effectively avoid the unreliability of the radiation transfer model under large solar zenith angles.

3. Results

Figure 8 shows a daily clear sky sea ice albedo in Arctic with the 4 km Ease Grid projection in 6 June 2017, for example. The average albedo retrieved is 0.76 for this day, which is in the starting of melt period in most parts of the region. The retrieval albedo is evaluated with aircraft measurements, the APP-x albedo product, and the OLCI MPF product in the following subsections.

3.1. Validation with Aircraft Measurements

Most parts of the Arctic region are covered by sea ice and surrounded by ocean throughout the year. Due to the harsh environment, field exploration is expensive and difficult to be conducted, and relevant measurements are quite few. The ACLOUD campaign and the AFLUX campaign measured the solar irradiance under different surface conditions in the melting season and in early spring, respectively.

The retrieved sea ice albedo of VIRR is matched with the aircraft measurements in temporal and spatial windows of 3 h and 4 km. Only the irradiance with the observation height below 100 m is considered to avoid atmospheric influences on the measurements. There were 20 research flights performed during the ACLOUD campaign and 14 research flights performed during AFLUX campaign. After matching and quality control, 21 research flights remain and 417 matching points are found. The space distribution of matching points and its corresponding sea ice concentration, which are in the north of Svalbard, Norway, are shown in Figure 9. Figure 10 shows the scatterplot and frequency distribution of matching points. The bias, standard deviation, and relative bias of difference between aircraft measurements and VIRR albedo are −0.040, 0.071, −4.68%, respectively.

Table 5. Comparison of the retrieved broadband albedo and measured broadband albedo.

Bias	Std	Median	Rsd	RMSE	Relative Bias	R2	Num.
−0.040	0.071	−0.039	0.071	0.081	−4.68%	0.83	391

presents more comparison results.

Since the ACLOUD measurement is performed in the melting season, the ice surface forms melt ponds and the sea ice concentration decreases, both of which can be reflected as the decrease in the sea ice albedo. From the end of May to the end of June, the average sea ice albedo of the ACLOUD matching dataset decreases from 0.80 to 0.64.

A part of negative bias may be caused by the different wavelength range between the aircraft measurements (0.2–3.6 μm) and simulated VIRR albedo (0.35–4 μm). A simulation based on the matching dataset is performed through the 6S model to evaluate how much will the error be. The discrepancy of simulated top of atmosphere broadband albedo is −0.008, which eventually causes an error of −0.010. Moreover, the great difference in match-up numbers in the 4 km grid of the aircraft measurements and VIRR albedo may cause a bias and standard deviation.

3.2. Comparison with APP-x Albedo Products

VIRR and AVHRR have similar spectral channel designs, both including one visible channel and one near-infrared channel. However, due to differences in spectral ranges and spectral response curves of the two sensors, different conversion coefficients of the narrowband to broadband conversion are bound to exist in the retrieval process. Therefore, 6S and FluxNet radiative transfer model are used in this study to determine the correction coefficients suitable for VIRR data.

To compare with the APP-x product, the VIRR retrieved albedo is projected to 25 km, which is the same as the APP-x spatial resolution. Then, the daily clear sky pixels at the same pixel from 2016 to 2019 are compared. The average deviation of the 4-year data is 0.040, and the average standard deviation is 0.077.

Table 6. Yearly average comparison of the retrieved broadband and APP-x/OLCI MPF V1.7 albedo product from 2016 to 2019.

Year	APP-x			OLCI MPF
	Bias	Std	R2	Bias	Std	R2
2016	0.052	0.070	0.85	\	\	\
2017	0.033	0.075	0.90	−0.009	0.083	0.91
2018	0.043	0.066	0.89	−0.013	0.097	0.89
2019	0.053	0.071	0.88	−0.015	0.092	0.88

presents the yearly comparison results. The deviation in the comparison every year is close, from which it can be seen that the two products have good consistency. Especially in April and May, the difference in average deviation is about 0.05, as shown in Figure 11. A bigger bias is shown in the melt season, caused possibly by the change in the surface type. Due to the limitation of the solar zenith angle, the data that can be used for comparison in March are few and the fluctuation is large. Compared with APP-x, the retrieved broadband albedo from VIRR is higher, but it is closer to the aircraft observations.

3.3. Comparison with OLCI MPF Products

The MPF albedo product is retrieved based on the MPD algorithm that different from the traditional algorithm. Since the time range of MERIS is not matched with VIRR, only OLCI data from 2017 to 2019 is compared. The VIRR retrieval daily albedo is projected to 12.5 km to match the OLCI MPF albedo and melt pond fraction.

The average deviation of the 3-year data is −0.011, and the average standard deviation is 0.091. Table 6 presents the yearly deviation and standard deviation. Since the OLCI albedo is retrieved based on melt ponds, the match-up dataset is distributed where there are melt ponds. A good consistency of these two products is shown in Figure 12. The correlation coefficient between OLCI albedo and melt pond fraction is −0.97, and the correlation coefficient between VIRR albedo and OLCI melt pond fraction is −0.93. An obvious negative correlation is found. The albedo is at the lowest value in July, and the number and fraction of melt ponds reach the highest relative values.

4. Discussion

The distribution of daily clear sky sea ice albedo is dispersive on account of the presence of clouds. The monthly composition is performed to analyze the spatial distribution pattern in the whole Arctic region. Figure 13 shows the monthly average sea ice albedo map from March to August, in the period 2016–2019. Limited by the high solar zenith angle in the high latitude region, the data at the center of the North Pole are missing.

Notably, in Figure 13, the albedo in the northern part of the East Siberian Islands is lower than that in other areas, with the sea ice completely melting in August. According to [38], this phenomenon can be explained as the Arctic oscillation. The offshore current in winter will lead to the fragmentation of fixed ice, and thin ice is easily generated in the fragmentation cracks. The thin ice melts earlier than other areas, with a significant decrease in sea ice concentration in May, as well as the albedo.

Due to the influence of snow on the ice surface, the average surface albedo of the Arctic in March and April can reach more than 0.8. Until then, the albedo of multi-year sea ice and first-year sea ice was indistinguishable due to the snow cover on the surface. However, compared with other regions, there is a more significant albedo difference in the Chukchi Sea and the Greenland Sea, which is caused by the difference in surface condition. Sea ice concentration and thickness are low near the open ocean. Starting in May, the ice and snow melt, forming melt ponds on the ice, and resulting in the albedo to drop rapidly to 0.5–0.6. Moving into June and July, it is evident that the albedo of the multi-year sea ice area is higher than that of the first-year sea ice area due to the melting of the snow on the surface. When the ice surface is covered by snow, the incident solar radiation is reflected by the covered snow. There is no apparent discrepancy whatever the characteristics of the underlying ice surface. When snow and ice melt, melt ponds are more likely to form and have a larger area on the surface of the first-year ice compared with the surface of the multi-year ice, and the difference between multi-year sea ice and first-year sea ice appears. In August, the melting ponds begin to freeze and the albedo increases.

There may be some errors in the retrieval process that influence the results. First, inadequate cloud detection is always an existing error resource. Second, simulation with radiative transfer models takes the atmosphere as a unity mixture, which in fact varies in water vapor, aerosol types, aerosol optical depth, ozone, and so on. For example, dust mixed in the surface or soot in the atmosphere could decrease the expected sea ice albedo [39]. This method can be improved in terms of these aspects in the future.

5. Conclusions

In this study, FY-3C/VIRR visible band and near-infrared data are used to estimate the clear sky sea ice albedo in the Arctic, developing the traditional algorithm appliable to the FY-3C/VIRR sensor. The 6S radiative transfer model is used to simulate reflectance in visible and near-infrared bands and shortwave broadband for narrowband to broadband reflectance conversion. The bidirectional reflectance distribution is corrected according to the anisotropy factors of snow and ice surface obtained from ERB and GOES data. The FluxNet radiative transfer model is introduced in the atmospheric correction to simulate the top of atmosphere downwelling irradiance and upwelling irradiance, as well as the surface. The results are compared with the aircraft measurements from ACLOUD and AFLUX campaign, and the average deviation of −0.040 and the standard deviation of 0.071 are obtained. Compared with the APP-x products, the average deviation of daily clear sky pixels is 0.040, and the standard deviation is 0.077. Compared with the OLCI MPF products, the average deviation of daily clear sky pixels is −0.011, and the standard deviation is 0.091. The suitable atmosphere component distribution range for the Arctic region is narrower than the range used in the APP-x albedo product, and the intervals for the lookup table are much smaller, which introduces a relative bias of 1.8%. Uncertainty is discussed in each step.

The albedo of sea ice changes with the season, with the changes in surface conditions, such as snow covering on the surface, snow melting, melt ponds formation, and melt ponds freezing. Albedo in the Arctic shows regular changes every year, and the surface albedo reaches above 0.8 with new snow covered. It begins to decrease in May and reaches a low of around 0.5 in July. The change in albedo is also related to the surface cover type, sea ice concentration, sea ice type, and other factors.

In the future, the retrieval method can be applied to FY series observations with visible and near-infrared bands to achieve longer time series of the sea ice albedo estimation in the Arctic, and can realize the research and application in climatology.