A method to derive satellite PAR albedo time series over first-year sea ice in the Arctic Ocean

Contrary to previous assumptions, under-ice phytoplankton blooms can be significant (Lowry et al., 2014), especially in the context of a thinner, more transparent Arctic sea ice cover (e.g., Horvat et al., 2017). To identify an under-ice bloom on a rather local scale, ice camps (Mundy et al., 2014; Oziel et al., 2019), oceanographic cruises (Arrigo et al., 2014; Randelhoff et al., 2019), ice-tethered buoys (Laney et al., 2014; Hill et al., 2018), bio-Argo floats (Mayot et al., 2018; Randelhoff et al., 2020) and sediment traps (Fortier et al., 2002) can be used. To recognize the extent of under-ice blooms at pan-Arctic scale, however, the in situ data from these cumbersome, expensive observational approaches is insufficient. Over synoptic spatial scales, a detailed description of abiotic factors controlling phytoplankton, such as the under-ice photosynthetically available radiation (PAR), could provide a better understanding of the changing Arctic Ocean ecosystem (Arrigo et al., 2012; Lewis et al., 2019).

From a radiative transfer perspective, without a full description of the above-surface radiance distribution and the inherent optical properties of sea ice, accurately estimating under-ice PAR is difficult. A solution envisioned by several researchers (e.g., Assmy et al., 2017; Horvat et al., 2017; Stroeve et al., 2021) is to estimate the above-surface plane irradiance, reduce it by the energy leaving the sea ice surface (conveyed quantitatively by the albedo) and apply a log-linear decrease of energy with depth, scaled by the thickness of the medium. Surface albedo is defined as the ratio between solar radiation reflected by Earth’s surface and solar radiation incident on the surface. Optically, it summarizes the interaction between the surface-incident radiation and everything at and beneath the atmospheric lower boundary. It is a key variable for determining downward irradiance above (Grenfell and Perovich, 2004; Laliberté et al., 2021) and below sea ice in the Arctic Ocean (Ebert et al., 1995; Nicolaus et al., 2010). Therefore, mapping PAR under sea ice and eventually assessing under-ice primary production at a large scale over the Arctic Ocean requires robust estimates of surface albedo based on remote sensing data.

A coarse solution to obtain PAR albedo over large scales is to use empirically fixed PAR albedo values for given surface categories (Castellani et al., 2017; Horvat et al., 2017; Lange et al., 2017). But this approach implies determining relative areal proportions of the different surface categories from grid cells much larger than the surface features (Rösel et al., 2012; Wright and Polashenski, 2020), leading to potentially important uncertainties. It also generally assumes too little spectral variation per surface category (Tschudi et al., 2008; Yackel et al., 2018). In parallel, dynamic semi-empirical approaches to derive seasonal patterns in sea ice evolution modulated by the timing of events (Perovich and Polashenski, 2012; Arndt and Nicolaus, 2014; Laliberté et al., 2021) can be used with satellite-retrieved information based on, for instance, polarized backscatter (Perovich et al., 2007) and a library of representative albedo values (e.g., De Abreu et al., 1995) to obtain the seasonal evolution of visible albedo. While suitable over large spatial and temporal scales, the spatial resolution of scatterometers or passive microwave radiometers (generally 25 km) would preclude the use of this approach near land where the signal would be contaminated. We also anticipate that, in many cases, 25 km may be too large to examine the relation between the icescape, phytoplankton drift and the dynamic of under-ice phytoplankton blooms. Finally, a more direct route is to estimate albedo using mature satellite-derived broadband albedo products, specifically available for the Polar Oceans, such as the Extended AVHRR Polar Pathfinder dataset, computed twice daily from AVHRR channels 1 and 2 (Key et al., 2016). These existing operational albedo products derived from remote sensing data have been developed to address heat flux at the Earth surface and related thermodynamic processes. Therefore, they are defined for the full spectrum of solar energy, between 300 and 3000 nm. The spectral range that is relevant to primary production by phytoplankton under sea ice, however, is the range of visible wavelengths (400–700 nm). Typically, when trying to assess visible light available under Arctic (Stroeve et al., 2021) or Antarctic (Pinkerton and Hayward, 2021) sea ice, the method consists of multiplying these broadband shortwave estimates by a conversion factor to reduce the wavelength range to the visible band. A conversion of a broadband albedo to a visible one is possible, but it is highly sensitive to atmospheric conditions and surface physical properties (Grenfell and Perovich, 2008), adding avoidable uncertainties to these estimates.

The use of visible spaceborne sensors should be a good alternative approach. Visible satellite imagery typically has high spatial resolution, sufficient radiometric signal-to-noise ratio, and a large swath to provide overlaps of adjacent orbits towards high latitudes to help overcome the prevalent cloud cover (Figure 1). As such, the mature and widely used moderate resolution imaging spectroradiometer (MODIS) land spectral albedo product (Schaaf et al., 2002) is built upon clear-sky atmospherically corrected surface reflectances from the combined Aqua and Terra observations. However, this approach cannot provide an estimate of sea ice albedo in the visible range at a relatively high spatial resolution because of the highly dynamic Arctic surface. Because of sea ice drift, changes in temperature, and precipitation, the extreme seasonal cycle of sea ice albedo over the Arctic Ocean is characterized by short periods of drastic change and is overall more variable than any other planetary surface. The continuous evolution of surface properties precludes the collection of sufficient multiangle observations over each pixel that would constrain a retrieval of albedo via inversion of a multidirectional anisotropy model (Lutch et al., 2000).

In this study, we used the satellite-derived surface reflectance product that is input to the operational MODIS spectral albedo algorithm. However, instead of generating albedo for a wide diversity of planetary surfaces and over a variety of spectral and broad bands, we focused on Arctic sea ice and the spectral range relevant to primary production modelling within the limited sun and viewing geometries of the polar latitudes. First, we present a brief theoretical background to describe the relation between surface reflectance and albedo. Second, we describe the model and the data used to develop it. Third, the model is evaluated with field measurements and is applied to a specific region in the Canadian Arctic. Finally, the model limitations are discussed.

We used the surface reflectance as a starting point to derive albedo. It is a directional quantity that is routinely derived by satellites. For any geophysical surface and given illumination, the radiance distribution leaving that surface is expected to emerge with some degree of directional irregularity (Goyens et al., 2018). Surfaces with such reflectance properties are said to be anisotropic. In contrast, under the same illumination, a surface with isotropic reflectance properties would reflect radiance equally in all directions. For our purpose, we defined the surface reflectance (SREF) as the ratio of sea-ice-leaving radiance $Lr(θi,ϕj)$ to sea-ice-incident irradiance E_d (W m^–2) normalized to a perfectly reflecting isotropic surface (⁠$1/π,sr−1$⁠):

In Equation 1, L_r (W m⁻² sr^–1) is the radiance in the upward direction at given zenith (θ_i) and azimuthal (ϕ_j) angles (Nicodemus et al., 1977). The subscripts i and j refer to a specific narrow direction for which the radiance is assumed homogeneous within a small finite field-of-view. Surface reflectance is a widely used measurable quantity (as opposed to conceptual; see Schaepman-Strub et al., 2006) often named hemispherical-conical reflectance factor (HCRF, dimensionless). Surface reflectance may refer to a particular satellite band, in which case

with $RSR(b,d,λ)$ as the relative spectral response for band b and sensor d. The relative spectral response indicates the quantum efficiency of the sensor at a particular wavelength.

Our goal was to achieve PAR albedo estimates from the satellite-derived surface reflectance. The albedo (α, dimensionless) is a compilation of the surface reflectances from all directions in a given band or spectral interval. If the albedo is computed for an isotropic surface, a surface reflectance of any direction is equal to its albedo. The albedo is defined as the ratio of upwelling planar irradiance E_u (W m^–2) to downwelling planar irradiance, and is computed as an angular weighted average of the surface reflectances over all discrete quadrilateral regions $ΔθiΔϕj$ of the upper hemisphere:

Here, radiances are approximated by non-equal regions of assumed constant fluxes spanning the entire upper hemisphere, and their conversion to a planar irradiance requires accounting for their projected area. It follows that the angular weighted sum of all surface reflectances approximates the albedo. Spectrally, the wavelength-integrated PAR albedo (⁠$αPAR$⁠) is defined as:

with $Ed(λ)$ as the spectral irradiance incident on the surface. To account for the variability of the reflectance from sea ice surfaces based on solar and viewing geometries (both known for any MODIS pixel), an anisotropy reflectance factor (ARF) can be defined by normalizing the surface reflectance (SREF in Equation 1) with the albedo (Equation 3). This normalization allows us to assess the angular and spectral dependence of contrasted melting surfaces using a common scale. The ARFs take on values below unity, equal to unity and above unity, depending on the viewing angle. Values equal to unity imply that the surface reflectance in that viewing direction is numerically equal to the albedo. An ARF value below (above) unity implies that the surface reflectance that is being recorded in that direction is numerically smaller (greater) than what would be recorded if the surface was isotropic; using the surface reflectance from that direction and simply assuming that the surface is isotropic would lead to an underestimation (overestimation) of the albedo.

Because reflectance from the sea ice surface is anisotropic, we needed a model to tie MODIS surface reflectance to the albedo. Each available satellite viewing geometry represents a specific observation of surface conditions and sample of anisotropic sea ice surfaces. Depending on the timing and the location, the surface could be composed of melting snow, bare ice and/or melt ponds. For a given reflectance measurement, we have no way of knowing the composition of the surface element mixture that was sampled, so we cannot apply an anisotropy correction based on a surface typology.

We solved this inverse problem using a well-defined and restricted approach. First, rather than deriving albedo for a wide range of planetary surfaces and spectral domains, we focused on a single surface cover class and the visible range. Second, we set a narrow frame of measurement conditions in which the departure from isotropy is small. In fact, our model only needs to represent a subset of surface anisotropy matching the sun and viewing geometry of the Arctic scanned by MODIS, in addition to the small spectral dependence of the surface over the visible range (discussed in this section).

When the MODIS albedo (code MCD43) is produced over land using the operational algorithm of Schaaf et al. (2002; 2011), an albedo retrieval is computed over sixteen days of satellite observation. On global scale, this was deemed the best tradeoff between the number of angular samples and the time for which a surface is not changing (Schaaf et al., 2011). The Arctic Ocean imposes exceptional constraints, mainly the prevalent cloud cover and the rapidly changing surface. Accordingly, we formulated an empirical model that converts each clear-sky remote sensing surface reflectance to albedo.

From a remote sensing perspective, to consider the orientation of surface roughness features was not realistic, especially as sea ice is potentially drifting around with winds and currents. Thus, Arctic sea ice surfaces, unlike ocean waves or wind-sculpted snow often found in Antarctica (Warren et al., 1998; Hudson et al., 2006), were simply assumed to be azimuthally homogeneous surfaces and, consequently, only the relative azimuth angle between the sun and the satellite was considered. Building on this assumption, we could rearrange the viewing geometries using only the relative azimuth angle with the Sun positioned at 180°. We observed that throughout the melt season at the northernmost latitudes both MODIS sensors have a limited domain of angular sampling (Barnsley et al., 1994; Privette et al., 1997). Over 66.6° N, between May and September, MODIS sweeps from across track scanning display only small variations in the azimuthal viewing geometries (Figure 2). This pattern of viewing geometries reduces the distribution of ARF values and constrains the conversion of surface reflectance to albedo.

The solar zenith angle in the Arctic (over 66.6° N) is always greater than 43°. The angular and spectral dependence of surface reflectance is expected to increase for large solar zenith and satellite viewing zenith angles (Perovich, 1994; Bourgeois et al., 2006). Additionally, at large solar zenith or satellite viewing zenith angles (towards 90°), we expect an increase in error due to a longer atmospheric path radiance and reduced at-sensor signal-to-noise ratio. For these reasons, we only examined situations with a solar zenith or satellite viewing zenithal angle less than 65° and 60°, respectively. The latitudinal limitation and avoidance of long optical paths thus also constrains the ARF values.

The Arctic exhibits a complex mixture of surfaces. For practical reasons, we decomposed Arctic sea ice into discrete surface categories. Perennial sea ice has a different physical structure than first-year ice because of the summer melt (Grenfell and Maykut, 1977). It is also thicker. In view of the distinctively complex topography and limited (and decreasing) spatial coverage (Kwok et al., 2009; Comiso, 2012), we decided to focus on first-year ice, again restricting the possible ARF values. Following Perovich et al. (2002), for a simplified representation of the seasonal evolution of first-year sea ice, four dominant categories were identified: dry snow, melting snow, bare ice and melt ponds.

The operational MODIS surface reflectance product is distributed for three visible bands not saturating over sea ice. MODIS also records in six visible ocean color sensor bands, but these bands have a high sensitivity that makes them reach digital saturation over sea ice surfaces. Figure 3 shows the spectral distribution of the selected three MODIS bands in the visible range i.e. bands 01 (613–682 nm), 03 (451–481 nm) and 04 (538–569 nm). These bands are well distributed to cover the spectral range that is of interest for studying photosynthesis. Aqua and Terra sensors have very similar spectral responses (Figure 3). Thus, we assume that Terra’s spectral response equals Aqua’s spectral response. Overlaid are the measured albedos of snow, bare ice and melt ponds surfaces measured at Resolute Bay, Canada, from Perovich (1994). These representative spectra are added to illustrate how smooth the albedos of sea ice surfaces are over the visible range (no sharp spectral features). This smoothness supports our premise that three sensor bands should capture most of the information contained in the visible range. In fact, both the magnitude of the signal and the spectral shape of the surface carry relevant information that is used by our model. For algorithm development, we had to represent the surface reflectances of the available satellite bands by mimicking surface reflectances from high spectral resolution data (in situ and modeled; see below) scaled to the satellite spectral response functions.

So far, we have stated that a model to convert from surface reflectance into albedo should yield an albedo estimate for every surface reflectance observation. This model should be developed to be valid under specific sun and viewing geometries (azimuthal and zenithal), applicable to a limited number of surface categories, and based on the current satellite band availability. As with any empirical model, our model will be limited by the amount of high-quality data used to develop it, and limited Arctic sea ice data exist. Even if in situ surface reflectance generally represents a valuable ground truth, it contains measurement uncertainties due to the intrinsic complexity involved, and, for the sake of algorithm development, covers too limited a range of variations in illumination conditions and in terms of surface categories (described in Section 2.3.1). Thus, physical descriptions of sea ice reported in the literature for our typical surface categories were used to compute a more diverse synthetic sea ice surface reflectance dataset simulated through a radiative transfer model (described in Section 2.3.2).

In this section, we first present the data used to develop our algorithm. We used in situ surface reflectance because, unlike models, in situ data have no chance of suffering from possible errors in computer programs or from the possible omission of the relevant physics. We used models because they are immune from possible uncertainties on calibration, temperature variations, gain or roll-off. A combination of both leverages their advantages and diminishes their weaknesses. We also present an independent in situ dataset that will be used to evaluate our model, as well as the remote sensing surface reflectance that will be input to the algorithm, namely the operational MODIS surface reflectance product. Finally, we present a set of ancillary data from a specific site that helps to verify if the PAR albedo presents realistic values when derived from remote sensing data.

The Green Edge field campaign took place at a fixed location on landfast first-year ice in Western Baffin Bay (67°28’ N, 63°48’ E; Figure 4 inset) near Qikiqtarjuaq along the east coast of Baffin Island, Canada. This area is where we will demonstrate the ability of the model to predict albedo. In consideration of the land surrounding the sampling site, we chose to define our study area as a 9-km² zone, which corresponds to 6 x 6 MODIS surface reflectance pixels as shown in Figure 4. The study period was from spring conditions (beginning of May, about –8°C) to ice break up (mid-July, about 10°C) for 2015 and 2016. Because the field sampling strategy was set to characterize the phytoplankton spring bloom dynamics, many light-related variables were collected on site, including key observations of in situ surface reflectance.

A dataset of highly resolved angular and spectral resolution field surface reflectance was derived from measurements taken during the 2015 Green Edge Campaign (Goyens et al., 2018). The complete dataset was generated using a regression analysis and includes:

1. measurements from a camera of hyperangular resolution (CamLum; Antoine et al., 2013) of 16,020 quadrilateral regions of uniform flux ratios at 1° resolution in θ and at 2° resolution in ϕ over the upper hemisphere;

2. measurements from a field spectroradiometer (Analytical Spectral Device) of hyperspectral resolution subsampled to give contiguous wavelength bands of 5 nm between 400 and 700 nm; and

3. three typical Arctic surface types measured under similar sun zenith angles, namely, dry snow (sun zenithal angle at 56°), bare ice (sun zenithal angle at 60°) and melt ponds (sun zenithal angle at 57°).

Full description of the method used to generate the dataset can be found in Goyens et al. (2018) and is available online (Massicotte et al., 2020). In brief, a multispectral circular fish-eye radiance camera (measuring radiances simultaneously in all directions) was combined to a hyperspectral field spectroradiometer guided by a metal arch (measuring surface leaving radiances successively with a goniometer for the angular increments). These in situ surface reflectances were scaled with the MODIS spectral response function (https://mcst.gsfc.nasa.gov/) to obtain three weighted band averages (Equation 2). To filter out part of the fine-scale irregular features emerging from a limited sample size (less relevant to a general solution for Arctic scale satellite remote sensing), the dataset was resampled to coarser 10° x 15° angular bins.

Viewing zenithal angles were limited to ≤ 60°. Towards larger angles, contamination from surrounding surfaces is hard to prevent. The lack of large viewing zenithal angles in the goniometer/camera-based dataset from Goyens et al. (2018) precluded the computation of the true albedo from surface reflectances (Equation 3). Given that the contribution of radiance at viewing zenith angles below 60° accounts for 75% of upward irradiance, the computed albedo from this dataset is in reality an approximation of the albedo. A total of three albedo values were computed from the in situ data.

Simulated surface reflectances from radiative transfer simulations were used to reproduce several environmental conditions. Complementing the in situ data, they were also used to train our algorithm as they simultaneously provide both the surface reflectances and the PAR albedo. Assuming that the atmosphere-sea ice-ocean system is sufficiently horizontally homogeneous, AccuRT (Hamre et al., 2017; Stamnes et al., 2018) should be a suitable radiative transfer tool for simulation of various sea ice surface reflectances, as it is built on the well tested and extensively used discrete-ordinate radiative transfer method DISORT (Stamnes et al., 1988), which has been developed to handle media of different refractive indices (Jin and Stamnes, 1994), with additional development of parameterizations of atmosphere, snow, sea ice, and ocean inherent optical properties (Hamre et al., 2004; Stamnes et al., 2011; 2018).

The fully coupled atmosphere-sea ice model was used to run 36 simulations (with 64 streams) over the visible range at 5-nm spectral resolution, with 10° angle resolution in θ, 15° resolution in ϕ, and all four surface types. The upwelling radiances are resolved at non-equal-area quadrilateral regions of constant fluxes over a hemisphere of 2π (surface and sky domes). We computed the downward and upward irradiances above the surface, and the upward radiances from 0 to 90° in θ and 0 to 360° in ϕ, to derive an equivalent physical quantity to what is provided by the MODIS surface reflectance product. These synthetic, hyperspectral simulated surface reflectance values were also scaled by the MODIS spectral response function.

The rationale was to represent these four distinct categories (dry snow, melting snow, bare ice, and melt ponds) under natural illuminations and typical physical variations. Three sun zenithal angles (45°, 55° and 65°) were chosen to be representative of relevant satellite observations for reasons established in the previous section. An oceanic aerosol model was used (Ahmad et al., 2010; IOCCG report, 2015). For each snow and ice medium, two layers were present. Realistic sea ice physical values were obtained from Light et al. (2003), Hamre et al. (2004), Marks and King (2014), Taskjell et al. (2017) and Oziel et al. (2019). The main features of the simulations are summarized as follows. For the dry and melting snow cases, the fine and coarse (respectively) geometrical thickness of the snow grain layer was changed from large to small (0.10 m, 0.05 m and 0.025 m) and overlaid by a constant 0.01-m layer of fresher snow, covering the transition at the melt onset. This approach resulted in a transitional decrease in the portion of light scattered upward (Warren, 2019). The main characteristic for the dry snow case is the presence of a snow with a density of 320 kg m^–3 and an effective snow grain radius of 300 µm. The snow in the melting snow case had the same thickness but a density of 400 kg m^–3 and an effective grain radius of 1000 µm. For the bare ice, the ice thickness was changed from thicker to thinner (0.95 m, 0.50 m and 0.30 m). On top, a surface scattering layer of 0.02 m was assumed. The surface scattering layer (here represented by old snow) had a density of 500 kg m^–3 and an effective grain radius of 1200 µm. This top layer was meant to dominate the scattering budget so that little variation in the reflected light was observed (Perovich, 2005). Finally, for the 0.10-m deep melt ponds, ice thickness was varied (0.95, 0.50 and 0.30 m) because the light reflectance depends very little on water thickness (Makshtas and Podgorny, 1996).

We used in situ and simulated surface reflectances to train our model, but it was desirable to have an additional collection of independently measured surface reflectances of sea ice to evaluate the performance of our model. Perovich (1994) examined five cases of first-year sea ice with a hyperspectral spectroradiometer guided by a goniometer at Resolute, in the Canadian Archipelago. The objective for making these field observations was to help with the interpretation and application of directional satellite measurements made over a variety of ice conditions. During the melting season, surface reflectance and corresponding albedo measurements were performed over dry snow, snow with a glazed surface (glazed snow), bare ice, blue ice and melt ponds. Each surface type was assessed by a single set of scans with a sun zenithal angle around 58° (see table 1 in Perovich, 1994). Measurements for the viewing geometry were made from 0 to 60° with 10° increments in θ and from 0 to 360° with 30° increments in ϕ. As formulated above, only the restricted zenithal viewing angles and MODIS-specific azimuthal range were of interest to us for the evaluation (Section 2.2; Figure 2). Because the irradiances and PAR albedo were not included in the dataset, the PAR albedo could not be computed with Equation 4, but it was approximated by the mean of the spectrally resolved albedo measurements of the visible interval (approximately 3-nm spectral resolution; shown in Figure 3). Using the simulated surface reflectances, we compared this approximation with the values derived from Equation 4 and found the approximation to be a reasonable assumption.

Since we are mostly interested in the performance of the model when applied to remote sensing data, we also used this evaluation dataset to assess the impact of atmospheric corrections. In theory, the remote sensing surface reflectance is corrected for atmospheric scattering and absorption and is supposed to represent what would be measured directly in the field (see next section), but this process is non-trivial. In fact, the level of accuracy of the atmospheric corrections for the surface remote sensing reflectance product is estimated to be typically ± 0.005 + 5% * surface reflectance (Vermote et al., 2015). Using that range as an imperfect estimate of uncertainty for surface reflectance, represented as a 99% confidence interval under the assumption that this error is normally distributed, we propagated the atmospheric correction uncertainties in the evaluation dataset of surface reflectance and estimated the resulting uncertainty in the albedo estimates.

We compared the predicted PAR albedo (M_i) to the observed PAR albedo (O_i) and computed the mean absolute error (MAE; Equation 5), the root mean square error (RMSE; Equation 6) and coefficient of determination (R²) metrics.

These metrics were chosen for their simplicity and different representativeness of the distribution of error in the model. Both Equations 5 and 6 present average error metrics in the same units as the PAR albedo, with the RMSE being more sensitive to outliers. The coefficient of determination represents the variation in the measured PAR albedo that is explained by our predictive model and is unitless.

The MODIS surface reflectance level 2G-Lite daily products (Aqua MYD09GA and Terra MOD09GA) include three 500-m spatial resolution visible bands (introduced in Figure 3). These geophysical products are corrected for the effects of atmospheric scattering (aerosols, thin cirrus clouds) and absorption (aerosols, ozone and water vapor) and thus aim to represent instantaneous estimates of surface reflectance that would be measured at the Arctic surface. In other words, the numerical values derived from remote sensing are meant to have the same physical meaning (SREF; Equation 1) than our in situ surface reflectance (Section 2.3.1), simulated surface reflectance (Section 2.3.2) and evaluation dataset of surface reflectance (Section 2.3.3).

We considered only the best quality satellite observations, defined as follows: 1) cloud masking was performed by MODIS MOD35 product in addition to the surface reflectance product internal cloud algorithm; 2) values flagged as none-high band quality (see Vermote et al., 2015, their table 10) were discarded; and 3) pixels with presumed cloud shadow or high aerosol loads, or having sun zenithal angles > 65° and viewing zenithal angles > 60°, were removed. Note that the in situ data (Section 2.3.1) were acquired with sun zenithal angles around 60° but the simulated data (Section 2.3.2) extended to sun zenithal angles of 65°, which makes us confident our algorithm can process satellite-derived surface reflectance with sun zenithal angles < 65°.

Combining Aqua and Terra sweeps across the viewing hemisphere regularly provides up to 14 observations of the Green Edge study area per day (Figure 1). In this investigation we used all satellite passes (layers), regularly yielding multiple individual viewing and illumination geometries in a single day (Figure 2). To derive area-averaged seasonal time series of PAR albedo over the study area, we used a 6 x 6 pixel spatial window of Aqua and Terra daily surface reflectance products centered on the Green Edge study area (Figures 4 and 5). Satellite passes were retained if, within the spatial window, they contained more than three valid pixels. Pixels were retained if, within an entire day, they had a coefficient of variation (standard deviation divided by the mean) less than or equal to 25% for the surface reflectance values of each band.

Other variables were compiled to provide a context for local PAR albedo fluctuations. For the Green Edge study area, spectral albedo data were acquired on site in the spring of both years (i.e., 2015 and 2016; Verin et al., 2019). SOLALB (a custom-built radiometer) measurements at 3-nm spectral resolution (Picard et al., 2016) were collected at random locations over the study area every second day using a metallic arm at 0.80 m above the surface. Integration times were automatically adjusted for the 10 spectra acquired in triplicate at each site (snow pits and transects, 5-m intervals over 100 m to 150 m). Daily averages of PAR albedo were computed from this dataset. Melt pond fraction was derived using a time-lapse of photographs taken from a high vantage point North of the study area (surrounding hill of elevation approximately 310 m; Figure 4). Melt pond fraction at the ice camp was derived based on a classification similar to that used by Perovich et al. (2002), Katlein et al. (2015), and Taskjelle et al. (2017). All photos (30-min temporal resolution, daylight) with a clear near-surface atmosphere were RGB-transformed and cropped to cover approximately the 9 km² area of interest. A fixed threshold of 0.5 in RGB intensity was chosen to distinguish melt ponds from bare ice. Melt pond fraction was averaged every second day and interpolated to obtain daily values. Finally, rain and snow precipitation data were summarized from hourly to daily from Qikiqtarjuaq airport (Figure 4) located 12 km away (provided by Environment and Climate Change Canada).

Altogether, the in situ (Section 2.3.1) and simulated datasets (Section 2.3.2) cover a large array of sea ice surface reflectances and corresponding PAR albedos. This dataset was used to train a model that converts remote sensing surface reflectances (Section 2.3.4) to PAR albedo. The technical details of the model are presented here.

First, a time series of n remote sensing data was retrieved and added to the training dataset to form a single matrix. This matrix has seven columns: surface reflectance in band 03 (⁠$SREF3$⁠), surface reflectance in band 04 (⁠$SREF4$⁠), surface reflectance in band 01 (⁠$SREF1$⁠), sun zenithal angle (θ_i), viewing zenithal angle (θ_v), relative azimuthal angle (Δϕ) and PAR albedo (α_PAR). It has a total of m rows, with k rows from the training data and n rows from remote sensing data $(m=k+n)$⁠. At this point, the PAR albedo is not available for the n rows from remote sensing data. Second, we selected only the six columns representing the three surface reflectance bands and their three geometrical configurations. This subset matrix has the dimensions of m x 6. We centered and scaled the subset matrix and computed its principal components. In the new coordinate system, the m rows are representing the scores of the six principal components. Third, the first k rows of the scores (associated to the training data) were regressed against the known PAR albedo to form a multilinear regression with seven coefficients (one intercept, B₀, and six coefficients, $Bw=1,…,6$⁠). To avoid overfitting the data, a subset of the coefficients was selected with a stepwise regression process based on the Akaike information criteria (Akaike, 1973). Fourth, the predictive model of the form

was used to map the unknown PAR albedo (⁠$αPARq$⁠) for the remaining scores (⁠$PCwq$⁠). A single observation $(q=1,…,n)$ can represent more than a single surface category, such as bare ice and melt ponds. The predictive model can be applied to linear mixtures of the represented cases used for the development of the model. However, the model will likely be making unreliable predictions for the cases outside of the scope for which it was developed. Accordingly, we rejected any inversion made from Equation 7 that would fall outside the range [0.2, 0.98] (i.e., possibility of open water and unrealistically high values).

All instantaneous visible albedo estimates $(αPARq)$ were then converted to a reference instantaneous albedo at a sun zenith angle of $θ0=65°$⁠, so that differences in albedos are ascribable to different atmospheric or surface conditions and not to different times and latitudes. Following Lindsay and Rothrock (1994), a normalization factor based on the function $h(θ)$ was used as

with

such that $h(θ0)≈0.67$⁠. In our work frame, 65° is the maximum sun zenith angle and the normalization factor equals one. However, this normalization factor can increase the albedo up to 7% of its original value, should the observation be made at the minimum sun zenith angle of 43°.

In this section we present details of the data that were used to build our predictive model.

We first consider the albedo spectra gathered to train our model to convert different types of Arctic first-year sea ice surface reflectances.

Albedo spectra derived from in situ measurements (dotted lines in Figure 6) differ slightly from those derived from radiative transfer simulations (solid lines in Figure 6). In the different groups of simulated albedos, variations of the sun zenithal angle and thickness (dry snow thickness, melting snow thickness, bare ice thickness and melt pond thickness) creates the spread between the different lines. Most simulated albedo spectra agree well with measurements found in the literature. For example, snow and melting snow are similar to those found in Grenfell and Maykut (1977) in their figure 1. With snow metamorphism, snow grain size increases, the optical path increases and the albedo decreases with time (e.g., Wiscombe and Warren, 1980). When the snow depth decreases, the underlying darker ice becomes apparent. Albedo increases with solar zenithal angle as the probability that a photon is scattered out increases when traveling at an angle close to that of the surface plane. Bare ice and melt ponds albedo are very similar to spectral albedos presented by Light et al. (2015) in their figure 4. Bare ice albedo decreases with decreasing underlying ice thickness, with the influence of underlying ice and seawater remaining minimal because of the ice surface scattering layer (Light et al., 2008). Even if the magnitudes of the bare ice spectra are generally comparable, there is a difference in spectral shape between in situ and modeled albedos, which is potentially due to a difference in amounts of impurities. Melt pond albedo decreases with decreasing underlying ice thickness (the ocean becomes increasingly apparent). Bare ice albedo and melt pond albedo have the smallest and the greatest albedo variability in the dataset, respectively. Overall, our simulated sea ice visible PAR albedo ranges between 0.95 and 0.2.

The ARF values derived from in situ surface reflectance for the three surface types and three MODIS spectral bands are shown in Figure 7. In these and all other polar plots presented in this paper, a common azimuthal reference has been adopted, with the viewing azimuthal angle pointing towards the Sun at 0°, and therefore the Sun is positioned at 180°. This azimuthal direction will be referred to as the principal plane, and the direction between 90° and 270° will be the perpendicular plane.

For all surface types, the largest in situ ARF values were observed in the upper part of the plots (viewing angle pointing towards the Sun). More generally, deviations from isotropy were more prominent along the principal plane while the view angles close to the perpendicular plane were more stable. The snow-covered surface showed smooth angular variations in ARF values for all bands. A notable feature was that within the circle circumscribing the 30° view angle (near-nadir viewing geometry), the variations were minimal and the snow ARF values were close to one. Despite more heterogeneous angular variations (speckles), which are typical of such surfaces, bare ice and melt ponds showed similar angular patterns, again with the maximum ARF values recorded when the viewing angle was pointing towards the Sun direction. A mixture of thawed ice and highly reflective ice grains caused this small-scale heterogeneity (see figure 4 in Goyens et al., 2018). The more extreme angular dependence found for bare ice and melt ponds resulted from fewer scattering events (Coackley, 2003) and internal reflections (Jiang et al., 2005; Smedley et al., 2020). For bare ice, the lower left and the lower right parts of the polar plot differ because of non-uniformly distributed fine-scale surface features. Also, because melt ponds were surrounded by a porous surface scattering layer, slightly higher ARF values were observed all around at the higher viewing zenith angles. All of these measurements are demonstrably more heterogeneous than what is typically computed from models with similar sea ice optical properties but with a surface that is perfectly homogeneous (Hudson et al., 2006). We suspect that this heterogeneity would be averaged out if the number of sites for data collection were increased.

For the ice medium, the spectral dependence of scattering in the visible wavelength interval is negligible and the absorption coefficient increases almost log-linearly with the wavelength (Warren and Brandt, 2008). As a consequence, for a given illumination condition and a given surface, the mean path length of the outgoing visible light decreases with the increasing wavelength. Therefore band04 (545–565 nm) is at an intermediate stage of anisotropy between band03 (459–479 nm) and band01 (620–670 nm), with band01 as the most anisotropic one.

Figure 8 presents the ARF values from our simulated dataset generated using the AccuRT model, showing four surface categories (dry snow, melting snow, bare ice and melt pond), for a single set of physical properties, and the medium solar zenith angle (55°).

Figure 8 shows a similar pattern as Figure 7, with a prominent hotspot in the viewing direction that opposes the azimuthal sun direction. However, this time, the surface reflectances are symmetrical with respect to the solar plane. As with the in situ data, the diverging ARF values increase with sun zenithal angle, wavelength and melt progression. Within this restricted zenithal range, the peaks in melt pond reflectance are weaker compared to that of in situ data. Although the values used for the simulations are not meant to mimic the conditions in which the in situ data were measured but rather to represent typical Arctic sea ice under additional environmental conditions, the ARF plots show comparable anisotropy patterns.

We now specifically address the range of variation in ARF that is actually relevant to the reflectance measurements made with the MODIS sensors over the Arctic Ocean. In Figures 7 and 8, we superimposed onto the ARF polar plots the MODIS viewing geometry observed at the Green Edge study area (largely representative of the pattern found elsewhere in Arctic; see Figure 2) and at any time when the sun zenith angle is less than or equal to 65°. Overall, the azimuthal angle ranges between 50° and 80° (230° and 260°) for Aqua and from 100° to 130° (280° and 310°) for Terra, lying closer to the perpendicular plane than the principal plane. From a practical viewpoint, this angular domain designates the range of variation in ARF that must be addressed in an effort to estimate the PAR albedo. It is much narrower than the full range of variation in ARF. The specific satellite geometric sampling avoids measuring the maximum reflection from the sun, which is centered on the forward direction (180° relative azimuth angle between sun and sensor, principal plane), as a result of the prevalent forward scattering within snow and ice, and to which specular reflection is added on melt ponds.

Figure 9a shows the frequency distribution of in situ ARF of the snow surface for the restricted zenithal range and for the specific azimuthal range of MODIS. We defined the restricted zenithal range as all ARF values having zenith angles below 60° (azimuthal range 0° to 360°). The specific azimuthal range is a subset of the restricted zenithal range corresponding only to the possible azimuthal viewing directions of MODIS-Aqua and MODIS-Terra in the Arctic. Figure 9b is an alternative way of viewing the slice of ARF data from the specific azimuthal range, shown as if the data overlaid by the viewing geometry in Figure 7 had been extracted. This slice represents only the MODIS Aqua viewing geometry (single azimuthal plane) observed at the Green Edge study area, with length of the arrows as the magnitude of the ARF values.

As can be seen in Figure 9a, the distribution of ARF values relevant to the measurements made with the MODIS sensors (i.e., specific azimuthal range) is narrower than when all azimuthal angles are included (i.e., restricted zenithal range). This difference is especially important in consideration of the long tail at the right of the ARF value of unity represented by the vertical dashed line in Figure 9a. Also, although this may not always be the case, the range of ARF values displayed in Figure 9b has little viewing zenithal dependence. For the snow-covered ice, the viewing geometries from MODIS are showing ARF values below unity for most of the viewing zenith angles, implying that if no correction were to be applied to the satellite surface reflectances (i.e., under the assumption of isotropy), the resulting albedo would be biased low. To generalize this observation, we extracted MODIS viewing geometries on all polar plots. In Table 1, we report the corresponding standard deviation of the ARF values as a measure of the spread of the distribution. We have included an additional category named all hemispheric values, which represents all ARF values within the 0° to 90° zenithal range and 0 to 360° azimuthal range.

To interpret the numbers shown in Table 1, consider the first block showing in situ dry snow. The middle column, band 04, reports on the example shown in Figure 9. The standard deviation for the distribution of ARF of the restricted zenithal range (white histogram in Figure 9) is 0.08 and that of the specific azimuthal range (black histogram in Figure 9) is 0.04. When showing all hemispheric values, the standard deviation from the in situ dataset was removed because of the artifacts in the large angles, but computing it for the simulations was possible. Table 1 illustrates how restricting the viewing zenith angle to 60° in conjunction with MODIS-specific azimuthal range greatly reduces the spread of the ARF values for all sea ice surface types. This reduction allows for the development of a simple and robust predictive model from satellite-derived surface reflectance.

Now that we have presented the spectral and angular considerations on which we base this study, we evaluate our model trained with the in situ and simulated data presented above. First, we evaluate our model directly with the 40 surface reflectances of relevant geometry and corresponding albedo from the field measurements of Perovich (1994). When compared with the corresponding field measurements of albedo, we found an MAE of 0.055, an RMSE of 0.072 and an R² value of 0.92. The predictions are gathered around the one-to-one line, representing a good agreement with the evaluation data (Figure 10a). Predicting the albedo for the glazed snow and blue ice surfaces appears somewhat less accurate than for other surface types, mainly because our model was not trained for those specific surface types. Because the error metrics are hard to interpret on their own, we used a simple model to compare the performance of our model. The averaging model does not use the angular information and computes the mean of the three bands of surface reflectances to estimate the albedo. This shortcut to derive PAR albedo from satellite data is tempting but, in principle, albedo should not be computed this way because the upward irradiances should be integrated first and then divided by the total downward irradiance. When the averaging model is used, we have an MAE of 0.085, an RMSE of 0.1, and an R² value of 0.82. Our model outperforms the averaging model, which yields greater scatter around the desired albedo (Figure 10b).

Second, to evaluate the accuracy of our model when used with remote sensing surface reflectance, we added the uncertainty emerging from the correction for the overlying atmosphere. The mean (standard deviation) MAE, RMSE and R² values were determined to be 0.057 (0.0026), 0.074 (0.004) and 0.91 (0.0071), respectively. These values represents a small decrease in performance compared with the evaluation of the dataset without the atmospheric uncertainty. The computed MAE and RMSE of the averaging model were twice as high as those of our model, and the R² was 0.77.

Next, we use Equations 7 to 9 to derive a time series of PAR albedo from spatially and temporally neighboring surface reflectance estimates at 500-m spatial resolution between cold spring (April) and ice break up (mid-July) conditions at the Green Edge study area.

Albedo time series were derived using the model developed above using MODIS surface reflectances for the Green Edge study area of 6 x 6 pixels (approximately 9 km²). To create the albedo time series, our method was to average the albedo values obtained for the entire area (36 pixels) and every satellite pass in a given day. In general, within a stack as presented in Figure 5, possible under- and over-estimations from our model are likely to counteract one another. This regression towards the mean is an advantage of using remote sensing over high latitudes where polar orbiting satellites often revisit a given location even over a relatively short time (Figures 1 and 5). The mean albedo time series with the coefficient of variation (over the 6 × 6 pixel area) are calculated daily and shown in Figure 11, along with related environmental conditions.

A comparison of satellite-derived albedo (black dots) with that derived from in situ data (open squares) is shown in Figure 11a and e. This comparison is not ideal for several reasons. The field data characterized an area of the order of one square meter, and the satellite data characterized the area-average properties of the non-uniform surface over several square kilometers. Moreover, no information on the time of acquisition for the in situ data was available, so normalizing to the sun position or matching data by a less-than-a-day time window was not possible. Nonetheless, even with the spatial heterogeneity of sea ice being potentially quite large at sub-pixel scale (Katlein et al., 2015) and inter-pixel scale (Figure 11b and f) and the potential influence of the diurnal albedo cycle (Lindsay and Rothrock, 1994), the satellite-derived values we obtained are consistent with the albedo derived from in situ data. The comparison with in situ measurements was made, however, for periods when the spatial heterogeneity was low (i.e., when in situ data represented the areal average albedo). Nine days of measurements occurred with valid satellite retrievals (matchups) and the comparison yielded an MAE of 0.040, an RMSE of 0.043 and an R² value of 0.96. This level of accuracy is better than what was obtained in the evaluation of the model (Section 3.2), but is only established from albedo values > 0.6, as shown in Figure 11a and e.

At the Green Edge study area over the course of the melt season, the albedo evolved similarly for years 2015 and 2016 overall (Figure 11a and e). The months of April and May exhibited relatively stable and high albedos with weak spatial heterogeneity, as expressed by the coefficient of variation calculated over the study area. In the first half of June 2015 (2016), when the surface was still very reflective to visible radiation, the PAR albedo dropped from > 0.95 to 0.90 (> 0.95 to 0.82) and the spatial coefficient of variation increased to 2% (3%). Then, during the second half of June, the areally averaged PAR albedo decreased to 0.67 (0.49). The spatial heterogeneity concomitantly increased reaching a coefficient of variation of 6% (8%) (Figure 11b and f). Finally, during the first half of July, the PAR albedo decreased to 0.10 for both years. The spatial heterogeneity was rather stable in July 2015 in comparison to the considerably fluctuating albedo in July 2016, with a coefficient of variation reaching 8% (13%) over the course of the month. This difference suggests that melt ponds were distributed more equally over the site in 2015. However, this suggestion is hard to verify because of the remarkable absence of valid satellite observations soon after the start of melt ponds due to the persistent cloud cover. For both years, the ice export due to landfast ice detachment (on July 18 in 2015; July 16 in 2016) created a drop in PAR albedo from approximately 0.40 to 0.08, producing the most drastic change (in absolute values) in underwater light availability for phytoplankton (Lewis et al., 2019). This transition to an open water albedo defines the limit of applicability for our model.

Context variables were chosen for their potential impact on the PAR albedo and because they were collected independently from satellite observations. There were no melt ponds prior to June 19, 2015, and June 14, 2016. The surface air temperature, measured on site by a Campbell Scientific HC2S3 temperature probe (not shown), was always below 0°C prior to May 24, 2015, and April 20, 2016, and the first day predominantly above 0°C was June 14, 2015, and May 15, 2016. Before the melt pond period, the above surface air temperature was impacting sea ice albedo, but whether variations of PAR albedo of that magnitude are ascribable to changes in temperature or simply to our inability to perfectly convert the remote sensing surface reflectances to the corresponding PAR albedo around the onset of surface melt is difficult to determine. In spite of this ambiguity, we clearly see how large variations in the independently measured melt pond fraction were captured successfully by the satellite PAR albedo time series (Figure 11c and g). The maximum recorded melt pond fractions were on June 25, 2015 (55%), and June 21, 2016 (70%). After drainage (Perovich and Polashenski, 2012), a late minimum in melt pond fraction was recorded on July 7, 2015 (8%) and July 4, 2016 (19%). Melt pond fraction increased again until ice break up (20%) in 2015. In 2016, it increased until July 7 (26%), then slowly decreased to 16% until ice break up. In the second half of June 2016, however, the intense fog prevented the camera from capturing melt pond coverage information on the surface for days. This situation is typical of the late spring and early summer in the Arctic, and clearly reveals the limitation of choosing visible imagery (including satellite remote sensing) for studying sea ice. Finally, precipitation (Figure 11d and h) should impact the PAR albedo time series enough to be captured by our method, at least when enough valid observations are available. Solid precipitation increases the PAR albedo with fresh snow having the highest PAR albedo occurring on Earth, and liquid precipitation decreases the PAR albedo either in fostering melt and promoting metamorphism or by accumulating in melt ponds. Important snowfall events were regularly observed before mid-June of both years, keeping the PAR albedo above 0.80. Snowfalls are also seen in the second half of June 2016, but because snow does not accumulate in ponds the impact was damped. The first rain events were recorded on June 4, 2015, and May 23, 2016. These short-lived and isolated episodes did not seem to affect the PAR albedo time series dramatically. However, it rained for five hours or more during three days in 2015 (June 25, July 1 and 15) and four days in 2016 (June 19, 25, 28 and July 9), and these events correspond with a period when the satellite PAR albedo time series was lower than 0.6.

We have shown here that achieving sea ice PAR albedo estimates over the Arctic is possible from satellite-derived surface reflectance. The main strength of our model is that it uses only the surface reflectance and associated geometry as a direct predictor of surface PAR albedo. This approach was justified by the persistent cloud cover and rapid changes occurring at the Arctic Ocean surface. Where scarcity in data availability would impede the identification of an appropriate sea ice anisotropy model, we established a conservative frame into which a predictive model can process the surface reflectances and provide a PAR albedo retrieval for every clear-sky satellite passage without using any other ancillary data.

Based on in situ measurements and the prescribed atmospheric correction uncertainty, the estimated PAR albedo had a mean absolute error of 0.057, a root mean square error of 0.074 and an R² value of 0.91. We recognize that this model would benefit from a more extensive error analysis in the future. The surface reflectance products are processed using a mature atmospheric correction model, and quality control is indicated by a large array of quality flags of the operational processing. Yet, uncertainty estimates from Vermote et al. (2015) that we used to refine the evaluation of our model were provided as a general measure of quantitative error. An evaluation of the surface reflectance product above large expanses of landfast or mobile sea ice in different clear-sky atmospheric conditions would be appropriate, especially because atmospheric corrections in the Arctic are known to be challenging (IOCCG, 2015). In acknowledgement that absolute values are sensitive to atmospheric corrections, we tried building a model from band ratios but, on the basis of ancillary data, found that it was not evolving dynamically enough over the course of the melt season. The band ratios would also be a strong function of impurities like biogenic material, dissolved organics or sediments. Absolute surface reflectance values, along with their associated geometries, were found to contain sufficient information to investigate the dynamic seasonal evolution of albedo. Our model to derive PAR albedo is also pending proper large-scale validation. Such an exercise would be possible by comparing satellite-retrieved values with airborne radiometric data collected near the surface. This exercise should at least be performed over sea ice at different melting stages to evaluate the full range of sea ice PAR albedo values (0.98 to 0.2).

At the moment, the model only achieves realistic PAR albedo estimates with similar sea ice conditions to those for which it was developed, that is, first-year sea ice surfaces such as dry snow, melting snow, bare ice and melt ponds (in situ and simulated surface reflectance described in Sections 2.3.1 and 2.3.2). Including more ground-based or airborne measurements of surface reflectance and corresponding albedo to train the model, providing a rigorous quality control of the data, would improve the model predictive power and broaden the scope of applicability to a wider range of sea ice conditions. In addition, having more data on sea ice physical properties and/or the inherent optical properties would facilitate the reproduction of cases of sea ice in radiative transfer models. In particular, the fragile upper boundary of sea ice is what matters most for albedo but is difficult to carry to the laboratory for optical measurements (Light et al., 2015). Efforts to study surface-most IOPs directly in the field would be useful. We did not attempt to compare modeled and measured cases because of the lack of information related to the field observations, but a rigorous comparison between AccuRT modeled reflectance and measured reflectance would be a valuable scientific inquiry, given complete sets of inherent and apparent sea ice optical properties are available. More generally, a better understanding of the combined optical properties from different Arctic sea ice types, including multi-year ice would help improve satellite-based sea ice PAR albedo monitoring.

Our model is convenient for anyone facing a tradeoff between a detailed description of the seasonal evolution of PAR sea ice albedo at high spatial resolution, and the need for efficiency when deriving an estimate of this geophysical variable from an operational standpoint over large spatial extents. Using a combination of aerial images and underwater remotely operated vehicle measurements on an ice floe in northern Fram Strait, Katlein et al. (2015) showed that under-ice light variability below the scale of the kilometer was driven more by surface albedo than by sea ice thickness, showing the need for reliable satellite-derived albedo products and time series of sea ice albedo with sub-kilometer resolution. Albedo is variable from year to year (Perovich et al., 2007) and place to place (Grenfell and Perovich, 2004), and given the spatial coefficient of variation depicted in Figure 11 for the Green Edge study area (approximately 9 km²), clearly there is high spatial heterogeneity locally (< 9 km²) and an expectedly high variability regionally (> 9 km²). Our model can be applied to examine the variations of PAR albedo at spatial scales of 500-m or greater during the melt period for different Arctic regions. We believe the differences in PAR albedo are key to a better understanding of in situ observed contrasted regimes of sympagic and pelagic biological production. In fact, deriving PAR albedo time series of melting sea ice to contextualize biological research conducted in ice camps, oceanographic cruises, or polar expeditions is possible. The combination of MODIS-Aqua and Terra surface reflectance product required to retrieve PAR albedo is consistent and continuous since the year 2002, offering an excellent potential to examine the light regime under sea ice (Stroeve et al., 2021). Applying our method over the remote Arctic Ocean will also support year-to-year comparisons of the rapidly evolving sea ice albedo and should thus provide new insights on the impact of the changing climate on Arctic phytoplankton phenology.

Data are accessible at the Green Edge database: http://www.obs-vlfr.fr/proof/php/GREENEDGE/greenedge.php.

The authors are thankful to Eric Brossier and Christian Haas for providing the time-lapse photography and to Environment and Climate Change Canada for the meteorological observations. They are also thankful to Sadashiva Devadiga for helping with the exploitation of the moderate resolution imaging spectroradiometer (MODIS) surface reflectance level 2G-Lite products and the MODIS Teams for the distribution of satellite-derived products and to Simon Bélanger and Robert Frouin for their feedback. They gratefully acknowledge scientific and financial support of Québec-Océan. Finally, they are grateful to the editor and reviewers for their constructive comments.

Julien Laliberté was supported by the FRQNT and the CREATE program (NSERC). Marcel Babin was funded NSERC (05175-2014). This work was funded by the Canadian Excellence Research Chair on Remote sensing of Canada’s new Arctic frontier, Canadian Space Agency Green Edge project, CNES Green Edge project (131425), Network of Centres of Excellence MEOPAR, and Observation core in remote sensing project (1-02-01-064.5). We gratefully acknowledge scientific and financial support of Québec-Océan.

The authors declare that there is no conflict of interest.

Designed and wrote the study: JL.

Helped with the satellite data and drafting of the text: ER.

Helped with the radiative transfer simulations and the writing: BH.

Collected and processed the in situ data and revised the text: GC.

Provided data and revised the text: DP.

Contributed to its conception and revise the text: MB.

Approved the submitted version for publication: All authors.

How to cite this article: Laliberté, J, Rehm, E, Hamre, B, Goyens, C, Perovich, DK, Babin, M. 2022. A method to derive satellite PAR albedo time series over first-year sea ice in the Arctic Ocean. Elementa: Science of the Anthropocene 10(1). DOI: https://doi.org/10.1525/elementa.2020.00080

Domain Editor-in-Chief: Jody W. Deming, University of Washington, Seattle, WA, USA

Associate Editor: Kevin R. Arrigo, Department of Earth System Science, Stanford University, Stanford, CA, USA

Knowledge Domain: Ocean Science

Part of an Elementa Special Feature: Green Edge