Introduction

The Greenland Ice Sheet has been losing mass since the 1990s and is currently the single largest cryospheric contributor to observed global sea-level rise1. This recent mass loss is primarily due to enhanced surface meltwater runoff, most of which derives from the narrow ablation zone where bare, dark glacial ice is exposed each summer2,3. There is a growing realization that clouds play an important role in determining the amount of meltwater runoff generated in the ablation zone4,5,6,7. Clouds are known to regulate surface energy balance through the competing effects of shortwave radiation shading (i.e. cooling) and longwave radiation trapping (i.e. warming)4,7,8,9. However, the precise radiative effects of clouds in the ablation zone remains uncertain due to a reliance on modeled radiative fluxes which are known to contain biases8,10,11. Therefore, it remains largely unknown how changes in cloudiness, such as those which occur due to atmospheric blocking events12, alter the rate of meltwater production from the ablation zone both currently and in the future.

While in situ observations of cloud radiative effects (CRE) have been acquired at research stations, they are sparsely distributed and mainly constrained to the accumulation zone (e.g.13,14,15,16. Satellite remote sensing therefore provides the only practicable technique for observing CRE across the whole ice sheet. A few studies have used CloudSat and Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO) to characterize cloud properties and their radiative effects over Greenland and the Arctic7,8,17,18. However, the active remote sensing instruments employed on these satellites have sparse spatiotemporal sampling, requiring aggregation into very coarse grid cells (e.g. 2 × 2° or 2 × 1°) to ensure a sufficient number of samples within each grid cell. Such products are therefore too coarse to resolve Greenland’s narrow ablation zone7,8,11, preventing characterization of CRE for the very part of the ice sheet that generates runoff. Optical remote sensing instruments such as the MODerate Imaging Spectroradiometer (MODIS) onboard Terra and Aqua offer much finer spatial and temporal resolution19. Yet visible-infrared wavelengths cannot provide a complete vertical profile of cloud structure, meaning that MODIS-derived cloud properties and associated radiative effects are inherently uncertain. The production of high-quality CRE observations at high spatial resolution therefore requires intelligent integration of cloud profiling with optical imagery, to combine the strengths of these different sensor technologies while overcoming their individual limitations.

Such a task is well-suited to supervised machine learning techniques that can identify complex, non-linear relationships in data to make predictions without explicit consideration of the intermediate physical processes. Machine learning has been particularly successful for integrating sparse, high-quality field measurements with more extensive remotely sensed observations for the prediction of phenomena across large spatial scales20. This approach has been widely used in the geosciences for predicting land cover21, evapotranspiration22, and CO2 emissions23. However, machine learning is relatively undeveloped in glaciology and, to our knowledge, has not been used to combine remote sensing data for investigations of Greenland Ice Sheet surface energy balance.

Here we use Random Forest Regression to fuse the respective strengths of CloudSat/CALIPSO (high accuracy vertical profiles of cloud properties) and MODIS (high spatial and temporal resolution visible-infrared imagery), enabling accurate, daily quantification of CRE over the Greenland Ice Sheet during summer for the 2003–2020 period (see Methods). We then exploit the unprecedented spatial resolution of our product (1 × 1 km grid cell size) to resolve fine spatial gradients and interannual variations of CRE within the narrow ice sheet ablation zone, where nearly all meltwater runoff is produced. We also use our product to investigate the impacts of anticyclonic circulation regimes (i.e. blocking events) on clouds and radiative fluxes in the ablation zone. Finally, we combine our product with projections from the CMIP6 experiment to assess the radiative effects of clouds and atmospheric blocking events on Greenland Ice Sheet surface energy balance in a future, warmer climate.

Results

Cloud radiative effects observed at high-resolution

We find that, averaged across the ice sheet, mean summer cloud cover is 54% for the 2003–2020 period, varying from 32% (5th percentile) in the Northeast Greenland to 71% (95th percentile) in South Greenland (Fig. 1a). Relative to clear skies, we find that clouds increase summer downward longwave radiation received at the surface by 16.4 ± 2.0% and decrease summer downward shortwave radiation received at the surface by 13.6 ± 2.6% (Fig. 1d, e). When these two competing effects are combined, we find that clouds exert a net warming effect on the Greenland Ice Sheet of +21.8 ± 5.9 W m−2 during the summer. These averaged values, however, mask strong spatial variations in CRE, which have a much larger warming effect in the accumulation zone (+23.5 ± 6.3 W m−2) than the ablation zone (+11.4 ± 3.1 W m−2), a finding that corroborates previous studies6,16 (Fig. 1f, S1). This can be attributed to lower albedo in the ablation zone (Fig. 1c), which amplifies the shortwave shading effect of clouds (Fig. 1d). Within the ablation zone, our high-resolution CRE product reveals that clouds have both a cooling and warming effect on average during our 18-year study period (Fig. 1f). CRE is negative near the lower elevation margins of the ablation zone (16% of the ablation zone area) but positive in the upper elevations of the ablation zone (84% of the ablation zone area; S1). These steep spatial gradients likely explain why previous studies have arrived at different conclusions about the role of clouds on ice sheet melt5,7.

Fig. 1: High-resolution (1 km × 1 km) predictions of summer cloudiness, albedo and cloud radiative effects (CRE) over the Greenland Ice Sheet derived from fusing different satellite datasets with machine learning.
figure 1

a cloud fraction, defined as the frequency of cloud cover (b) cloudiness, defined as the percentage reduction of shortwave radiative fluxes at the surface due to clouds; (c) surface albedo; (d) cloud radiative effect on shortwave radiation received at the ice surface; (e) cloud radiative effect on longwave radiation received at the surface; (f) cloud radiative effect on net all-wave radiation at the surface. This product was derived by fusing CloudSat/CALIPSO radiative flux and MODIS cloud and surface albedo products.

Full size image

Although we find no statistically significant trends (either regionally or at the ice sheet scale) during our 2003–2020 study period, we find substantial interannual variability (std. dev. = ±1.6%) in summer “cloudiness” (Fig. S2a, b; defined as the percentage reduction of shortwave radiative fluxes at the surface due to clouds, see Methods). Cloudiness ranges from as high as 16.2% in the summer of 2018 to as low as 10.8% in the summer of 2015. Temporal variability in summer cloudiness is somewhat higher in the ablation zone (std. dev. = ±2.3%) than in the accumulation zone (std. dev. = ±1.5%). We find no significant correlation between summer cloudiness and mean summer air temperatures in the ablation zone (Fig. S3a and S4), indicating that cloudiness likely will not change substantially in a warmer climate. We also observe substantial interannual variability in summer albedo which ranges by 4.0% with much more variation in the ablation zone (std. dev. = ±3.4%) than the accumulation zone (std. dev. = ±0.7%) (Fig. S2c, d). Variations in albedo can be attributed to processes such as snowline fluctuations, ice algae growth, meltwater ponding, surface roughening, and cryoconite hole expansion, all of which are captured by the high-resolution remotely sensed albedo products used in this study (see Methods2,24. Unlike cloudiness, we find that mean summer albedo is strongly linearly correlated with mean summer air temperatures over our 2003–2020 study period (R2 = 0.76), indicating that the ablation zone will darken in a warmer climate (Fig. S3c).

As expected, we find a strong positive correlation between cloudiness and downward longwave radiative fluxes and a strong negative correlation between cloudiness and downward shortwave radiative fluxes in our high-resolution CRE product (Fig. 2). In the high-albedo accumulation zone, the net effect is that cloud longwave trapping overwhelms cloud shortwave shading. For every 1% change in summer cloudiness, downward longwave radiation changes by 2.67 W m−2 while net shortwave radiation changes by only 0.83 W m−2 (Fig. 2a, c). The ~1.8 W m−2 difference between the two competing effects signifies that variations in cloudiness exert an important control on the surface energy balance of the accumulation zone. However, since only a small amount of meltwater runoff is produced from the accumulation zone3, this finding has little or no effect on Greenland’s contribution to global sea-level.

Fig. 2: Correlative relationships between summer cloudiness and radiative fluxes at the ice sheet surface.
figure 2

a Net shortwave response to cloudiness in the accumulation zone, (b) Net shortwave response to cloudiness in the ablation zone, (c) Downward longwave response to cloudiness in the accumulation zone, (d) Downward longwave response to cloudiness in the ablation zone. Note that x and y axes are scaled to the same range. The slopes of these relationships signify the sensitivity of radiative fluxes to changes in summer cloudiness.

Full size image

Interestingly, we also find that clouds have relatively little impact on meltwater production in the ablation zone, but for different reasons than the accumulation zone. During our study period, the shortwave shading response to changing cloudiness is approximately equal to the corresponding response of longwave trapping (±1.39 W m−2 and ±1.55 W m−2, respectively) (Fig. 2b, d). The compensating effects of clouds (when averaged across the entire ablation zone) signify that variations in summer cloudiness currently only have a small effect (±0.16 W m−2) on surface energy balance and associated meltwater production from the Greenland Ice Sheet. This finding is supported by the fact that intense melt years (e.g. 2010, 2012, and 2019) were not significantly cloudier/clearer than low melt years (e.g. 2013, 2015) (Fig. S3a).

Cloud radiative effects during atmospheric blocking events

Since MODIS has daily temporal resolution, our high-resolution product enables investigation of the radiative effects of short-lived atmospheric circulation features occurring for days-to-weeks over Greenland. One of the most important of these for ice sheet surface melt is anticyclonic blocking events, which enhance ice sheet melt through cloud suppression and warm air advection25. Blocking events have been attributed to negative phases of the North Atlantic Oscillation (NAO) and have significantly increased in frequency over the past twenty years26,27. Ward et al. (2020)25 find that summer blocking occurred on 205 days (15%) during the 2003–2018 period using used from Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA-2) reanalysis data. The net radiative effects of summer blocking in the ablation zone, however, have yet to be directly observed.

We investigate the radiative effects of atmospheric blocking events by combining our high-resolution CRE product with blocking days identified by Ward et al. (2020)25 who locate blocking events into four quadrants (NW, NE, SE, SW) based on local maximum 500 hPa geopotential height anomalies. Depending on where the blocking is centered, we find that blocking events reduce cloudiness in the ablation zone by 0.2–6.4% percentage points (Table S1). When blocking events are centered over Northwest, Northeast, and Southwest Greenland, net radiative fluxes in the ablation zone increase by +11.5 ± 3.1 W m−2, +4.9 ± 1.3 W m−2 and +4.0 ± 1.1 W m−2, respectively (Table S1; Fig. 3), relative to mean summer conditions. These increased radiative fluxes can be attributed to a relatively large increase in downward shortwave radiation combined with a smaller decrease in downward longwave radiation due to warmer near-surface air temperatures25 (Fig. 3). When blocking events are centered over Southeast Greenland, net radiative fluxes are actually reduced by −1.6 ± 0.4 W m−2, because the loss of shortwave shading over the narrow Southeast Greenland ablation zone is countered by enhanced cloudiness over the larger Southwest and Northwest Greenland ablation zones (Fig. S5). Generally, these findings signify that the radiative effects of blocking events enhance meltwater production in the ablation zone, with their net effect varying depending on location.

Fig. 3: Radiative effects of all summer atmospheric blocking events on the surface of the Greenland Ice Sheet relative to mean (2003–2020) conditions.
figure 3

(a) cloudiness; (b) net shortwave radiative fluxes; (c) downward longwave radiative fluxes; (d) allwave radiative fluxes. Blocking days, which were provided by Ward et al.25, have a net positive (warming) effect on radiation receipt at the surface of the ablation zone.

Full size image

Cloud radiative effects in a warmer climate

In a warmer climate, the shortwave shading effect of clouds will strengthen in the ablation zone due to decreased surface albedo. To forecast future changes in ice sheet albedo, we used observed air temperature-albedo relationships to develop an empirical model that predicts albedo from future air temperature forecasts (see the section “Cloud radiative effects in the future”). By applying our model to a range of projected future air temperatures from CMIP6 global climate model simulations, we find that summer surface albedo in the ablation zone will darken by ~0.02–0.14 depending on the emissions scenario (Fig. S6). Interannual variations in cloudiness will therefore have a much greater effect on radiative fluxes on the ablation zone surface than they do presently. Under a high emissions scenario (SSP5-8.5), for example, we project that by 2100, a 1% change in cloudiness will alter net shortwave radiation by ±1.94 W m−2 (Fig. 4), which is much greater than its effect on downward longwave radiation (±1.55 W m−2; Fig. 2d). The difference between cloud shortwave shading and longwave trapping effects on the ablation zone (i.e. ±0.39 W m−2 in 2100 [1.55–1.94] versus ±0.16 W m−2 during our 2003–2020 study period [from Fig. 2; 1.39–1.55]) suggests that the sensitivity of ice sheet meltwater production to cloudiness will more than double in a future, warmer climate.

Fig. 4: Sensitivity of net shortwave radiation to changes in summer cloudiness for a range of climate scenarios.
figure 4

Black line represents the current observed sensitivity derived from our high-resolution radiative flux products with shaded areas representing uncertainty (see Methods). Red, orange, and yellow lines represent projected changes in sensitivity for SSP1-2.6, SSP2-4.5, SSP5-8.5 scenarios based on the ensemble mean summer air temperature of 26 models from the CMIP6 experiment. Shaded areas represent the 25th and 75th percentile projections of all models.

Full size image

This increasing sensitivity of Greenland ablation zone energy balance to cloudiness implies that, in a warmer climate, atmospheric blocking events will enhance ice sheet melt even more than they do currently. By 2100, in an SSP5-8.5 scenario, we compute that the same blocking events centered over Northwest, Northeast, and Southwest Greenland will increase net radiation fluxes in the ablation zone by +15.7 ± 4.2, +6.6 ± 1.8, and +5.3 ± 1.4 W m−2 respectively, relative to mean 2100 conditions (Table S1). An overall ~34% increased radiative impact of atmospheric blocking events compared with the current study period (+4.7 for 2003–2020 period vs. +6.3 W m−2 by 2100 in SSP5-8.5), together with enhanced sensible heat fluxes from warm air advection25, will make these anomalous circulation events extremely potent for ice sheet melt in the future. It is therefore critical that global climate models used to forecast future ice mass loss accurately represent the frequency and duration of atmospheric blocking events over the Greenland Ice Sheet.

To conclude, we used machine learning to merge accurate Cloudsat/CALIPSO cloud profiling retrievals with high temporal and spatial resolution MODIS visible-infrared imagery, enabling retrieval and investigation of fine-scale CRE variations across the Greenland Ice Sheet surface from 2003 to 2020. The high spatial resolution (1 × 1 km) of our CRE product represents a major advance for understanding interactions between clouds and specific areas of the ice sheet surface. Previous observational studies of CRE were limited to either very coarsely gridded satellite products (e.g. 2 × 2°) or sparsely located in situ instruments, most of which are located in the accumulation zone. Our high-resolution spatially-continuous product overcomes these limitations, thus enabling systematic assessment of radiative fluxes over the entire Greenland Ice Sheet, including the ablation zone where the vast majority of ice sheet meltwater runoff is produced.

Analysis of our fine-scale CRE product reveals an important emerging role for clouds on future meltwater production on the Greenland Ice Sheet surface. During our 2003−2020 study period, the competing radiative impacts of clouds (i.e. longwave warming vs. shortwave cooling) approximately compensated over the ablation zone, signifying that clouds currently have a small regulatory role on meltwater runoff production. But as the ablation zone darkens2, the shortwave shading effect of clouds will become increasingly critical for determining meltwater production and associated ice sheet contributions to global sea-level rise. By the century’s end, we project that surface energy balance in the Greenland Ice Sheet ablation zone will be twice as sensitive to clouds than it is today and that atmospheric circulation patterns that suppress clouds (i.e. blocking events) will enhance meltwater production ~34% more than they already do.

Methods

General approach

Our general approach for investigating the impacts of clouds on Greenland Ice Sheet mass balance from satellite remote sensing observations and machine learning is outlined here with further details provided in the section “Prediction of cloud enhancement factors”. Central to our approach is the computation of CRE which describes the instantaneous effect of clouds on the surface energy budget relative to clear skies28:

$${{{{{{\rm{CRE}}}}}}}_{{{{{{\rm{SW}}}}}}}=({{{{{{\rm{SW}}}}}}}_{\downarrow {{{{{\rm{all}}}}}}-{{{{{\rm{sky}}}}}}}-{{{{{{\rm{SW}}}}}}}_{\downarrow {{{{{\rm{clear}}}}}}-{{{{{\rm{sky}}}}}}}).\,(1-\alpha )$$
(1)
$${{{{{{\rm{CRE}}}}}}}_{{{{{{\rm{LW}}}}}}}=({{{{{{\rm{LW}}}}}}}_{\downarrow {{{{{\rm{all}}}}}}-{{{{{\rm{sky}}}}}}}-{{{{{{\rm{LW}}}}}}}_{\downarrow {{{{{\rm{clear}}}}}}-{{{{{\rm{sky}}}}}}})$$
(2)
$${{{{{{\rm{CRE}}}}}}}_{{{{{{\rm{NET}}}}}}}={{{{{{\rm{CRE}}}}}}}_{{{{{{\rm{SW}}}}}}}+{{{{{{\rm{CRE}}}}}}}_{{{{{{\rm{LW}}}}}}}$$
(3)

where SW↓all−sky and LW↓all−sky are the downward shortwave and longwave fluxes at the surface during all-sky conditions, SW↓clear−sky and LW↓clear−sky are the downward shortwave and longwave fluxes that would occur at the surface in the absence of clouds, and α is the surface albedo. A positive CRENET (in W m−2) indicates net cloud warming at the surface, whereas a negative CRENET indicates net cloud cooling.

We first derive CRENET using high-quality all-sky and clear-sky radiative flux observations from CloudSat/CALIPSO and surface albedo observations from MODIS. Yet, while this approach may be accurate, CloudSat/CALIPSO has a relatively narrow ground footprint (~1.4 km), long repeat interval (16 days), and short timespan (2006−2010) which limits its usefulness of for investigating CRE at fine spatiotemporal resolution over long time periods16,18. To overcome this, we train a machine learning model, using CloudSat/CALIPSO observations, to predict CRE from passive visible-infrared satellite observations (i.e. from MODIS). This approach involves computation of a cloud enhancement factor (F) from CloudSat/CALIPSO given as:

$${{{\rm{F}}}}_{{{\rm{SW}}}} = {\frac {{{{\rm{SW}}}}_{\downarrow {{{\rm{all}}}}-{{{\rm{sky}}}}}} {{{{\rm{SW}}}}_{\downarrow {{{\rm{clear}}}}-{{{\rm{sky}}}}}}}$$
(4)

and

$${{{\rm{F}}}}_{{{\rm{LW}}}} = {\frac {{{{\rm{LW}}}}_{\downarrow {{{\rm{all}}}}-{{{\rm{sky}}}}}} {{{{\rm{LW}}}}_{\downarrow {{{\rm{clear}}}}-{{{\rm{sky}}}}}}}$$
(5)

The cloud enhancement factor represents the change in downward shortwave or longwave radiation due to clouds relative to clear skies7. This value depends upon cloud optical thickness, phase, temperature, and height, as well as solar zenith angle for shortwave fluxes and near-surface air temperature for longwave fluxes7,29. Our supervised machine learning algorithm is trained to predict FSW and FLW from cloud properties derived from passive visible-infrared sensors such as MODIS. Since MODIS has much wider spatial coverage and shorter repeat intervals (i.e. sub-daily) than CloudSat/CALIPSO, this approach enables computation of observational based CRE at finer spatial resolution and longer time-scales (2003−2020 vs. 2006–2010) than was previously possible. Further details on our approach are provided in the section “Prediction of cloud enhancement factors”.

Data

We derived high-quality downward all-sky and clear-sky radiative fluxes at the surface from the CloudSat/CALIPSO Level-2 ‘Fluxes and Heating Rates’ (2B-FLXHR-LIDAR) product which was the first to use active remotely sensed cloud observations to retrieve surface radiative fluxes at a global scale30,31. This product incorporates (1) cloud and aerosol observations from CloudSat/CALIPSO; (2) ancillary temperature and humidity profiles from ECMWF atmospheric re-analyses; and (3) International Geosphere-Biosphere Programme surface albedo/emissivity data to constrain broadband radiative fluxes using the two-stream radiative transfer model “BugsRad”31. Since our study is focused on the summer, we only use 2B-FLXHR-LIDAR data acquired during June, July, and August. We use Release 05 version of the 2B-FLXHR-LIDAR product which provides a more robust partitioning of ice-only and liquid-bearing clouds than Release 047. The radiative fluxes of the 2B-FLXHR-LIDAR algorithm have been extensively evaluated7,31 and previously used to study cloud impacts on surface energy budgets in the polar regions7,11,18.

We derived summer cloud properties from the MODIS Level-2 Cloud product “MYD06_L2” Collection 6.1. This product uses remotely sensed infrared, visible and near-infrared reflected radiances to derive cloud properties such as cloud top temperature, cloud top height, effective emissivity, cloud phase (ice vs. water), cloud optical thickness, and effective particle radius19. These parameters are available at a spatial resolution of either 1 km or 5 km resolution (at nadir). Since NASA’s Aqua satellite orbits within two to three minutes of CloudSat/CALIPSO in the A-Train constellation, cloud properties from the MYD06_L2 product are near-coincident with the 2B-FLXHR-LIDAR product. We also derived surface albedo from the MODIS daily snow cover product, MOD10A1 Collection 6.1.2. These data were averaged to provide mean summer albedo for the Greenland Ice Sheet at 1 × 1 km grid cell resolution. Since satellite measurements of albedo are acquired under clear-sky conditions and snow/ice albedo is known to be higher under cloud cover, we applied a correction to our clear-sky albedo climatologies following Key et al.32. To this, we multiplied the mean cloud effect (+0.05) with the cloud fraction of each grid cell (Fig. 1a) and added this value to the clear-sky albedo climatology. Finally, hourly near-surface (2 m) air temperatures were derived from ERA5 climate reanalysis33.

Prediction of cloud enhancement factors

We predicted cloud enhancement factors for each MODIS pixel using a Random Forests regression algorithm trained from CloudSat/CALIPSO radiative flux observations. The approach was introduced in the section “General approach” but the specific steps are described here:

Response variable

We first computed cloud enhancement factors (FSW and FLW) from all CloudSat/CALIPSO downward all-sky and clear-sky radiative flux observations acquired over the Greenland Ice Sheet between 13 June 2006 and 18 August 2010 (Eqs. 4 and 5).

Predictive variables

For each CloudSat/CALIPSO cloud enhancement factor (n = 317,614), we acquired coincident cloud properties from the MYD06_L2 product and near-surface temperatures from ERA5 climate reanalysis using a nearest neighbors approach (k = 1). Cloud properties used to predict FSW include cloud optical thickness, cloud top pressure, cloud effective radius, cloud top temperature, cloud top height, and cloud water path from MODIS, and near-surface air temperature from ERA5. Cloud properties used to predict FLW include cloud optical thickness, cloud top pressure, cloud effective radius, cloud top temperature, cloud top height, cloud water path, and cloud phase from MODIS and near-surface air temperature from ERA5. Since downward shortwave all-sky radiative fluxes are highly dependent on solar zenith angle, we grouped the training data into 5° bins (between 35 and 80°) before statistical modeling.

Statistical modeling

For each group of training data (ten groups for FSW and one group for FLW), we fitted a Random Forests regression model with 100 trees to a random sample containing 50% of the training data (n = 158,807 for FLW and n = ~15,000 for each FSW group). We then used the regression models (which depend on solar zenith angle for the FSW model) to predict FSW and FLW for every valid pixel in every MYD06_L2 product acquired over the Greenland Ice Sheet during the 2003–2020 period. The importance of each variable for predicting FSW and FLW is provided in Table S2 and S3. We evaluated our models using the other 50% of the data that was not used for training (see the section “Error budget”).

Computation of metrics

For each predicted FSW and FLW, we derived corresponding downward clear-sky radiative fluxes from ERA5 reanalysis. To account for possible bias in the ERA5 reanalysis, we first performed linear regression between 2B-FLXHR-LIDAR and ERA5 reanalysis clear-sky radiative fluxes and calibrated ERA5 clear-sky fluxes using the slope and intercept of this function (Fig. S7). We then combined clear-sky fluxes with the aforementioned cloud enhancement factors (FSW and FLW) to compute downward all-sky fluxes (using Eqs. 4 and 5).

Next, we combined these downward all-sky fluxes with downward clear-sky fluxes and mean summer surface albedo derived in Ryan et al. (2019)2 to compute summer CRESW, CRELW and CRENET (using Eqs. 1, 2 and 3) for every pixel in every MYD06_L2 product acquired over the Greenland Ice Sheet during the 2003–2020 period. We also derive a “cloudiness” metric which was defined as the percentage reduction in shortwave radiation at the surface due to clouds (Eq. 6).

$${{{{{\rm{Cloudiness}}}}}} = \frac {\left({{{{{{\rm{SW}}}}}}}_{\downarrow {{{{{\rm{clear}}}}}}-{{{{{\rm{sky}}}}}}}-{{{{{{\rm{SW}}}}}}}_{\downarrow {{{{{\rm{all}}}}}}-{{{{{\rm{sky}}}}}}}\right)}{{{{{{{\rm{SW}}}}}}}_{\downarrow {{{{{\rm{clear}}}}}}-{{{{{\rm{sky}}}}}}}}$$
(6)

Radiative flux climatologies

Our satellite-derived CRE observations reflect cloud conditions at the time of satellite overpass which for parts of South Greenland is only once or twice per day (Fig. S8). This biases our radiative flux climatologies high because all-sky radiative fluxes are derived during daytime. To overcome this issue, we computed the mean difference between hourly and mean daily shortwave radiative fluxes in ERA5 between 2003 and 2020 for every ERA5 grid cell on the Greenland Ice Sheet. We then subtracted this correction from our instantaneous all-sky and clear-sky shortwave radiative fluxes. After calibration, our satellite-derived all-sky and clear-sky shortwave radiative fluxes represent daily means. Downward longwave radiation mostly depends on cloudiness and does not exhibit any predictable diurnal cycle in AWS data or ERA5 reanalysis, so we assumed our instantaneous satellite-derived all-sky and clear-sky longwave radiative fluxes represent daily means which introduces an uncertainty of ±11.6 W m−2 (Fig. S9).

We bilinearly resampled our products to a standard 1 × 1 km study grid used by the Greenland Surface Mass Balance Intercomparison Project (GrSMBMIP)34. We separated our results into the ablation and accumulation zones based on a bare ice presence index threshold described by Ryan et al. (2019)2. The ablation zone was defined as all areas of the ice sheet where the bare ice is present for more than 10% of the record-setting melt summer of 2012 with the rest of the ice sheet classified as the accumulation zone. The ablation zone represents 14.8% of the ice sheet area. We average our sub-daily CRE products for each summer (June, July, and August) during the 2003–2020 period, providing an 18-year record of cloud conditions and radiative effects over the Greenland Ice Sheet.

Error budget

The uncertainties in our MODIS-derived all-sky radiative fluxes are a combination of uncertainties in predicted cloud enhancement factors (FSW and FLW) and clear-sky radiative fluxes (Eqs. 4 and 5). We evaluated our predicted cloud enhancement factors against high-quality cloud enhancement factors derived from CloudSat/CALIPSO observations (n = 158,807 for FLW and n = ~15,000 for FSW) that were not used to train the machine learning algorithm. We found a mean absolute error (MAE) of ±0.075 (±11.7%) for predicted FSW (Fig. S10a) and ±0.036 (±2.6%) for predicted FLW (Fig. S10b).

We evaluated our clear-sky radiative fluxes by comparing the clear-sky fluxes from ERA5 reanalysis with the 2B-FLXHR-LIDAR product. We find a MAE of ±20.8 W m−2 (±5.6%) between clear-sky shortwave fluxes (Fig. S7a) and a MAE ± 13.7 W m−2 (±6.7%) for the clear-sky longwave fluxes (Fig. S7b and S9). We assigned an uncertainty of ±0.03 (±5.8%) to our MODIS-derived summer albedo product based on previous comparisons with surface albedo data from PROMICE AWS35. When these sources of error are combined, we compute an uncertainty of ±17.3% for downward SWall-sky, ±9.3% for downward Lwall-sky, ±23.9% for CRESW, ± 11.9% for CRELW, and ±26.7% for CRENET.

We carried out an independent check on our satellite-derived product by comparing it to high-quality CRE observations acquired by in situ instruments at Summit Station (72.68 °N, −38.58 °W)13. Between the period 2011–2013, Miller et al. (2015)13 observed a CRENET of +35.3 W m−2 during the summer which is very similar to what we find over the same time period (+32.3 ± 8.7 W m−2).

Cloud radiative effects in the future

While we find no evidence that cloudiness in the ablation zone will substantially change as the climate warms, the albedo of the ablation zone will likely reduce in a warmer climate due to migrating snowlines, algal growth, and enhanced meltwater production. To forecast future albedo in the ablation zone, we developed an empirical model that predicts albedo for each grid cell based on a reference albedo grid (2015–2020) and near-surface air temperatures from ERA5 climate reanalysis and 26 models from the CMIP6 climate model experiment. We first linearly regressed mean summer albedo from MODIS with downscaled air temperature from ERA5 for each 1 × 1 km grid cell for the 2003–2020 period. ERA5 summer air temperatures were downscaled onto the target 1 km DEM by using the Local Inverse Distance Weighting (IDW) Lapse Algorithm from the MeteoIO library36. Regional lapse rates were calculated from all ERA5 grid points within 100 km distance from the target grid point. The obtained lapse rates were used to remove the elevation gradient from the ERA5 grid points. We then interpolated the residuals over the 1 km target DEM by inverse distance weighting using a scale of 1 km, and the power value alpha set to 1. Finally, the locally calculated lapse rate is used again to add the elevation gradient to the 1 km DEM to the interpolated grid.

After downscaling, non-significant (>0.05) and unrealistic relationships (i.e. positive correlations) were removed from analysis. We then computed the mean slope of these linear relationships within nine temperature bands (Fig. S11). This means that the sensitivity of grid cell albedo to air temperature is dynamic in response to a warming climate. For example, the albedo of a grid cell that is permanently snow-covered during the 2015–2020 reference period might not be initially very sensitive to warmer summer air temperatures, because the only increasing snow grain size can modify the grid cell albedo. However, if sufficient warming (and melting) occurs, bare ice could be exposed in the same grid cell, dramatically reducing the albedo. Our empirical albedo scheme would capture the addition of this process by placing the grid cell into the next, more sensitive, temperature band that represent grid cells whose albedos respond to both snow grain size evolution and bare ice exposure. Our empirical model therefore captures expected non-linearity between air temperature and albedo. We evaluated our model using temperature data from before our reference period (i.e. 2003–2014). We find a strong correlation between observed and predicted albedo (R2 = 0.83) with an accuracy of ±0.04 (Fig. S12).

We forecasted the albedo of the ice sheet to 2100 using near-surface air temperatures provided by the CMIP6 climate model experiment for three Shared Socioeconomic Pathways (SSP) (SSP1-2.6, SSP2-4.5, SSP3-7.0 and SSP5-8.5). We used 26 models, namely: ACCESS-CM2, ACCESS-ESM1-5, CESM2-WACCM, CESM2, CMCC-CM2-SR5, CNRM-CM6-1-HR, CNRM-CM6-1, CNRM-ESM2-1, CanESM5-CanOE, CanESM5, EC-Earth3-Veg, EC-Earth3, GFDL-ESM4, GISS-E2-1-G, INM-CM4-8, INM-CM5-0, IPSL-CM6A-LR, KACE-1-0-G, MIROC-ES2L, MIROC6, MPI-ESM1−2-HR, MPI-ESM1-2-LR, MRI-ESM2-0, NorESM2-LM, NorESM2-MM, and UKESM1-0-LL. This ensemble has a size and breadth of models that is comparable to the ensembles used in parts of the IPCC AR6 WG1 report. For each model we resampled near-surface air temperatures onto a common study grid using bilinear interpolation. For each model we also normalized air temperature by subtracting the mean summer ice sheet air temperature to account for any temperature bias. We then derived albedo for each grid cell based on the temperature deviation from a 2015–2020 reference period and the empirical albedo scheme described above. During this process, grid cell albedo is prevented from going below 0.3, below which would be an unrealistically low albedo for bare ice (Fig. S6). Finally, we computed the linear relationship between cloudiness (see the section “ Computation of metrics”) and SWnet [i.e. SW↓all−sky * (1 – α)] based on the forecasted albedo in the ablation zone. This enabled us to assess the sensitivity of ice sheet meltwater production to changing cloudiness in a future, warmer climate (Fig. 4).