Satellite Remote Sensing of the Greenland Ice Sheet Ablation Zone: A Review

The Greenland Ice Sheet is now the largest land ice contributor to global sea level rise, largely driven by increased surface meltwater runoff from the ablation zone, i.e., areas of the ice sheet where annual mass losses exceed gains. This small but critically important area of the ice sheet has expanded in size by ~50% since the early 1960s, and satellite remote sensing is a powerful tool for monitoring the physical processes that influence its surface mass balance. This review synthesizes key remote sensing methods and scientific findings from satellite remote sensing of the Greenland Ice Sheet ablation zone, covering progress in (1) radar altimetry, (2) laser (lidar) altimetry, (3) gravimetry, (4) multispectral optical imagery, and (5) microwave and thermal imagery. Physical characteristics and quantities examined include surface elevation change, gravimetric mass balance, reflectance, albedo, and mapping of surface melt extent and glaciological facies and zones. The review concludes that future progress will benefit most from methods that combine multi-sensor, multi-wavelength, and cross-platform datasets designed to discriminate the widely varying surface processes in the ablation zone. Specific examples include fusing laser altimetry, radar altimetry, and optical stereophotogrammetry to enhance spatial measurement density, cross-validate surface elevation change, and diagnose radar elevation bias; employing dual-frequency radar, microwave scatterometry, or combining radar and laser altimetry to map seasonal snow depth; fusing optical imagery, radar imagery, and microwave scatterometry to discriminate between snow, liquid water, refrozen meltwater, and bare ice near the equilibrium line altitude; combining optical reflectance with laser altimetry to map supraglacial lake, stream, and crevasse bathymetry; and monitoring the inland migration of snowlines, surface melt extent, and supraglacial hydrologic features.


Introduction
The Greenland Ice Sheet (GrIS) is the second-largest ice mass on earth. If the entire ice sheet were to melt, global sea level would rise by about seven meters [1]. At present, the GrIS is losing mass, consistent with its expected response to anthropogenic climate warming [2][3][4][5][6][7][8][9]. Its mass loss averaged -171 Gt yr −1 between 1991 and 2015, equivalent to 0.47 ± 0.23 mm yr −1 global eustatic sea level rise (SLR) [1]. The rate of SLR contribution accelerated to >1 mm yr −1 during the latter half of this period with a maximum 1.2 mm contribution in the record melting year 2012 [1,[10][11][12][13]. Since 2013, satellite gravimetry data suggest mass loss acceleration has stalled but remains negative with an average 0.69 mm yr −1 SLR contribution during the 2012-2016 period [10,14,15]. Consequently, there is an urgent

Ice Sheet Mass Balance, Surface Mass Balance, and Energy Balance
Before reviewing satellite remote sensing of the GrIS ablation zone, we briefly define key concepts and terminology that are discussed throughout the review. Ice sheets and glaciers gain mass through snowfall and lose mass through sublimation, meltwater runoff, and solid ice discharge. Following Lenaerts et al. [21] and van den Broeke et al. [60], the mass balance MB of an ice sheet is usually written as: where SMB is surface mass balance and D is solid ice discharged to surrounding oceans. The components of MB are commonly expressed in units of mass per time (e.g., kg yr −1 or Gt yr −1 ), mass per unit area per time (e.g., kg m −2 yr −1 ), and mass per unit area per time normalized by density of liquid water (e.g., m water equivalent, m w.e. yr −1 ) [21]. The SMB is the sum of mass inputs and outputs at the ice surface: where precipitation P is the sum of liquid precipitation P liquid (e.g., rain) and solid precipitation P solid (e.g., snow), E is the sum of evaporation and sublimation, ER ds is net snow erosion, and runoff R is the sum of the liquid water balance: where M is surface meltwater production, CO is condensation, RF is refreezing, and RT is liquid water retention, e.g., water stored in lakes, aquifers, or under capillary tension in snow and ice. The GrIS ablation zone SMB is dominated by P solid and M with the remaining terms representing smaller but uncertain components [1,[82][83][84][85]. The accumulation and ablation zones of the ice sheet are defined as those areas where the local SMB > 0 and the local SMB < 0, respectively. The equilibrium line altitude (ELA) separates the two zones where local SMB = 0 [60]. The total change in ice sheet surface elevation over time (the quantity observed by spaceborne altimeters) is usually defined as: where dH/dt is the total change in surface elevation (m yr −1 ), SMB (kg m −2 yr −1 ) is surface mass balance, ρ s (kg m −3 ) is surface snow, firn, and/or ice density, BMB (kg m −2 yr −1 ) is the change in basal mass balance, ρ b (kg m −3 ) is basal ice density, D (kg m −2 yr −1 ) is the local change in ice mass due to ice dynamic motion, ρ i is solid ice density, v c (m yr −1 ) is vertical velocity due to snow and/or firn compaction, and v br (m yr −1 ) is vertical bedrock velocity (e.g., glacial isostatic adjustment) [86,87]. The vertical velocity terms contribute to dH/dt, but are not associated with change in MB. Equation (4) is used to estimate volumetric changes in the polar ice sheets, and also to estimate MB. This requires correction for the vertical velocity terms and an estimate of ρ s , which for the ablation zone is typically taken as a constant value between 830 kg m −3 and 917 kg m −3 [88], whereas for snow and firn above the ELA, measured or modeled values of ρ s are necessary to convert dH/dt to SMB [86]. In some cases, all but the first term on the right-hand side of Equation (4) can be neglected and the local SMB can be inferred from remotely sensed dH/dt [89].
The energy available for melting snow or ice is determined by the sum of positive and negative heat fluxes at the surface, referred to as the surface energy balance SEB (W m −2 ): s f c + SHF + LHF + G S (5) where SW and LW are shortwave (solar) and longwave (terrestrial) radiation fluxes, SHF and LHF are the turbulent surface fluxes of sensible and latent heat (which are proportional to the aerodynamic surface roughness length, ζ (m −1 )), G S is subsurface (conductive) heat flux, α is the broadband surface albedo, or the ratio of upwelling (reflected) solar radiation to downwelling (incoming) solar radiation (SW up /SW in ), σ is the Stefan Boltzmann constant (5.67 × 10 −8 W m −2 K −2 ), and T sfc is the ice surface kinetic temperature. In Equation (4), all heat fluxes have units W m −2 and are defined as positive toward the surface [60]. If T sfc reaches the melting temperature of ice, positive SEB provides melt energy ME (W m −2 ) and liquid meltwater M is produced. The dominant source of ME for the GrIS is SW NET [90]. Typical values of α are 0.9 for freshly fallen snow, 0.6 for melting snow, and 0.4 for bare glacier ice [73]. This wide variability in α demonstrates the critical role of albedo in modulating the SEB and consequently M in the GrIS ablation zone [26,29,30].

Ice Surface Elevation Change
The magnitude of GrIS dH/dt and its spatio-temporal variability has been quantified for nearly four decades using spaceborne radar (Radio Detection and Ranging) altimetry [37,38,[91][92][93][94][95], and since 2001, using spaceborne laser altimetry [96][97][98][99][100][101]. Recent advances in laser and radar altimeter precision and spatial resolution allow discrimination of dH/dt in slower-flowing, land-terminating sectors of the GrIS ablation zone, principally caused by a change in SMB rather than D [102,103]. For the accumulation zone, regional climate models and firn compaction models are used to partition the component of dH/dt caused by a change in SMB from those caused by a change in firn density and D [104,105]. The following subsections summarize spaceborne radar altimetry of the GrIS, providing a historical perspective that emphasizes major challenges and advances in the field, and identifies future opportunities for radar altimetry observations of the GrIS ablation zone.

Radar Altimetry
Conventional radar altimeters are nadir pointing single-beam radars that transmit and receive tens to thousands of microwave pulses per second, yielding ground measurement footprint diameters on the order of 1-10 km posted at 0.3-0.6 km along-track spacing [39,106,107]. An estimate of dH/dt is obtained by comparing measurements at crossover points, where an earlier satellite ground track intersects a later one [92]. Two critical factors affect the accuracy of radar altimeter elevation retrievals over ice sheets: (1) variation in topographic slope, which controls the average distance between the altimeter and the measured ground footprint surface [39,108,109], and (2) changes in the electrical permittivity of the ice sheet surface, which controls the reflection, transmission, and absorption of the microwave signal [110][111][112][113].
For typical radar altimeter frequencies (Ku-band~13 GHz), solid ice and liquid water are highly reflective, whereas dry snow and firn are semi-transparent, with Ku-band penetration depths on the order of several meters [114,115]. Spatial and temporal variations in signal penetration depth caused by seasonal bare ice exposure, snow and firn densification, surface meltwater presence, and refreezing of meltwater introduce spurious height change signals that are poorly quantified for the GrIS but are estimated to exceed >0.50 m in some cases [116][117][118]. Variations in surface slope dominate radar altimeter accuracy in the topographically-complex bare-ice ablation zone, with slope-induced errors that exceed >20 m for slopes >1 • [39,109,113]. The recent CryoSat-2 radar altimeter mission employs synthetic aperture radar (SAR) and interferometric synthetic aperture radar (InSAR) technologies, which decrease the effective along-track footprint diameter to the order 0.1-0.3 km and provide the cross-track slope from InSAR processing [93]. Together with denser orbital-track spacing and advances in waveform retracking, topographic mapping of ice sheet ablation zone areas is accurate to within a few meters on average, and dH/dt to within sub-meter uncertainty relative to laser altimetry [93,119,120].

Radar Altimetry Sensors and Datasets
The first satellite altimeter measurements of surface topography for the Greenland Ice Sheet were made by the GEOS-3, Geosat, and Seasat oceanographic Ku-band altimeters (Table 1) [37,91,106,121,122]. Their geographical coverage was limited to ±72 • , but the data they collected showed it was possible to calculate volumetric changes in the polar ice sheets from space, providing proof of the concept for future polar-orbiting altimeters [92]. A first demonstration using GEOS-3, Geosat, and Seasat data found that for areas south of <72 • N, the GrIS thickened by 23 ± 6 cm yr −1 between 1978-1986, with enhanced snowfall in a warmer climate hypothesized to explain the observed thickening [123]. These early missions supported important developments in waveform retracking, slope correction, and physical and empirical models for subsurface volume scattering [106,108,111,124]. These corrections are critical for accurate elevation retrieval and reliable change detection [125,126]. For example, Davis et al. [127] used an improved geodetic model and an improved retracking model to reinterpret the earlier findings al meters [114,115]. Spatial and temporal variations in signal penetration depth caused are ice exposure, snow and firn densification, surface meltwater presence, and eltwater introduce spurious height change signals that are poorly quantified for the estimated to exceed >0.50 m in some cases [116][117][118]. Variations in surface slope r altimeter accuracy in the topographically-complex bare-ice ablation zone, with slopes that exceed >20 m for slopes >1° [39,109,113]. The recent CryoSat-2 radar altimeter ys synthetic aperture radar (SAR) and interferometric synthetic aperture radar (InSAR) hich decrease the effective along-track footprint diameter to the order 0.1-0.3 km and track slope from InSAR processing [93]. Together with denser orbital-track spacing and aveform retracking, topographic mapping of ice sheet ablation zone areas is accurate w meters on average, and / to within sub-meter uncertainty relative to laser 19,120].
ltimetry Sensors and Datasets satellite altimeter measurements of surface topography for the Greenland Ice Sheet the GEOS-3, Geosat, and Seasat oceanographic Ku-band altimeters (Error! Reference nd.) [37,91,106,121,122]. Their geographical coverage was limited to ±72 o but the data showed it was possible to calculate volumetric changes in the polar ice sheets from ng proof of concept for future polar-orbiting altimeters [92]. A first demonstration , Geosat, and Seasat data found that for areas south of <72 o N, the GrIS thickened by 23 een 1978-1986, with enhanced snowfall in a warmer climate hypothesized to explain hickening [123]. These early missions supported important developments in waveform pe correction, and physical and empirical models for subsurface volume scattering 4]. These corrections are critical for accurate elevation retrieval and reliable change ,126]. For example, Davis et al. [127] used an improved geodetic model and an acking model to reinterpret the earlier findings of Zwally et al. [92], finding a smaller r -1 thickening for the period 1978-1986 that highlighted the importance of slope-and tion-correction. ummary of radar and laser altimeter remote sensing platforms, instruments, temporal observed wavelength, and managing agencies.  (Table A1) and may not reflect joint collaborations.
The ESA European Remote-sensing Satellite (ERS-1) carried the first polar-orbiting radar altimeter, extending coverage of the polar ice sheets to ±81.5 o [95]. Together with the follow-on ERS-2 and EnviSat missions, these satellites provided near complete GrIS coverage and 21 years of continuous data acquisition [38]. The first digital elevation model (DEM) spanning the entire GrIS surface was produced from combined ERS-1 and Geosat altimetry data, supplemented by stereophotogrammetric and cartographic data sets, with pan-GrIS average accuracy -0.33 ± 6.97 m [128,129]. ERS-1/2 altimetry data were used to discriminate thinning rates in the ablation zone (defined as <1500 m a.s.l.) of -2.0 ± 0.9 cm yr -1 from accumulation zone (>1500 m a.s.l.) thickening of 6.4 ± 0.2 cm yr -1 with an overall net thickening of 5.4 ± 0.2 cm yr -1 for the period 1992-2003 [130]. To extend the temporal coverage and increase spatial measurement density, Khvorostovsky and Johannessen [131] developed an algorithm to merge ERS-1/2 and EnviSat datasets and detect and eliminate intersatellite biases. For the period 1992-2008, their dataset suggests net thinning of the GrIS initiated around 2000, and accelerated between 2006-2008, with a thickening of 4.0 ± 0.2 cm yr -1 above 1500 m a.s.l. and thinning of -7.0 ± 1.0 cm yr -1 below 1500 m a.s.l. [94].
Conventional radar altimeters such as the instruments included on ERS-1/2 and EnviSat were optimized for detecting elevation over open oceans and low-gradient polar ice sheet interiors with mean accuracies typically <0.2 m for ice sheet surface slopes <0.2 o [39,109]. Elevation accuracy and / uncertainty is much higher for the topographically complex ablation zone due to their large spatial measurement footprint and wide orbital track spacing, with slope errors >20 m for slopes >1 o and systematic biases over rough surfaces [39,95,109,132]. Careful processing has been used to discriminate ablation zone thinning rates using ERS-1/2 and EnviSat altimetry data, including repeattrack analysis, but the accuracy of these data and their process-based interpretation is ambiguous owing to the aforementioned slope errors and sparse spatial measurement density in the ablation zone [118,[133][134][135].
The CryoSat-2 satellite (launched April 2010) provides continuity with ERS-1/2 and EnviSat and employs SAR delay-Doppler along-track processing and interferometric cross-track processing [120,136]. CryoSat-2 is the first ESA radar altimeter dedicated to polar studies, with orbital coverage to ±88 o and 1.6 km cross-track spacing at 70 o . CryoSat-2 carries two Ku-band SAR/Interferometric Radar Altimeter (SIRAL) systems. SIRAL operates in three distinct modes depending on location: 1) low resolution mode (LRM) over open ocean and ice sheet interiors, 2) synthetic aperture radar mode (SAR) over sea ice and coastal regions, and 3) SAR interferometric mode (SARIn) over ice sheet margins. The SARIn mode employs two SAR instruments oriented across the satellite track. Interferometric phase processing of the dual waveforms provides cross-track slope at ~0.3 km alongtrack resolution [137,138].
Over mild slopes (<1 o ), phase unwrapping [136] of CryoSat-2 interferometric data provides 5 km wide-swath elevation retrievals with two orders of magnitude more individual measurements per See Abbreviations at end of article for expanded acronyms.  (Table A1) and may not reflect joint collaborations.
The ESA European Remote-sensing Satellite (ERS-1) carried the first polar-orbiting radar altimeter, extending coverage of the polar ice sheets to ±81.5 o [95]. Together with the follow-on ERS-2 and EnviSat missions, these satellites provided near complete GrIS coverage and 21 years of continuous data acquisition [38]. The first digital elevation model (DEM) spanning the entire GrIS surface was produced from combined ERS-1 and Geosat altimetry data, supplemented by stereophotogrammetric and cartographic data sets, with pan-GrIS average accuracy -0.33 ± 6.97 m [128,129]. ERS-1/2 altimetry data were used to discriminate thinning rates in the ablation zone (defined as <1500 m a.s.l.) of -2.0 ± 0.9 cm yr -1 from accumulation zone (>1500 m a.s.l.) thickening of 6.4 ± 0.2 cm yr -1 with an overall net thickening of 5.4 ± 0.2 cm yr -1 for the period 1992-2003 [130]. To extend the temporal coverage and increase spatial measurement density, Khvorostovsky and Johannessen [131] developed an algorithm to merge ERS-1/2 and EnviSat datasets and detect and eliminate intersatellite biases. For the period 1992-2008, their dataset suggests net thinning of the GrIS initiated around 2000, and accelerated between 2006-2008, with a thickening of 4.0 ± 0.2 cm yr -1 above 1500 m a.s.l. and thinning of -7.0 ± 1.0 cm yr -1 below 1500 m a.s.l. [94].
Conventional radar altimeters such as the instruments included on ERS-1/2 and EnviSat were optimized for detecting elevation over open oceans and low-gradient polar ice sheet interiors with mean accuracies typically <0.2 m for ice sheet surface slopes <0.2 o [39,109]. Elevation accuracy and / uncertainty is much higher for the topographically complex ablation zone due to their large spatial measurement footprint and wide orbital track spacing, with slope errors >20 m for slopes >1 o and systematic biases over rough surfaces [39,95,109,132]. Careful processing has been used to discriminate ablation zone thinning rates using ERS-1/2 and EnviSat altimetry data, including repeattrack analysis, but the accuracy of these data and their process-based interpretation is ambiguous owing to the aforementioned slope errors and sparse spatial measurement density in the ablation zone [118,[133][134][135].
The CryoSat-2 satellite (launched April 2010) provides continuity with ERS-1/2 and EnviSat and employs SAR delay-Doppler along-track processing and interferometric cross-track processing [120,136]. CryoSat-2 is the first ESA radar altimeter dedicated to polar studies, with orbital coverage to ±88 o and 1.6 km cross-track spacing at 70 o . CryoSat-2 carries two Ku-band SAR/Interferometric Radar Altimeter (SIRAL) systems. SIRAL operates in three distinct modes depending on location: 1) low resolution mode (LRM) over open ocean and ice sheet interiors, 2) synthetic aperture radar mode (SAR) over sea ice and coastal regions, and 3) SAR interferometric mode (SARIn) over ice sheet margins. The SARIn mode employs two SAR instruments oriented across the satellite track. Interferometric phase processing of the dual waveforms provides cross-track slope at ~0.3 km alongtrack resolution [137,138].
Over mild slopes (<1 o ), phase unwrapping [136] of CryoSat-2 interferometric data provides 5 km wide-swath elevation retrievals with two orders of magnitude more individual measurements per Managing agencies are identified by the WMO OSCAR database (Table A1) and may not reflect joint collaborations.
The ESA European Remote-Sensing Satellite (ERS-1) carried the first polar-orbiting radar altimeter, extending coverage of the polar ice sheets to ±81.5 • [95]. Together with the follow-on ERS-2 and EnviSat missions, these satellites provided near-complete GrIS coverage and 21 years of continuous data acquisition [38]. The first digital elevation model (DEM) spanning the entire GrIS surface was produced from combined ERS-1 and Geosat altimetry data, supplemented by stereo-photogrammetric and cartographic data sets, with pan-GrIS average accuracy of −0.33 ± 6.97 m [128,129]. ERS-1/2 altimetry data were used to discriminate thinning rates in the ablation zone (defined as <1500 m a.s.l.) of -2.0 ± 0.9 cm yr −1 from accumulation zone (>1500 m a.s.l.) thickening of 6.4 ± 0.2 cm yr −1 with an overall net thickening of 5.4 ± 0.2 cm yr −1 for the period 1992-2003 [130]. To extend the temporal coverage and increase spatial measurement density, Khvorostovsky and Johannessen [131] developed an algorithm to merge ERS-1/2 and EnviSat datasets and detect and eliminate inter-satellite biases. For the period 1992-2008, their dataset suggests net thinning of the GrIS initiated around 2000, and accelerated between 2006-2008, with a thickening of 4.0 ± 0.2 cm yr −1 above 1500 m a.s.l. and thinning of -7.0 ± 1.0 cm yr −1 below 1500 m a.s.l. [94].
Conventional radar altimeters, such as the instruments included on ERS-1/2 and EnviSat, were optimized for detecting elevation over open oceans and low-gradient polar ice sheet interiors with mean accuracies typically of <0.2 m for ice sheet surface slopes <0.2 • [39,109]. Elevation accuracy and dH/dt uncertainty are much higher for the topographically complex ablation zone due to their large spatial measurement footprint and wide orbital track spacing, with slope errors >20 m for slopes >1 o and systematic biases over rough surfaces [39,95,109,132]. Careful processing has been used to discriminate ablation zone thinning rates using ERS-1/2 and EnviSat altimetry data, including repeat-track analysis, but the accuracy of these data and their process-based interpretation is ambiguous owing to the aforementioned slope errors and sparse spatial measurement density in the ablation zone [118,[133][134][135].
The CryoSat-2 satellite (launched April 2010) provides continuity with ERS-1/2 and EnviSat, and employs SAR delay-Doppler along-track processing and interferometric cross-track processing [120,136]. CryoSat-2 is the first ESA radar altimeter dedicated to polar studies, with orbital coverage to ±88 • and 1.6 km cross-track spacing at 70 • . CryoSat-2 carries two Ku-band SAR/Interferometric Radar Altimeter (SIRAL) systems. SIRAL operates in three distinct modes depending on location: (1) low resolution mode (LRM) over open ocean and ice sheet interiors, (2) synthetic aperture radar mode (SAR) over sea ice and coastal regions, and (3) SAR interferometric mode (SARIn) over ice sheet margins. The SARIn mode employs two SAR instruments oriented across the satellite track. Interferometric phase processing of the dual waveforms provides a cross-track slope at~0.3 km along-track resolution [137,138].
Over mild slopes (<1 • ), phase unwrapping [136] of CryoSat-2 interferometric data provides 5 km wide-swath elevation retrievals with two orders of magnitude more individual measurements per surface area than any prior radar altimeter [119,120,137]. The first spatially continuous swath DEMs were generated from CryoSat-2 data for areas of the Devon Ice Cap and western GrIS with ±3 m precision [119,137]. Swath processing provides information about within-swath topography that is not provided by traditional radar processing techniques. For example, Gray et al. [119] demonstrated a novel method for supraglacial lake water surface elevation retrievals in the ablation zone from CryoSat-2 swath retrievals and provided a detailed technical description of swath processing.
The SIRAL system nominally resolves dH/dt to 3.3 cm yr −1 near the ice sheet margins and 0.7 cm yr −1 in the ice sheet interior, at a 104 km 2 and 106 km 2 spatial scale, respectively [138]. Using three years (2011-2014) of CryoSat-2 standard retrievals and an updated threshold retracking algorithm, a new pan-GrIS DEM was generated with a 3 ± 15 m elevation accuracy relative to ICESat over 80% of the GrIS and 5 ± 65 m accuracy over 90% of the GrIS (Figure 1) [93]. These data suggest reliable CryoSat-2 elevation retrievals remain limited to areas of the ablation zone with a surface slope <1.5 o (Figure 1). Elevation change from these data suggests a 2.5 factor increase in pan-GrIS volume loss (-375 ± 24 km 3 yr −1 ) compared with the ICESat (2003-2009) period, with large losses concentrated in the west and southeast marginal ablation zones [93]. A separate analysis of the same data combined with a firn density model found an equivalent mass loss of 269 ± 51 Gt yr −1 [105].
Currently operating beyond its design lifetime, CryoSat-2 will be succeeded by the SAR Radar Altimeter (SRAL) onboard Sentinel-3A/B, and by the AltiKa Ka-band radar altimeter onboard the Satellite with ARgos and ALtiKa (SARAL) ( Table 1). SARAL is a joint French Space Agency (CNES) and Indian Space Research Organisation (ISRO) oceanographic mission launched in 2013 with the secondary goal of monitoring polar ice sheet surface elevations [139]. SARAL has the same 35-day repeat orbit as EnviSat and ERS-1/2, but the AltiKa instrument operates at Ka-band (~36 GHz), thereby offering unique opportunities for cross-platform validation of prior Ku-band radar altimeters [140]. In particular, the Ka-band should improve diagnosis of Ku-band signal penetration bias, owing to its higher (~36 GHz) signal frequency with~10 × lower theoretical penetration depth into dry snow relative to Ku-band, higher spatial resolution (~8 km footprint and 175 m along-track spacing), and higher pulse-repetition frequency [135,[141][142][143]. SARAL also provides an independent record of ice sheet dH/dt. Using the 1 km ICESat DEM for the period 2003-2005 [144] as a baseline, data from SARAL-AltiKa for the period 2014-2016 was used to estimate a pan-GrIS volume loss rate of 247 km 3 yr −1 , with the largest basin-scale volumetric decreases found in the north and northwest [145].  [93], with slope indicated by shaded relief; (b) median error (difference between CryoSat-2 DEM and ICESat elevations) vs. surface slope, with elevation indicated by colormap in (a); (c) standard deviation of error vs. surface slope, with elevation indicated by colormap in (a). The median and standard deviation of error are calculated from the CryoSat-2 DEM error grid binned by a slope with a 0.01 o bin size, following Helm et al. [93] (Figure 9). The median and standard deviation of error increase with slopes with a step shift toward higher error at slopes >1.5 o . Elevation and slope values calculated from the DEM indicate that 16% of the ice sheet area has a surface slope >1.5 o , and nearly all (95%) such areas have elevation <2000 m a.s.l., which is generally representative of the ablation zone (calculations performed by the first author). The error grid was produced by first calculating elevation differences between the CryoSat-2 DEM and individual ICESat elevations (campaign 3F, 3G, and 3H), corrected for elevation change between the individual ICESat observation and the DEM reference time (1 Jul 2012) and then calculating a weighted error as a function of surface roughness, surface slope, and number of cross-validation data points [93]. The DEM and associated error grid are publicly available at https://doi.pangaea.de/10.1594/PANGAEA.831394 [146].

Current Challenges and Future Opportunities
Progress in waveform retracking and SAR/InSAR technology has improved radar altimetry performance over complex terrain, but retracking errors over rough surfaces and spatiotemporal variations in signal penetration depth remain important sources of uncertainty. Reliable performance in the ablation zone appears limited to areas with surface slope <1.5° (Figure 1), and signal penetration biases may exceed 0.5 m in the lower accumulation zone [116,118,119,147]. Nevertheless, radar altimetry provides the longest record of GrIS / , and two challenges stand out for leveraging this unique dataset: (1) homogenization of the inter-mission and cross-platform methods and datasets used to detect / [133,148], and (2) better understanding of spatial and temporal changes in snow, firn, and ice permittivity caused by surface meltwater presence, meltwater percolation, and ice lens formation [116,147]. These challenges are closely related. For example, methods such as waveform smoothing, waveform model fitting, and optimal thresholding increase absolute accuracy, but at the expense of the homogenization required for statistical change detection [125,126,149]. In other cases, waveform retracking may mask real / signals [150]. Similarly, real changes in snow and firn permittivity may introduce false / signals, such as the apparent 56 ± 26 cm elevation increase attributed to ice lens formation in the percolation zone following the 2012 melt event [116]. Here, the SARAL-AltiKa Ka-band altimeter may provide new opportunities for diagnosing radar penetration bias and waveform interpretation. In addition, its replication of the  (Figure 9). The median and standard deviation of error increase with slopes with a step shift toward higher error at slopes >1.5 o . Elevation and slope values calculated from the DEM indicate that 16% of the ice sheet area has a surface slope >1.5 o , and nearly all (95%) such areas have elevation <2000 m a.s.l., which is generally representative of the ablation zone (calculations performed by the first author). The error grid was produced by first calculating elevation differences between the CryoSat-2 DEM and individual ICESat elevations (campaign 3F, 3G, and 3H), corrected for elevation change between the individual ICESat observation and the DEM reference time (1 July 2012) and then calculating a weighted error as a function of surface roughness, surface slope, and number of cross-validation data points [93]. The DEM and associated error grid are publicly available at https://doi.pangaea.de/10.1594/PANGAEA.831394 [146].

Current Challenges and Future Opportunities
Progress in waveform retracking and SAR/InSAR technology has improved radar altimetry performance over complex terrain, but retracking errors over rough surfaces and spatiotemporal variations in signal penetration depth remain important sources of uncertainty. Reliable performance in the ablation zone appears limited to areas with surface slope <1.5 • (Figure 1), and signal penetration biases may exceed 0.5 m in the lower accumulation zone [116,118,119,147]. Nevertheless, radar altimetry provides the longest record of GrIS dH/dt, and two challenges stand out for leveraging this unique dataset: (1) homogenization of the inter-mission and cross-platform methods and datasets used to detect dH/dt [133,148], and (2) better understanding of spatial and temporal changes in snow, firn, and ice permittivity caused by surface meltwater presence, meltwater percolation, and ice lens formation [116,147]. These challenges are closely related. For example, methods such as waveform smoothing, waveform model fitting, and optimal thresholding increase absolute accuracy, but at the expense of the homogenization required for statistical change detection [125,126,149]. In other cases, waveform retracking may mask real dH/dt signals [150]. Similarly, real changes in snow and firn permittivity may introduce false dH/dt signals, such as the apparent 56 ± 26 cm elevation increase attributed to ice lens formation in the percolation zone following the 2012 melt event [116]. Here, the SARAL-AltiKa Ka-band altimeter may provide new opportunities for diagnosing radar penetration bias and waveform interpretation. In addition, its replication of the ERS-1/2 and EnviSat repeat orbit will improve the statistical reliability of along-track change detection [140].
In contrast to the percolation zone, the ablation zone is typically conceptualized as homogenous solid ice with uniform electromagnetic properties. Ku-band backscatter over the bare ice ablation zone is dominated by surface scattering from the rough (and seasonally saturated) ice surface, but volume scattering has received little direct study and may be important over superimposed ice, seasonally-snow covered surfaces, or areas with multi-modal surface roughness distributions [151]. Penetration depths at C-and L-band exhibit an order of magnitude range over bare ice surfaces on the GrIS, possibly reflecting variations in ice thermal structure and surface roughness [152]. Similar comparisons at the Ku-band do not exist to our knowledge but could be facilitated by dual-frequency radar altimeters such as SRAL or RA-2, or by comparison with SARAL-AltiKa Ka-band altimetry ( Table 1). Factors that affect bare ice microwave permittivity include its grain size, temperature, porosity, water content, crystal structure, and chemical and physical impurity content [153][154][155][156]. Although slope errors dominate elevation retrieval uncertainty in the ablation zone, seasonal and spatial variations in ablation zone surface properties and their effect on radar backscatter may be an overlooked source of uncertainty [117,135]. Conversely, this variation and its effect on waveform shape could provide new opportunities for understanding the ablation zone, such as detecting water surface elevation change in supraglacial lakes, inferring changes in near-surface ice structure related to physical weathering of solid ice, or estimating seasonal snow thickness using a dual-frequency radar, as was recently demonstrated for Arctic sea ice using combined CryoSat-2 and SARAL-AltiKa retrievals [119,156,157].
At the ice-sheet scale, knowledge of snow, firn, and ice density is the main source of uncertainty for conversion between ice sheet volume change and mass change [36]. Density changes are caused by snow and firn compaction and by meltwater refreezing, which both affect waveform interpretation by changing the scattering properties of the snow and firn [105,116]. Field investigations suggest ice lenses are widespread in the percolation zone and have thickened in recent decades [20,158]. The vertical and horizontal distribution of meltwater percolation and refreezing is difficult to model and may not be accurately represented by regional climate models [158,159]. Although these features present challenges for dH/dt detection, they also provide unique opportunities for characterizing snow and firn processes, including detection of extreme melt events, ice lens formation, and snow accumulation rates following ice lens formation [115].

Laser Altimetry
Laser altimeters use lidar (light detection and ranging) technology to measure the two-way travel time of narrow-beam monochromatic laser pulses transmitted between the altimeter and the earth's surface. In contrast to radar altimetry, signal penetration in ice and snow is minimal at common lidar wavelengths (visible and near-infrared), and the narrow laser beam illuminates a much smaller ground area, which reduces slope errors over complex terrain [39,160]. The small laser footprint increases absolute accuracy but may also introduce uncertainty because interpolation between narrow footprints is needed to obtain spatially continuous elevations [99]. Laser altimeters are sensitive to atmospheric scattering and lack the all-weather capability of radar altimeters but are uniquely adept at measuring local surface elevation in the topographically complex ablation zone, allowing resolution of thinning rates at the scale of individual outlet glaciers and providing benchmark datasets for altimeter accuracy [97,161]. The following subsections summarize spaceborne laser altimetry of the GrIS, emphasizing the operational performance of ICESat, and future opportunities for laser altimetry observations of the GrIS ablation zone from the recently launched ICESat-2.

Laser Altimetry Sensors, Methods, and Datasets
The NASA Lidar In-Space Technology Experiment (LITE), flown on the Discovery Space Shuttle in September 1994, was the first spaceborne lidar [162]. The LITE was focused on the vertical structure of clouds and aerosols and provided operational proof of the concept. The NASA Geoscience Laser Altimeter System (GLAS) onboard the ICESat satellite (2003)(2004)(2005)(2006)(2007)(2008)(2009) (Table 1) was the first spaceborne lidar designed for polar science. The GLAS transmitted short pulses of near infrared (1064 nm) light for surface altimetry and visible green light (532 nm) for vertical distribution of clouds and aerosols [101]. The ICESat mission was marked by several important innovations in polar altimetry, including unprecedented spatial resolution (70 m footprint, 172 m along-track spacing, and 3 cm vertical resolution), geographic coverage to ±86 • , and the use of a narrow-beam 1064 nm laser which reduced sensitivity to both slope errors and signal penetration into snow and firn [163].
An important design objective of ICESat was to measure ice sheet surface elevation change over regions of high surface slope and complex topography, with primary mission objectives to measure spatially-averaged (104 km 2 ) ice sheet surface elevation to <15 cm absolute accuracy and elevation change to <1.5 cm yr −1 accuracy [36,101]. Field validation suggests operational absolute accuracy achieved 2 ± 3 cm for optimal conditions [164]. At-a-point measurement precision is ± 3 cm over the ice sheets and the repeat crossover measurement uncertainty is typically 10-15 cm [163,165], and up to 59 cm in regions of high surface slope [39]. Performance is affected by forward scattering of the return signal by intervening clouds, detector saturation, uncertain laser pointing angle, and laser transmit power decline [39,164,166].
The ICESat payload included three separate laser systems. Technical issues with the lasers prevented planned continuous operation, and instead, data collection took place during 18 separate campaigns [166,167]. Campaigns were designated by the laser system used and the campaign-specific orbital repeat frequency. For example, the first laser system operated for 38 days in an 8-day calibration and validation mode repeat frequency before it failed [167]. The second and third laser systems operated in campaign mode with a 91-day repeat frequency constituting individual campaigns~33 days in length separated by 4-6 months. Inter-campaign range biases up to~20 cm have been found [165,168]. The biases owe in part to a coding error in the ICESat signal processing algorithm that has since been corrected [168]. On-orbit and post-processing bias corrections were developed to compensate for orbital drift and systematic laser orientation (pointing) errors, but residual inter-campaign biases are an important measurement uncertainty in the ICESat data [86,165,169].
As with all spaceborne altimeters, the ICESat orbit followed ascending and descending orbital reference tracks with crossover reference points at the ascending and descending orbit intersections ( Figure 2). Several factors affect the geolocation of ICESat ground footprints at repeat-track and crossover locations. These include orbital vibrations, orbital drift, random errors in ICESat's laser orientation determination system, and declines in laser transmit power through time that modify the illuminated footprint diameter [164,167,170]. These factors produce both real and apparent deviations between illuminated ground tracks and reference tracks up to a few hundred meters [163,171]. Consequently, exact co-located repeat track and crossover point measurements are extremely rare, and spatial interpolation is required to recover dH/dt at-a-point, which has become a central challenge for ICESat data interpretation [171]. To estimate pan-GrIS dH/dt, an additional spatial interpolation of the point dH/dt estimates is typically performed, which is complicated by temporal differences among the point dH/dt estimates e.g., [86].
The XO method calculates dH/dt at the ascending and descending orbital crossover points. The method was originally developed for radar altimetry and was expected to be the primary method for interpreting ICESat data [101]. The method was used to estimate a relative elevation accuracy for ICESat of 0.14-0.59 m depending on the surface slope [39]. Using ICESat XO points as reference values, absolute accuracies for ERS-1/2 and EnviSat radar retrievals are~0.10-0.56 m for surface slopes <0.1 • , but are 2.27 m on average and up to 30 m for surface slopes exceeding 0.7 • [39]. These differences demonstrate the utility of laser altimetry for validating radar altimetry but also suggest caution when combining retrievals at discrete locations with high surface slopes and/or topographic variability c.f. [109]. As detailed by Felikson et al. [99], four methods have been developed to recover / from ICESat data: (1) repeat tracks (RTs), (2) crossovers (XOs), (3) overlapping footprints (OFPs), and (4) surface fitting methods, such as triangulated irregular networks (TINs). The Surface Elevation Reconstruction and Change (SERAC) method [171] combines elements of RT and XO [99]. Alberti and Biscaro [172] summarize the ICESat orbital parameters in detail and describe a flexible algorithm for determining RT and XO points from them.
The XO method calculates / at the ascending and descending orbital crossover points. The method was originally developed for radar altimetry and was expected to be the primary method for interpreting ICESat data [101]. The method was used to estimate a relative elevation accuracy for ICESat of 0.14-0.59 m depending on the surface slope [39]. Using ICESat XO points as reference values, absolute accuracies for ERS-1/2 and EnviSat radar retrievals are ~0.10-0.56 m for surface slopes <0.1°, but are 2.27 m on average and up to 30 m for surface slopes exceeding 0.7° [39]. These differences demonstrate the utility of laser altimetry for validating radar altimetry but also suggest caution when combining retrievals at discrete locations with high surface slopes and/or topographic variability c.f. [109].
The RT method calculates / from repeat measurements at reference points along the reference tracks, which increases spatial measurement density and coverage relative to XO points ( Figure 3) [173]. However, RT points are typically further apart than XO points, which may increase slope errors [99]. Owing to higher measurement density, the RT method has seen greater operational use than the XO method and is more appropriate for local or regional scale analysis [99]. The RT method was combined with Advanced Spaceborne Thermal Emission and Reflection Radiometer The RT method calculates dH/dt from repeat measurements at reference points along the reference tracks, which increases spatial measurement density and coverage relative to XO points ( Figure 3) [173]. However, RT points are typically further apart than XO points, which may increase slope errors [99]. Owing to higher measurement density, the RT method has seen greater operational use than the XO method and is more appropriate for local or regional scale analysis [99]. The RT method was combined with Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) digital elevation model differencing, which together suggest an average volume loss in southeast Greenland of~108 km 3 yr −1 for the period 2003-2005 [161]. Concentrated thinning of narrow outlet glaciers contributed less volume loss than dispersed, but smaller thinning over larger interior areas, highlighting the enhanced capability of laser altimetry for observing local-scale processes in the ablation zone [161].
The OFP method uses both XO and RT points but requires that the laser footprints of repeat measurements overlap. ICESat footprints are ellipses with a 50-90 m major axis length [174]. The overlapping criterion is determined by a predefined semi-major axis distance between laser footprints [100]. In principle, the method does not account for the slope between footprints, but slope corrections have been applied, for example, to validate ICESat with NASA Airborne Topographic Mapper (ATM) laser altimeter data [175]. Using a similar method, an average 0.07 m offset was found between NASA airborne LVIS lidar and ICESat elevation retrievals over the GrIS [176]. The OFP method has greater spatial coverage owing to the combined use of RT and XO points. Consequently, dH/dt uncertainty appears to be the lowest among methods [99]. The method was used to estimate an average thickening of 0.02 m yr −1 (equivalent to 21 km 3 yr −1 volumetric change) for areas above 2000 m a.s.l., and an average thinning of −0.24 m yr −1 (168 km 3 yr −1 ) below 2000 m a.s.l. for the period 2003-2007 [100]. Thinning rates near the margin were found to be slightly larger than estimates from other methods due to greater sample size and coverage in the ablation zone.
Remote Sens. 2019, 11, x FOR PEER REVIEW 12 of 52 (ASTER) digital elevation model differencing, which together suggest an average volume loss in southeast Greenland of ~108 km 3 yr −1 for the period 2003-2005 [161]. Concentrated thinning of narrow outlet glaciers contributed less volume loss than dispersed, but smaller thinning over larger interior areas, highlighting the enhanced capability of laser altimetry for observing local-scale processes in the ablation zone [161].

Figure 3.
Illustration of altimeter (e.g., ICESat) reference track and repeat tracks. Repeat tracks are parallel to reference tracks but separated in the cross-track direction. An area of interest in the alongtrack direction is used to collect ground footprint measurements that are spatially interpolated to recover the area-average surface elevation and subsequent change in area-average surface elevation / . In practice, repeat tracks are often not parallel owing to orbital variations, and in the case of ICESat laser altimetry, footprint diameter may change owing to laser transmit power variation.
The OFP method uses both XO and RT points but requires that the laser footprints of repeat measurements overlap. ICESat footprints are ellipses with a 50-90 m major axis length [174]. The overlapping criterion is determined by a predefined semi-major axis distance between laser footprints [100]. In principle, the method does not account for the slope between footprints, but slope corrections have been applied, for example, to validate ICESat with NASA Airborne Topographic Mapper (ATM) laser altimeter data [175]. Using a similar method, an average 0.07 m offset was found between NASA airborne LVIS lidar and ICESat elevation retrievals over the GrIS [176]. The OFP method has greater spatial coverage owing to the combined use of RT and XO points. Consequently, / uncertainty appears to be the lowest among methods [99]. The method was used to estimate an average thickening of 0.02 m yr −1 (equivalent to 21 km 3 yr −1 volumetric change) for areas above 2000 m a.s.l., and an average thinning of −0.24 m yr −1 (168 km 3 yr −1 ) below 2000 m a.s.l. for the period 2003-2007 [100]. Thinning rates near the margin were found to be slightly larger than estimates from other methods due to greater sample size and coverage in the ablation zone.
The TIN method [103] is similar to RT as it aggregates all available along-track measurements within a predefined distance around each RT reference point. Unlike the standard RT method, a TIN surface is fitted to a subset of points at each RT reference point collected within a reference two-year time period.
/ is estimated from all combinations of surface elevation measurements outside the reference period relative to the TIN surface. The TIN surface implicitly accounts for cross-track and along-track slope and sacrifices temporal resolution for both spatial resolution and reference period consistency. Using the TIN method, ICESat / uncertainty was estimated to be ±0.07 m yr −1 at 1σ level [103]. However, the method may systematically under-sample some areas of the ice sheet margin with insufficient points to fit the TIN, leading to an underestimated total ice sheet volume change [99]. Qualitatively, each method shows a similar pattern of increasing ice thinning Figure 3. Illustration of altimeter (e.g., ICESat) reference track and repeat tracks. Repeat tracks are parallel to reference tracks but separated in the cross-track direction. An area of interest in the along-track direction is used to collect ground footprint measurements that are spatially interpolated to recover the area-average surface elevation and subsequent change in area-average surface elevation dH/dt. In practice, repeat tracks are often not parallel owing to orbital variations, and in the case of ICESat laser altimetry, footprint diameter may change owing to laser transmit power variation.
The TIN method [103] is similar to RT as it aggregates all available along-track measurements within a predefined distance around each RT reference point. Unlike the standard RT method, a TIN surface is fitted to a subset of points at each RT reference point collected within a reference two-year time period. dH/dt is estimated from all combinations of surface elevation measurements outside the reference period relative to the TIN surface. The TIN surface implicitly accounts for cross-track and along-track slope and sacrifices temporal resolution for both spatial resolution and reference period consistency. Using the TIN method, ICESat dH/dt uncertainty was estimated to be ±0.07 m yr −1 at 1σ level [103]. However, the method may systematically under-sample some areas of the ice sheet margin with insufficient points to fit the TIN, leading to an underestimated total ice sheet volume change [99]. Qualitatively, each method shows a similar pattern of increasing ice thinning along the ice sheet margins with the highest thinning rates along the southern coasts, especially in the marine-terminating Jakobshavn Isbrae region and in southeast Greenland. Across methods, pan-GrIS volume loss ranges from~200-275 km 3 yr −1 during the ICESat era [99].

ICESat-2 and Future Opportunities
The Advanced Topographic Laser Altimeter System (ATLAS) was launched in September 2018 onboard the Ice, Cloud, and land Elevation Satellite-2 (ICESat-2) ( Table 1). ATLAS is a dual-beam single-photon counting laser altimeter marked by several improvements to the ICESat/GLAS measurement strategy [40,177]. These include its dual-beam laser, organized into three pairs with a 90 m beam separation and 3.3 km pair separation. Each beam has a nominal 17 m spatial footprint diameter and 0.7 m along-track measurement spacing ( Figure 4). This configuration permits the determination of cross-track slope, increases spatial measurement density, and increases local measurement accuracy [177]. ICESat-2 mission requirements include the determination of ice sheet dH/dt to 0.4 cm yr −1 accuracy (0.25 cm yr −1 for outlet glaciers) when averaged over 100 km 2 areas [177]. Data collection began on Oct 14, 2018, with initial ICESat-2 data products released on 7 June 2019, including an initial release of the ATLAS/ICESat-2 L3A Land Ice Height (ATL06), Version 1 [178], which provides mean ice sheet surface elevation averaged along 40 m segments of ground track posted at 20 m along-track spacing ( Figure 5) (https://nsidc.org/data/ATL06/). showing the ICESat-2 multi-beam configuration (two beam-pairs shown, the third beam-pair was not included in ATL06 Version 1 at this time and location). Single-beam altimeters such as ICESat give one crossover measurement at each crossover location, and the cross-track slope cannot be determined in the along-track direction. The ICESat-2 multi-beam configuration allows determination of cross-track slope and gives multiple unique crossover measurements at each crossover location (up to nine possible, four shown in this example).
In addition to an improved multi-beam configuration, the ATLAS instrument will transmit laser pulses at 532 nm wavelength, which will enhance both photon-return density and optical penetration capability relative to the GLAS instrument on ICESat. Water, both in liquid and solid form, is nearly transparent to 532 nm light. For example, the path length in pure water required to attenuate 532 nm light to 50% of its incident intensity is ~16 m, whereas this path length is ~0.01 m for 1064 nm light, owing to an ~1400x respective increase in the imaginary index of refraction of water [179,180]. Consequently, ICESat-2 laser pulses at 532 nm will effectively "see through" liquid water and thereby facilitate spaceborne measurement of supraglacial lake, stream, and water-filled crevasse bathymetry [181]. Conversely, the reflectivity of glacier ice and snow is near maximum at 532 nm owing to the extremely low absorption coefficient of ice and extremely high scattering efficiency from ice grains and air bubbles [182]. This enhanced reflectivity should produce higher photon-return density relative to the 1064 nm laser but will also include a non-zero volume scattering component that could introduce a range bias similar to radar penetration into snow and firn [183,184]. (b) example crossover location (red box inset in (a)) showing the ICESat-2 multi-beam configuration (two beam-pairs shown, the third beam-pair was not included in ATL06 Version 1 at this time and location). Single-beam altimeters such as ICESat give one crossover measurement at each crossover location, and the cross-track slope cannot be determined in the along-track direction. The ICESat-2 multi-beam configuration allows determination of cross-track slope and gives multiple unique crossover measurements at each crossover location (up to nine possible, four shown in this example).
In addition to an improved multi-beam configuration, the ATLAS instrument will transmit laser pulses at 532 nm wavelength, which will enhance both photon-return density and optical penetration capability relative to the GLAS instrument on ICESat. Water, both in liquid and solid form, is nearly transparent to 532 nm light. For example, the path length in pure water required to attenuate 532 nm light to 50% of its incident intensity is~16 m, whereas this path length is~0.01 m for 1064 nm light, owing to an~1400x respective increase in the imaginary index of refraction of water [179,180]. Consequently, ICESat-2 laser pulses at 532 nm will effectively "see through" liquid water and thereby facilitate spaceborne measurement of supraglacial lake, stream, and water-filled crevasse bathymetry [181]. Conversely, the reflectivity of glacier ice and snow is near maximum at 532 nm owing to the extremely low absorption coefficient of ice and extremely high scattering efficiency from ice grains and air bubbles [182]. This enhanced reflectivity should produce higher photon-return density relative to the 1064 nm laser but will also include a non-zero volume scattering component that could introduce a range bias similar to radar penetration into snow and firn [183,184]. ICESat-2 research to date has focused on pre-mission proof of concept, airborne and ground validation, and sensor calibration [186][187][188]. In addition to its unique ability to map surface water bathymetry, other novel applications of laser altimetry to the ablation zone that will benefit from ICESat-2 continuity include assimilation of / observations into transient ice flow simulations [189], fusion of multi-sensor (e.g., stereophotogrammetry or radar altimetry) datasets to increase the spatial density of surface elevation measurements [174], and application of lidar-based snow grain size and surface roughness retrievals independent of or in combination with optical sensors [190,191]. Comparison between ice sheet surface elevation change estimates by the laser altimetry and regional climate model (RCM) SMB is another area that is underexplored [89]. Such comparisons require accurate knowledge of the surface density and change in surface elevation due to ice dynamic motion, and, therefore, are most appropriate for bare ice areas during summer when surface melting dominates the elevation change signal [89]. Finally, it was recently reported that the CryoSat-2 science team is considering adjusting the satellite's orbit to overlap its ground tracks with those of ICESat-2 once every ~2 days [21]. If so, these overlapping data would provide cross-platform / validation and could be used to infer changes in the CryoSat-2 Ku-band penetration depth into snow and firn. ICESat-2 research to date has focused on pre-mission proof of concept, airborne and ground validation, and sensor calibration [186][187][188]. In addition to its unique ability to map surface water bathymetry, other novel applications of laser altimetry to the ablation zone that will benefit from ICESat-2 continuity include assimilation of dH/dt observations into transient ice flow simulations [189], fusion of multi-sensor (e.g., stereophotogrammetry or radar altimetry) datasets to increase the spatial density of surface elevation measurements [174], and application of lidar-based snow grain size and surface roughness retrievals independent of or in combination with optical sensors [190,191]. Comparison between ice sheet surface elevation change estimates by the laser altimetry and regional climate model (RCM) SMB is another area that is underexplored [89]. Such comparisons require accurate knowledge of the surface density and change in surface elevation due to ice dynamic motion, and, therefore, are most appropriate for bare ice areas during summer when surface melting dominates the elevation change signal [89]. Finally, it was recently reported that the CryoSat-2 science team is considering adjusting the satellite's orbit to overlap its ground tracks with those of ICESat-2 once everỹ 2 days [21]. If so, these overlapping data would provide cross-platform dH/dt validation and could be used to infer changes in the CryoSat-2 Ku-band penetration depth into snow and firn.

Remote Sensing of Mass Balance
Three methods are used to estimate MB from satellite remote sensing: (1) changes in ice sheet surface elevation from altimetry, (2) the mass budget or input-output method, and (3) time-variable gravimetry [192]. Satellite and airborne altimeters measure dH/dt, which is converted to MB using modeled firn density and modeled SMB to partition D. The input-output method uses remotely-sensed ice sheet surface velocity and ice sheet thickness to calculate D for outlet glacier drainage basins and uses modeled SMB to calculate MB. Satellite gravimetry measures changes in the earth's gravity field, to directly estimate MB. The methods are complementary and largely independent, and together provide a comprehensive view of spatial and temporal patterns in the GrIS MB, since 2002 for gravimetry, 1992 for altimetry, and 1972 for IOM [192]. We review each method below and compare estimates of MB from each.

Converting Ice Surface Elevation Change to Mass Change
Satellite altimeters observe the total change in ice sheet surface elevation, which includes both real and apparent changes in ice sheet thickness dH. Glaciers and ice sheets thicken from snow accumulation and thin by melting, sublimation, horizontal ice flow, snow redistribution, and increases in ice, snow, and firn density. Apparent changes in ice sheet thickness are caused by the vertical motion of underlying bedrock. Converting altimetric measurements of dH/dt to change in ice volume and mass balance requires an estimate of each individual contribution to dH/dt, and a correction for each term in Equation (4) that is not associated with a change in ice volume and/or change in mass balance.
The typical procedure for estimating each term on the right-hand side of Equation involves obtaining SMB from a climate model and defining the accumulation and ablation zones as those areas where SMB is positive and negative, respectively, separated by the ELA where modeled SMB = 0 [86]. Above the ELA, positive dH/dt is caused by net accumulation and ρ s is estimated with a firn density model, whereas negative dH/dt is caused by ice dynamics and ρ s is solid ice density (typically taken as 917 kg m −3 ). Below the ELA, dH/dt is caused by both negative SMB and ice dynamics, and solid ice density is used to convert both to mass c.f. [86,97,105,193]. The v c term causes an apparent change in mass that is estimated with a firn compaction model and removed from dH/dt [194]. BMB is negligible in the accumulation zone and very small (~0.015-0.020 m yr −1 ) in the ablation zone, and v br contributes about 1 GT yr −1 in apparent mass change that can be corrected or ignored [86,195].
When corrected for the terms in Equation, the ICESat altimetry data suggest mass loss rates of -191 ± 23 Gt yr −1 to -240 ± 28 Gt yr −1 for the period 2003-2008 [86]. The spread arises from the choice of dH/dt interpolation method (Section 3.2.1) and, more substantially, the choice of firn density model. Removal of v c reduces mass loss rates by 33-36 Gt yr −1 , almost all of which is due to compaction in the zone between the ELA and 2000 m a.s.l. [86]. Despite the uncertainties associated with firn densification, a good agreement is found between altimetry and mass balance estimates from GRACE, which suggest a mass loss of -230 ± 33 Gt yr −1 during the period 2002-2009 [196], -179 ± 25 Gt yr −1 for 2002-2007 [197], and -237 ± 20 Gt yr −1 for 2003-2008 [17]. Applying the same methods as Sørensen et al. [86], Greenland's peripheral ice caps lost mass at -28 ± 11 Gt yr −1 (0.08 ± 0.03 mm yr −1 SLR equivalent), representing 20% of total GrIS mass loss during the ICESat era [96].
A combined analysis of ICESat data and airborne laser altimeter data from the NASA Airborne Thematic Mapper (ATM) suggests a mass loss rate of -243 ± 18 Gt yr −1 (0.68 ± 0.05 mm yr −1 SLR equivalent) for the 2003-2009 period [97]. The SERAC method was used to combine ICESat and ATM XO and RT points, and an SMB model was used to diagnose spatial differences in thinning, likely representing the most comprehensive analysis of spatiotemporal variability in dH/dt and its attribution to SMB vs. D for the ICESat era. Negative SMB accounted for 52% of total mass loss and accelerated during the study period. Persistent SMB decreases were concentrated in the southwest and southeast marginal ablation zones, but an acceleration in thinning due to SMB was also observed in the northwest. On average, dynamic thinning of outlet glaciers decelerated during the study period, but large spatiotemporal variability was found with some outlets accelerating, owing to local scale controls on ice dynamics (e.g., bed geometry). Regionally, the largest contribution to D was from southeast outlet glaciers, with the largest thinning rates observed on Helheim and Kangerdlugssuaq glaciers and on Sermeq Kujalleq (Jakobshavn Isbrae) in the west. As with prior studies, e.g., [198], a slight thickening of the interior ice sheet above the ELA was found.
As with altimetric surveys of dH/dt, the IOM reveals a complex spatial pattern of ice sheet mass loss concentrated in narrow outlet glaciers along the coastal margins, with longer term (e.g., 1990-present) mass losses concentrated in the southeast and the Jakobshavn basin in the west, and recent (post-2005) increases concentrated in the northwest [18,192,200,220,221]. In general, a small number of isolated glaciers drive the majority of D with just four (Sermeq Kujalleq, Kangerdlugssuaq, Køge Bugt, and Ikertivaq South) accounting for~50% of total D acceleration during the period 2000-2012 [18,192,222]. The IOM provides the longest record of D and, therefore, uniquely places recent MB trends in a long term context. For example, although negative SMB accounts for approximately 60% of MB for the period 1990-present, the proportion is reversed for the period 1972-2018 [192,220,221].
The recent availability of digital bed elevation models, ice surface velocity mosaics, and digital surface elevation models for the entire ice sheet provide new insight into the spatial and temporal drivers of GrIS MB at the scale of individual outlet glaciers and their upstream drainage basins [185,206,212,218]. Mouginot et al. [192] reconstruct 46 years of MB for 260 individual outlet glaciers for the period 1972-2018, using surface velocity constructed from all available SAR (ALOS/PALSAR, ENVISAT/ASAR, RADARSAT-1/2, and Sentinel-1b) and Landsat-8 optical imagery [212], ice thickness from BedMachine v3 [218], and surface topography from the Greenland Ice Mapping Project DEM [185], WorldView DEMs, and historical orthophotographs [223]. Their results show that MB switched from +47 ± 21 Gt yr −1 in 1972-1980 to -187 ± 17 Gt yr −1 in 2000-2008, and -286 ± 20 Gt yr −1 in 2010-2018, within the bounds of most estimates from altimetry and GRACE for the latter two periods. As in Enderlin et al. [18], they find four glaciers control half the mass loss during 2000-2012, but for the 1972-2018 study period, Ikertivaq South gained mass, while Steenstrup-Dietrichson in northwest, Humboldt in north, and Midgårdgletscher in southeast join the three others as dominant drivers of mass loss. The largest single historical contributor is Sermeq Kujalleq (Jakobshavn Isbrae), which has experienced long term retreat and episodic rapid acceleration [224], including the disintegration of its floating ice tongue between 2000-2003 and rapid retreat of its calving front ( Figure 6) [224][225][226][227][228][229]. Whereas D has historically been dominated by the southeast, northwest, and western sectors, they conclude that the north and northeast sectors are of greatest importance to future sea level rise owing to their present-day slow velocities and potential for large increases in D. The annual average ice surface velocity for the period 2011-2016 is shown for upstream areas, calculated from Landsat-8 (optical), Sentinel-1, and RADARSAT-2 (interferometric SAR) image processing [212]. The maximum surface velocity value is ~12 km yr −1 . Calving fronts, which mark the terminus position of the outlet glacier where ice is discharged to the ocean, are mapped from Landsat-7, Landat-8, and Sentinel-2B optical imagery [229]. Sermeq Kujalleq has lost 137 km 2 of surface area between 1998-2018 (compare thick red line to thick black line), resulting in near-complete disintegration of its floating ice tongue [224,227]. Its flow speed has nearly doubled since the early 1990s, with upstream thinning rates exceeding 15 m yr −1 [226]. Inset is the Greenland Ice Mapping Project (GIMP) ice mask [185] with the red box showing the extent of the image area.

Time Variable Gravimetry and the Twin-GRACE Mission
The NASA/German Aerospace Center twin Gravity Recovery and Climate Experiment (GRACE) satellites (2002-2017) are conceptually unique from all other remote sensing platforms because they do not measure the interaction between the earth and electromagnetic radiation [230]. As the twin GRACE satellites orbit, their proximal distance varies with earth's gravitational pull. These small changes in their acceleration are used to measure variations in the density of earth at ~400-500 km spatial resolution. Consequently, GRACE provides the only direct measurement of ice sheet and independent validation of estimates made from satellite altimetry and RCMs (Figure 7) [231][232][233]. The annual average ice surface velocity for the period 2011-2016 is shown for upstream areas, calculated from Landsat-8 (optical), Sentinel-1, and RADARSAT-2 (interferometric SAR) image processing [212]. The maximum surface velocity value is~12 km yr −1 . Calving fronts, which mark the terminus position of the outlet glacier where ice is discharged to the ocean, are mapped from Landsat-7, Landat-8, and Sentinel-2B optical imagery [229]. Sermeq Kujalleq has lost 137 km 2 of surface area between 1998-2018 (compare thick red line to thick black line), resulting in near-complete disintegration of its floating ice tongue [224,227]. Its flow speed has nearly doubled since the early 1990s, with upstream thinning rates exceeding 15 m yr −1 [226]. Inset is the Greenland Ice Mapping Project (GIMP) ice mask [185] with the red box showing the extent of the image area.

Time Variable Gravimetry and the Twin-GRACE Mission
The NASA/German Aerospace Center twin Gravity Recovery and Climate Experiment (GRACE) satellites (2002-2017) are conceptually unique from all other remote sensing platforms because they do not measure the interaction between the earth and electromagnetic radiation [230]. As the twin GRACE satellites orbit, their proximal distance varies with earth's gravitational pull. These small changes in their acceleration are used to measure variations in the density of earth at~400-500 km spatial resolution. Consequently, GRACE provides the only direct measurement of ice sheet MB and independent validation of MB estimates made from satellite altimetry and RCMs (Figure 7) [231][232][233].  [234]. The GRACE mass balance product is produced by the Danish Institute of Space (DTU Space) and is provided as a five-year average for the period 2012-2016 with a 500 km nominal spatial resolution [235]. Both datasets are Essential Climate Variables provided by the ESA Climate Change Initiative and are available online at http://products.esa-icesheets-cci.org.
As an independent dataset, many studies have combined or compared GRACE with satellite altimetry or other traditional remote sensing techniques [236]. GRACE has been compared to ERS radar altimetry [237], to surface melt extent from MODIS thermal imagery [238], and to from combined RCM and InSAR derived [239]. These comparisons confirmed earlier findings from both satellite altimetry and RCMs that the GrIS declined at an accelerating rate during the post-2002 period [196,237,240]. GRACE estimates range from −191 ± 21 Gt yr −1 between 2002-2009 [98] to −244 ± 6 Gt yr −1 between 2002-2016, with an acceleration of −28 ± 9 Gt yr −2 during this period [241]. An updated analysis of GRACE data suggests the negative acceleration halted in 2013 owing to the suppression of surface runoff production in west Greenland caused by a shift in the North Atlantic Oscillation to its cool phase [14,93].
The accuracy of GRACE data is inherently limited by its spatial resolution, with accuracy decreasing as spatial resolution increases [197,232,242]. Sophisticated signal processing methods are required to determine the optimal spatial resolution and to quantify the spatial dependence of the signal processing error [240,243,244]. For example, north-south oriented spatially correlated errors ("striping") are a persistent issue [245]. Most early methods used the Mascon approach to extract the  [234]. The GRACE mass balance product is produced by the Danish Institute of Space (DTU Space) and is provided as a five-year average for the period 2012-2016 with a 500 km nominal spatial resolution [235]. Both datasets are Essential Climate Variables provided by the ESA Climate Change Initiative and are available online at http://products.esa-icesheets-cci.org.
As an independent dataset, many studies have combined or compared GRACE MB with satellite altimetry or other traditional remote sensing techniques [236]. GRACE MB has been compared to ERS radar altimetry [237], to surface melt extent from MODIS thermal imagery [238], and to MB from combined RCM SMB and InSAR derived D [239]. These comparisons confirmed earlier findings from both satellite altimetry and RCMs that the GrIS MB declined at an accelerating rate during the post-2002 period [196,237,240]. GRACE MB estimates range from −191 ± 21 Gt yr −1 between 2002-2009 [98] to −244 ± 6 Gt yr −1 between 2002-2016, with an acceleration of −28 ± 9 Gt yr −2 during this period [241]. An updated analysis of GRACE data suggests the negative MB acceleration halted in 2013 owing to the suppression of surface runoff production in west Greenland caused by a shift in the North Atlantic Oscillation to its cool phase [14,93].
The accuracy of GRACE data is inherently limited by its spatial resolution, with accuracy decreasing as spatial resolution increases [197,232,242]. Sophisticated signal processing methods are required to determine the optimal spatial resolution and to quantify the spatial dependence of the signal processing error [240,243,244]. For example, north-south oriented spatially correlated errors ("striping") are a persistent issue [245]. Most early methods used the Mascon approach to extract the GRACE signal onto a geometric grid [233,246]. Recent work has shown spherical Slepian functions theoretically maximize the spatial resolution of the GRACE signal [241]. The Slepian solutions were used to extract MB trends at a previously unresolvable scale, showing nearby Baffin Island and Ellesmere Island MB trends of −22 ± 2 and −38 ± 2 Gt yr −1 , respectively [241]. In addition to signal processing errors, the glacial isostatic adjustment effect on the GRACE signal is~5-20 Gt yr −1 for the GrIS, but is uncertain [235,247].
The GRACE Follow-On (GRACE-FO) mission (launched 22 May 2018) is the successor to the GRACE mission [248]. GRACE-FO uses the same orbital acceleration measurement technology as the original GRACE mission and will continue the GRACE measurement time series for a planned ten years. In addition to providing gravity time series continuity, GRACE-FO will test a laser ranging interferometer that is expected to improve the satellite-to-satellite distance measurement relative to the GRACE microwave ranging system [42,43]. Following a period of in-orbit checks, GRACE-FO entered the science phase of its mission in January 2019 [249]. The first GRACE-FO Level-1 data products were released on 24 May 2019, and Level-2 gravity field products were released on 10 June 2019, with planned monthly release updates thereafter [250].

Remote Sensing of Ice Surface Reflectance and Albedo
Surface albedo modulates the absorption of shortwave radiation by ice and snow surfaces [75]. Shortwave radiation is the dominant contributor to melt energy in the GrIS ablation zone [90]. Consequently, albedo is an important control on the production of surface meltwater and, by extension, the SMB [93]. Satellite remote sensing instruments measure multispectral reflectance, which is used to estimate albedo [251]. Remotely sensed albedo is used to understand spatial patterns in surface melting, as input to land surface models and regional climate models that simulate the SMB, and as a diagnostic marker of the changing ice surface [26,30,[252][253][254]. The following subsections summarize key terminology relevant to spaceborne measurements of surface reflectance and albedo, the sensors and datasets used to retrieve the GrIS surface albedo, observed changes in GrIS albedo and their diagnosis, and current challenges and future opportunities for understanding the GrIS ablation zone albedo.

Definition of Reflectance, BRDF, and Albedo
Standardized nomenclature for reflectance quantities were defined in terms of incident and reflected beam geometry by Nicodemus et al. [255] and later adapted for common optical remote sensing measurement configurations [256][257][258][259]. Here we briefly review key definitions of reflectance quantities, following Schaepman-Strub et al. [259].  [259]. For passive optical remote sensing measurements, the incident beam geometry is hemispherical and is composed of both direct-beam and diffuse solar radiation (Figure 8). The exitent (reflected) beam geometry is conical, defined by the sensor instantaneous field of view (IFOV), and is composed of both direct-beam and diffuse reflected solar radiation, corresponding to the hemispherical-conical reflectance function (HCRF) described by Schaepman-Strub et al. [259].
The bidirectional reflectance distribution function BRDF [sr −1 ] describes the scattering of an infinitesimal beam of incident light from one direction in a hemisphere into another direction, and is where θ i and θ r are the incident and reflected viewing angle, respectively, and φ i and φ r are the incident and reflected azimuth angles, respectively. As a ratio of infinitesimal quantities, the BRDF is a theoretical construct that cannot be directly measured but describes the intrinsic reflectance properties of a surface from which measurable quantities can be derived [259]. For example, integration of the BRDF across θ r (0 → π/2) and φ r (0 → 2π) yields bihemispherical reflectance, or what is commonly called albedo.
Remote Sens. 2019, 11, x FOR PEER REVIEW 21 of 52 a theoretical construct that cannot be directly measured but describes the intrinsic reflectance properties of a surface from which measurable quantities can be derived [259]. For example, integration of the BRDF across r (0 → π/2) and r (0 → 2 ) yields bihemispherical reflectance, or what is commonly called albedo. Figure 8. Diagram of the incident and reflected viewing angles for a directional light source with direct and diffuse (hemispherical) incoming radiance and conical reflected radiance, corresponding to the hemispherical-conical reflectance factor (HCRF) described by Schaepman-Strub et al. [259]. The incident viewing angle ( i ), reflected viewing angle ( r ), and the azimuth angle ( ) together define the angular coordinates of both the directional illumination and the light reflected toward the sensor, with respect to the surface normal, z.
Viewing geometry is important because every natural earth surface is an anisotropic reflector. For example, snow and ice are strongly forward scattering at the particle scale, especially at visible wavelengths, but ice surfaces can be strongly backward scattering due to surface roughness [260]. Consequently, an estimate of the BRDF is required to convert multi-angular reflectance measured by satellite remote sensing instruments to albedo. To achieve this, the BRDF for each satellite ground footprint is estimated (not measured) from repeat multi-angular reflectance measurements collected over time and, therefore, under varying illumination conditions [251] or collected instantaneously by spaceborne multi-angular instruments [261]. A complete technical discussion of the BRDF is available from Schaepman-Strub et al. [259].

Optical Reflectance and Albedo Sensors and Datasets
Spatially-and temporally-consistent repeat measurements of surface reflectance, and estimates of BRDF and albedo for the Greenland Ice Sheet are provided by satellite-borne imaging spectrometers, including the NASA Moderate Resolution Imaging Spectroradiometer (MODIS), the Figure 8. Diagram of the incident and reflected viewing angles for a directional light source with direct and diffuse (hemispherical) incoming radiance and conical reflected radiance, corresponding to the hemispherical-conical reflectance factor (HCRF) described by Schaepman-Strub et al. [259]. The incident viewing angle (θ i ), reflected viewing angle (θ r ), and the azimuth angle (φ r ) together define the angular coordinates of both the directional illumination and the light reflected toward the sensor, with respect to the surface normal, z.
Viewing geometry is important because every natural earth surface is an anisotropic reflector. For example, snow and ice are strongly forward scattering at the particle scale, especially at visible wavelengths, but ice surfaces can be strongly backward scattering due to surface roughness [260]. Consequently, an estimate of the BRDF is required to convert multi-angular reflectance measured by satellite remote sensing instruments to albedo. To achieve this, the BRDF for each satellite ground footprint is estimated (not measured) from repeat multi-angular reflectance measurements collected over time and, therefore, under varying illumination conditions [251] or collected instantaneously by spaceborne multi-angular instruments [261]. A complete technical discussion of the BRDF is available from Schaepman-Strub et al. [259].

Optical Reflectance and Albedo Sensors and Datasets
Spatially-and temporally-consistent repeat measurements of surface reflectance, and estimates of BRDF and albedo for the Greenland Ice Sheet are provided by satellite-borne imaging spectrometers, including the NASA Moderate Resolution Imaging Spectroradiometer (MODIS), the NASA Multi-angle Imaging SpectroRadiometer (MISR), and the National Oceanic and Atmospheric Administration's Advanced Very High Resolution Radiometers (AVHRR) ( Table 2) [261][262][263][264]. Other notable spaceborne sensors providing surface reflectance that have been used to study the GrIS ablation zone include the Enhanced Thematic Mapper Plus (ETM+) onboard Landsat 7, the Operational Land Imager (OLI) onboard Landsat 8, the ASTER onboard Terra, the High Resolution imagers onboard the Satellite Pour l'Observation de la Terre (SPOT 1-7), the Global Land Imager (GLI) onboard the Advanced Earth Observing Satellite II (ADEOS-2), the Multispectral Imager (MSI) onboard Sentinel-2, and the Ocean Colour and Land Instrument (OLCI) onboard Sentinel-3A/B (Table 2) [265][266][267][268][269][270][271]. Table 2. Summary of optical and near-infrared remote sensing platforms, instruments, temporal coverage, observed wavelengths, and managing agencies.

Platform * Instrument
GrIS but are estimated to exceed >0.50 m in some cases [116][117][118]. Variations in surface slope dominate radar altimeter accuracy in the topographically-complex bare-ice ablation zone, with slopeinduced errors that exceed >20 m for slopes >1° [39,109,113]. The recent CryoSat-2 radar altimeter mission employs synthetic aperture radar (SAR) and interferometric synthetic aperture radar (InSAR) technologies which decrease the effective along-track footprint diameter to the order 0.1-0.3 km and provide cross-track slope from InSAR processing [93]. Together with denser orbital-track spacing and advances in waveform retracking, topographic mapping of ice sheet ablation zone areas is accurate to within a few meters on average, and / to within sub-meter uncertainty relative to laser altimetry [93,119,120].

Radar Altimetry Sensors and Datasets
The first satellite altimeter measurements of surface topography for the Greenland Ice Sheet were made by the GEOS-3, Geosat, and Seasat oceanographic Ku-band altimeters (Error! Reference source not found.) [37,91,106,121,122]. Their geographical coverage was limited to ±72 o but the data they collected showed it was possible to calculate volumetric changes in the polar ice sheets from space, providing proof of concept for future polar-orbiting altimeters [92]. A first demonstration using GEOS-3, Geosat, and Seasat data found that for areas south of <72 o N, the GrIS thickened by 23 ± 6 cm yr -1 between 1978-1986, with enhanced snowfall in a warmer climate hypothesized to explain the observed thickening [123]. These early missions supported important developments in waveform retracking, slope correction, and physical and empirical models for subsurface volume scattering [106,108,111,124]. These corrections are critical for accurate elevation retrieval and reliable change detection [125,126]. For example, Davis et al. [127] used an improved geodetic model and an improved retracking model to reinterpret the earlier findings of Zwally et al. [92], finding a smaller 1.5 ± 0.5 cm yr -1 thickening for the period 1978-1986 that highlighted the importance of slope-and signal penetration-correction.

Temporal Coverage
Observed Wavelength Agency dominate radar altimeter accuracy in the topographically-complex bare-ice ablation zone, with slopeinduced errors that exceed >20 m for slopes >1° [39,109,113]. The recent CryoSat-2 radar altimeter mission employs synthetic aperture radar (SAR) and interferometric synthetic aperture radar (InSAR) technologies which decrease the effective along-track footprint diameter to the order 0.1-0.3 km and provide cross-track slope from InSAR processing [93]. Together with denser orbital-track spacing and advances in waveform retracking, topographic mapping of ice sheet ablation zone areas is accurate to within a few meters on average, and / to within sub-meter uncertainty relative to laser altimetry [93,119,120].

Radar Altimetry Sensors and Datasets
The first satellite altimeter measurements of surface topography for the Greenland Ice Sheet were made by the GEOS-3, Geosat, and Seasat oceanographic Ku-band altimeters (Error! Reference source not found.) [37,91,106,121,122]. Their geographical coverage was limited to ±72 o but the data they collected showed it was possible to calculate volumetric changes in the polar ice sheets from space, providing proof of concept for future polar-orbiting altimeters [92]. A first demonstration using GEOS-3, Geosat, and Seasat data found that for areas south of <72 o N, the GrIS thickened by 23 ± 6 cm yr -1 between 1978-1986, with enhanced snowfall in a warmer climate hypothesized to explain the observed thickening [123]. These early missions supported important developments in waveform retracking, slope correction, and physical and empirical models for subsurface volume scattering [106,108,111,124]. These corrections are critical for accurate elevation retrieval and reliable change detection [125,126]. For example, Davis et al. [127] used an improved geodetic model and an improved retracking model to reinterpret the earlier findings of Zwally et al. [92], finding a smaller 1.5 ± 0.5 cm yr -1 thickening for the period 1978-1986 that highlighted the importance of slope-and signal penetration-correction. ion zone, regional climate models and firn compaction models are used to partition the / caused by change in from those caused by change in firn density and following subsections summarize spaceborne radar altimetry of the GrIS, providing an pective that emphasizes major challenges and advances in the field, and identifies nities for radar altimetry observations of the GrIS ablation zone. etry onal radar altimeters are nadir pointing single-beam radars that transmit and receive ands of microwave pulses per second, yielding ground measurement footprint he order 1-10 km posted at 0.3-0.6 km along-track spacing [39,106,107]. An estimate of ned by comparing measurements at crossover points, where an earlier satellite ground s a later one [92]. Two critical factors affect the accuracy of radar altimeter elevation r ice sheets: 1) variation in topographic slope, which controls the average distance timeter and the measured ground footprint surface [39,108,109], and 2) changes in the ittivity of the ice sheet surface, which controls the reflection, transmission, and the microwave signal [110][111][112][113].
l radar altimeter frequencies (Ku-band ~13 GHz), solid ice and liquid water are highly reas dry snow and firn are semi-transparent, with Ku-band penetration depths on the al meters [114,115]. Spatial and temporal variations in signal penetration depth caused are ice exposure, snow and firn densification, surface meltwater presence, and eltwater introduce spurious height change signals that are poorly quantified for the estimated to exceed >0.50 m in some cases [116][117][118]. Variations in surface slope r altimeter accuracy in the topographically-complex bare-ice ablation zone, with slopes that exceed >20 m for slopes >1° [39,109,113]. The recent CryoSat-2 radar altimeter ys synthetic aperture radar (SAR) and interferometric synthetic aperture radar (InSAR) hich decrease the effective along-track footprint diameter to the order 0.1-0.3 km and track slope from InSAR processing [93]. Together with denser orbital-track spacing and aveform retracking, topographic mapping of ice sheet ablation zone areas is accurate w meters on average, and / to within sub-meter uncertainty relative to laser 19,120].
ltimetry Sensors and Datasets satellite altimeter measurements of surface topography for the Greenland Ice Sheet the GEOS-3, Geosat, and Seasat oceanographic Ku-band altimeters (Error! Reference nd.) [37,91,106,121,122]. Their geographical coverage was limited to ±72 o but the data showed it was possible to calculate volumetric changes in the polar ice sheets from ng proof of concept for future polar-orbiting altimeters [92]. A first demonstration , Geosat, and Seasat data found that for areas south of <72 o N, the GrIS thickened by 23 een 1978-1986, with enhanced snowfall in a warmer climate hypothesized to explain hickening [123]. These early missions supported important developments in waveform pe correction, and physical and empirical models for subsurface volume scattering 4]. These corrections are critical for accurate elevation retrieval and reliable change ,126]. For example, Davis et al. [127] used an improved geodetic model and an acking model to reinterpret the earlier findings of Zwally et al. [92], finding a smaller r -1 thickening for the period 1978-1986 that highlighted the importance of slope-and tion-correction. ummary of radar and laser altimeter remote sensing platforms, instruments, temporal observed wavelength, and managing agencies. The ESA European Remote-sensing Satellite (ERS-1) carried the first polar-orbiting radar altimeter, extending coverage of the polar ice sheets to ±81.5 o [95]. Together with the follow-on ERS-2 and EnviSat missions, these satellites provided near complete GrIS coverage and 21 years of continuous data acquisition [38]. The first digital elevation model (DEM) spanning the entire GrIS surface was produced from combined ERS-1 and Geosat altimetry data, supplemented by stereophotogrammetric and cartographic data sets, with pan-GrIS average accuracy -0.33 ± 6.97 m [128,129]. ERS-1/2 altimetry data were used to discriminate thinning rates in the ablation zone (defined as <1500 m a.s.l.) of -2.0 ± 0.9 cm yr -1 from accumulation zone (>1500 m a.s.l.) thickening of 6 (Table A1) and may not reflect joint collaborations.
The ESA European Remote-sensing Satellite (ERS-1) carried the first polar-orbiting radar altimeter, extending coverage of the polar ice sheets to ±81.5 o [95]. Together with the follow-on ERS-2 and EnviSat missions, these satellites provided near complete GrIS coverage and 21 years of continuous data acquisition [38]. The first digital elevation model (DEM) spanning the entire GrIS surface was produced from combined ERS-1 and Geosat altimetry data, supplemented by stereophotogrammetric and cartographic data sets, with pan-GrIS average accuracy -0.33 ± 6.97 m [128,129]. ERS-1/2 altimetry data were used to discriminate thinning rates in the ablation zone (defined as <1500 m a.s.l.) of -2.0 ± 0.9 cm yr -1 from accumulation zone (>1500 m a.s.l.) thickening of 6.4 ± 0.2 cm yr -1 with an overall net thickening of 5.4 ± 0.2 cm yr -1 for the period 1992-2003 [130]. To extend the Managing agencies are identified by the WMO OSCAR database (Table A1)  The NOAA Extended AVHRR Polar Pathfinder (App-X) product provides surface broadband all-sky albedo and a suite of surface radiative flux and cloud property variables on a 25 km equal area grid twice-daily for the period 1982 to present for the Arctic and Antarctic regions [272]. The AVHRR App-X product provides a continuous record of GrIS surface albedo and an important climatological dataset [253]. An early demonstration of AVHRR surface albedo retrieval revealed a region of anomalously dark ice in the southwestern GrIS ablation zone later termed the "Dark Zone" (Section 5.3) [273]. Although rigorous quality control procedures are applied during dataset production [49], its utility for investigations of GrIS ablation zone albedo may be limited by geolocation errors, four-channel spectral resolution, and inter-mission biases [274,275]. Spatial and spectral resolution are both improved with MODIS, which measures surface reflectance in 30 narrow spectral bands at 250-500 m spatial resolution.
Higher spectral and/or spatial resolution reflectance products are provided by Landsat, MISR, SPOT, Sentinel, and other spaceborne spectrometers ( Table 2). Various albedo estimations exist for most sensors [276]. Many rely on the MODIS BRDF or empirical conversions from narrowband reflectance to broadband albedo, and for others (e.g., MISR), a unique BRDF and physically-based albedo estimation exists [261,277]. Most studies of the GrIS albedo have used MODIS albedo owing to the development of operational, publicly available data products spanning nearly two decades. Two MODIS albedo products are currently available: (1) the 500-m Terra/Aqua MODIS 8-day BRDF/Albedo product (MCD43) and, (2) the 500-m Terra daily snow/ice cover and BRDF/Albedo product (MOD10A) [251,263,278]. The MCD43 algorithm uses cloud-free, multi-date, atmospherically corrected input data from combined Terra (MOD09) and Aqua (MYD09) multi-angular surface reflectance retrievals to estimate a unique BRDF for each MODIS pixel in the seven MODIS 'Land' bands (1-7) every eight days. Albedo is calculated using this BRDF estimation for each of bands 1-7, and also for three broadbands (0.4-0.7, 0.7-3.0, and 0.4-3.0 um) [251,278].
The MOD10A algorithm is a daily product that uses the discrete-ordinate radiative transfer (DISORT) model to convert MOD09 daily surface reflectance to spectral albedo [279]. The daily frequency of the MOD10A product makes it particularly useful for studies of the GrIS ablation zone, where day-to-day albedo variability is substantial [26,280]. Both the MCD43 and MOD10A retrieval algorithms require cloud-free images, as identified by the MOD35 cloud mask. For the MCD43 product, data gaps within the 16-day period window may result in under-sampling of the BRDF, requiring the use of a backup BRDF algorithm. Albedo values estimated with the backup algorithm are generally considered less reliable, and quality flags are provided to distinguish between them [281].
The accuracy of MODIS albedo over the GrIS has been evaluated by comparison with surface-based measurements of albedo at Greenland Climate Network (GC-Net) automatic weather stations (AWS) [282]. The MCD43 product has an average root mean squared difference (RMSD) of ±0.07 and mean bias +0.022 relative to AWS albedo, which has ±0.035 RMSD uncertainty relative to precision pyranometer measurements [278,283]. The MCD43 RMSD is reduced to ±0.04 when only the highest-quality flagged values are used [283]. MODIS albedo accuracy is reduced at a high solar zenith angle (SZA) [263,283]. For example, Wang and Zender [284] found a large MCD43 albedo bias relative to GC-Net AWS albedo in the dry-snow zone for SZA > 55 • . However, these biases are mostly eliminated if only the highest-quality flagged values are used as recommended [281]. In general, the MODIS albedo products are not recommended for comparison with surface measurements for SZA > 70 • [281,283].
Relative to three AWS stations in the southwestern GrIS ablation zone, the MOD10A albedo is lower by up to 0.10 and by 0.06 on average during July for the period 2009-2016 [280]. The low bias is attributed to the inadequate sampling of surface heterogeneity by the AWS, which are preferentially located on flat areas of bare ice rather than areas of lower albedo, such as meltwater channels or crevasses that influence the satellite albedo [280]. Using one day of high quality in situ spectral albedo measurements (N = 232) collected within two MODIS satellite pixel ground footprints in the southwest GrIS ablation zone, Moustafa et al. [285] showed that the MODIS Collection V006 daily blue-sky albedo was accurate to within −0.04 to +0.07 for the homogeneous surface of one pixel and was biased low by −0.04 to −0.07 for the heterogeneous surface of the second pixel. As in Ryan et al. [280], the low bias was attributed to under sampling of dark meltwater channels and shadowed crevasses. In addition to spatial heterogeneity, comparisons between satellite albedo and surface measurements are complicated by factors including sensor spectral sensitivity and whether hemispherical or angular reflectance is measured [281,283]. Owing to these varied processes, ground validation of satellite albedo products for the GrIS ablation zone is challenging, and satellite retrieval uncertainty is enhanced relative to the interior snow-covered accumulation-zone surface, especially at sub-pixel scale [280,285,286].

Dark Ice in the Ablation Zone: Albedo Trends and Drivers
Satellite albedo retrievals provide an important marker of the changing ice surface. The MOD43 data suggest the GrIS area-averaged albedo for June, July, and August (JJA) decreased by −0.044 between 2000-2012 [278]. The JJA albedo anomaly during the record 2012 pan-GrIS melt event [287] was more than two standard deviations below the period mean. The most negative albedo trends are observed for the southern and western ice sheet ablation area where July albedo trends are −0.12 to −0.24 per decade, attributed to both reduced seasonal snow cover duration and increased bare ice exposure during this period [29,30,278]. The MOD10A product indicates a similar JJA albedo trend of −0.083 averaged over the GrIS ablation zone for the period 2000-2010 [30] and −0.078 for the period 2000-2013 [288]. Although MOD10A reflectance data suggest a modest brightening of the ablation zone since 2013 that likely reflects anomalously cold and snowier spring and summer conditions during this period [26,51], there is strong potential for melt-albedo feedbacks to accelerate negative SMB in the coming years [30,32].
For the ablation zone, spatial and temporal albedo variability is driven foremost by the seasonal snowline position and its control on bare ice exposure (Figure 9) [26]. Where snow is present, snow grain growth is the strongest modifier of surface albedo [252,253]. Bare ice albedo is scale dependent but is principally controlled by air bubble size and shape distribution, shadowing by cracks and crevasses [285,289], and by ice grain metamorphism, surface meltwater presence, and exposure and deposition of mineral and biological light-absorbing impurities (LAI) [51,52,58,253,[290][291][292]. The darkening effect of LAI is particularly apparent in satellite images of the GrIS ablation zone that reveal foliated bands of "dark ice" with anomalously low albedo relative to surrounding ice ( Figure 10) [52,273,293,294]. The wavy appearance of these foliated bands and their static position indicate their provenance as outcropping ice layers rich in dust deposited during past millennia [52,[294][295][296].
The spatial extent of bare ice and albedo variability within the bare ice zone exerts a primary control on ice sheet albedo and surface meltwater production in the ablation zone [26,30]. Shimada et al. [297] used MODIS/Terra surface reflectance to quantify the regional distribution of bare ice and dark ice extent during July for the period 2000-2014. The spatial extent of bare ice ranged from 5-16% of the GrIS surface area, dark ice ranged from 4-10%, and both exhibited a positive trend (4.4% yr −1 for bare ice and 7.6% yr −1 for dark ice) with the greatest expansion in the southwest ablation zone. Bare ice extent was strongly correlated (r = 0.66) with air temperature, whereas dark ice extent was weakly correlated with air temperature and was negatively correlated with solar radiation, suggesting that bare ice weathering by solar radiation may reduce dark ice extent and increase surface albedo in the GrIS ablation zone (Figure 11).
Recent work highlights the darkening effect of biological LAI on bare ice albedo, including assemblages of biologically-active dust termed "cryoconite" that melt quasi-cylindrical holes into the ice (Figure 11) [298] and distributed communities of algae and cyanobacteria that inhabit the ice surface [299]. For example, Tedstone et al. [51] applied the bare ice and dark ice detection algorithm from Shimada et al. [297] to MODIS surface reflectance imagery of the southwest ablation zone during June-July-August for the period 2000-2016. In contrast to Shimada et al. [297], they conclude that distributed algae blooms, rather than bare ice weathering and cryoconite hole growth, likely explains dark ice dynamics in their study region, owing to the synchronous and abrupt timing of dark ice exposure across the study area and the progressive rather than episodic increase in dark ice extent from June to August. In general, the deepening of cryoconite holes into the ice surface appears to limit their regional effect on albedo, especially at a low sun angle (Figure 11d), whereas distributed LAI that accumulate on the ice surface are considered stronger agents of bare ice darkening in the GrIS ablation zone [58,300]. Understanding seasonal and interannual relationships between bare ice structure, mineral and biological LAI, and bare ice albedo driven by cryoconite hole deepening and removal, weathering crust growth and decay, rotting and removal of superimposed ice, and variations in ice grain and air bubble geometry remain important areas of research that currently elude confident detection by satellite remote sensing [51,85,297,301,302].  [26] (courtesy Johnathan Ryan, Brown University). Daily reflectance maps for June, July, and August were classified into bare ice, snowcovered, and water-covered pixels using supervised random forest classification. The bare ice presence index is an exposure frequency representing the fraction of total days classified as bare ice for each pixel. The average end-of-summer snowline elevation is 1520 ± 113 m in this sector, with interannual variation of ±385 m. Interannual snowline variability explains 53% of MOD10A albedo variability [26].
The spatial extent of bare ice and albedo variability within the bare ice zone exerts a primary control on ice sheet albedo and surface meltwater production in the ablation zone [26,30]. Shimada et al. [297] used MODIS/Terra surface reflectance to quantify the regional distribution of bare ice and Daily reflectance maps for June, July, and August were classified into bare ice, snow-covered, and water-covered pixels using supervised random forest classification. The bare ice presence index is an exposure frequency representing the fraction of total days classified as bare ice for each pixel. The average end-of-summer snowline elevation is 1520 ± 113 m in this sector, with interannual variation of ±385 m. Interannual snowline variability explains 53% of MOD10A albedo variability [26]. Recent work highlights the darkening effect of biological LAI on bare ice albedo, including assemblages of biologically-active dust termed "cryoconite" that melt quasi-cylindrical holes into the ice (Figure 11) [298] and distributed communities of algae and cyanobacteria that inhabit the ice surface [299]. For example, Tedstone et al. [51] applied the bare ice and dark ice detection algorithm from Shimada et al. [297] to MODIS surface reflectance imagery of the southwest ablation zone during June-July-August for the period 2000-2016. In contrast to Shimada et al. [297], they conclude that distributed algae blooms, rather than bare ice weathering and cryoconite hole growth, likely explains dark ice dynamics in their study region, owing to the synchronous and abrupt timing of dark ice exposure across the study area and the progressive rather than episodic increase in dark ice extent from June to August. In general, the deepening of cryoconite holes into the ice surface appears to limit their regional effect on albedo, especially at a low sun angle (Figure 11d), whereas distributed LAI that accumulate on the ice surface are considered stronger agents of bare ice darkening in the GrIS ablation zone [58,300]. Understanding seasonal and interannual relationships between bare ice structure, mineral and biological LAI, and bare ice albedo driven by cryoconite hole deepening and removal, weathering crust growth and decay, rotting and removal of superimposed ice, and variations in ice grain and air bubble geometry remain important areas of research that currently elude confident detection by satellite remote sensing [51,85,297,301,302]. At the low and mid-elevation sites, the ice surface is glazed and smooth, and the water table within the cryoconite holes is nearly coincident with the ice surface. At the high-elevation site (c), the surface is rougher, and evidence of nocturnal refreezing is visible. (d) Same location as (c), showing the rough, weathered ice surface at low sun angle, reprinted from Cooper et al. [85]. Quadrat shown in (b) is 3 m wide. Cryoconite holes in (c) are on the order of 1-5 cm wide. Seasonal weathering of the ice surface, including cryoconite hole deepening and removal, exerts a primary control on ice roughness and grain morphology but its effect on bare ice albedo has received little direct study [297].

Current Challenges and Future Opportunities
Contemporary research on GrIS ablation zone surface reflectance and albedo is focused on understanding and diagnosing agents of albedo change [253]. Important questions remaining unresolved include whether observed reductions in GrIS snow albedo are caused by enhanced snow grain metamorphism due to a warming climate [252], or by deposition of LAI such as dust and black carbon from industrial emissions and/or forest fire [292,303]. These questions cannot currently be answered owing to the spectral and radiometric limitations of existing spaceborne optical sensors Figure 11. Images of cryoconite hole-studded ice surface in the western Greenland Ice Sheet ablation zone, collected at (a)~850 m a.s.l.; (b)~950 m a.s.l.; and (c)~1200 m a.s.l. along an elevation transect outside Kangerlussuaq (Søndre Strømfjord). At the low and mid-elevation sites, the ice surface is glazed and smooth, and the water table within the cryoconite holes is nearly coincident with the ice surface. At the high-elevation site (c), the surface is rougher, and evidence of nocturnal refreezing is visible. (d) Same location as (c), showing the rough, weathered ice surface at low sun angle, reprinted from Cooper et al. [85]. Quadrat shown in (b) is 3 m wide. Cryoconite holes in (c) are on the order of 1-5 cm wide. Seasonal weathering of the ice surface, including cryoconite hole deepening and removal, exerts a primary control on ice roughness and grain morphology but its effect on bare ice albedo has received little direct study [297].

Current Challenges and Future Opportunities
Contemporary research on GrIS ablation zone surface reflectance and albedo is focused on understanding and diagnosing agents of albedo change [253]. Important questions remaining unresolved include whether observed reductions in GrIS snow albedo are caused by enhanced snow grain metamorphism due to a warming climate [252], or by deposition of LAI such as dust and black carbon from industrial emissions and/or forest fire [292,303]. These questions cannot currently be answered owing to the spectral and radiometric limitations of existing spaceborne optical sensors and inadequate surface validation measurements [280,304]. For example, it is unlikely that the albedo-reducing effect of industrial black carbon emissions on relatively clean accumulation zone snow can be detected from space [305,306]. Similarly, it is unknown if bare ice albedo reduction due to deposition or emergence of inorganic mineral LAI can be distinguished using satellite imagery from albedo reduction due to biological LAI [292,301].
To date, efforts to remotely sense changing bare ice albedo have focused on quasi-physical proxies such as "darkening" in broad spectral bands [51,58]. Hyperspectral reflectance is used to map the albedo-reducing effect of LAI, including snow algae and dust in mountain environments and valley glaciers elsewhere [307,308]. However, ice algae are taxonomically distinct from snow algae [291], and optical methods developed by the snow albedo community need to be adapted to detect ice algae or separate its effect on ice albedo from inorganic impurities or variations in ice grain size and structure [292,301]. At present, no spaceborne hyperspectral imager operates over the polar regions, although several hyperspectral missions are planned for the coming years [309]. Recently, multispectral reflectance imagery from the OLCI onboard Sentinel-3A was used to infer the spatiotemporal distribution of ice algal blooms in the southwest GrIS ablation zone, providing insight into the capability of enhanced spectral resolution for diagnosing albedo change in the polar regions [269]. While this represents an important first step toward bioalbedo detection from space, additional ground validation of cellular optical properties specific to the microbial communities on the GrIS surface is needed to discriminate the individual drivers of bare ice albedo reduction using satellite remote sensing [301,310].

Mapping Surface Melt and Glaciological Zones
Mapping of surface zones or facies, such as dry snow, wet snow, and bare ice, was pioneered by Benson [35], and later mapped across the GrIS using both ERS SAR imagery [50] and Seasat-A scatterometer backscatter [311]. In contrast to the climatological facies of Benson [35], the concept of radar glacier zones was introduced to distinguish dynamic zones driven by the seasonal cycle of surface melt onset, surface roughening, and snowline migration, thereby revealing the seasonal extent of the bare ice ablation zone [312][313][314]. Multi-angular optical reflectance, together with derived surface roughness, provides an independent method to map glacier zones on the GrIS with particular relevance to ablation zone studies owing to the unique ability of angular reflectance to detect crevasse fields and the lower limit of superimposed ice [53]. These classification approaches are complemented by an extensive legacy of surface melt detection studies using active microwave backscatter, and thermal and passive microwave brightness temperature, which together provide a comprehensive view of the ablation zone surface and changes in ice sheet melt extent with time [55,315-319].

Active Microwave Detection of Surface Melt and Glacier Zones
SAR imagers measure radar backscatter (image brightness) at (typically) C-band (5 GHz) and Ku-band (14 GHz) frequencies (Table 3). Backscatter strength, or the normalized radar cross-section (σ • ), is principally controlled by the complex electrical permittivity of the ice surface and by geometric properties, including ice surface roughness, grain size, incidence angle, and ice thermal and physical structure [50,319]. Variations in σ • associated with these controls are used to map zones of unique snow and ice composition on glacier and ice sheet surfaces [320]. An early demonstration of the method used variations in C-band SAR image brightness to map zones of dry snow, percolation, wet snow, and bare ice, extending the early facies work of Benson [35] into the satellite era [50]. Cold, dry, fine-grained snow in the accumulation zone appears dark in C-band imagery, whereas the percolation zone and wet snow zone appear bright, owing to reflections from refrozen ice lenses in the snow/firn column [316]. The bare ice zone is distinguished from the wet snow zone using differential brightness between winter and summer images and the bright reflections of rough, crevassed surfaces [50,321]. Table 3. Summary of active microwave (synthetic aperture radar and microwave scatterometer) remote sensing platforms, instruments, temporal coverage, observed frequencies, and managing agencies.

Radar Altimetry Sensors and Datasets
The first satellite altimeter measurements of surface topography for the Greenland Ice Sheet were made by the GEOS-3, Geosat, and Seasat oceanographic Ku-band altimeters (Error! Reference source not found.) [37,91,106,121,122]. Their geographical coverage was limited to ±72 o but the data they collected showed it was possible to calculate volumetric changes in the polar ice sheets from space, providing proof of concept for future polar-orbiting altimeters [92]. A first demonstration using GEOS-3, Geosat, and Seasat data found that for areas south of <72 o N, the GrIS thickened by 23 ± 6 cm yr -1 between 1978-1986, with enhanced snowfall in a warmer climate hypothesized to explain the observed thickening [123]. These early missions supported important developments in waveform retracking, slope correction, and physical and empirical models for subsurface volume scattering [106,108,111,124]. These corrections are critical for accurate elevation retrieval and reliable change detection [125,126]. For example, Davis et al. [127] used an improved geodetic model and an improved retracking model to reinterpret the earlier findings of Zwally et al. [92], finding a smaller 1.5 ± 0.5 cm yr -1 thickening for the period 1978-1986 that highlighted the importance of slope-and signal penetration-correction.  (Table A1) an The ESA European Remote-sensing S altimeter, extending coverage of the polar ice 2 and EnviSat missions, these satellites pro continuous data acquisition [38]. The first di surface was produced from combined ERS-1 photogrammetric and cartographic data sets, ERS-1/2 altimetry data were used to discrimin m a.s.l.) of -2.0 ± 0.9 cm yr -1 from accumulat with an overall net thickening of 5.4 ± 0.2 temporal coverage and increase spatial measu developed an algorithm to merge ERS-1/2 satellite biases. For the period 1992-2008, th around 2000, and accelerated between 2006-2 a.s.l. and thinning of -7.0 ± 1.0 cm yr -1 below 1 Conventional radar altimeters such as t optimized for detecting elevation over open mean accuracies typically <0.2 m for ice shee / uncertainty is much higher for the top spatial measurement footprint and wide orbi and systematic biases over rough surfaces discriminate ablation zone thinning rates usin track analysis, but the accuracy of these dat owing to the aforementioned slope errors an zone [118,[133][134][135].
The CryoSat-2 satellite (launched April 2 employs SAR delay-Doppler along-track p [120,136]. CryoSat-2 is the first ESA radar alti to ±88 o and 1.6 km cross-track spacing at 70 Radar Altimeter (SIRAL) systems. SIRAL ope low resolution mode (LRM) over open ocean (SAR) over sea ice and coastal regions, and margins. The SARIn mode employs two Interferometric phase processing of the dual w track resolution [137,138].
Over mild slopes (<1 o ), phase unwrappin wide-swath elevation retrievals with two ord surface area than any prior radar altimeter the microwave signal [110][111][112][113]. l radar altimeter frequencies (Ku-band ~13 GHz), solid ice and liquid water are highly reas dry snow and firn are semi-transparent, with Ku-band penetration depths on the al meters [114,115]. Spatial and temporal variations in signal penetration depth caused are ice exposure, snow and firn densification, surface meltwater presence, and eltwater introduce spurious height change signals that are poorly quantified for the estimated to exceed >0.50 m in some cases [116][117][118]. Variations in surface slope r altimeter accuracy in the topographically-complex bare-ice ablation zone, with slopes that exceed >20 m for slopes >1° [39,109,113]. The recent CryoSat-2 radar altimeter ys synthetic aperture radar (SAR) and interferometric synthetic aperture radar (InSAR) hich decrease the effective along-track footprint diameter to the order 0.1-0.3 km and track slope from InSAR processing [93]. Together with denser orbital-track spacing and aveform retracking, topographic mapping of ice sheet ablation zone areas is accurate w meters on average, and / to within sub-meter uncertainty relative to laser 19,120].
ltimetry Sensors and Datasets satellite altimeter measurements of surface topography for the Greenland Ice Sheet the GEOS-3, Geosat, and Seasat oceanographic Ku-band altimeters (Error! Reference nd.) [37,91,106,121,122]. Their geographical coverage was limited to ±72 o but the data showed it was possible to calculate volumetric changes in the polar ice sheets from ng proof of concept for future polar-orbiting altimeters [92]. A first demonstration , Geosat, and Seasat data found that for areas south of <72 o N, the GrIS thickened by 23 een 1978-1986, with enhanced snowfall in a warmer climate hypothesized to explain hickening [123]. These early missions supported important developments in waveform pe correction, and physical and empirical models for subsurface volume scattering 4]. These corrections are critical for accurate elevation retrieval and reliable change ,126]. For example, Davis et al. [127] used an improved geodetic model and an acking model to reinterpret the earlier findings of Zwally et al. [92], finding a smaller r -1 thickening for the period 1978-1986 that highlighted the importance of slope-and tion-correction. ummary of radar and laser altimeter remote sensing platforms, instruments, temporal observed wavelength, and managing agencies.  (Table A1) and may not reflect joint collaborations.
The ESA European Remote-sensing Satellite (ERS-1) carried the first polar-orbiting radar altimeter, extending coverage of the polar ice sheets to ±81.5 o [95]. Together with the follow-on ERS-2 and EnviSat missions, these satellites provided near complete GrIS coverage and 21 years of continuous data acquisition [38]. The first digital elevation model (DEM) spanning the entire GrIS surface was produced from combined ERS-1 and Geosat altimetry data, supplemented by stereophotogrammetric and cartographic data sets, with pan-GrIS average accuracy -0.33 ± 6.97 m [128,129]. ERS-1/2 altimetry data were used to discriminate thinning rates in the ablation zone (defined as <1500 m a.s.l.) of -2.0 ± 0.9 cm yr -1 from accumulation zone (>1500 m a.s.l.) thickening of 6.4 ± 0.2 cm yr -1 with an overall net thickening of 5.4 ± 0.2 cm yr -1 for the period 1992-2003 [130]. To extend the temporal coverage and increase spatial measurement density, Khvorostovsky and Johannessen [131] developed an algorithm to merge ERS-1/2 and EnviSat datasets and detect and eliminate intersatellite biases. For the period 1992-2008, their dataset suggests net thinning of the GrIS initiated around 2000, and accelerated between 2006-2008, with a thickening of 4.0 ± 0.2 cm yr -1 above 1500 m a.s.l. and thinning of -7.0 ± 1.0 cm yr -1 below 1500 m a.s.l. [94]. Conventional radar altimeters such as the instruments included on ERS-1/2 and EnviSat were optimized for detecting elevation over open oceans and low-gradient polar ice sheet interiors with mean accuracies typically <0.2 m for ice sheet surface slopes <0.2 o [39,109]. Elevation accuracy and / uncertainty is much higher for the topographically complex ablation zone due to their large spatial measurement footprint and wide orbital track spacing, with slope errors >20 m for slopes >1 o and systematic biases over rough surfaces [39,95,109,132]. Careful processing has been used to discriminate ablation zone thinning rates using ERS-1/2 and EnviSat altimetry data, including repeattrack analysis, but the accuracy of these data and their process-based interpretation is ambiguous owing to the aforementioned slope errors and sparse spatial measurement density in the ablation zone [118,[133][134][135].
The CryoSat-2 satellite (launched April 2010) provides continuity with ERS-1/2 and EnviSat and employs SAR delay-Doppler along-track processing and interferometric cross-track processing [120,136] (Table A1) and may not reflect joint collaborations.
The ESA European Remote-sensing Satellite (ERS-1) carried the first polar-orbiting radar altimeter, extending coverage of the polar ice sheets to ±81.5 o [95]. Together with the follow-on ERS-2 and EnviSat missions, these satellites provided near complete GrIS coverage and 21 years of continuous data acquisition [38]. The first digital elevation model (DEM) spanning the entire GrIS surface was produced from combined ERS-1 and Geosat altimetry data, supplemented by stereophotogrammetric and cartographic data sets, with pan-GrIS average accuracy -0.33 ± 6.97 m [128,129]. ERS-1/2 altimetry data were used to discriminate thinning rates in the ablation zone (defined as <1500 m a.s.l.) of -2.0 ± 0.9 cm yr -1 from accumulation zone (>1500 m a.s.l.) thickening of 6.4 ± 0.2 cm yr -1 with an overall net thickening of 5.4 ± 0.2 cm yr -1 for the period 1992-2003 [130]. To extend the temporal coverage and increase spatial measurement density, Khvorostovsky and Johannessen [131] developed an algorithm to merge ERS-1/2 and EnviSat datasets and detect and eliminate intersatellite biases. For the period 1992-2008, their dataset suggests net thinning of the GrIS initiated around 2000, and accelerated between 2006-2008, with a thickening of 4.0 ± 0.2 cm yr -1 above 1500 m a.s.l. and thinning of -7.0 ± 1.0 cm yr -1 below 1500 m a.s.l. [94].
Conventional radar altimeters such as the instruments included on ERS-1/2 and EnviSat were optimized for detecting elevation over open oceans and low-gradient polar ice sheet interiors with mean accuracies typically <0.2 m for ice sheet surface slopes <0.2 o [39,109]. Elevation accuracy and / uncertainty is much higher for the topographically complex ablation zone due to their large spatial measurement footprint and wide orbital track spacing, with slope errors >20 m for slopes >1 o and systematic biases over rough surfaces [39,95,109,132]. Careful processing has been used to discriminate ablation zone thinning rates using ERS-1/2 and EnviSat altimetry data, including repeattrack analysis, but the accuracy of these data and their process-based interpretation is ambiguous owing to the aforementioned slope errors and sparse spatial measurement density in the ablation zone [118,[133][134][135].
The CryoSat-2 satellite (launched April 2010) provides continuity with ERS-1/2 and EnviSat and employs SAR delay-Doppler along-track processing and interferometric cross-track processing [120,136]. CryoSat-2 is the first ESA radar altimeter dedicated to polar studies, with orbital coverage to ±88 o and 1.6 km cross-track spacing at 70 o . CryoSat-2 carries two Ku-band SAR/Interferometric Radar Altimeter (SIRAL) systems. SIRAL operates in three distinct modes depending on location: 1) low resolution mode (LRM) over open ocean and ice sheet interiors, 2) synthetic aperture radar mode (SAR) over sea ice and coastal regions, and 3) SAR interferometric mode (SARIn) over ice sheet margins. The SARIn mode employs two SAR instruments oriented across the satellite track. Interferometric phase processing of the dual waveforms provides cross-track slope at ~0.3 km along-Managing agencies are identified by the WMO OSCAR database (Table A1) and may not reflect joint collaborations. C-and Ku-band wind scatterometers measure σ • at both vertical (VV) and horizontal (HH) polarization and are used to map the timing of melt onset and the seasonal progression of surface melt extent [313,322]. As with SAR, the method exploits the strong reduction in σ • caused by the presence of liquid water at or near the ice surface. Relative to SAR imagers, scatterometers provide σ • at coarse spatial resolution but high temporal resolution, for example, a 25 km grid size twice-daily for the SeaWinds on NASA's Quick Scatterometer (QuikSCAT) [322]. Resolution enhancement techniques increase the effective resolution to~8-10 km [323].
Surface melt on the GrIS and peripheral ice caps has been detected by the Seasat-A scatterometer (SASS), the Advanced Microwave Instrument-Scatterometer (AMI-SCAT) on ERS-1/2, and the SeaWinds scatterometer on both NASA's Quick Scatterometer (QuikSCAT) and ADEOS-2 (Table 3) [311,320,322,[324][325][326]. Methods to detect surface melt include single-channel absolute σ • thresholds [324], diurnal σ • variability (relative thresholds) [322], physical and statistical model-based methods [326], and dual-frequency/polarization thresholds that also use the diurnal difference between ascending and descending orbital passes [327]. The diurnal method exploits contrasts in σ • caused by diurnal melt-freeze cycles. The use of relative thresholds reduces errors from sensor drift, cross-mission biases, or step changes in surface properties, such as ice lens formation that affect absolute thresholds [115]. In addition to discrete melt onset, time-integrated σ • reduction is used to infer seasonal melt intensity [324][325][326].
In addition to mapping surface melt onset and extent, seasonal changes in σ • are used to infer the timing and spatial extent of ice layer formation in snow and firn [115,316,328]. Understanding ice layer formation is important because meltwater refreezing increases firn density without reducing mass and raises the effective backscattering surface detected by radar altimeters, leading to errors in radar mass balance estimates (Section 3.1.2). Nghiem et al. [115] developed a field-validated method that relates threshold increases in QuickSCAT HH-polarized σ o before and after melt seasons to ice layer formation in the GrIS percolation zone. The method also provides a basis for estimating snow accumulation by integrating σ o reduction following ice layer formation. The method assumes threshold increases in σ o are caused by enhanced reflections from newly formed ice layers, whereas the gradual attenuation of σ • following ice layer formation is caused by new snow accumulation. Wang et al. [328] applied the method to five years of enhanced-resolution QuickSCAT imagery [329] and found extensive increases in ice layer formation following a short three-day melt event in 2002, highlighting the disproportionate impact of extreme melt events on ice layer formation.

Passive Microwave and Thermal Radiometry
At thermal and microwave wavelengths (beyond about 3 um) spectral radiance (converted via the Planck function to brightness temperature, T b ) is approximately linear with T sfc : T b = ε * T sfc , where ε is the material emissivity [330]. Whereas σ • is dramatically reduced by liquid water in snow, T b is dramatically increased, forming the basis for surface melt detection threshold algorithms, typically using a threshold value below freezing to indicate melting [315,331,332]. Passive microwave radiometers measure background microwave emission at (primarily) K-band (~19 GHz) and Ka-band (~37 GHz) frequencies and various combinations of VV and HH polarization (Table 4). In contrast to the higher spatial resolution (~10 km 2 ) but the lower temporal resolution (weeks to months) of spaceborne active microwave sensors, passive microwave sensors provide T b twice-daily at 25-50 km 2 spatial resolution and complete ice sheet coverage. Melt detection methods using T b include single-channel thresholds [332,333], dual frequency/polarization combinations (cross polarized gradient ratio XPGR) [315,334], and diurnal amplitude variations (DAV) on ascending and descending passes [335,336]. As with σ • , microwave T b is strongly modified by liquid meltwater presence at the ice sheet surface but typically does not provide information about the internal snow or firn structure [324,337]. Table 4. Summary of passive microwave remote sensing platforms, instruments, temporal coverage, observed frequencies, and managing agencies. [104,105]. The following subsections summarize spaceborne radar altimetry of the GrIS, providing an historical perspective that emphasizes major challenges and advances in the field, and identifies future opportunities for radar altimetry observations of the GrIS ablation zone.

Radar Altimetry
Conventional radar altimeters are nadir pointing single-beam radars that transmit and receive tens to thousands of microwave pulses per second, yielding ground measurement footprint diameters on the order 1-10 km posted at 0.3-0.6 km along-track spacing [39,106,107]. An estimate of / is obtained by comparing measurements at crossover points, where an earlier satellite ground track intersects a later one [92]. Two critical factors affect the accuracy of radar altimeter elevation retrievals over ice sheets: 1) variation in topographic slope, which controls the average distance between the altimeter and the measured ground footprint surface [39,108,109], and 2) changes in the electrical permittivity of the ice sheet surface, which controls the reflection, transmission, and absorption of the microwave signal [110][111][112][113].
For typical radar altimeter frequencies (Ku-band ~13 GHz), solid ice and liquid water are highly reflective whereas dry snow and firn are semi-transparent, with Ku-band penetration depths on the order of several meters [114,115]. Spatial and temporal variations in signal penetration depth caused by seasonal bare ice exposure, snow and firn densification, surface meltwater presence, and refreezing of meltwater introduce spurious height change signals that are poorly quantified for the GrIS but are estimated to exceed >0.50 m in some cases [116][117][118]. Variations in surface slope dominate radar altimeter accuracy in the topographically-complex bare-ice ablation zone, with slopeinduced errors that exceed >20 m for slopes >1° [39,109,113]. The recent CryoSat-2 radar altimeter mission employs synthetic aperture radar (SAR) and interferometric synthetic aperture radar (InSAR) technologies which decrease the effective along-track footprint diameter to the order 0.1-0.3 km and provide cross-track slope from InSAR processing [93]. Together with denser orbital-track spacing and advances in waveform retracking, topographic mapping of ice sheet ablation zone areas is accurate to within a few meters on average, and / to within sub-meter uncertainty relative to laser altimetry [93,119,120].

Radar Altimetry Sensors and Datasets
The first satellite altimeter measurements of surface topography for the Greenland Ice Sheet were made by the GEOS-3, Geosat, and Seasat oceanographic Ku-band altimeters (Error! Reference source not found.) [37,91,106,121,122]. Their geographical coverage was limited to ±72 o but the data they collected showed it was possible to calculate volumetric changes in the polar ice sheets from space, providing proof of concept for future polar-orbiting altimeters [92]. A first demonstration using GEOS-3, Geosat, and Seasat data found that for areas south of <72 o N, the GrIS thickened by 23 ± 6 cm yr -1 between 1978-1986, with enhanced snowfall in a warmer climate hypothesized to explain the observed thickening [123]. These early missions supported important developments in waveform retracking, slope correction, and physical and empirical models for subsurface volume scattering [106,108,111,124]. These corrections are critical for accurate elevation retrieval and reliable change detection [125,126]. For example, Davis et al. [127] used an improved geodetic model and an improved retracking model to reinterpret the earlier findings of Zwally et al. [92], finding a smaller 1.5 ± 0.5 cm yr -1 thickening for the period 1978-1986 that highlighted the importance of slope-and signal penetration-correction.  (Table A1) and may no The ESA European Remote-sensing Satellite ( altimeter, extending coverage of the polar ice sheets t 2 and EnviSat missions, these satellites provided n continuous data acquisition [38]. The first digital ele surface was produced from combined ERS-1 and G photogrammetric and cartographic data sets, with pan ERS-1/2 altimetry data were used to discriminate thinn m a.s.l.) of -2.0 ± 0.9 cm yr -1 from accumulation zone with an overall net thickening of 5.4 ± 0.2 cm yr -1 temporal coverage and increase spatial measurement developed an algorithm to merge ERS-1/2 and Env satellite biases. For the period 1992-2008, their data around 2000, and accelerated between 2006-2008, wit a.s.l. and thinning of -7.0 ± 1.0 cm yr -1 below 1500 m a Conventional radar altimeters such as the instru optimized for detecting elevation over open oceans a mean accuracies typically <0.2 m for ice sheet surfac / uncertainty is much higher for the topograph spatial measurement footprint and wide orbital track and systematic biases over rough surfaces [39,95,1 discriminate ablation zone thinning rates using ERS-1 track analysis, but the accuracy of these data and th owing to the aforementioned slope errors and spars zone [118,[133][134][135].
The CryoSat-2 satellite (launched April 2010) pro employs SAR delay-Doppler along-track processin [120,136] following subsections summarize spaceborne radar altimetry of the GrIS, providing an pective that emphasizes major challenges and advances in the field, and identifies nities for radar altimetry observations of the GrIS ablation zone.
etry onal radar altimeters are nadir pointing single-beam radars that transmit and receive ands of microwave pulses per second, yielding ground measurement footprint he order 1-10 km posted at 0.3-0.6 km along-track spacing [39,106,107]. An estimate of ned by comparing measurements at crossover points, where an earlier satellite ground s a later one [92]. Two critical factors affect the accuracy of radar altimeter elevation r ice sheets: 1) variation in topographic slope, which controls the average distance timeter and the measured ground footprint surface [39,108,109], and 2) changes in the ittivity of the ice sheet surface, which controls the reflection, transmission, and the microwave signal [110][111][112][113].
l radar altimeter frequencies (Ku-band ~13 GHz), solid ice and liquid water are highly reas dry snow and firn are semi-transparent, with Ku-band penetration depths on the al meters [114,115]. Spatial and temporal variations in signal penetration depth caused are ice exposure, snow and firn densification, surface meltwater presence, and eltwater introduce spurious height change signals that are poorly quantified for the estimated to exceed >0.50 m in some cases [116][117][118]. Variations in surface slope r altimeter accuracy in the topographically-complex bare-ice ablation zone, with slopes that exceed >20 m for slopes >1° [39,109,113]. The recent CryoSat-2 radar altimeter ys synthetic aperture radar (SAR) and interferometric synthetic aperture radar (InSAR) hich decrease the effective along-track footprint diameter to the order 0.1-0.3 km and track slope from InSAR processing [93]. Together with denser orbital-track spacing and aveform retracking, topographic mapping of ice sheet ablation zone areas is accurate w meters on average, and / to within sub-meter uncertainty relative to laser 19,120].
ltimetry Sensors and Datasets satellite altimeter measurements of surface topography for the Greenland Ice Sheet the GEOS-3, Geosat, and Seasat oceanographic Ku-band altimeters (Error! Reference nd.) [37,91,106,121,122]. Their geographical coverage was limited to ±72 o but the data showed it was possible to calculate volumetric changes in the polar ice sheets from ng proof of concept for future polar-orbiting altimeters [92]. A first demonstration , Geosat, and Seasat data found that for areas south of <72 o N, the GrIS thickened by 23 een 1978-1986, with enhanced snowfall in a warmer climate hypothesized to explain hickening [123]. These early missions supported important developments in waveform pe correction, and physical and empirical models for subsurface volume scattering 4]. These corrections are critical for accurate elevation retrieval and reliable change ,126]. For example, Davis et al. [127] used an improved geodetic model and an acking model to reinterpret the earlier findings of Zwally et al. [92], finding a smaller r -1 thickening for the period 1978-1986 that highlighted the importance of slope-and tion-correction. ummary of radar and laser altimeter remote sensing platforms, instruments, temporal observed wavelength, and managing agencies.  (Table A1) and may not reflect joint collaborations.
The ESA European Remote-sensing Satellite (ERS-1) carried the first polar-orbiting radar altimeter, extending coverage of the polar ice sheets to ±81.5 o [95]. Together with the follow-on ERS-2 and EnviSat missions, these satellites provided near complete GrIS coverage and 21 years of continuous data acquisition [38]. The first digital elevation model (DEM) spanning the entire GrIS surface was produced from combined ERS-1 and Geosat altimetry data, supplemented by stereophotogrammetric and cartographic data sets, with pan-GrIS average accuracy -0.33 ± 6.97 m [128,129]. ERS-1/2 altimetry data were used to discriminate thinning rates in the ablation zone (defined as <1500 m a.s.l.) of -2.0 ± 0.9 cm yr -1 from accumulation zone (>1500 m a.s.l.) thickening of 6.4 ± 0.2 cm yr -1 with an overall net thickening of 5.4 ± 0.2 cm yr -1 for the period 1992-2003 [130]. To extend the temporal coverage and increase spatial measurement density, Khvorostovsky and Johannessen [131] developed an algorithm to merge ERS-1/2 and EnviSat datasets and detect and eliminate intersatellite biases. For the period 1992-2008, their dataset suggests net thinning of the GrIS initiated around 2000, and accelerated between 2006-2008, with a thickening of 4.0 ± 0.2 cm yr -1 above 1500 m a.s.l. and thinning of -7.0 ± 1.0 cm yr -1 below 1500 m a.s.l. [94]. Conventional radar altimeters such as the instruments included on ERS-1/2 and EnviSat were See Abbreviations at end of article for expanded acronyms.  (Table A1) and may not reflect joint collaborations.
The ESA European Remote-sensing Satellite (ERS-1) carried the first polar-orbiting radar altimeter, extending coverage of the polar ice sheets to ±81.5 o [95]. Together with the follow-on ERS-2 and EnviSat missions, these satellites provided near complete GrIS coverage and 21 years of continuous data acquisition [38]. The first digital elevation model (DEM) spanning the entire GrIS surface was produced from combined ERS-1 and Geosat altimetry data, supplemented by stereophotogrammetric and cartographic data sets, with pan-GrIS average accuracy -0.33 ± 6.97 m [128,129]. ERS-1/2 altimetry data were used to discriminate thinning rates in the ablation zone (defined as <1500 m a.s.l.) of -2.0 ± 0.9 cm yr -1 from accumulation zone (>1500 m a.s.l.) thickening of 6.4 ± 0.2 cm yr -1 with an overall net thickening of 5.4 ± 0.2 cm yr -1 for the period 1992-2003 [130]. To extend the temporal coverage and increase spatial measurement density, Khvorostovsky and Johannessen [131] developed an algorithm to merge ERS-1/2 and EnviSat datasets and detect and eliminate intersatellite biases. For the period 1992-2008, their dataset suggests net thinning of the GrIS initiated around 2000, and accelerated between 2006-2008, with a thickening of 4.0 ± 0.2 cm yr -1 above 1500 m a.s.l. and thinning of -7.0 ± 1.0 cm yr -1 below 1500 m a.s.l. [94].
Conventional radar altimeters such as the instruments included on ERS-1/2 and EnviSat were Managing agencies are identified by the WMO OSCAR database (Table A1) and may not reflect joint collaborations.
With more than 30 years of continuous data collection, pan-GrIS spatial coverage, and all-weather capability, spaceborne passive microwave radiometers provide a unique insight into the climatic drivers of ice sheet surface mass balance processes, including changes in the location and extent of surface melting ( Figure 12) [315,317,318,334,336,338,339]. Data from the Special Sensor Microwave/Imager (SSM/I) suggests the area of active surface melt over the GrIS has approximately doubled since the early 1990s, with a 40,000 km 2 yr −1 trend during this period [318,336]. Data from the Scanning Multichannel Microwave Radiometer (SMMR, 1978(SMMR, -1987 and SSM/I were used to characterize the effects of the Mt. Pinatubo eruption on GrIS surface melt patterns [340]. Surface melt extent from SSM/I has been correlated with downstream sediment plume concentrations in Greenland fjords driven by ice sheet meltwater discharge [341], to validate surface melt extent calculated from surface energy balance models [342,343], and to quantify extreme events such as the record July 2012 melt event when 98.6% of the GrIS surface was actively melting [287,339].  [287]. Surface melt frequency is calculated by the authors from the NASA MEaSUREs Greenland Surface Melt Daily 25 km EASE-Grid 2.0 data set [344]. Thermal radiance (~3-14 μm) is used to map sfc and provides an additional method for surface melt detection [345]. Variations in thermal and optical radiance form the basis for mapping "reflectance zones" to characterize melt presence and changes in the thermal structure on glacier surfaces [346]. Whereas σ° and microwave B are primarily diagnostic of melt presence, both at or near the surface, thermal radiance is diagnostic of the surface "skin" temperature, and provides little to no information about subsurface temperature, meltwater presence, or snow and firn structure. Consequently, thermal sfc is strictly an indicator of surface melt and, together with the higher spatial and radiometric resolution of thermal sensors such as MODIS, provides an independent method for validating active and passive microwave melt presence products [347,348]. When combined with microwave melt presence and remotely sensed albedo, thermal sfc may improve the discrimination of surface and subsurface melt areas and diurnal variations in melt-freeze cycles [347]. As with microwave surface melt detection, thermal radiance melt detection has been used to validate modeled surface melt [318,342], and to study the spatiotemporal variation of surface melt extent and duration on the GrIS and its relation to climatic variability [55,238,[348][349][350].
The primary spaceborne thermal sensors used to obtain sfc are MODIS and AVHRR. The NOAA Extended AVHRR Polar Pathfinder (App-X) product provides sfc on a 25 km polar equalarea grid twice-daily for the period 1982 to present [272]. As with the AVHRR albedo product (Section 5.2), there are spectral, radiometric, and inter-mission homogenization issues that limit its utility [351]. The MODIS land surface temperature product (MOD11A1) uses the split-window technique  [287]. Surface melt frequency is calculated by the authors from the NASA MEaSUREs Greenland Surface Melt Daily 25 km EASE-Grid 2.0 data set [344].
Thermal radiance (~3-14 µm) is used to map T sfc and provides an additional method for surface melt detection [345]. Variations in thermal and optical radiance form the basis for mapping "reflectance zones" to characterize melt presence and changes in the thermal structure on glacier surfaces [346]. Whereas σ • and microwave T B are primarily diagnostic of melt presence, both at or near the surface, thermal radiance is diagnostic of the surface "skin" temperature, and provides little to no information about subsurface temperature, meltwater presence, or snow and firn structure. Consequently, thermal T sfc is strictly an indicator of surface melt and, together with the higher spatial and radiometric resolution of thermal sensors such as MODIS, provides an independent method for validating active and passive microwave melt presence products [347,348]. When combined with microwave melt presence and remotely sensed albedo, thermal T sfc may improve the discrimination of surface and subsurface melt areas and diurnal variations in melt-freeze cycles [347]. As with microwave surface melt detection, thermal radiance melt detection has been used to validate modeled surface melt [318,342], and to study the spatiotemporal variation of surface melt extent and duration on the GrIS and its relation to climatic variability [55,238,[348][349][350].
The primary spaceborne thermal sensors used to obtain T sfc are MODIS and AVHRR. The NOAA Extended AVHRR Polar Pathfinder (App-X) product provides T sfc on a 25 km polar equal-area grid twice-daily for the period 1982 to present [272]. As with the AVHRR albedo product (Section 5.2), there are spectral, radiometric, and inter-mission homogenization issues that limit its utility [351]. The MODIS land surface temperature product (MOD11A1) uses the split-window technique developed for AVHRR to calculate T sfc from radiance at 10.78 µm and 11.77 µm [352]. The MOD11A1 data are provided on a 1 km grid globally for clear-sky conditions as discriminated by the MOD35 cloud mask. The MOD11A1 product is accurate to ± 1 • C on average over snow and ice surfaces and agrees to within ± 0.5 • C with T sfc calculated from ASTER and Landsat ETM+ thermal radiance over the GrIS [350]. Surface melt extent from MOD11A1 corresponds closely to surface melt inferred from QuickSCAT using the diurnal σ o method [347].
The Greenland Ice Surface Temperature (IST) product is an enhanced version of MOD11A1 that provides IST and binary melt absence/presence for the period 2000-2014 on a 1.25 km polar stereographic grid for the GrIS [349,351]. For values of T sfc near 0 • C, the IST data are −0.5 • C cooler than surface-based T sfc measurements collected at the Summit Station [353,354]. The bias increases to −5.0 • C for values of T sfc near −60 • C and increases with SZA, suggesting the bias may be related to reduced accuracy of the MOD35 cloud mask at high SZA. Under-sampling of T sfc during warm inversions with dense cloud coverage may also introduce bias but this effect has not been systematically evaluated [353,354]. In general, spaceborne thermal T sfc retrievals are effective for detecting surface melt during clear-sky conditions but are sensitive to the atmospheric aerosol and water vapor profile [351,355]. To overcome cloud-cover data gaps, Välisuo et al. [348] gap-filled the IST product with modeled values of T sfc from the European Centre for Medium-Range Weather Forecasts ERA-Interim reanalysis.

Multi-Angular Reflectance and Surface Roughness
The Multi-angle Imaging SpectroRadiometer (MISR) measures coincident-in-time bi-directional radiance in four spectral bands between 450-850 nm at 0 • (nadir), 26.1 • , 45.6 • , 60.0 • , and 70.5 • fore and aft of nadir. Reflectance is provided on a 275 m grid and nine-day global coverage to ±83 • (Table 2) [260]. The normalized difference angular index (NDAI) was developed to estimate ice surface roughness from MISR angular reflectance [260]. The method uses an empirical relationship between MISR red band reflectance at 60 o fore and aft of nadir and ATM lidar-derived surface roughness to develop spatially-continuous maps of surface roughness over ice sheets and sea ice [356]. The MISR NDAI and near-infrared albedo (Section 5.2) were used to map unique signatures of surface glaciological features in the western GrIS ablation zone, including crevasse fields, wet snow, and bare glacier ice [53]. With this approach, MISR angular reflectance appears to provide a unique method for detecting the superimposed ice zone and may improve the identification of changes in crevasse field roughness relative to SAR and optical sensors [53,260].
Surface roughness controls ice-atmosphere interactions via the aerodynamic roughness length and the net vapor flux, an understudied component of the GrIS SMB [60,82,357]. The MISR surface roughness product [356] was used to define the aerodynamic roughness length of the BMF13 vapor flux model to improve the spatial realism of net vapor flux in the ablation zone where roughness is highly variable [357]. The average annual modeled vapor flux was 14.6 ± 3.6 Gt yr −1 , or 6 ± 2% of annual SMB for the period 2003-2014. The average annual difference between modeled vapor flux with and without MISR roughness was 30 ± 15%. In addition to angular reflectance, surface roughness has been quantified from ICESat and ATM laser altimetry waveforms [173,358,359]. The forthcoming ICESat-2 laser altimeter may provide additional capability for measuring surface roughness on the GrIS, for example using the multi-sensor lidar-angular reflectance approach of Nolin and Mar [356] and Nolin et al. [260], which may also be useful for radar altimetry waveform interpretation of the leading edge of the beam footprint [93].

Future Opportunities for Mapping the Changing GrIS Ablation Zone Surface
The GrIS surface is undergoing rapid change, driven by increased surface meltwater production in the ablation zone. Spaceborne SAR imagers, wind scatterometers, and passive microwave, thermal, and angular radiometers are used to map and monitor diagnostic features of change on the ablation zone surface, including surface melt presence and extent, subsurface ice layer formation, ice surface temperature, and ice surface roughness. These characteristic features are used to define and map the dry snow zone, wet snow zone, percolation zone, superimposed ice zone, and the bare ice zone. Mapping of supraglacial hydrologic features, including meltwater lakes, rivers, and moulins stands out as an additional research priority c.f. [68], in particular, their expansion and inland migration toward sensitive (e.g., high elevation) areas of the ice sheet [56,360,361]. Other diagnostic markers of change detectable in satellite imagery include the end of summer snowline position (Figure 9) [26], the lower limit of superimposed ice [53], and the lower limit of the slush zone [293]. Multi-sensor methods, for example, optical imagery combined with SAR imagery or scatterometry, shows promise for detecting dynamic regions characterized by changing snow, firn, and ice surface types, including supraglacial lakes and slush fields obscured by snow or clouds, and may reduce detection bias caused by cloud cover [362][363][364].

Conclusions
For over forty years, earth-observing satellites sensitive to visible, infrared, and microwave electromagnetic radiation, together with gravimetry, have documented the patterns of change on the GrIS ablation zone surface. Satellite remote sensing data show an ablation zone expanded in size, its albedo and surface elevation lower in response to enhanced melting and ice discharge, and an ice sheet transition from steady state to negative mass balance that now represents the largest land ice contributor to global sea level rise. Already, some 78-85% of the total liquid runoff produced from surface melting is generated in the bare ice ablation zone, despite it covering~22% of the ice sheet's total surface area, up from~15% in the early 1960s [20,30,[33][34][35]. Although often conceptualized as a uniform surface of solid ice, the ablation zone is a dynamic region with widely varying electromagnetic properties controlled by diverse physical, biological, and hydrologic processes. Future progress in remote sensing the ablation zone will likely benefit most from methods that directly address this complexity, for example, using multi-sensor, multi-wavelength, and cross-platform datasets. Examples include fusing radar and laser altimetry with optical stereophotogrammetry to discriminate and diagnose causes of surface elevation change [174], or fusing radar and laser backscatter with optical imagery to discriminate snow, ice, liquid water, and refrozen meltwater in sensitive areas near the equilibrium line altitude [362,363]. Other areas of opportunity recommended for future research include spaceborne detection of subsurface refrozen meltwater and its effects on radar backscatter, which requires additional in-situ validation [115,116], the partitioning of ablation zone thinning into ice-dynamic and surface mass balance components [97], cross-validation of ice surface elevation change from altimetry with modeled surface mass balance [89] and modeled ice dynamic motion [189], spaceborne diagnosis of changing bare ice albedo [269] and grain size [191], and monitoring the inland migration of snowlines, surface melt extent, and surface hydrologic features including lakes, streams, and moulins [26,56]