**1. Introduction**

Recently, more evidence reveals the shrinking trend of Arctic sea ice [1,2]. The alteration from high-albedo sea ice to a low-albedo ocean would increase the amount of absorbed solar radiation, leading to a warmer effect and further accelerating the ice melting. Therefore, sea-ice albedo variation has attracted more attention when studying Arctic and global climate change. Satellite observations are

essential for providing sustained, consistent, and near real-time albedo over large, remote, and sparsely populated areas such as sea ice [3,4].

Despite the unprecedented demand for authoritative information on sea-ice albedo, local data resources are limited since the polar region is one of the most under-sampled domains in the climate system. Satellite observations are essential for providing sustained, consistent, and near real-time albedo estimates over large, remote, and sparsely populated areas. Since the last century, based on the reliable operational global imagery from Advanced Very High Resolution Radiometer (AVHRR) data, several studies have discussed the algorithms for mapping broadband albedo of sea ice [5–9]. The available sea-ice albedo products include the AVHRR Polar Pathfinder (APP) albedo product [10], the APP-extended (APP-x) albedo product [11,12], and the Satellite Application Facility on Climate Monitoring (CM-SAF) surface albedo (SAL) product [13]. However, only two spectral visible/near-infrared bands of AVHRR exist, which has limited its accuracy and sensitivity of broadband albedo [14].

With its higher spectral/spatial characteristics and measurement precision as compared to AVHRR, Moderate Resolution Imaging Spectroradiometer (MODIS) produces an operational 500 m/1 km daily albedo product. However, the operational MODIS albedo product left the sea-ice pixels blank due to the high cloud coverage in the Arctic region [12] and there is increasing uncertainty of the MODIS atmospheric correction algorithm on surfaces without dense vegetation coverage [15,16]. Moreover, the traditional bidirectional reflectance distribution function (BRDF) integration algorithm applied for MODIS albedo production assumes a 16-day stable earth surface BRDF pattern to accumulate enough multi-angle observations for BRDF retrieval. This requirement may be easily violated due to the dynamic coverage of sea ice, especially as the Arctic Ocean is experiencing a shift from perennial to seasonal sea ice [17,18]. Although some studies are trying to combine the observations from different satellites to form the multi-angle dataset [19,20], the rolling window length is still too long to grasp the dynamic sea-ice change.

Alternatively, Liang et al. [21] developed a direct estimation algorithm for land surface albedo estimation over the Greenland ice sheet from MODIS top-of-atmosphere (TOA) reflectance data. Instead of the atmospheric correction and BRDF inversion approaches, the direct retrieval algorithm uses an explicit multiple linear regression method for the albedo estimation. Later, it was further improved by applying surface BRDF models in the regression computation. To grasp the influence of BRDF anisotropy, the incident/reflected hemisphere was divided into regular grids, and linear regression relationships were developed for each grid to represent specific solar-viewing geometries. This prototype of the direct regression algorithm was operationally applied for the land surface albedo productions in the mission of Global LAnd Surface Satellite (GLASS) land surface products from MODIS observations and in the mission of S-NPP (Suomi National Polar-orbiting Partnership) from the VIIRS (Visible Infrared Imaging Radiometer Suite) dataset [15,22]. However, the sea-ice albedo production was still a problem due to the lack of sea-ice surface BRDF observations available for training the dataset built-up in the regression approach. No sea-ice BRDF records were available from satellite observations due to the limited solar zenith angles and large cloud fraction in polar regions. Thus, the sea-ice albedo of VIIRS was not released simultaneously. Currently, no operational sea-ice albedo at a 1 km scale is available to the user community.

Qu et al. [23] simulated sea-ice surface BRDF from physical models, which has been proved to be able to represent different sea-ice types and mixed pixels of ice/snow/seawater with different fractions. In this study, their BRDF simulation method was applied to the VIIRS instrument to produce a VIIRS sea-ice albedo (VSIA) product.

As the successor of MODIS, VIIRS started its observation from October of 2011. VIIRS has improved its design upon the capabilities of the operational AVHRR and provides observation continuity with MODIS. VSIA will service weather prediction applications and climate model calibrations. Moreover, combined with the historical observations from MODIS and from the future Joint Polar Satellite System (JPSS), VSIA will dramatically contribute to providing high-accuracy routine sea-ice albedo products and irreplaceable records for monitoring the long-term sea-ice albedo for climate studies. Recently, VSIA has passed its operational readiness review by the Satellite Products and Services Review Board at National Oceanic and Atmospheric Administration's (NOAA) National Environmental Satellite, Data, and Information Service (NESDIS). It will be produced operationally through the NOAA S-NPP data Exploration (NDE) system.

In this study, we aim to evaluate the VSIA accuracy using ground measurements, regarded as the "true value" of the ice/snow albedo. The accuracy assessment attempts to provide the community and users with information about the VSIA performance and its reliability in applications. The rest of the paper is organized as follows. Section 2 introduces the VIIRS albedo product, VSIA algorithm, and the operational production framework. Ground measurements and the validation results are presented in Section 3. More information about the VSIA algorithm is shown in Section 4, including the illustration of large-scale time-series retrieval instances from the VSIA algorithm, probing of look up table (LUT) profiles at specific angles, and highlights of limitations in this validation attempt. Section 5 summarizes this study.

## **2. Method and Dataset**

#### *2.1. VIIRS Albedo Product*

The VIIRS sensor is a component of the Suomi National Polar-orbiting Partnership (S-NPP) satellite and the Joint Polar Satellite System (JPSS). S-NPP was launched on 28 October 2011, and JPSS-1 (NOAA-20) was launched on 18 November 2017. The VIIRS was designed to improve the series of measurements initiated by the Advanced Very High Resolution Radiometer (AVHRR) and the Moderate Resolution Imaging Spectroradiometer (MODIS) [24].

The VIIRS albedo product has been designed as an environmental data record (EDR) generated in a granule-based format and is currently processed in the NOAA near real-time Interface Data Processing Segment (IDPS); archived and distributed by NOAA's Comprehensive Large Array-data Stewardship System (CLASS). The next version of the VIIRS albedo product will be generated by the S-NPP Data Exploration (NDE) system, in which sea-ice albedo retrievals will be included.

The VIIRS albedo deploys a direct estimation algorithm from single-date/angular observations, which is capable of grasping the dynamic variation of surface BRDF change and is critical for the accumulation and melting seasons of sea ice. Due to the high spatial resolution (750 m at nadir) and high temporal resolution (crosses the equator about 14 times daily in the afternoon orbit), the VIIRS can repeatedly observe the polar regions in several paths. The overlapped observations have a similar field of view owing to the design features of VIIRS; by its controlling of pixel growth rate at the edge of the scan to minimize the bow-tie effect. The VIIRS pixel size at the edge of the scan is 2.1 times that of the nadir pixel, while the MODIS pixel size at the edge of the scan is 4.8 times of the nadir pixel [25]. Note that the VIIRS albedo product is defined as the blue-sky albedo, which can be directly compared with the simultaneous ground measurements.

#### *2.2. The BRDF-Based Direct Regression LUT Generation*

The broadband instantaneous blue-sky sea-ice albedo is calculated from multispectral top-of-atmosphere (TOA) reflectance using an angular-specific linear regression relationship.

$$\mathfrak{a}(\theta\_{\mathfrak{s}}) = \mathfrak{c}\_{0}(\theta\_{\mathfrak{s}\prime}\theta\_{\mathfrak{v}\prime}\varphi\_{\mathfrak{s}}) + \sum\_{\mathfrak{b}} \mathfrak{c}\_{\mathfrak{i}}(\theta\_{\mathfrak{s}\prime}\theta\_{\mathfrak{v}\prime}\varphi\_{\mathfrak{s}})\rho\_{\mathfrak{i}}(\theta\_{\mathfrak{s}\prime}\theta\_{\mathfrak{v}\prime}\varphi\_{\mathfrak{s}})\_{\mathfrak{r}} \tag{1}$$

where *<sup>α</sup>*(*<sup>θ</sup>s*) is the broadband blue-sky albedo; *θs* is the solar zenith angle (SZA); *θv* is the view zenith angle (VZA); and *ϕs* is the relative azimuth angle (RAA). *i* (1,2,3,4,5,7,8,10,11) represent the nine VIIRS moderate resolution bands used in sea-ice albedo retrieval. *ρi* is the TOA reflectance from sensor data records (SDRs). *c*0 and *ci* are the retrieval coefficients; *c*0 is the constant term. The coefficients are stored in a pre-defined look up table (LUT) for evenly spaced angular bins in SZA, VZA, and RAA.

The pre-defined LUT trained from the representative dataset renders the algorithm highly efficient and accurate. The algorithm will ge<sup>t</sup> the coefficients for an actual (*<sup>θ</sup><sup>s</sup>*, *θ<sup>v</sup>*, *ϕs*) combination through linear interpolation in the surrounding angular bin, which runs fast in operational practices and avoids the discontinuity in neighboring albedo values. The LUT configuration covers SZA from 0◦ to 80◦ with an increment of 2◦, VZA from 0◦ to 64◦ with an increment of 2◦, and RAA from 0◦ to 180◦ with an increment of 5◦.

Figure 1 shows the data flow used to generate the VSIA LUT. The BRDF database is the basis for deriving TOA reflectance and broadband albedo, respectively. This database can be retrieved from satellite data for land surface pixels but has to be simulated from physical models for sea-ice pixels due to the lack of clear and low-SZA satellite observations in polar regions. A satellite pixel in the sea-ice zone was simplified as a combination of snow, ice, pond, and seawater [23]. The sea-ice BRDF can be regarded as a linear composition of the component BRDFs. The contributions of some distributed surface impurities containing dust and soot are expressed in the snow/ice BRDF model [23,26]. The BRDFs of the snow/ice component are simulated using the asymptotic radiative transfer (ART) model, with the input of inherent optical properties (IOPs) calculated from a variety of snow/ice physical parameters [26]. For ocean water, each BRDF is a linear combination of its three components (sun glint, whitecaps, and water-leaving reflectance from just beneath the air-water interface) [27]. The pond BRDF simulation in VSIA LUT deploys the analytical model proposed by Zege et al. [28] with the optical characteristic values referred to by Morassutti and LeDrew [29].

Sea-ice BRDF is considered as the linear mixing of these components' BRDFs. The fractions of different components in sea-ice pixels vary through time, thus the inherent heterogeneity of sea-ice BRDF should be considered in the simulation process. The Monte Carlo simulation method is used to generate samples of fractions in assembling the sea-ice BRDF for efficiency. The fraction of the first three components is determined by a uniform random number within 0~1. The fractions of the four components used sum to 1 in each BRDF item.

We generated a sea-ice BRDF dataset consisting of 120,000 simulated sea-ice BRDF items. The next key step is to generate the surface albedos. For each sea-ice BRDF item, the surface broadband albedo, including the black-sky albedos (BSAs) and white-sky albedo (WSA), are derived through an angular integration and narrowband-to-broadband conversion. Another method to achieve the sea-ice albedos is by aggregating the BSAs and WSA of each component using the simulated fractions. Its convenience embodies in the direct acquisition of snow/ice/seawater albedo from their BRDF models. The results of the two methods are consistent. The second one was used in our calculation.

The direct estimation algorithm was to directly infer surface albedo from TOA reflectance. Another key step is to simulate the TOA reflectance and diffuse skylight factor in each angular bin from sea-ice BRDF through atmospheric simulation using the 6S (Second Simulation of the Satellite Signal in the Solar Spectrum) tool. To eliminate the uncertainty resulting from atmospheric effects, multiple possible atmospheric conditions have been considered in the training data setup. Sea ice is mainly distributed in the Arctic, Southern Ocean, and Antarctic. The possible atmospheric influence in these regions can be represented in 6S using three predefined atmospheric models including sub-Arctic winter (SAW), mid-latitude winter (MLW), and sub-Arctic summer (SAS), while the typical aerosol types can be described using rural and maritime [30]. In practice, the atmospheric parameters are pre-calculated and stored in an atmospheric LUT. We transferred the surface reflectance spectra under each atmospheric condition type, which is the combination of an aerosol model and an atmospheric model. Then the number of simulated TOA spectra is much expanded than that of the surface BRDF spectra.

For the convenience of users, our retrieval object parameter was set as blue-sky albedo, which is defined as the ratio of up-welling radiation fluxes to down-welling radiation fluxes in a given wavelength range and is comparable with the in situ albedo observations. The blue-sky albedo is estimated using a linear relationship of black-sky albedo (directional-hemispherical surface reflectance) and the white-sky albedo (bi-hemispherical surface reflectance). The diffuse skylight factor has been

recorded in the atmospheric simulation and is used as the regulatory factor in combining the BSA and WSA into blue-sky albedo [31] as the regulatory factor (2).

$$
\alpha\_{bluc-sky} = (1 - \beta)\alpha\_{black-sky} + \beta\alpha\_{wlite-sky} \tag{2}
$$

It is assumed that the BRDF dataset is representative to all possible sea-ice pixels in the polar region. Based on the calculation results above, we built linear relationships between TOA reflectance and surface blue-sky broadband albedo for different solar/viewing angular bins using a least squares method. The assumption is that the regression relationship in each angular bin is applicable to various sea-ice surface types with a minimum overall square error over the BRDF database. The regressed coefficients form the sea-ice albedo LUT.

**Figure 1.** The data flow of VIIRS sea-ice albedo (VSIA) look up table (LUT) development. VIIRS: Visible Infrared Imaging Radiometer Suite; BRDF: bidirectional reflectance distribution function; TOA: top-of-atmosphere; BSAs: black-sky albedos; WSA: white-sky albedo; 6S: Second Simulation of the Satellite Signal in the Solar Spectrum.

#### *2.3. The Operational VSIA Algorithm*

The NDE albedo process consists of two components, as shown in Figure 2. The granule albedo, i.e., the published VIIRS albedo product, is computed online from a combination of the directly estimated albedo and a historical temporally filtered gap-free albedo; the historical albedo is computed offline from the granule albedo of previous days.

The instant retrieval of sea-ice albedo is calculated online from VIIRS TOA reflectance stored in the sensor data record (SDR), with latitude/longitude, solar zenith/azimuth angle, and view zenith/azimuth angle that are stored in a geolocation file. VIIRS NDE cloud mask EDR is used to remove the cloud contaminated pixels. The algorithm also needs the VIIRS Ice Concentration intermediate product (IP) [32] to distinguish the ocean pixels covered by ice.

Instantly-retrieved albedo values are only available for clear-sky pixels which results in gaps in the online generated albedo due to cloud contamination. Therefore, an offline process was designed to provide noise-eliminated albedo for online gap-filling. The offline branch deploys a temporal filtering algorithm [33], which deduces the albedo from a time series of instant retrievals (including the previous 8 day retrieval plus the current day) based on maximum likelihood estimation. Albedo climatology acts as a critical static input to provide statistics about inter-annual trends and variation of albedo for weight calculation and backup value estimation.

To solve the geometric registration problem between multi-day historical albedo images, temporal filtering is conducted on a fixed grid of sinusoidal projection. The global grid is evenly divided into 72 (horizontal) by 72 (vertical) tiles. After online instant retrieval, the albedo granules will go through a forward granule-to-tile transformation prior to offline filtering and a backward tile-to-granule transformation subsequently. The online albedo after gap-filling with the filtered albedo can provide continuous and reliable sea-ice albedo for users.

**Figure 2.** The flowchart of operational S-NPP Data Exploration (NDE) VIIRS sea-ice albedo product (legends shown at the bottom of the figure). S-NPP: Suomi National Polar-orbiting Partnership; EDR: environmental data record; IP: intermediate product; LSA: VIIRS Surface Albedo Product.

#### **3. Validation of VSIA Using Ground Measurements**

### *3.1. Ground Measurement Dataset*

In situ albedo observation of sea ice is expected to verify the algorithm. However, this type of measurement is rare and its match-up with VIIRS observations is more so. Therefore, we firstly used Programme for Monitoring of the Greenland Ice Sheet (PROMICE) and Greenland Climate Network (GC-NET) measurements as substitutes. Greenland is the world's second-largest ice mass. Its ice sheet has existed for more than 2.4 million years and covers more than 80% of the island area. Considering the comparability of samples from different seasons, the match-ups around local solar noon (13:30~15:30 UTC time) were used for the comparison. Only the clear-sky observations were retained. The term "clear-sky" is used to distinguish a sky without cloud influence. The cloud effect on the radiation budget, referred to as 'cloud forcing', is one of the largest uncertainties in radiative transfer models. It is acceptable to screen out the cloudy-sky days before the comparison.

We also compared the VSIA with some in situ datasets from existing studies. Although the sample size is limited and the measurement footprint is not enough to catch the spatial heterogeneities within the satellite pixel, the comparison would embody the performance of VSIA in actual sea-ice surfaces.
