Consequences of the 2019 Greenland Ice Sheet Melt Episode on Albedo

In mid-June 2019, the Greenland ice sheet (GrIS) experienced an extreme early-season melt event. This, coupled with an earlier-than-average melt onset and low prior winter snowfall over western Greenland, led to a rapid decrease in surface albedo and greater solar energy absorption over the melt season. The 2019 melt season resulted in significantly more melt than other recent years, even compared to exceptional melt years previously identified in the moderate-resolution imaging spectroradiometer (MODIS) record. The increased solar radiation absorbance in 2019 warmed the surface and increased the rate of meltwater production. We use two decades of satellite-derived albedo from the MODIS MCD43 record to show a significant and extended decrease in albedo in Greenland during 2019. This decrease, early in the melt season and continuing during peak summer insolation, caused increased radiative forcing of the ice sheet of 2.33 Wm−2 for 2019. Radiative forcing is strongly influenced by the dramatic seasonal differences in surface albedo experienced by any location experiencing persistent and seasonal snow-cover. We also illustrate the utility of the newly developed Landsat-8 albedo product for better capturing the detailed spatial heterogeneity of the landscape, leading to a more refined representation of the surface energy budget. While the MCD43 data accurately capture the albedo for a given 500 m pixel, the higher spatial resolution 30 m Landsat-8 albedos more fully represent the detailed landscape variations.


Introduction
The Greenland ice sheet (GrIS) is currently experiencing a decadal-scale trend of intensified snow and ice melt [1][2][3][4][5] corresponding to a positive trend in melt season intensity and duration [3]. Increased temperatures are magnified by the snow-albedo feedback, in which progressively earlier melt onset and increased seasonal melt intensity serve to lower the albedo, thus resulting in reduced reflection of solar irradiance and further warming [1]. Increased warming and subsequent melting then lead to a rapid and sequential decrease in albedo, as snow grain growth darkens the upper snowpack. Incipient melting reduces albedo still further, and ultimately melt, runoff, and evaporation expose bare ice in the ablation zone of the ice sheet [1,6]. While fresh, fine-grained snow typically has an albedo in the range of~0.8-0.9, coarse-grained older snow has an albedo of~0.7-~0.8, and wet snow and bare clean ice of~0.4-~0.7 [7][8][9]. This reduced reflectance of solar radiation alters the surface energy budget, with a positive feedback of increasing solar energy absorbance causing an increase in the rate of meltwater production [1,6].
Previous research has shown that the rate of summer albedo decrease has accelerated during the 2000s [10]. During the 2019 melt season, the GrIS experienced a record early-season melt event, with an earlier-than-average melt onset over much of the coastal ice sheet, and a record melt area (relative to 1978-2018) on June 12, encompassing 45 to 55% of the ice sheet [11,12]. This intense early melt start followed exceptionally low prior winter snowfall totals over western Greenland [13]. Here, we investigate the impacts of this melt episode on albedo using the Terra/Aqua moderate-resolution imaging spectroradiometer (MODIS) bidirectional reflectance distribution function (BRDF) and albedo product (MCD43) and provide a preliminary estimate of the albedo-driven radiative forcing dynamics for 2019 relative to the multi-decade 2000-2020 time period. Additionally, we explore the utility of the newly produced higher spatial resolution (30 m) albedo product derived from the Landsat-8 operational land imager (OLI) instrument [14], to capture the finer-spatial-scale albedo patterns, and thus the full range of albedo variability in the landscape.
Previous research has established that the GrIS is experiencing both long-term mass loss and exceptional seasonal melt episodes [1,3,15]. MODIS-derived surface temperature has been used to examine GrIS melt-season dynamics, finding that over the period from 2000 to 2005, years 2002 and 2005 showed unusually extensive melt, and that surface temperature variability increased over the time period [15]. A positive trend has also been found in ice surface temperature between 2000 and 2012, with the greatest increases in northwestern Greenland during that study period [3]. For the GrIS as a whole, the surface temperatures in both the summer and winter seasons have increased. Major melt events in western Greenland occurred in 2002,2003,2007,2010,2012, and 2019 [2,3,12]. Extensive and persistent melt has been documented in most summer seasons compared to a period extending between 1981-2010, accompanying the highest ice-sheet-wide summer average temperatures in the MODIS data record.
Several studies have employed satellite and in situ sensors to evaluate the albedo data displaying these trends [2,16,17]. It has been shown that the Terra-only MOD43 product, the antecedent to the modern MCD43 product (Terra and Aqua), agreed closely with an in situ sensor network of 16 Greenland Climate Network (GC-Net) automatic weather stations distributed over homogenous and semi homogenous ice-covered surface throughout the GrIS over the years 2000 to 2003 (root mean square error (RMSE) = 0.04 for high quality retrievals) [16]. Stroeve et al. (2013) extended this work by examining MCD43 data from 2000 to 2012, confirming the agreement between field and satellite sources (RMSE = 0.067) [2]. These authors found a negative trend of summer (June, July, and August, hereafter: JJA) albedo across the 2000 to 2012 period, and highlighted 2012 as an extreme melt year. This negative albedo trend was found to be largest during July along western Greenland. Moustafa et al. (2017) also found good agreement between MCD43 and in situ data during a campaign in the GrIS ablation zone, with error varying from −4% to +7%, within the range of the standard error of the field spectrometers used [17].
In this investigation, our research questions are as follows: what were the albedo consequences for 2019 in relationship to the entire multi-decadal MCD43 record? Do fine-scale albedos provided by the Landsat-8 provide additional information? What are the radiative forcing (RF) consequences of the 2019 melt? Figure 1 shows the study area, which comprises the entirety of Greenland, as defined by the ice portion of the ocean/land/ice mask [18]; the ice mask is indicated by the blue line just inland of the coast and comprises 1,754,779 km 2 . While the 2019 melt episode affected a large portion of the entire island, the impacts were particularly pronounced along the western coast [2,3], and therefore a subset area comprising all catchments draining to the west coast is shown as an orange outline (756,773 km 2 ). Additionally, the catchments of the Helheim and Russell Glaciers were examined in particular, based on their importance in Greenland's seasonal ice melt dynamics [19][20][21]. Therefore, we show several areas of interest (AOIs) on the western and southeastern portions of Greenland ( Figure 1). The labeled light blue areas are the drainage basins for the Russell (134,974 km 2 ) and Helheim (64,223 km 2 ) glaciers, as defined previously from Ice, Cloud, and land Elevation Satellite (ICESat) data [22]. The Helheim Glacier is a large marine-terminating outlet glacier on the southeastern coast, which has been characterized by fluctuating speeds and advancement/retreat associated with dynamic discharge behavior [23,24]. The Russell Glacier is a large land-terminating outlet glacier on the western margin of the GrIS, where extreme inter-annual variations in meltwater availability have been shown to have complex implications for glacier velocity [19]. The orange line indicates the union of all drainage basins along the west coast of Greenland, constituting another AOI that we refer to as the "west coast catchments AOI" (731,410 km 2 ). Within each AOI, we explore albedo trajectories of each year from 2000-2020. The green and blue boxes and red transect represent subset areas used to further explore the spatiotemporal pattern of the albedo and are explored in greater detail in Figures 2-8.

Study Area Description
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 17 importance in Greenland's seasonal ice melt dynamics [19][20][21]. Therefore, we show several areas of interest (AOIs) on the western and southeastern portions of Greenland (Figure 1). The labeled light blue areas are the drainage basins for the Russell (134,974 km 2 ) and Helheim (64,223 km 2 ) glaciers, as defined previously from Ice, Cloud, and land Elevation Satellite (ICESat) data [22]. The Helheim Glacier is a large marine-terminating outlet glacier on the southeastern coast, which has been characterized by fluctuating speeds and advancement/retreat associated with dynamic discharge behavior [23,24]. The Russell Glacier is a large land-terminating outlet glacier on the western margin of the GrIS, where extreme inter-annual variations in meltwater availability have been shown to have complex implications for glacier velocity [19]. The orange line indicates the union of all drainage basins along the west coast of Greenland, constituting another AOI that we refer to as the "west coast catchments AOI" (731,410 km 2 ). Within each AOI, we explore albedo trajectories of each year from 2000-2020. The green and blue boxes and red transect represent subset areas used to further explore the spatiotemporal pattern of the albedo and are explored in greater detail in Figures 2-8. This research relies heavily on the satellite-based albedo product MCD43 V006, which is generated using data from both the Aqua and Terra MODIS sensors. This daily product uses a kernel-based, semi-empirical approach, with a combination of RossThick and LiSparse reciprocal kernels to retrieve a BRDF model for a given gridded 500 m 2 pixel, and subsequently to determine diffuse white-sky bi-hemispherical and direct black-sky directional-hemispherical components of albedo [25][26][27][28][29]. The MODIS product makes use of all high quality, cloud-free observations acquired over a 16-day observation window, with emphasis on the day of interest in the center of the window, to accumulate a sufficiently diverse angular sampling of the target pixel, such that the surface anisotropy can be accurately modeled [26]. The actual (blue-sky) surface albedo for a particular time and date can be calculated using the aerosol optical depth and solar zenith angle for the time of interest to infer the observed proportions of diffuse and direct sunlight. Here we use the method described by Román et al. to calculate the blue-sky albedo, which considers multiple scattering between the atmosphere and land surface [30]. We use the solar zenith angle associated with local solar since this is the maximum solar illumination, and we This research relies heavily on the satellite-based albedo product MCD43 V006, which is generated using data from both the Aqua and Terra MODIS sensors. This daily product uses a kernel-based, semi-empirical approach, with a combination of RossThick and LiSparse reciprocal kernels to retrieve a BRDF model for a given gridded 500 m 2 pixel, and subsequently to determine diffuse white-sky bi-hemispherical and direct black-sky directional-hemispherical components of albedo [25][26][27][28][29]. The MODIS product makes use of all high quality, cloud-free observations acquired over a 16-day observation window, with emphasis on the day of interest in the center of the window, to accumulate a sufficiently diverse angular sampling of the target pixel, such that the surface anisotropy can be accurately modeled [26]. The actual (blue-sky) surface albedo for a particular time and date can be calculated using the aerosol optical depth and solar zenith angle for the time of interest to infer the observed proportions of diffuse and direct sunlight. Here we use the method described by Román et al. to calculate the blue-sky albedo, which considers multiple scattering between the atmosphere and land surface [30]. We use the solar zenith angle associated with local solar since this is the maximum solar illumination, and we obtained aerosol optical depth from the MOD08 product. Note that albedos for solar zenith angles less than 70 degrees are considered lower quality and thus suspect, and as such have been omitted from this analysis. Validation at stage 3 has been achieved for the MCD43 albedo products, meaning that the product's uncertainties are well understood and quantified relative to numerous in situ reference data over a broad range of environments, including snow-covered pixels and Greenland's ablation region [2,16,17,26,31,32].
Changes in surface albedo alter the surface energy balance, driving radiative forcing [33]. Radiative forcing measures perturbations to the solar energy balance at the planetary surface or at the top of the atmosphere. Changes to radiative forcing are particularly relevant in high-latitude regions dominated by snow/ice melt dynamics, which are subject to melt-albedo radiative forcing feedbacks [1,34,35]. Albedo changes Earth's energy budget, thus its effect on global climate is evaluated in terms of the changes to this energy balance. This methodology allows the impact of changing albedo to be quantified and compared to the climate impacts of greenhouse gases. Using all available MCD43 blue-sky albedo observations (i.e., between 2000-02-24 and 2020-08-07) at local solar noon, we describe the albedo-driven radiative forcing implications of the 2019 melt episode, relative to the 2000-2020 baseline.

Satellite Albedo Retrievals
The MCD43 product performs a high-quality full model inversion, as described above, when sufficient angularly distributed cloud-free surface reflectance observations are available within the 16-day window. Otherwise, a backup algorithm is used for each pixel, performing a lower quality magnitude inversion using available observations and archetypal BRDF parameters from the most recent high-quality retrieval from the backup database. The MCD43A2 product provides extensive quality assurance (QA) flags. Here, we used only the highest quality (full inversion), clear sky albedo retrievals (no cloudcontaminated pixels are used), and only included observations with solar zenith angle (SZA) ≤ 70 • . This constraint, in addition to the lack of sunlight above the Arctic Circle during winter, limited our observation window to late February to late October. The Landsat-8 albedo is based on the BRDF parameters from the spatiotemporally coincident MCD43 data, and is produced with QA flags. Therefore, the same QA standards were applied to Landsat-8 data. Additionally, it should be noted that snow and snow-free BRDFs are processed separately by the algorithm, with snow status determined as that of the day of interest [26].
MCD43 provides daily retrievals of the directional and diffuse components of albedo, referred to as black-sky and white-sky, respectively. The black-sky albedo was retrieved at the SZA for local solar noon (the time of greatest solar illumination). To calculate the actual blue-sky albedo observed on a particular pixel and date/time, we used the method presented by Róman et al., 2010, which calculates blue-sky albedo as a function of the diffuse and direct albedo components, aerosol optical depth, and solar zenith angle at the time of interest. Aerosol optical depth was retrieved from the coarse resolution MOD08 daily aerosol product.
The MCD43 albedo provided the basis for the majority of the analysis presented in this work, since it is available at a daily frequency from 2000 to present. Figures presenting average values were aggregated spatially using the AOIs shown in the manuscript to calculate mean and standard deviation. We removed any dates when the number of valid observations was lower than the 25th percentile of valid pixels calculated across the entire data range; this served to screen out dates with excessive cloud cover or otherwise low QA values, which would otherwise potentially bias the results. Seasonal averages were calculated for the December, January, February (DJF), March, April, May (MAM), JJA, and September, October, November (SON) periods.
Landsat-8 has a revisit time of 16 days, which when combined with interference from cloud cover, limits the number of observations available for analysis. For this work, we Remote Sens. 2021, 13, 227 5 of 16 calculated Landsat-8 albedo for the spatial subsets detailed in Figure 1 and utilized images with low cloud coverage.
To determine the statistical significance of each year's aggregated albedo value from the 2000-2020 mean and to determine statistical separability of seasonal mean albedos, a series of t-tests were conducted. For annual tests, as shown in Figure 4, each year in the time series was compared to the 2000-2020 mean, using a 99% confidence level. Similarly, for seasonal analysis, a given season within a year was compared to the same season in each other year, using a 99% confidence level.

Spatial Autocorrelation Analysis
Spatial autocorrelation, measured with Moran's I, was measured for spatiotemporally coincident Landsat-8 and MCD43 shortwave blue-sky albedo measurements. Moran's I is a global metric used to measure spatial autocorrelation, with −1 representing perfectly negatively autocorrelated data with dissimilar values close together (e.g., a checkerboard) and +1 representing perfectly spatially clustered data, in which similar values are close together [36]. This metric of spatial autocorrelation was calculated for each available pair of images from the two datasets for 2019. All available images were used, yielding 18 image pairs after cloud-covered images were removed. The area used for this analysis is shown in Figure 6.

Radiative Forcing
The climate impact of changing albedo was expressed in terms of radiative forcing (Wm −2 ). Calculations were made on a monthly basis and normalized per square kilometer of land conversion. The radiative forcing of surface albedo at the top of the atmosphere (RF TOA ∆a SFC ) was calculated from the average monthly difference between the annual albedo and mean albedo from years 2000-2020 (∆a SFC ), from the monthly incoming shortwave radiation at the surface (SW SFC↓ ), and from the local clearness factor (T) (Equation (1)) [33,37]: where A earth represents the total surface area of the Earth and A local represents the ice-only portion of the compared site n, y represents the comparison year, n represents the compared site, and m represents month. T represents the local clearness index, which is the fraction of shortwave radiation received at the ground over the shortwave radiation received at the top of the atmosphere [33,[37][38][39][40]. Radiative forcing at the surface (RF SFC ∆a SFC ) was calculated similarly (Equation (2)) [33,37]: ∆a SFC describes the average monthly difference between the annual albedo and mean 2000-2020 albedo; and SW SFC↓ represents the monthly incoming shortwave radiation at the surface [33,[37][38][39][40][41].
Monthly SW TOA↓ for the study focal location was obtained from the Clouds and the Earth's Radiant Energy System (CERES) Science Team: CERES EBAF-TOA Edition 4.0. NASA Atmospheric Science and Data Center (ASDC) [42]. Monthly SW SFC↓ for the study focal location was obtained from the ICECAPS Summit Station, Greenland (72 • 36 N, 38 • 25 W, 3250 m) [43]. Only cloud-free data were included in the analysis. Due to the low number of albedo retrievals during the November to February period of each year, we omitted these months from the radiative forcing analysis. anomalies during the melt season dates for 2019, with corresponding average albedo for those dates between 2000 and 2020. Here, the extensive early melt in 2019 is most clearly visible along the western coast. Figure 3 shows the transect AOI outlined in red in Figure 1 at ten-day time steps between April 30 and July 29 of 2019, illustrating the extensive melt, particularly on the western and southeastern coasts. The bottom panel of the figure shows MCD43 albedo and MODGRNLD ice surface temperature (IST) [44,45] as a function of longitude along a transect spanning from 65 • N, 54 • W to 65 • N, 38 • W. The IST tends to decrease as the surface albedo increases, but note that the IST does not exceed 0 • C for the snow/ice portion of the transect, which helps to explain the relatively low R 2 value of 0.26. Note that the transitions between land and ice occur at approximately 49.7 • W and 40.5 • W, between which albedo can be seen to increase rapidly along the transect.

MODIS Albedo
Because the most intense albedo decreases occurred in the western coast region, we used the west coast catchments AOI (see Figure 1) to aggregate the albedo trajectories for the 2000 to 2020 period, as shown in Figure 4a Figure 3 shows the transect AOI outlined in red in Figure 1 at ten-day time steps between April 30 and July 29 of 2019, illustrating the extensive melt, particularly on the western and southeastern coasts. The bottom panel of the figure shows MCD43 albedo and MODGRNLD ice surface temperature (IST) [44,45] as a function of longitude along a transect spanning from 65° N, 54° W to 65° N, 38° W. The IST tends to decrease as the surface albedo increases, but note that the IST does not exceed 0 °C for the snow/ice portion of the transect, which helps to explain the relatively low R 2 value of 0.26. Note that the transitions between land and ice occur at approximately 49.7° W and 40.5° W, between which albedo can be seen to increase rapidly along the transect.   Because the most intense albedo decreases occurred in the western coast region, we

Landsat-8 Albedo
As Landsat-8 is a near-nadir instrument, unlike MODIS, its observations alone cannot provide the multi-angular observations required to retrieve an accurate surface BRDF and subsequent measures of albedo. Rather, retrieval of albedo from Landsat-8 relies on the use of spatiotemporally coincident MCD43 BRDF parameters to characterize the BRDF for the representative near-nadir spectral reflectances associated with each Landsat acquisition [14,46]. Using within the Landsat-8 image, the MCD43 BRDF parameters are used to calculate the albedo-to-nadir ratio for the viewing geometry of the given Landsat-8 acquisition, which in turn is applied across the Landsat-8 surface reflectance scene to yield 30 m gridded measures of surface albedo. Figure 5 illustrates the additional spatial detail that can be captured by the 30 m Landsat-8 albedo product, again showing the 2019 melt season with all available clear retrievals. The location of this subset is shown by the green box in Figure 1 and is situated in the ablation zone. In this figure, the surface topography is much more evident, and melt ponds and the interfaces between snow, bare ice, and substrate are clearly delineated.
We conducted a spatial autocorrelation analysis for the 30 m resolution Landsat-8 and 500 m MCD43 albedos in the area shown in Figure 6 (the blue inset box in Figure 1). Moran's I of +1 represents perfectly spatially clustered data, in which similar values are close together [36]. Using global Moran's I to indicate spatially clustered patterns, results indicate that Landsat-8 albedo images produce scores very close to 1 (0.988 < i < 0.997), implying high spatial autocorrelation, whereas the same analysis over the same area using the MODIS albedo data showed markedly lower values (0.578 < i <0.831). This suggests that while MCD43 data accurately represent albedo for a given gridded 500 m pixel, accurate mapping of the intrinsic landscape pattern requires the higher spatial resolution afforded by Landsat-8. Detailed spatial pattern information is particularly important in areas with rough surface topography such as western Greenland, where albedo values may quite widely. that can be captured by the 30 m Landsat-8 albedo product, again showing the 20 season with all available clear retrievals. The location of this subset is shown by th box in Figure 1 and is situated in the ablation zone. In this figure, the surface topo is much more evident, and melt ponds and the interfaces between snow, bare substrate are clearly delineated. We conducted a spatial autocorrelation analysis for the 30 m resolution La and 500 m MCD43 albedos in the area shown in Figure 6 (the blue inset box in Fi Moran's I of +1 represents perfectly spatially clustered data, in which similar va close together [36]. Using global Moran's I to indicate spatially clustered patterns indicate that Landsat-8 albedo images produce scores very close to 1 (0.988 < i < implying high spatial autocorrelation, whereas the same analysis over the same ar the MODIS albedo data showed markedly lower values (0.578 < i <0.831). This s that while MCD43 data accurately represent albedo for a given gridded 500 m p curate mapping of the intrinsic landscape pattern requires the higher spatial re afforded by Landsat-8. Detailed spatial pattern information is particularly impo areas with rough surface topography such as western Greenland, where albedo may quite widely.

Evaluation Against in Situ Measurements
To ensure that these satellite products accurately capture the surface albed the greater study area, an evaluation of the MCD43 blue-sky albedo retrievals w pared to in situ albedo data farther inland from the National Oceanic and Atm Administration (NOAA) Greenland Environmental Observatory (GEOSummit, 38.48° W). This station continuously measures both upwelling and down shortwave radiation using Kipp & Zonen CM22 pyranometers, providing app evaluation data. Using all available satellite retrievals and in situ daily data for 2018, we found a high correspondence of albedo (RMSE = 0.08, Figure 7), indica applicability of the satellite record for monitoring the broader region. Additiona screened Landsat-8-based blue-sky albedo was also compared to the sensor m ments for all available pixels from May to July of 2015-2017, showing an RMSE To ensure that the tower-based in situ data are appropriate to evaluate the 500 m Although the overall pattern is consistent between the two sensors, the increased spatial variability afforded by Landsat-8 provides considerably more detail.

Evaluation against In Situ Measurements
To ensure that these satellite products accurately capture the surface albedo within the greater study area, an evaluation of the MCD43 blue-sky albedo retrievals was compared to in situ albedo data farther inland from the National Oceanic and Atmospheric Administration (NOAA) Greenland Environmental Observatory (GEOSummit, 72.58 • N 38.48 • W). This station continuously measures both upwelling and downwelling shortwave radiation using Kipp & Zonen CM22 pyranometers, providing appropriate evaluation data. Using all available satellite retrievals and in situ daily data for 2013 to 2018, we found a high correspondence of albedo (RMSE = 0.08, Figure 7), indicating the applicability of the satellite record for monitoring the broader region. Additionally, QA-screened Landsat-8-based blue-sky albedo was also compared to the sensor measurements for all available pixels from May to July of 2015-2017, showing an RMSE of 0.053. To ensure that the tower-based in situ data are appropriate to evaluate the 500 m gridded MCD43 data, we performed a spatial representativeness analysis, also shown in Figure 7. This analysis, which uses variography to assess spatial dependence within 1.0, 1.5, and 2.0 km boxes around the tower location, indicated adequate spatial homogeneity around this tower location, within a MODIS 500 m 2 gridded pixel [30,31,47,48].
2018, we found a high correspondence of albedo (RMSE = 0.08, Figure 7), indicating the applicability of the satellite record for monitoring the broader region. Additionally, QAscreened Landsat-8-based blue-sky albedo was also compared to the sensor measurements for all available pixels from May to July of 2015-2017, showing an RMSE of 0.053. To ensure that the tower-based in situ data are appropriate to evaluate the 500 m gridded MCD43 data, we performed a spatial representativeness analysis, also shown in Figure 7. This analysis, which uses variography to assess spatial dependence within 1.0, 1.5, and 2.0 km boxes around the tower location, indicated adequate spatial homogeneity around this tower location, within a MODIS 500 m 2 gridded pixel [30,31,47,48]. Figure 7. Evaluation results at the Greenland Environmental Observatory (GEOSummit). The tower data were obtained from the National Oceanic and Atmospheric Administration's broadband radiation pyranometers. The lower two panels show the spatial representativeness assessment, based on corresponding 30 m Landsat-8 OLI data.

Radiative Forcing
Using the MODIS MCD43 BRDF/albedo product, we investigated the climate impact of changing albedo during the 2019 melt season by establishing estimates of radiative forcing [30,37]. As indicated by Figure 8, our results show that in 2019 the GrIS experienced a positive surface radiative forcing of producing of 2.33 Wm −2 , relative to the 2000-2020 mean (Equation (2)). Additionally, we found that radiative forcing at the surface specifically for the west coast, Helheim, and Russell Glacier AOIs was 2.69 Wm −2 y −1 , 7.99 Wm −2 y −1 , and 5.24 Wm −2 y −1 , respectively. When considering these impacts at a global scale, the albedo-driven radiative forcing contribution of the GrIS for 2019 was 6.22 × 10 −3 Wm −2 (Equation (1)).

Discussion
Analysis of the spatial patterns of cloud-free, blue-sky albedo during the major melt episode of 2019 on the Greenland ice sheet shows that the areas of lower albedo cause positive radiative forcing feedback due to greater absorption of solar insolation in the western part of the ice sheet. Compared to the entire multi-decadal MODIS record, 2019 shows significantly reduced albedo for the entire GrIS, and particularly in the west coast catchments AOI. These results are shown both at a daily time step (Figures 3a,c) and when averaged annually (Figures 3b,d). As previously noted in the literature, the year 2012 also experienced extreme melt [1,3], evidenced here by markedly low albedos shown in Figure  4. However, 2019's mean albedo was even lower than during the early melt season in 2012, when solar insolation is near maximum, increasing its impact on the surface energy balance. The severity of 2019 is highlighted by the two previous years, which showed above average albedo for the entire study area and west coast catchments AOI (Figure 4), corresponding to cooler and more snow-rich conditions [13].
The severity of the 2019 melt season is even clearer when observing individual catchments, such as the Russell

Discussion
Analysis of the spatial patterns of cloud-free, blue-sky albedo during the major melt episode of 2019 on the Greenland ice sheet shows that the areas of lower albedo cause positive radiative forcing feedback due to greater absorption of solar insolation in the western part of the ice sheet. Compared to the entire multi-decadal MODIS record, 2019 shows significantly reduced albedo for the entire GrIS, and particularly in the west coast catchments AOI. These results are shown both at a daily time step (Figure 3a,c) and when averaged annually (Figure 3b,d). As previously noted in the literature, the year 2012 also experienced extreme melt [1,3], evidenced here by markedly low albedos shown in Figure 4. However, 2019's mean albedo was even lower than during the early melt season in 2012, when solar insolation is near maximum, increasing its impact on the surface energy balance. The severity of 2019 is highlighted by the two previous years, which showed above average albedo for the entire study area and west coast catchments AOI (Figure 4), corresponding to cooler and more snow-rich conditions [13].
The severity of the 2019 melt season is even clearer when observing individual catchments, such as the Russell (Figure 4c). The Helheim catchment also showed a strong early season darkening for MAM of 2020 (0.773), although the albedo tended towards the long-term mean after approximately DOY 150.
Radiative forcing results show that the 2019 melt season had a substantial impact on the surface energy balance for the GrIS, with an additional 2.33 Wm −2 for 2019 relative to the 2000-2020 baseline period (Figure 7). When considering the catchment-based AOIs, the radiative forcing impacts show increases of 2.69 Wm −2 , 7.99 Wm −2 , and 5.24 Wm −2 , for the west coast catchments AOI, Helheim Glacier, and Russell Glacier, compared to the 2000-2020 mean for each region, respectively. In the global context, 2019 contributed 6.22 × 10 −3 Wm −2 to Earth's radiative forcing. While this forcing constituted a rather small proportion of global radiative forcing, the regional ramifications of the melt include intensified sea level rise and North Atlantic freshening. These findings emphasize both the global and regional importance of the GrIS albedo as a part of the global climate system [1,49]. Albedo decreases during peak insolation (near the solstice), as occurred in 2019, are particularly impactful on radiative forcing in this ice-and snow-covered region.
The MODIS sensors aboard the Aqua and Terra satellites provide daily coverage of the study area, providing an invaluable data record spanning from 2000 to present. However, at a gridded 500 m spatial resolution, the MCD43 BRDF/albedo products cannot capture the full landscape pattern, particularly along the GrIS margins where melt dynamics create a complex mosaic of snow and ice of different grain size, melt pools, and exposed land surface (Figures 4 and 5). The OLI instrument aboard Landsat-8 provides 30 m spatial resolution reflectances and can therefore resolve much more of this spatial detail. The sequence of Landsat-8 images shown in Figure 4 demonstrates the spatial detail available in this product. Here, the land/ice boundary is quite clear, as are melt ponds and snow/ice with locally variable albedo. We show a direct side-by-side comparison of the observed spatial pattern of MCD43 and the Landsat-8 albedo product in Figure 5. What is evident to the human interpreter is confirmed by a measurement of global spatial autocorrelation, Moran's I: the intricate patterning of melt ponds, striated snow, and ice darkness are captured much more clearly by Landsat-8. Moran's I results showed that all Landsat-8 albedo scenes had values very close to 1 (0.988 < I < 0.997), which is the maximum possible value and indicates strong positive spatial autocorrelation. Conversely, the Moran's I values for the same date and area for MCD43 were much more moderate (0.578 < I < 0.831), indicating a pattern closer to spatial randomness (which is greatest at 0). The Landsat-8 albedo product, which incorporates BRDF information from the commensurate MCD43 retrievals and adds value to the Landsat-8 surface reflectance product by providing whitesky and black-sky albedo, thus allowing for a complete surface energy balance accounting at 30 m, which will be further explored in future research. Additionally, information pertaining to the location and extents of ice-dammed lakes visible at the Landsat-8 spatial scale may prove important for future research on rapid drainage events [20].
The inverse relationship between IST and albedo, measured by MODGRNLD [45] and MCD43, respectively, is quite clear in Figure 3b. The transect spans the width of Greenland at 65 • N, and therefore shows land area at the extremes, particularly on the west end, and ice in the middle. The agreement between these independently created MODIS data products helps verify the findings by corroborating the physical conditions under observation.
Note that due to the solar zenith angle-based constraints of the albedo product, the vast majority of valid 2020 data has already been incorporated, and we do not expect the results to change dramatically by including the remaining retrievals from 2020, since very few additional valid observations will become available (typically none beyond early September). Future work will focus on a production of the entire Landsat-8 albedo archive, in order to fully capture the fine spatial resolution albedo patterns within the study region. An effort is also underway to retrieve albedo from Sentinel-2A/B imagery [50], launched in 2015 and 2017, which, when combined with Landsat-8 and the future Landsat-9, will dramatically increase the frequency of observations, potentially allowing for a much finer spatial resolution analysis of recent GrIS melt trajectories.

Conclusions
This research extends earlier work on Greenland's albedo characteristics and dynamics, and in particular helps better characterize the albedo trajectories of the extreme 2019 melt season. We demonstrate that 2019 exhibited particularly low albedo values, with only 2012 showing a lower JJA mean albedo. Mean albedo during the early season (MAM) 2019 was notably low, and statistically similar to the extreme melt years of 2003 and 2010. When considering the west coast catchments in particular, 2019 showed the lowest MAM albedo mean, even considering other extreme melt years of 2003, 2010, and 2012. Additionally, in this region, the JJA mean was also notably low, and statistically similar to 2010 and 2012. Within the Russell and Helheim glacier catchments, the albedo effects of the 2019 melt episode were even more clearly apparent, with Russell showing lower MAM albedo means than any other year, and Helheim showing lower JJA albedo means than any other year. The 2019 melt season saw particularly extensive melt, even compared to other exceptional years in the recent MODIS record. We show that the albedo decreases associated with this melt season led to increased local radiative forcing, which has implications for regional and global climate change. The GrIS represents a fragile and critically important element of the global climate system, with rapidly decreasing mean albedos hastening ice melt and thus increasing absorbed radiation. It is critical to continue to monitor albedo of the GrIS and other high latitude areas, in order to more fully understand climate impacts on and from these environments. Funding: This research was funded by the following grants: USGS grant number 140G0118C0010 and 140G0118C005 for Landsat Science Team Support; NASA grant 80NSSC18K0642 for MODIS Science Team support; and NASA grant 80NCSSC18K0479 for Land Cover and Land Use Change Circumpolar Albedo support.

Conflicts of Interest:
The authors declare no conflict of interest.