Surface Melting Drives Fluctuations in Airborne Radar Penetration in West Central Greenland

Greenland Ice Sheet surface melting has increased since the 1990s, affecting the rheology and scattering properties of the near‐surface firn. We combine firn cores and modeled firn densities with 7 years of CryoVEx airborne Ku‐band (13.5 GHz) radar profiles to quantify the impact of melting on microwave radar penetration in West Central Greenland. Although annual layers are present in the Ku‐band radar profiles to depths up to 15 m below the ice sheet surface, fluctuations in summer melting strongly affect the degree of radar penetration. The extreme melting in 2012, for example, caused an abrupt 6.2 ± 2.4 m decrease in Ku‐band radar penetration. Nevertheless, retracking the radar echoes mitigates this effect, producing surface heights that agree to within 13.9 cm of coincident airborne laser measurements. We also examine 2 years of Ka‐band (34.5 GHz) airborne radar data and show that the degree of penetration is half that of coincident Ku‐band.


Introduction
In recent decades, increased melting at the surface of the Greenland Ice Sheet (van den Broeke et al., 2016) has had a marked impact on rates of runoff (Enderlin et al., 2014; and glacier flow (van de Wal et al., 2008), and has also affected the structure of the near-surface firn owing to the redistribution and refreezing of surface meltwater (de la Peña et al., 2015;Machguth et al., 2016). These processes, and the associated changes in firn properties, present challenges for satellite radar altimetry surveys of the ice sheet mass balance when converting observations of volume change to mass change (McMillan et al., 2016). Surface melting has a large impact on the firn stratigraphy and density as meltwater can refreeze at the surface, or percolate into the snowpack and refreeze to form ice lenses (Benson, 1996), or refreeze in between already existing ice layers and form thicker ice slabs . Ice cores provide records of density and stratigraphy at point locations (Mosley-Thompson et al., 2001) and firn densification models provide estimates of densities across the ice sheet (e.g., Kuipers Munneke et al., 2015). Radars have also been widely used over glaciers and ice sheets to map their structure (MacGregor et al., 2015), calculate accumulation rates (Miège et al., 2017), and track changes in their elevation .
Radar systems transmit electromagnetic pulses and record the amplitude and time delay of the waves scattered back from discontinuities in the dielectric properties. Their echoes are sensitive to density variations in the firn column, and the firn structure can reveal continuous internal scattering horizons, corresponding to isochrones (Hawley et al., 2006). Both ground-based (Brown et al., 2012) and airborne radar, such as the snow radar flown during NASA Operation IceBridge (Koenig et al., 2016;Medley et al., 2013;Montgomery et al., 2020) or the European Space Agency's Airborne SAR/Interferometric Radar Altimeter System (ASIRAS) operating at Ku-band (de la Peña et al., 2010;Helm et al., 2007;Overly et al., 2016;Simonsen et al., 2013), have been used to track isochrones and derive accumulation rates. Unlike ground-based and airborne radar, satellite radar measurements lack the vertical resolution to resolve the internal structure of the firn column due to the smaller bandwidth and coarser spatial resolution of the radar footprint. Nonetheless, satellite radar altimeters are sensitive to variations in the firn properties as the radar signal penetrates into the snowpack. Additionally, the radar signal penetration is frequency dependent. At present, two frequencies are used by satellite radar altimeters: Ku-band (13.5 GHz) is used by CryoSat-2 and Sentinel 3A/B, and Ka-band (37 GHz) is used by AltiKa. Studies have shown that at Ku-band, the radar signals can penetrate up to~15 m into firn, while the penetration depth of higher-frequency Ka-band radars is reduced to~0.5 m (Rémy et al., 2015). The altimeter echo recorded is therefore a combination of surface and volume scattering (Ridley & Partington, 1988). The ratio between surface and volume scattering varies spatially and temporally according to changes in the surface and subsurface properties and may impact the height retrieval from radar altimeters (Simonsen & Sørensen, 2017).
In this study, we use airborne Ku-band radar data acquired using the ASIRAS instrument, Ka-band radar data acquired using the KAREN radar (the MetaSensing Ka-band altimeter), airborne laser data, shallow firn cores (< 6 m), and firn density models to characterize and assess the impact of spatial and temporal fluctuations in the properties of the near-surface firn in West Central Greenland. This study is performed along the glaciological transect established during the Expéditions Glaciologiques Internationales au Groenland (EGIG) in 1958 (Renaud et al., 1963). The EGIG line extends from the ablation zone at the Western margin of the ice sheet, across the percolation and dry snow zones to the Summit, and further toward the Eastern margin and is therefore a representative location of density variations across the Greenland Ice Sheet .

Data and Methods
The EGIG line ( Figure 1a) has been surveyed for more than a decade as part of ESA's Cryosat Validation Experiment (CryoVEx). According to previous in situ investigations of snow density and stratigraphy, EGIG line sites T3 to T21 lie in the percolation zone with site T21 marking the start of the dry snow zone (Morris & Wingham, 2011;Scott et al., 2006). In this study, we use data collected between 2006 and 2017 over a~675 km transect of the EGIG line, starting about 14 km from Ilullissat airport at an elevation of 157 m above sea level (m.a.s.l.), 115 km before site T1, and ending 149 km beyond site T41 at an altitude of 2,956 m.a.s.l.
Shallow firn cores were collected in October 2016 (T1, T4, and T5) and April 2017 (T5, T9, T12, T19, T30, and T41) (Table S1 in the supporting information). The stratigraphy was recorded by illuminating the cores to identify ice lenses. Firn density was measured by weighing the different stratigraphic layers (see Supporting Information S1). These in situ density measurements are used to evaluate two firn densification models: (1) the stand-alone Institute for Marine and Atmospheric Research Utrecht firn densification model (IMAU-FDM) (Ligtenberg et al., 2011(Ligtenberg et al., , 2018, forced at the surface by RACMO . IMAU-FDM simulates the density and temperature in a vertical, one-dimensional firn column through time at a vertical resolution of 5 cm and (2) the CROCUS snow model (Brun et al., 1992) embedded in the Modèle Atmosphérique Régional (MAR) (version 3.10, Fettweis et al., 2017Fettweis et al., , 2020 to derive the MAR firn densification model (MAR-FDM). CROCUS simulates the transfer of mass and energy between a fixed number of layers of snow, firn, or ice with a vertical resolution varying from 5 to 500 cm.
Airborne radar data were collected in 2004, 2006, 2008, 2011, 2012, 2014, 2016, and 2017, along segments of the EGIG line of varying length. For all CryoVEx airborne surveys, the same aircraft and instrumental setup were used, which includes an airborne laser scanner (ALS), the ASIRAS radar, and the KAREN radar in two surveys in 2016 and 2017. The ALS is a Riegl LMS-Q140-i60 laser scanner operating at 904 nm (red), which gives surface heights in 0.7 m intervals across a 300-m-wide swath with an accuracy of 0.1 m (Skourup et al., 2019). ASIRAS is a Ku-band radar designed as a prototype for the SIRAL altimeter on board CryoSat-2 operating at a central frequency of 13.5 GHz with a bandwidth of 1.0 GHz, range resolution of 0.109 m in air, and a nominal footprint size of 10 m across-track and 3 m along-track at a flight elevation of~300 m a.g.l. (Cullen, 2010). KAREN is a Ka-band radar altimeter operating at a central frequency of 34.5 GHz-the same frequency used by AltiKa-with a bandwidth of~0.5 GHz, range resolution of 0.165 m in air, and a footprint size of 10 m across-track and 5 m along-track (version "levc"). Data acquired when the aircraft roll angle exceeded ±1.5°were discarded. We applied three different retracking algorithms to locate the ice sheet surface in the radar echoes as different processing strategies have been shown to affect the elevation measurements when the radar penetration depth varies . We used (1) an offset center of gravity (OCOG) retracker (Wingham et al., 1986), (2) a 30% threshold center of gravity (TCOG) retracker (Davis, 1997) similar to the retracker implemented in CryoSat-2 ground segment, and (3) a 50% threshold first maximum retracking algorithm (TFMRA, Helm et al., 2014). Each retracker's performance was evaluated by comparing their elevations to those recorded by the ALS, after calibrating ASIRAS and KAREN relative to the ALS along runway over-flights (Tables S2 and S3). We aligned each radar waveform to the ice sheet surface and took the mean radar waveforms in 1 km segments along-track to reduce noise. We converted the radar two-way travel time to depth using a depth-density profile from the MAR output. Finally, internal layers present within each ASIRAS profile were traced in an automated way with their chronology resolved relative to the surface (de la Peña et al., 2010;Hawley et al., 2006) (see Supporting Information S1).

Results
We evaluated the modeled firn densities by comparison to those measured in the shallow cores ( Figure 1b). The IMAU-and MAR-FDM densities are highly correlated with the cores (r ¼ 0.93 and 0.89, respectively) and show good overall agreement, with root-mean-square differences (RMSD) of 64 and 104 kg/m 3 , respectively ( Figure 2a). However, we note a spatial pattern in the difference between the in situ and modeled densities. At sites below 2,000 m (T1 to T5), IMAU-FDM overestimates firn density by 10% on average and underestimates firn density by 11% at higher elevation sites. MAR-FDM exhibits a similar bias, with an overestimation of 21% on average at sites below 2,000 m, and an 11% underestimation of firn density at higher elevation sites. The largest departure from the firn cores is recorded for both models at site T1, located at an elevation of 1,698 m in the low percolation zone. At this site, we measured a total of 40 cm of ice from the firn core while IMAUand MAR-FDM simulated a total of 188 and 294 cm of ice in the corresponding firn column, indicating that the firn ice content is overestimated in the lower section of the EGIG line ( Figure S1a).
The ASIRAS radar profiles across the EGIG line ( Figure 1c) show a clear sequence of internal layers starting at an elevation of~2,200 m, which we attribute to melt and refreezing in the percolation zone and to autumn hoar in the dry snow zone. We compare the distribution and sequence of internal layers present within the 2017 ASIRAS profile to isochrones derived from the 2017 IMAU-FDM firn density and chronology. We identify annual isochrones in the IMAU-FDM profile by locating the maximum firn density in data of the same age. The ASIRAS internal layers and IMAU-FDM chronology are highly correlated (r ¼ 0.99, Figure 2b) with a robust dispersion estimate (RDE) of 71 cm. However, compared to the radar layers, IMAU-FDM annual isochrones are consistently found deeper in the firn column, with a systematic bias of 17.1%. The bias between the modeled isochrones and the observed layers accumulates with depth, shifting the distribution of annual layers compared to the radar observations. After adjusting the IMAU-FDM isochrones for this systematic bias, the radar layers and modeled isochrones are well-aligned with an overall RDE of 26 cm and  Figure S5). The close agreement between the sequence and depth of the ASIRAS internal layers with the IMAU-FDM chronology leads us to conclude they are recording the same physical features.
The layers recorded in the ASIRAS profiles show marked interannual variability, with a clear transition after 2012 ( Figure S2). Until 2012, the top two isochrones are the strongest peaks in the radar return. Afterwards in 2014, however, the strongest peaks correspond to the surface layer and the 2012 isochrone. In 2016 and 2017, even though the 2012 isochrone is located deeper in the snowpack, the associated waveform peak remains relatively high-at 33% and 29% of the maximum peak respectively. The strong dielectric contrast of the 2012 melt layer-reducing the energy transmitted to the deeper firn column-is linked to the formation of an ice lens following the intense melt event of that year (Nghiem et al., 2012). This attenuation of the radar backscatter is seen in both the percolation and the dry snow zone. In the percolation zone, we recorded a 1-cm-thick ice lens at a depth of 4.4 m in the firn core collected at site T12 in 2017 (at an elevation of 2,352 m). This layer is aligned with the 2012 isochrone measured at a depth of 4.4 ± 0.2 m in ASIRAS data collected in the same year. At elevations above 2,500 m and prior to the melt event in 2012, the Ku-band radar records significant power (above 1% of the maximum surface return) to a depth of 11.5 ± 1.3 m below the ice sheet surface (Figure 4a). After the melt event, the attenuation in the firn column is mainly driven by the strong reflection at the 2012 melt layer, and as a result, the degree of radar penetration is reduced to 5.3 ± 2.0 m in 2014 and 7.5 ± 2.0 m in 2017. The strong reduction in radar penetration coincides with a peak in the density anomaly recorded by IMAU-FDM of 32.9 kg/m 3 and by MAR-FDM of 63.0 kg/m 3 , which demonstrates that the near-surface firn densities and degree of the radar penetration are linked (Figure 4b).
We observe strong spatial variations in the degree of radar penetration into the near-surface firn along the EGIG line (Figure 3). In all years, few or no internal layers are present in the ASIRAS data in the ablation and percolation zones below~2,200 m, after which their abundance begins to increase with surface elevation as the firn density falls. The number of layers with significant power (above 10% of the maximum surface

Geophysical Research Letters
OTOSAKA ET AL. return) varies in a similar manner to the OCOG retracker width, indicating that strong near-surface scattering has masked scattering at depth. In the ablation and percolation zones, water percolates into the winter snow and new ice lenses or layers are formed each year, preventing the radar signal from penetrating deep into the firn. This process leads to a reduction of the OCOG width, because the main scattered energy is concentrated nearer the ice sheet surface. In all years, the OCOG retracker width and the number of layers show a tendency to increase with elevation and reach maxima at the highest elevation of the transect. However, the maximum number of layers visible and the maximum OCOG retracker width vary from year to year; in 2014, for example, the maxima of both parameters above 2,500 m are more than three times lower than in 2012. Furthermore, compared to 2012, the range of variations in OCOG width is reduced by 86% in 2014 over the same section, which shows that spatial variations in volume scattering are also less prominent after the 2012 melt event.
We evaluate the impact of volume scattering fluctuations on the performance of three alternative radar retrackers-OCOG, TCOG, and TFMRA. Both the TCOG and TFMRA retrackers track the very first peak recorded by the radar altimeter and identify the surface at similar locations with a mean difference of 7.1 cm and standard deviation of 1.9 cm. On the other hand, the OCOG retracker follows the center of gravity of the radar echoes. At elevations above 3,000 m, the scattering horizon is shifted toward the snow surface after the melt event in 2014 compared to 2012, resulting in a 73-cm increase in the altimetry range measurement using the OCOG retracker. We compare the heights from each of the retrackers to the ALS heights (Table S4). Over a total of 2,496 km of flights acquired on five different years, the mean difference between the ALS and OCOG heights was 107 cm, with a standard deviation of 55 cm. By comparison, the TCOG and TFMRA retracked heights are in far better agreement with the ALS data, with mean differences of 14 and 20 cm, and standard deviations of 20 and 21 cm, respectively. Despite the large fluctuations in volume scattering that have occurred along the EGIG line as a consequence of changes in the snowpack structure, the TCOG and TFMRA retracking of radar heights is stable, demonstrating that these algorithms are effective methods of mitigating the impact of radar penetration variations.

Discussion
The two firn densification models we have tested at the EGIG line are able to reproduce the mean density of the shallow firn column with typical differences of 10% (IMAU) and 15% (MAR-FDM) by comparison to in situ measurements. At site T1 in particular, the firn ice content is largely overestimated by the models with more than 59% and 92% of the total length of the firn core simulated as ice by IMAU-and MAR-FDM, respectively, compared to only 13% from the field observations. This indicates that surface melt and refreezing might not be quantified properly in the lower percolation zone of the EGIG line. However, using the density measured from the firn cores or the density outputs from either of these models to convert the radar travel time to depth leads to mean differences within 13 cm. Nevertheless, small biases in modeled firn density do accumulate with depth, offsetting the vertical distribution of annual ice layers (Figure 2b). We also investigate to which extent firn densification models are able to capture the distribution of annual melt layers within the firn column compared to radar-derived layers. This requires a firn model with a fine vertical resolution, such as the IMAU-FDM, to resolve the different layers, which are typically spaced between 30 and 100 cm apart. Although the chronology and spatial distribution of isochrones derived from the IMAU-FDM show good agreement with the internal layers detected by ASIRAS, there is a systematic bias of 17.1% between the two data sets, which we suspect is due to an underestimation of the firn densification rate. This shows that capturing the density variability in the firn column is more challenging than simulating the mean density of the column as suggested by firn model intercomparison studies (Vandecrux et al., 2018;Verjans et al., 2019).

Geophysical Research Letters
We link the sharp reduction in the degree of radar penetration depth after 2012 to the 2012 melt event. Over elevations of 2,500 m, the 2012 isochrone is recorded in the Ku-band altimetry echoes at a depth of 1.5 ± 0.2 m with a significant power of 82% of the maximum surface return and in 2014 at a depth of 3.7 ± 0.2 m with a significant power of 60% (Figure 4a). Over the same part of the transect, the melt layer is also captured by the IMAU-and MAR-FDM with a high density peak after the 2012 summer, which is double the density of the previous summer's peak. In general, when surface melting occurs, the degree of radar penetration into the near-surface firn reduces sharply. The density fluctuations recorded above 2,500 m in 2012 are the largest since 2006 (Figure 4b) and coincide with the fluctuations in radar penetration recorded in the same year. This is shown by the 63% decrease in the OCOG width, the 68% reduction in the number of layers with a high power return, and the 6.2 ± 2.4 m decrease in the radar penetration depth. Such density fluctuations lead to instantaneous upwards (toward the surface) shifts in the distribution of power within the radar echo, followed by gradual downwards (away from the surface) returns to pre-melt conditions. These saw-tooth variations in radar penetration lead to aliased fluctuations in the elevation of the scattering surface, explaining effects that have been highlighted (and corrected for) in analyses of satellite radar altimetry (McMillan et al., 2016;Nilsson et al., 2015;Slater et al., 2019).
We explored two approaches to mitigate the impact of these density fluctuations on surface elevation derived from radar altimetry. First, we found that the application of threshold waveform echo retracking is able to provide estimates that agreed to within 20 cm of coincident laser altimetry. We also examined the sensitivity of higher-frequency Ka-band radar data to volume scattering, as an alternative means of mitigating firn density fluctuations. At higher frequencies, surface scattering is relatively stronger than volume scattering, and we find that the penetration depth is smaller at Ka-band than Ku-band. For example, in 2016 and 2017, although significant power was recorded in ASIRAS Ku-band radar at depths of up to 6.8 ± 0.3 m below the surface along the high altitude section of the EGIG line (300 to 600 km), the corresponding KAREN Ka-band radar penetration was 3.4 ± 0.3 m. We also compared Ka-band surface elevation estimates to those recorded by the ALS (Table S5). Over a total of 782 km of KAREN flight tracks, the mean difference between the KAREN and laser data, when using OCOG, TCOG, and TFMRA retracking algorithms, was 16.0, 12.5, and 15.6 cm, respectively, with standard deviations of 10.8, 10.9, and 11.3 cm. A more detailed assessment of the Ka-band penetration depth was not possible due to the reduced bandwidth of KAREN compared to ASIRAS, which prevents internal layering from being resolved (Figure 1d).

Conclusions
We present an extensive and coincident set of near-surface firn density and airborne radar and laser measurements acquired between 2006 and 2017 along the EGIG line in West Central Greenland. Using these data, we examine the impacts of firn density fluctuations on spatial and temporal variations in the scattering of airborne ASIRAS Ku-band radar waveforms. The largest fluctuations in radar penetration over this period are recorded after 2012 with an abrupt decrease of 6.2 ± 2.4 m in the Ku-band radar penetration. We link this decrease in radar penetration to the density fluctuations associated with the 2012 extreme melt event. As the frequency and extent of extreme melt events is likely to increase in the coming decades (Hall et al., 2013), the effects of fluctuations in radar penetration are an important consideration for satellite radar altimetry. We find that simple methods of threshold retracking are efficient at mitigating this effect on Ku-band airborne radar data and the impact of such events are likely to last for a shorter period of time on Ka-band data due to its reduced penetration depth.

Data Availability Statement
The ASIRAS, KAREN, and ALS raw data are freely available from the European Space Agency at https:// earth.esa.int/web/guest/campaigns. Ice core data, IMAU-FDM, MAR-FDM, ASIRAS, and KAREN profiles are available on PANGAEA at https://doi.pangaea.de/10.1594/PANGAEA.921673.