An empirical algorithm to map perennial ﬁrn aquifers and ice slabs within the Greenland Ice Sheet using satellite L-band microwave radiometry

. Perennial ﬁrn aquifers are subsurface meltwater reservoirs consisting of a meters-thick water-saturated ﬁrn layer that can form on spatial scales as large as tens of kilometers. They have been observed within the percolation facies of glaciated regions experiencing intense seasonal 5 surface melting and high snow accumulation. Widespread perennial ﬁrn aquifers have been identiﬁed within the Greenland Ice Sheet (GrIS) via ﬁeld expeditions, airborne ice-penetrating radar surveys, and satellite microwave sensors. In contrast, ice slabs are nearly continuous ice layers that can 10 also form on spatial scales as large as tens of kilometers as a result of surface and subsurface water-saturated snow and ﬁrn layers sequentially refreezing following and Quantifying the possible rapid expansion of these sub-facies using satellite L-band microwave radiometry has signiﬁcant for understanding ice-sheet-wide variability in hydrology that may drive meltwater-induced hydrofracturing and ice ﬂow as well as 5 high-elevation meltwater runoff that can impact the mass balance and stability of the GrIS.

Abstract. Perennial firn aquifers are subsurface meltwater reservoirs consisting of a meters-thick water-saturated firn layer that can form on spatial scales as large as tens of kilometers. They have been observed within the percolation facies of glaciated regions experiencing intense seasonal 5 surface melting and high snow accumulation. Widespread perennial firn aquifers have been identified within the Greenland Ice Sheet (GrIS) via field expeditions, airborne icepenetrating radar surveys, and satellite microwave sensors. In contrast, ice slabs are nearly continuous ice layers that can 10 also form on spatial scales as large as tens of kilometers as a result of surface and subsurface water-saturated snow and firn layers sequentially refreezing following multiple melting seasons. They have been observed within the percolation facies of glaciated regions experiencing intense seasonal 15 surface melting but in areas where snow accumulation is at least 25 % lower as compared to perennial firn aquifer areas. Widespread ice slabs have recently been identified within the GrIS via field expeditions and airborne ice-penetrating radar surveys, specifically in areas where perennial firn aquifers 20 typically do not form. However, ice slabs have yet to be identified from space. Together, these two ice sheet features represent distinct, but related, sub-facies within the broader percolation facies of the GrIS that can be defined primarily by differences in snow accumulation, which influences the 25 englacial hydrology and thermal characteristics of firn layers at depth.
Here, for the first time, we use enhanced-resolution vertically polarized L-band brightness temperature (T B V ) imagery (2015-2019) generated using observations collected 30 over the GrIS by NASA's Soil Moisture Active Passive (SMAP) satellite to map perennial firn aquifer and ice slab areas together as a continuous englacial hydrological system. We use an empirical algorithm previously developed to map the extent of Greenland's perennial firn aquifers 35 via fitting exponentially decreasing temporal L-band signatures to a set of sigmoidal curves. This algorithm is recalibrated to also map the extent of ice slab areas using airborne ice-penetrating radar surveys collected by NASA's Operation IceBridge (OIB) campaigns (2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017). Our SMAP- 40 derived maps show that between 2015 and 2019, perennial firn aquifer areas extended over 64 000 km 2 , and ice slab areas extended over 76 000 km 2 . Combined together, these sub-facies are the equivalent of 24 % of the percolation facies of the GrIS. As Greenland's climate continues to warm, 45 seasonal surface melting will increase in extent, intensity,

Introduction
The recent launches of several satellite L-band microwave radiometry missions by NASA (Aquarius mission, Le Vine 10 et al., 2007; Soil Moisture Active Passive (SMAP) mission, Entekhabi et al., 2010) and ESA (Soil Moisture and Ocean Salinity (SMOS), Kerr et al., 2001) have provided a new Earth-observation tool capable of detecting meltwater stored tens of meters to kilometers beneath the ice sheet surface. 15 Jezek et al. (2015) recently demonstrated that in the highelevation (3500 m a.s.l.) dry snow facies of the Antarctic Ice Sheet, meltwater stored in subglacial Lake Vostok can be detected as deep as 4 km beneath the ice sheet surface. Subglacial lakes represent radiometrically cold subsurface melt- 20 water reservoirs. Upwelling L-band emission from the radiometrically warm bedrock underlying the subglacial lakes is effectively blocked by high reflectivity and attenuation at the interface between the bedrock and the overlying lake bottom. This results in a lower observed microwave brightness tem- 25 perature (T B ) at the ice sheet surface as compared to other dry snow facies areas where bedrock contributes to L-band emission depth-integrated over the entire ice sheet thickness.
Similar to subglacial lakes, perennial firn aquifers also represent radiometrically cold subsurface meltwater reser- 30 voirs (Miller et al., 2020) consisting of a 4-25 m thick watersaturated firn layer Montgomery et al., 2017;Chu et al., 2018) that can form on spatial scales as large as tens of kilometers . Perennial firn aquifers have been identified via field expeditions (Forster 35 et al., 2014), airborne ice-penetrating radar surveys (Miège et al., 2016), and satellite microwave sensors (Brangers et al., 2020;Miller et al., 2020) in the lower-elevation (< 2000 m a.s.l.) percolation facies of the Greenland Ice Sheet (GrIS) at depths from between 1 and 40 m beneath the ice 40 sheet surface. They exist in areas that experience intense seasonal surface melting and rain (> 650 mm w.e. yr −1 ) during the melting season and high snow accumulation (> 800 mm w.e. yr −1 ) during the freezing season . High snow accumulation in perennial firn aquifer ar- 45 eas thermally insulates water-saturated firn layers from the cold atmosphere, allowing seasonal meltwater to be stored in liquid form year-round if the overlying seasonal snow layer is sufficiently thick (Kuipers Munneke et al., 2014).  estimated that the volumetric fraction of meltwa- 50 ter stored within the pore space of Greenland's perennial firn aquifers just prior to melt onset ranges from between 10 % and 25 %, which limits the upward propagation of electromagnetic energy from greater depths within the ice sheet. Large volumetric fractions of meltwater within the firn pore 55 space result in high reflectivity and attenuation at the interface between water-saturated firn layers and the overlying refrozen firn layers as well as between glacial ice or an impermeable layer and the overlying water-saturated firn layers. Upwelling L-band emission from deeper glacial ice and the 60 underlying bedrock is effectively blocked.
While perennial firn aquifers are radiometrically cold, the slow refreezing of deeper firn layers saturated with large volumetric fractions of meltwater represents a significant source of latent heat that is continuously released throughout the 65 freezing season. Refreezing of seasonal meltwater by the descending winter cold wave (Pfeffer et al., 1991), and the subsequent formation of embedded ice structures (i.e., horizontally oriented ice layers and ice lenses as well as vertically oriented ice pipes; Benson et al., 1960;Humphrey et 70 al., 2012;Harper et al., 2012) within the upper snow and firn layers, represents a secondary source of latent heat. These heat sources help maintain meltwater at depth. Perennial firn aquifer areas are radiometrically warmer than other percolation facies areas where the single source of latent heat is via 75 refreezing of seasonal meltwater. This results in a higher observed T B at the ice sheet surface during the freezing season as compared to other percolation facies areas where seasonal meltwater is fully refrozen and stored exclusively as embedded ice.

80
Recently, mapping the extent of Greenland's perennial firn aquifers from space was demonstrated using satellite L-band microwave radiometry (Miller et al., 2020). Exponentially decreasing temporal L-band signatures observed in enhanced-resolution vertically polarized L-band bright-85 ness temperature (T B V ) imagery (2015-2016) generated using observations collected over the GrIS by the microwave radiometer on NASA's SMAP satellite (Long et al., 2019) were correlated with a single year of perennial firn aquifer detections (Miège et al., 2016). These detections were identified 90 via the Center for Remote Sensing of Ice Sheets (CReSIS) Multi-Channel Coherent Radar Depth Sounder (MCoRDS) flown by NASA's Operation IceBridge (OIB) campaigns (Rodriguez-Morales et al., 2014). An empirical algorithm to map extent was developed by fitting temporal L-band signa-95 tures to a set of sigmoidal curves derived from the continuous logistic model.
The relationship between the radiometric, and thus the physical, temperature of perennial firn aquifer areas, as compared to other percolation facies areas, forms the basis of the 100 empirical algorithm. Miller et al. (2020) hypothesized that the dominant control on the relatively slow exponential rate of T B decrease over perennial firn aquifer areas is physical temperature versus depth. L-band emission from the radiometrically warm upper snow and firn layers decreases dur-105 ing the freezing season as embedded ice structures slowly re-freeze at increased depths below the ice sheet surface. In the percolation facies, refreezing of seasonal meltwater results in the formation of an intricate network of embedded ice structures that are large (10-100 cm long, 10-20 cm wide; Jezek et al., 1994) relative to the L-band wavelength (21 cm). Em-5 bedded ice structures induce strong volume scattering (Rignot et al., 1993;Rignot, 1995) that decreases T B (Zwally, 1977;Swift et al., 1985;Jezek et al., 2018).
Ice slabs are 1-16 m thick nearly continuous ice layers that can form on spatial scales as large as tens of kilometers as 10 a result of surface and subsurface water-saturated snow and firn layers sequentially refreezing following multiple melting seasons (Machguth et al., 2016;MacFerrin et al., 2019). Over time, they become dense low-permeability solid-ice layers overlying deeper permeable firn layers. Ice slabs have been 15 identified via field expeditions and airborne ice-penetrating radar surveys in the lower-elevation (< 2000 m a.s.l.) percolation facies of the GrIS at depths from between 1 and 20 m beneath the ice sheet surface (MacFerrin et al., 2019). They exist in areas that experience intense seasonal surface melt-20 ing and rain (excess melt of 266-573 mm w.e. yr −1 ; see Mac-Ferrin et al., 2019, for a description) during the melting season and lower snow accumulation (< 572±32 mm w.e. yr −1 ) during the freezing season as compared to perennial firn aquifer areas (MacFerrin et al., 2019). Lower snow accumu-25 lation in ice slab areas results in a seasonal snow layer that is insufficiently thick to thermally insulate water-saturated firn layers and seasonal meltwater is instead stored as embedded ice.
Refreezing of seasonal meltwater by the descending win- 30 ter cold wave, and the subsequent formation of ice slabs as well as other embedded ice structures within the upper snow and firn layers, is the single source of latent heat. While ice slab areas are radiometrically warmer than other percolation facies areas with a lower volumetric fraction of embedded 35 ice, they are radiometrically colder than perennial firn aquifer areas. This results in typically higher observed T B at the ice sheet surface during the freezing season in ice slab areas, as compared to other percolation facies areas, but typically lower observed T B as compared to perennial firn aquifer ar-40 eas. Similar to temporal L-band signatures over perennial firn aquifer areas, temporal L-band signatures over ice slab areas are exponentially decreasing during the freezing season; however, the rate of T B decrease is slightly more rapid.
In this study, we exploit the observed sensitivity of L-band 45 emission to variability in the depth-and time-integrated dielectric and geophysical properties of the percolation facies of the GrIS to map perennial firn aquifer and ice slab areas together as a continuous englacial hydrological system using satellite L-band microwave radiometry. 50

Methods
We adapt our previously developed empirical algorithm to map the extent of Greenland's perennial firn aquifers (Miller et al., 2020) using a multi-year calibration technique. We use enhanced-resolution L-band T B V imagery (2015-2019) 55 generated using observations collected over the GrIS by the microwave radiometer on NASA's SMAP satellite (Long et al., 2019) and airborne ice-penetrating radar surveys collected by NASA's OIB campaigns (Rodriguez-Morales et al., 2014). First, we correlate (1) a "firn saturation" parameter de-60 rived from a simple two-layer L-band brightness temperature model; (2) maximum and (3) minimum T B V values; and (4) exponentially decreasing temporal L-band signatures, with 5 years of perennial firn aquifer detections (2010-2014) identified via the CReSIS Accumulation Radar (AR) (Miège et 65 al., 2016) and 3 years of additional detections (2015-2017) more recently identified via MCoRDS (Miller et al., 2020). Next, we extend our empirical algorithm to map the extent of ice slab areas. We correlate the SMAP-derived parameters with 5 years of ice slab detections (2010-2014) recently 70 identified via AR (MacFerrin et al., 2019). Finally, we recalibrate our empirical model to map the extent of perennial firn aquifer and ice slab areas over the percolation facies. Interannual variability in extent is not resolved in this study; however, it will be explored further in future work. 75

SMAP enhanced-resolution L-band T B imagery
The key science objectives of NASA's SMAP mission (https: //smap.jpl.nasa.gov/, last access: 4 January 2022) are to map terrestrial soil moisture and freeze/thaw state over Earth's land surfaces from space. However, the global L-band T B 80 observations collected by the SMAP satellite also have cryospheric applications. Mapping perennial firn aquifer and ice slab areas over Earth's polar ice sheets represents an interesting analog and an innovative extension of the SMAP mission's science objectives. The SMAP satellite was launched 85 on 31 January 2015 and carries a microwave radiometer that operates at an L-band frequency of 1.41 GHz (Enkentabi et al., 2010). It is currently collecting observations of vertically and horizontally polarized T B over Greenland. The surface incidence angle is 40 • , and the radiometric accuracy is ap-90 proximately 1.3 K (Piepmeier et al., 2017).
The Scatterometer Image Reconstruction (SIR) algorithm was originally developed to reconstruct coarse-resolution satellite radar scatterometry imagery on a higher-spatialresolution grid (Long et al., 1993;Early and Long, 2001). 95 The SIR algorithm has been adapted for coarse-resolution satellite microwave radiometry imagery (Long and Daum, 1998;Long and Brodzik, 2016;Long et al., 2019). The microwave radiometer form of the SIR algorithm (rSIR) uses the measurement response function (MRF) for each observa-100 tion, which is a smeared version of the antenna pattern. Using the overlapping MRFs, the rSIR algorithm reconstructs T B from the spatially filtered low-resolution sampling provided by the observations. In effect, it generates an MRFdeconvolved T B image. Combining multiple orbital passes increases the sampling density, which improves both the accuracy and resolution of the SMAP enhanced-resolution T B imagery (Long et al., 2019).
Over Greenland, the rSIR algorithm combines satellite orbital passes that occur between 08:00 and 16:00 local time of day to reconstruct SMAP enhanced-resolution T B imagery twice-daily (i.e., morning and evening orbital pass in-10 terval, respectively). T B imagery is projected on a Northern Hemisphere (NH) Equal-Area Scalable Earth Grid (EASE-Grid 2.0; Brodzik et al., 2012) at a 3.125 km rSIR grid cell spacing (e.g., Fig. 1). The effective resolution for each grid cell is dependent on the number of observations used in 15 the rSIR reconstruction and is coarser than the rSIR grid cell spacing. While the effective resolution of conventionally processed SMAP T B imagery posted on a 25 km grid is approximately 30 km (e.g., Fig. 1a), the effective resolution of SMAP enhanced-resolution T B imagery posted on a 20 3.125 km grid is approximately 18 km (e.g., Fig. 1b), an improvement of 60 % (Long et al., 2019).
As previously noted, for our analysis of the percolation facies we use SMAP enhanced-resolution T B V imagery over the GrIS. Compared to the horizontally polarized channel, 25 the vertically polarized channel exhibits decreased sensitivity to variability in the volumetric fraction of meltwater, which is attributed to reflection coefficient differences between channels (Miller et al., 2020). Using the vertically polarized channel also results in a reduced chi-squared error 30 statistic when fitting T B V time series to the sigmoid function (Sect. 2.3.4). We construct T B V imagery that alternates morning and evening orbital pass observations annually, beginning and ending just prior to melt onset. The Greenland Ice Mapping Project (GIMP) Land Ice and Ocean Classification 35 Mask and Digital Elevation Model  are projected on the NH EASE-Grid 2.0 at a 3.125 km rSIR grid cell spacing. The derived ice mask includes the Greenland Ice Sheet and the peripheral ice caps, including Maniitsoq and Flade Isblink. T B V imagery between 1 April 2015 and 31 40 March 2019 is ice-masked, and an elevation for each rSIR grid cell is calculated.  Miège et al. (2016). Bright lower reflectors that undulate with the local topographic gradient underneath which reflectors are absent in the percolation facies are interpreted as the upper sur-75 face of meltwater stored within perennial firn aquifers (e.g., Fig. 3a). The large dielectric contrast between refrozen and water-saturated firn layers results in high reflectivity at the interface. However, the presence of meltwater increases attenuation, limiting the downward propagation of electromag-80 netic energy through the water-saturated firn layer. The total number of AR derived perennial firn aquifer detections is 325 000, corresponding to a total extent of 98 km 2 . The analysis assumes a smooth surface, which is typical of much of the percolation facies, and a grid cell size of 15 m×20 m. The 85 total number of MCoRDS-derived perennial firn aquifer detections is 142 000, corresponding to a total extent of 80 km 2 . This analysis also assumes a smooth surface and a grid cell size of 14 m × 40 m. The combined total number of grid cells (467 000) and total extent (178 km 2 ) is significantly 90 larger than the total number of MCoRDS-derived grid cells (78 000) and total extent (44 km 2 ) calculated for 2016 (Miller et al., 2020). Perennial firn aquifer detections are mapped in northwestern, southern, and south and central eastern Greenland as well as the Maniitsoq and Flade Isblink ice caps 95 (Figs. 1c and 2a).

Airborne ice-penetrating radar surveys
We project AR-and MCoRDS-derived perennial firn aquifer detections on the NH EASE-Grid 2.0 at an rSIR grid cell spacing of 3.125 km. Each rSIR grid cell has an extent of approximately 10 km 2 . The total number of rSIR grid cells 100 with at least one perennial firn aquifer detection is 800, corresponding to a total extent of 8000 km 2 . However, given the limited AR and MCoRDS grid cell coverage, less than 1 % of the rSIR grid cell extent has airborne ice-penetrating radar survey coverage. As compared to the total number of 105 MCoRDS-derived perennial firn aquifer detections (780) calculated for 2016 (Miller et al., 2020), the total number of rSIR grid cells with at least one detection is only increased multi-year calibration technique. Thick dark surface-parallel regions of low reflectivity in the percolation facies are interpreted as ice slabs (e.g., Fig. 3b). The large dielectric contrast between ice slabs and the overlying and underlying snow and firn layers results in high reflectivity at the interfaces.

10
However, electromagnetic energy is not scattered or absorbed within the homogeneous ice slab; it instead propagates downward through the layer and into the deeper firn layers. The total number of AR-derived ice slab detections is 505 000, corresponding to a total extent of 283 km 2 . Ice slab detections 15 are mapped in western, central and northeastern, and northern Greenland as well as the Flade Isblink Ice Cap (Figs. 1c and 2b).
We project the AR-derived ice slab detections on the NH EASE-Grid 2.0 at an rSIR grid cell spacing of 3.125 km. The 20 total number of rSIR grid cells with at least one ice slab detection is 2000, corresponding to a total extent of 20 000 km 2 . However, less than 2 % of the rSIR grid cell extent has airborne ice-penetrating radar survey coverage.
An advantage of the multi-year calibration technique as 25 compared to the single-coincident year calibration technique (Miller et al., 2020) is that it increases the number of rSIR grid cells that can be assessed. It also provides repeat targets that can account for variability in the depth-and timeintegrated dielectric and geophysical properties that influ-30 ence the radiometric temperature in stable perennial firn aquifer and ice slab areas. Uncertainty is introduced by correlating the SMAP-derived parameters with AR-and MCoRDS-derived detections that are not coincident in time. The multi-year calibration technique assumes the extent of 35 each area remains stable, which is not necessarily the case as climate extremes (Cullather et al., 2020) can influence each of these sub-facies. The assumption of stability neglects boundary transitions in the extent of perennial firn aquifer areas associated with refreezing of shallow water-saturated firn layers, englacial drainage of meltwater into crevasses at the periphery (Poinar et al., 2017(Poinar et al., , 2019, and transient upslope 5 expansion (Montgomery et al., 2017). Once formed, ice slabs are essentially permanent features within the upper snow and firn layers of the percolation facies until they are compressed into glacial ice. However, they may transition into superimposed ice at the lower boundary of ice slab areas or rapidly 10 expand upslope, particularly following extreme melting seasons (MacFerrin et al., 2019). Thus, we simply consider our mapped extent a high-probability area for the preferential formation of each of these sub-facies, with continued presence dependent on seasonal surface melting and snow accu-15 mulation in subsequent years.
Annual perennial firn aquifer and ice slab detections that may introduce significant uncertainty into the multi-year calibration technique include those following the 2010 melting season, which was exceptionally long (Tedesco et al., 2011); 20 the anomalous 2012 melting season, during which seasonal surface melting extended across 99 % of the GrIS (Nghiem et al., 2012); and the 2015 melting season, which was especially intense in western and northern Greenland . Following these extreme melting seasons, sig-25 nificant changes in the dielectric and geophysical properties likely occurred across large portions of the GrIS, including perennial firn aquifer recharging resulting in increases in meltwater volume and decreases in the depth to the upper surface of stored meltwater. The upper snow and firn layers of the dry snow facies and percolation facies were also sat-5 urated with relatively large volumetric fractions of meltwater as compared to the negligible-to-limited volumetric fractions of meltwater that percolates during more typical seasonal surface melting over the GrIS.
Seasonal meltwater was refrozen into spatially coherent 10 melt layers following the 2010 and 2012 melting seasons (Culberg et al., 2021) as well as more recently following the 2015 and 2018 melting seasons identified as part of the temporal L-band signature analysis in this study (Sect. 2.3.1). As compared to ice slabs, which are dense low-permeability 15 solid-ice layers, spatially coherent melt layers are a network of embedded ice structures primarily consisting of discontinuous horizontally oriented ice layers and ice lenses sparsely connected via vertically oriented ice pipes (Culberg et al., 2021). Spatially coherent melt layers are relatively thin 20 (0.2 cm-2 m) and can rapidly form across the high-elevation (up to 3200 m a.s.l.) dry snow facies at depths of less than 1 m beneath the ice sheet surface following a single extreme melting season. They can further merge together into thicker solid-ice layers following multiple extreme melting seasons. 25 Spatially coherent melt layers are exceptionally bright in AR radargrams (e.g., Fig. 3a). The large dielectric contrast between the spatially coherent melt layer and the overlying, underlying, and interior snow and firn layers results in high reflectivity at the interfaces. However, electromagnetic en-30 ergy still propagates downward through the high-reflectivity layer into the deeper firn layers. Culberg et al. (2021) recently demonstrated mapping the extent of spatially coherent melt layers formed following the 2012 melting season (Nghiem et al., 2012) via AR (Figs. 1c and 2). integrated dielectric and geophysical properties of the ice sheet (Ulaby et al., 2014). The most significant geophysical property influencing T B is the volumetric fraction of meltwater within the snow and firn pore space (Mätzler and Hüppi, 1989). During the melting season, the upper snow 10 and firn layers of the percolation facies are saturated with large volumetric fractions of meltwater that percolates vertically into the deeper firn layers (Benson, 1960;Humphrey et al., 2012). Increases in the volumetric fraction of meltwater result in rapid relative increases in the imaginary part 15 of the complex dielectric constant (Tiuiri et al., 1984). This typically increases T B and decreases volume scattering and penetration depth. The L-band penetration depth can rapidly decrease from tens to hundreds of meters to less than a meter, dependent on the local snow and firn conditions. During the 20 freezing season, surface and subsurface water-saturated snow and firn layers and embedded ice structures subsequently refreeze. Decreases in the volumetric fraction of meltwater result in rapid relative decreases in the imaginary part of the complex dielectric constant. This decreases T B and increases 25 volume scattering and penetration depth. The L-band penetration depth increases back to tens to hundreds of meters on variable timescales. We analyze melting and freezing seasons in temporal Lband signatures exhibited in T B V time series over and near 30 the AR-and MCoRDS-derived perennial firn aquifer and ice slab detections projected on the NH EASE-Grid 2.0 ( Fig. 4 and Table 1). We project ice surface temperature observations calculated using thermal infrared brightness temperature collected by the Moderate Resolution Imaging Spectroradiome-35 ter (MODIS) on the Terra and Aqua satellites (Hall et al., 2012) on the NH EASE-Grid 2.0 at a 3.125 km rSIR grid cell spacing. We then derive melt onset and surface freezeup dates for each rSIR grid cell using the methodology described in Miller et al. (2020). We set a threshold of ice sur-40 face temperature > −1 • C for meltwater detection (Nghiem et al., 2012), consistent with the ±1 • C accuracy of the ice surface temperature observations. For temperatures that are close to 0 • C, ice surface temperatures are closely compatible with contemporaneous NOAA near-surface air temperature 45 observations (Shuman et al., 2014). Melt onset and surface freeze-up dates are overlaid on T B V time series to partition the melting and freezing seasons. Melt onset dates typically occur between April and July, and surface freeze-up dates typically occur between July and September. The melting season 50 increases in duration moving downslope from the dry snow facies and ranges from a single day in the highest elevations (> 2500 m) of the percolation facies to 150 d in the ablation facies. Similarly, the freezing season decreases in duration moving downslope and ranges between 215 and 365 d.

55
Over perennial firn aquifer areas (e.g., Fig. 4a (Fountain and Walder, 1998). Minimum T B V (T B V,min ) values remain radiometrically warm during the freezing season as a result of latent heat continuously released by the slow refreezing of the deeper firn layers that are saturated with large volumetric fractions of meltwater (Miller 65 et al., 2020). Temporal L-band signatures exhibit slow exponential decreases and approach, and sometimes achieve, stable T B V values. T B V can decrease by more than 50 K during the freezing season, which represents the descent of the upper surface of stored meltwater by depths of meters to tens of 70 meters beneath the ice sheet surface (Miège et al., 2016).
Over ice slab areas (e.g., Fig. 4b, SMAP Test Site B: 66.8850 • N, 42.7765 • W; 1817 m a.s.l.), T B V,max values are typically radiometrically colder than over perennial firn aquifer areas during the melting season. The presence of 75 dense low-permeability solid-ice layers reduces the snow and firn pore space available to store seasonal meltwater at depth. Meltwater may alternatively run off ice slabs downslope towards the wet snow facies. T B V,min values are also typically radiometrically colder than over perennial firn aquifer areas 80 during the freezing season as a result of the absence of meltwater stored at depth. Temporal L-band signatures exhibit exponential decreases that are slightly more rapid than over perennial firn aquifer areas and often achieve stable T B V values.

85
Over other percolation facies areas (e.g., Fig. 4c, SMAP Test Site C: 66.9024 • N, 44.7528 • W; 2350 m a.s.l.), where seasonal meltwater is fully refrozen and stored exclusively as embedded ice, T B V,max values are typically radiometrically colder than over perennial firn aquifer and ice slab areas 90 during the melting season. T B V,min values are also typically radiometrically cold during the freezing season. Temporal L-band signatures exhibit rapid exponential decreases and achieve stable T B V values. However, over the highest elevations (> 2500 m a.s.l.) of the percolation facies approaching 95 the dry snow line, where seasonal surface melting and the formation of embedded ice structures is limited, T B V,min values remain radiometrically warm during the freezing season. T B V decreases, often step responses exceeding 10 K, are a result of an increase in volume scattering from newly formed 100 embedded ice structures within a spatially coherent melt layer. Temporal L-band signatures that increase several K on timescales of years indicate the burial of spatially coherent melt layers formed following the 2010, 2012, 2015, and 2018 melting seasons by snow accumulation. 105 Exponentially decreasing temporal L-band signatures transition smoothly between perennial firn aquifer, ice slab, and other percolation facies areas -there are no distinct tem-  lation facies, where temporal L-band signatures exhibit rapid increases following melt onset, temporal L-band signatures reverse and exhibit rapid decreases. These reversals are a result of high reflectivity and attenuation at the fully water saturated snow layer and/or at the wet rough superimposed ice-10 air interface. Meltwater runs off superimposed ice downslope towards the ablation facies. T B V,min values remain radiometrically warm during the freezing season. Temporal L-band signatures exhibit rapid increases and achieve stable T B V values.
2.3.2 Two-layer L-band brightness temperature model 15 Based on our analysis of T B V,max and T B V,min in temporal Lband signatures over the percolation facies (Sect. 2.3.1), we derive a firn saturation parameter using a simple two-layer L-band brightness temperature model (Ashcraft and Long, 2006). The firn saturation parameter is similar to the "melt 20 intensity" parameter derived in Hicks and Long (2011) that uses enhanced-resolution vertically polarized Ku-band radar backscatter imagery (2003) collected by the SeaWinds radar scatterometer that was flown in tandem on NASA's Quick Scatterometer (QuikSCAT) satellite (Tsai et al., 2000) and 25 JAXA's Advanced Earth Observing Satellite 2 (ADEOS-II) (Freilich et al., 1994). We use the firn saturation parameter to estimate the maximum seasonal volumetric fraction of meltwater within the saturated upper snow and firn layers of the percolation facies using T B V,max and T B V,min values extracted 30 from T B V time series. We calculate the firn saturation parameter for each rSIR grid cell within the ice-masked extent of the GrIS as part of our adapted empirical algorithm (Sect. 2.3.4).
We assume a base layer underlying a water-saturated firn layer with a given depth and volumetric fraction of meltwa- 35 ter. Each of the layers is homogenous. The ice sheet is discretely layered to calculate T B V at an oblique incidence angle (Eq. 1). Emission from the base layer is a function of both the macroscopic roughness and the dielectric properties of the layer. It occurs in conjunction with volume scattering at 40 depth and is locally dependent on embedded ice structures, spatially coherent melt layers, ice slabs, and perennial firn aquifers. Reflectivity at depth (i.e., at the base layer-watersaturated firn layer interface) and at the ice sheet surface (i.e., at the water-saturated firn layer-air interface) is neglected. 45 The contribution from each layer is individually calculated.
The two-layer L-band brightness temperature model is represented analytically by where T B V,max is the maximum vertically polarized L-band 50 brightness temperature at the ice sheet surface and represents emission from the maximum seasonal volumetric fraction of meltwater stored within the water-saturated firn layer. T B V,min is the minimum vertically polarized L-band brightness temperature emitted from the base layer. T is the physical tem-55 perature of the water-saturated firn layer, θ is the transmission angle, κ e is the extinction coefficient, and d is depth.
We invert Eq. (1) and solve for the firn saturation parameter (ξ ) 60 where ξ = κ e d. The maximum vertically polarized L-band brightness temperature asymptotically approaches the physical temperature of the water-saturated firn layer as the extinction coefficient and the depth of the water-saturated firn layer increases. For simplicity, we follow Jezek et al. (2015) 65 and define the extinction coefficient as the sum of the Raleigh scattering coefficient (κ s ) and the absorption coefficient (κ a ). This assumes scattering from snow grains, which are small (millimeter scale) relative to the L-band wavelength (21 cm), and neglects Mie scattering from large (centimeter scale) em-70 bedded ice structures. However, for water-saturated firn, absorption dominates over scattering, and increases in the extinction coefficient are controlled by the volumetric fraction of meltwater (m v ). We assume that thicker water-saturated firn layers with 75 larger volumetric fractions of meltwater generate higher firn saturation parameter values. However, the thickness of the Figure 5. Theoretical L-band penetration depths for a uniform layer of (a) refrozen and (b) water-saturated firn. Penetration depths (1/(κ s + κ a )) are calculated as a function of the Raleigh scattering coefficient (κ s ; Eq. 8) and the absorption coefficient (κ a ; Eq. 10). The complex dielectric constant is calculated using the empirically derived models described in Tiuri et al. (1984). Refrozen firn penetration depths are calculated as a function of firn density (ρ firn ), and the curves are plotted for snow grain radii (r) set to r = 0.5 mm (upper curve) and r = 4 mm (lower curve). Water-saturated firn penetration depths are calculated as a function of the volumetric fraction of meltwater (m v ), and the curves are plotted for firn density set to ρ firn = 400 kg m −3 (upper curve) and ρ firn = 917 kg m −3 (lower curve). Given the complexity of modeling embedded ice structures, they are excluded from the penetration depth calculation. Increases in the volumetric fraction of embedded ice within the firn will result in an increase in volume scattering, which will decrease and compress the distance between the penetration depth curves for both refrozen and water-saturated firn.
water-saturated firn layer is limited by the L-band penetration depth. Theoretical L-band penetration depths calculated for a water-saturated firn layer range from between 10 m for small volumetric fractions of meltwater (m v < 1 %) and 1 cm for large volumetric fractions of meltwater (m v = 20 %) (Fig. 5). 5 Large volumetric fractions of meltwater result in high reflectivity and attenuation at the water-saturated firn layer-air interface and a radiometrically cold firn layer.

Continuous logistic model
We adapt our previously developed empirical algorithm to 10 map the extent of Greenland's perennial firn aquifers (Miller et al., 2020) to also map the extent of ice slab areas. The empirical algorithm is derived from the continuous logistic model, which is based on a differential equation that models the decrease in physical systems as a function of time using 15 a set of sigmoidal curves. These curves begin at a maximum value with an initial interval of decrease that is approximately exponential. Then, as the function approaches its minimum value, the decrease slows to approximately linear. Finally, as the function asymptotically reaches its minimum value, the 20 decrease exponentially tails off and achieves stable values. We use the continuous logistic model to parametrize the refreezing rate within the water-saturated upper snow and firn layers of the percolation facies using T B V time series that are partitioned using T B V,max and T B V,min values. We calculate the 25 refreezing rate for each rSIR grid cell within the percolation facies extent as part of our adapted empirical algorithm (Sect. 2.3.4).
The continuous logistic model is described by a differential equation known as the logistic equation that has the solution where x o is the function's initial value, ζ is the function's exponential rate of decrease, and t is time. The function x(t) 35 is also known as the sigmoid function. We use the sigmoid function to model the exponentially decreasing temporal Lband signatures observed over the percolation facies as a set of decreasing sigmoidal curves. We first normalize T B V time series for each rSIR grid cell 40 where T B V,min is the minimum vertically polarized L-band brightness temperature, and T B V,max is the maximum vertically polarized L-band brightness temperature. We then apply the sigmoid fit T B V,N (t ∈ [t max , t min ]) is the normalized vertically polarized L-band brightness temperature on the time interval t ∈ [t max , t min ], where t max is the time the function achieves a maximum value, and t min is the time the function achieves a 50 minimum value. The initial normalized vertically polarized L-band brightness temperature (T B V,N (t max )) is the function's maximum value. The final normalized vertically polarized L-band brightness temperature (T B V,N (t min )) is the function's minimum value. The function's exponential rate of decrease 55 represents the refreezing rate parameter (ζ ). An example set of simulated sigmoidal curves is shown in Fig. 6.

SMAP-derived extent mapping
Our adapted empirical algorithm is implemented in two steps: (1) mapping the extent of the percolation facies using the firn saturation parameter derived from the simple two-layer L-band brightness temperature model (Sect. 2.3.2) 5 and (2) mapping the extent of perennial firn aquifer and ice slab areas over the percolation facies using the continuous logistic model (Sect. 2.3.3) we calibrate using airborne icepenetrating radar surveys (Sect. 2.2). Using Eq.
(2), we first set a threshold for the firn saturation 10 parameter (ξ T ) defined by the relationship We calculate the Raleigh scattering coefficient (κ s ) in Eq. (7) using 15 where N d is the particle density, k o is the wave number of the background medium of air, r is the snow grain radius set to r = 2 mm, and ε r is the complex dielectric constant. The particle density is defined by 20 where ρ firn is firn density set to ρ firn = 400 kg m −3 , and ρ ice is ice density set to ρ ice = 917 kg m −3 . Our grain radius and firn density estimates are consistent with measurements within the upper snow and firn layers of the percolation facies of southeastern Greenland at the Helheim Glacier field site 25 (Fig. 2a, blue circle), where in situ perennial firn aquifer measurements have recently been collected (Miller et al., 2017). We calculate the absorption coefficient (κ a ) in Eq. (7) using 30 where I{ } represents the imaginary part. We calculate the complex dielectric constant of the water-saturated firn layer in Eqs. (8) and (10) using the empirically derived models described in Tiuri et al. (1984). We set the volumetric fraction of meltwater to m v = 1 %. We set the depth of the 35 water-saturated firn layer in Eq. (7) to d = 1 m. These values are consistent with typical lower frequency (e.g., 37, 13.4, 19 GHz) passive (e.g., Mote et al., 1995;Abdalati and Steffen, 1997;Ashcraft and Long, 2006) and active (e.g., Hicks and Long, 2011) microwave algorithms used to de-40 tect seasonal surface melting over the GrIS. Using the results of Eqs. (7-10), we calculate the firn saturation parameter threshold to be ξ T = 0.1. The first step in our adapted empirical algorithm is to map the extent of the percolation facies. For each rSIR grid 45 cell within the ice-masked extent of the GrIS, we smooth the corresponding T B V time series using a 14-observation (1 week) moving window. We extract the minimum vertically polarized L-band brightness temperature (T B V,min ) and the maximum vertically polarized L-band brightness temperature (T B V,max ). We set the physical temperature of the watersaturated firn layer to T = 273.15 K and the transmission angle to θ = 40 • . We then calculate the firn saturation parameter (ξ ) using Eq. (2). If the calculated firn saturation parameter exceeds the firn saturation parameter threshold, the rSIR grid cell is converted to a binary parameter to map the total extent of the percolation facies.
We note that smoothing T B V time series will mask brief low-intensity seasonal surface melting that occurs in the 10 high-elevation (> 2500 m) percolation facies, where seasonal meltwater is rapidly refrozen within the colder snow and firn layers (e.g., Fig. 4d). Thus, the calculated firn saturation parameter will not exceed the firn saturation parameter threshold, and these rSIR grid cells are excluded from 15 the algorithm. The exclusion of rSIR grid cells in the highelevation percolation facies is not expected to have a significant impact on our results as our algorithm targets rSIR grid cells in areas that experience intense seasonal surface melting. The exclusion of rSIR grid cells may slightly underesti-20 mate the mapped percolation facies extent.
The second step in our adapted empirical algorithm is to map the extent of perennial firn aquifer and ice slab areas over the percolation facies. For each rSIR grid cell within the mapped percolation facies extent, we normalize 25 the corresponding T B V time series (T B V,N (t)) using Eq. (5). We then extract the initial normalized vertically polarized L-band brightness temperature (T B V,N (t max )) and the final normalized vertically polarized L-band brightness temperature (T B V,N (t min )) and partition T B V,N (t) on the time interval t ∈ [t max , t min ]. We smooth T B V,N (t ∈ [t max , t min ]) using a 56observation (4 week) moving window. The sigmoid fit is then iteratively applied using Eq. (6). Smoothing reduces the chisquared error statistic when fitting T B V,N (t ∈ [t max , t min ]) to the sigmoid function. We fix the initial normalized vertically 35 polarized L-band brightness temperature at T B V,N (t max ) = 0.99, which provides a uniform parameter space in which the refreezing rate parameter (ζ ) can be analyzed. Variability in T B V,N (t max ) is controlled by the volumetric fraction of meltwater within the upper snow and firn layers of the percolation 40 facies and is accounted for in the firn saturation parameter (ξ ), which is analyzed separately. T B V,N (t ∈ [t max , t min ]) values iteratively fit to the sigmoid function converge quickly (i.e., algorithm iterations I ∈ [5,15]), and observations are a good fit (i.e., chi-squared error statistic is χ 2 ∈ [0, 0.1]). 45 Using the SMAP-derived T B V,N (t max ) and T B V,N (t min ), rather than the MODIS-derived initial normalized vertically polarized L-band brightness temperature at the surface freezeup date (T B V,N (t sfu )), and final normalized vertically polarized L-band brightness temperature at the melt onset date 50 (T B V,N (t mo )) that were used in the empirical algorithm described in Miller et al. (2020), has several advantages. The key advantage of this approach is that maps can be generated using T B imagery collected from a single satellite, which simplifies our adapted empirical algorithm. Another 55 advantage is that unlike T B collected at shorter-wavelength thermal infrared frequencies (e.g., MODIS), T B collected at longer-wavelength microwave frequencies (e.g., SMAP) is not sensitive to clouds, which eliminates observational gaps and cloud contamination and provides more accurate time 60 series partitioning and more robust curve fitting.
We calibrate our adapted empirical algorithm using the AR-and MCoRDS-derived perennial firn aquifer and ice slab detections projected on the NH EASE-Grid 2.0. For each rSIR grid cell with at least one detection, we extract the 65 correlated maximum vertically polarized L-band brightness temperature (T B V,max ), the minimum vertically polarized Lband brightness temperature (T B V,min ), the firn saturation parameter (ξ ), and the refreezing rate parameter (ζ ). For each of the extracted calibration parameters, we calculate the stan-70 dard deviation (σ ). Thresholds of ±2σ are set in an attempt to eliminate peripheral rSIR grid cells near the ice sheet edge and near the boundaries of each sub-facies, where L-band emission can be influenced by morphological features, such as crevasses, and superimposed and glacial ice, and spatially 75 integrated with emission from rock, land, the ocean, and adjacent percolation facies and wet snow facies areas. The calibration parameter intervals are given in Table 2. We apply the calibration to each rSIR grid cell within the percolation facies extent. If the extracted calibration parameters are within the 80 intervals, the rSIR grid cell is converted to a binary parameter to map the total extent of each of these sub-facies. Miller et al. (2020) cited significant uncertainty in the SMAP-derived perennial firn aquifer extent as a result of the lack of a distinct temporal L-band signature delineating 85 the boundary between perennial firn aquifer areas and adjacent percolation facies areas. In this study, similar uncertainty exists in the SMAP-derived perennial firn aquifer and ice slab extents. This uncertainty could, at least in part, be a result of the rSIR algorithm. An rSIR grid cell corresponds to 90 the weighted average of T B over SMAP's antenna footprint (Long et al., 2019). The weighting is the grid cell's spatial response function (SRF), which is approximately 18 km (i.e., the effective resolution) in diameter. The SRF is centered on the rSIR grid cell. Since the effective resolution (i.e., the size 95 of the 3 dB contour of the SRF) is greater than the rSIR grid cell spacing, the rSIR grid cell SRF's overlap and the T B values are not statistically independent. This uncertainty, however, could also have a geophysical basis, as it is unlikely that the boundaries between sub-facies as well as between facies 100 are distinct. The thickness of the water-saturated firn layer or ice slab may thin and taper off at the periphery, and subfacies and facies may become spatially scattered and merge together.
The limited extent (AR, 15 m × 20 m; MCoRDS, 14 m × 105 40 m) of the airborne ice-penetrating radar surveys as compared to the rSIR grid cell extent (3.125 km) and the effective resolution of the SMAP enhanced-resolution T B V imagery is also cited in Miller et al. (2020) as a source of uncertainty in the empirical algorithm. In this study, similar uncertainty exists in our adapted empirical algorithm. The total rSIR grid cell extent with airborne ice-penetrating radar survey coverage is less than 2 %. Thus, 98 % of the total rSIR grid cell 5 extent from which the SMAP-derived calibration parameter intervals are extracted is unknown. Calculating the total rSIR grid cell extent where detections are absent along OIB flight lines and statistically integrating this calculation into the multi-year calibration technique may help reduce the un-10 certainty, particularly the significant uncertainty in the interannual variability in extent, which we have yet to resolve. A sensitivity analysis suggests that even small changes in the SMAP-derived calibration parameter intervals (i.e., several K for T B V,min and T B V,max , several tenths of a percentage point 15 for ξ , and several hundredths of a percentage point for ζ ) can result in variability in the mapped extents of hundreds of square kilometers and boundary transitions between perennial firn aquifer and ice slab areas. Thus, the mapped extent of each of these sub-facies should simply be considered an 20 initial result demonstrating the potential of our adapted empirical algorithm for future work.

Results and discussion
The SMAP-derived maximum vertically polarized L-band brightness temperature values generated by our adapted  (Table 1). Firn saturation parameter values range from between ξ = 0.1 and 4.0. Refreezing rate parameter values range from between ζ = −0.09 and −0.01. The observed lower bound (ζ = −0.09) of the refreezing rate parameter is significantly higher than the predicted lower bound 35 (ζ = −1) in our example set of simulated sigmoidal curves (black line, Fig. 6). The SMAP-derived perennial firn aquifer, ice slab, and percolation facies extents are shown in Figs. 7a-9a. The percolation facies extent (5.8×10 5 km 2 ) is mapped at elevations 40 between 500 and 3000 m a.s.l. and extends over 32 % of the GrIS extent (1.8 × 10 6 km 2 ). The perennial firn aquifer extent (64 000 km 2 ) is mapped at elevations between 600 and 2600 m a.s.l. and extends over 11 % of the percolation facies extent and 4 % of the GrIS extent. Predominately high T B V,max , 45 T B V,min , ξ , and ζ values mapped within the perennial firn aquifer extent indicate the widespread presence of thicker water-saturated firn layers with larger volumetric fractions of meltwater that are radiometrically warm during both the melting and freezing seasons and have extended refreezing 50 rates. The ice slab extent (76 000 km 2 ) is mapped at elevations between 800 and 2700 m a.s.l. and extends over 13 % of the percolation facies extent and 4 % of the GrIS extent. As compared to perennial firn aquifer areas, decreased T B V,max , T B V,min , ξ , and ζ values in ice slab areas indicate the presence 55 of thinner water-saturated firn layers with lower volumetric fractions of meltwater that are radiometrically colder and have slightly more rapid refreezing rates. Combined together, the total extent (140 000 km 2 ) is the equivalent of 24 % of the percolation facies extent and 10 % of the GrIS extent. The 60 extents of these sub-facies are generally isolated and somewhat scattered within the percolation facies. However, in several areas in south, south and central eastern, and northern Greenland, the sequential formation of sub-facies and facies (dry snow facies-percolation facies-ice slab-perennial firn 65 aquifer-ablation facies) are mapped. Figures 7b-9b show perennial firn aquifers, ice slabs, and spatially coherent melt layers detected by airborne icepenetrating radar surveys overlaid on the SMAP-derived percolation facies extent. The SMAP-derived perennial firn 70 aquifer extent mapped in southern as well as south and central eastern Greenland is consistent with the AR-and MCoRDS-derived perennial firn aquifer detections. Additional smaller perennial firn aquifer areas are mapped in northern Greenland. The SMAP-derived ice slab extent 75 mapped in southwestern and central eastern Greenland is generally consistent with the spatial patterns of the ARderived ice slab detections but is significantly expanded upslope in each of these areas. In northern Greenland, perennial firn aquifers areas are alternatively mapped, and additional 80 expansive ice slab areas are mapped upslope of perennial firn aquifer areas. Additional smaller ice slab areas are mapped in south and southeastern Greenland. We note that the AR-and MCoRDS-derived perennial firn aquifer and ice slab detections are limited in space and time, particularly in northern 85 Greenland, with a time interval as large as 9 years between the airborne ice-penetrating radar surveys and the SMAP enhanced-resolution T B V imagery we use in our adapted empirical algorithm. In western and northern Greenland, the 2015 melting season was especially intense . And, in northern Greenland, the ablation facies have recently (2010-2019) increased in extent (Noël et al., 2019), and supraglacial lakes have recently (2014-2019) advanced inland (Turton et al., 2021), indicating a possible geophys-5 ical basis for the observed formation, boundary transitions, and expansion. Neither perennial firn aquifer nor ice slab areas are mapped on the Maniitsoq and Flade Isblink ice caps, where spatially integrated L-band emission results in calibration parameter values outside the defined intervals for each of these sub-facies.
Although the AR-derived spatially coherent melt layer detections are often observed to be adjacent to perennial firn aquifer and ice slab areas, these sub-facies were masked in the original airborne ice-penetrating radar survey analysis by 15 Culberg et al. (2021). Spatially coherent melt layers often overlay perennial firn aquifers (e.g., Fig. 3a) and merge with ice slabs (Culberg et al., 2021;Fig. 4).
Shallow buried supraglacial lakes have recently been identified within the percolation facies of western, northern, 20 and north and central eastern Greenland using airborne icepenetrating radar surveys (Koenig et al., 2015) and satellite synthetic aperture radar imagery (Miles et al., 2017;Schröder et al., 2020;Dunmire et al., 2021). These buried supraglacial lakes are within the SMAP-derived perennial firn aquifer and 25 ice slab extents but are not expected to significantly influence L-band emission in these areas for two reasons. (1) As compared to SMAP's 18 km effective resolution, the mean extent of buried supraglacial lakes is limited (less than 1 km 2 ), and they are sparsely distributed in perennial firn aquifer 30 and ice slab areas (Dunmire et al., 2021). (2) Supraglacial lakes form during the melting season as a result of meltwater storage within topographic depressions at the ice sheet surface (Echelmeyer et al., 1991). Similar to subglacial lakes (Jezek et al., 2015) and perennial firn aquifers (Miller et 35 al., 2020), supraglacial lakes represent radiometrically cold near-surface meltwater reservoirs. Upwelling L-band emission from deeper firn layers, superimposed and/or glacial ice, and the underlying bedrock are effectively blocked by high reflectivity and attenuation at the interface between the lake bottom and the underlying impermeable layer. This re-5 sults in low observed T B at the upper surface of meltwater stored within supraglacial lakes. During the freezing season, the upper surface of meltwater refreezes and forms a partial or solid-ice cap that is sometimes buried by snow accumulation (Koenig et al., 2015). Airborne ice-penetrating radar surveys in April and May between 2009 and 2012 suggest the mean depth to the upper surface of meltwater stored within buried supraglacial lakes is approximately 2 m (Koenig et al., 2015). Over buried supraglacial lakes, L-band emission from the refreezing partial or solid-ice cap, which is smooth rela- 15 tive to the L-band wavelength (21 cm), likely induces surface scattering. As a result, T B V decreases over buried supraglacial lakes are likely negligible. Thus, over SMAP's 18 km effec-tive resolution, we postulate water-saturated firn layers dominate L-band emission over the percolation facies of the GrIS. 20 The SMAP-derived perennial firn aquifer extent (64 000 km 2 ) generated by our adapted empirical algorithm and the multi-year calibration technique (2015-2019) is consistent with the extent (66 000 km 2 ) generated by the previously developed empirical algorithm and the single-25 coincident year calibration technique (2016) described in Miller et al. (2020). The SMAP-derived perennial firn aquifer extent is generally consistent with previous C-band (5.3 GHz) satellite-radar-scatterometer-derived perennial firn aquifer extents mapped using the Advanced SCATterometer   (Nghiem et al., 2012) in which significant changes in the dielectric and geophysical properties that influence radar backscatter likely occurred. The unreasonably expansive (i.e., more than twice the mean) mapped 10 extent is a result of ASCAT's shallow (several meters) C-band penetration depth (Jezek et al., 1994) and the simple threshold-based algorithm, which was not calibrated for an extreme melting season that included saturation of the upper snow and firn layers of the dry snow facies and percolation 15 facies with relatively large volumetric fractions of meltwater (Miller et al., 2019). Water-saturated firn layers had extended refreezing rates; however, seasonal meltwater was not stored at depth. Widespread spatially coherent melt layers were alternatively formed in many of the mapped areas (Culberg 20 et al., 2021). The SMAP-derived ice slab extent (76 000 km 2 ) is also consistent with previous AR-derived ice slab extents (2010-2014, 64 800-69 400 km 2 ; MacFerrin et al., 2019).
Although we simply consider our mapped extents a highprobability area for preferential formation, the maps gener-25 ated by our adapted empirical algorithm and the multi-year calibration technique for individual years suggest there is reasonable interannual variability in perennial firn aquifer and ice slab extents (Table 3). Our results demonstrate sensitivity to variability in the depth-and time-integrated dielec-30 tric and geophysical properties of the percolation facies that influence the radiometric temperature, even during the 2015 melting season .

Implications
Seasonal surface melting over the GrIS has increased in extent, intensity, and duration since early in the satellite era (Steffen et al., 2004;Tedesco et al., 2008Tedesco et al., , 2011Tedesco et al., , 2016Nghiem et al., 2012;Tedesco and Fettweis, 2020;Cullather et al., 2020). Consistent with recent seasonal surface melting trends, meltwater runoff has accelerated to become the dominant mass loss mechanism over the GrIS (van den Broeke et al., 2016). Meltwater storage in both solid (i.e., embedded ice structures, including ice slabs, and spatially coherent 10 melt layers) and liquid (i.e., perennial firn aquifers) form can buffer meltwater runoff in the percolation facies and delay its eventual release into the ocean . However, significant uncertainty remains in meltwater runoff estimates as a result of the lack of knowledge of heterogeneous 15 infiltration and refreezing processes within the snow and firn layers (Pfeffer and Humphrey, 1996) as well as the depths to which meltwater can descend beneath the ice sheet surface . If the increasing seasonal surface melting trend contin-20 ues (Franco et al., 2013;Noël et al., 2021), perennial firn aquifer formation and expansion may increase the possibility of crevasse deepening via meltwater-induced hydrofracturing (Alley et al., 2005;van der Veen, 2007), especially if crevasse fields expand into perennial firn aquifer areas 25 as a result of accelerated ice flow . Meltwater-induced hydrofracturing is an important component of supraglacial lake drainage during the melting season (Das et al., 2008;Stevens et al., 2015), leading to at least temporary localized accelerated ice flow velocities (Zwally 30 et al., 2002;Joughin et al., 2013;Moon et al., 2014) as well as ice discharge from outlet glaciers (Chudley et al., 2019) and mass balance changes (Joughin et al., 2008). Perennial firn aquifers may also support meltwater-induced hydrofracturing, even during the freezing season (Poinar et al., 2017, 35 2019).
The formation and expansion of ice slabs reduces permeability within the upper snow and firn layers and facilitates lateral meltwater flow with minimum vertical percolation into the deeper firn layers, thereby enhancing melt- 40 water runoff and mass loss at the periphery (Machguth et al., 2016;MacFerrin et al., 2019). Lateral meltwater flow across ice layers overlying deeper permeable firn layers was first postulated by Müller (1962). The theory was then further developed by Pfeffer et al. (1991) as an end-member 45 case for meltwater runoff in the percolation facies, with the other end-member case being lateral meltwater flow across superimposed ice. Lateral meltwater flow and high-elevation (1850 m a.s.l.) meltwater runoff across ice slabs in the percolation facies were first observed in visible satellite imagery 50 collected by the NASA-USGS Landsat 7 mission during the 2012 melting season (Machguth et al., 2016).
Spatially coherent melt layers represent a recently identified refreezing mechanism in the dry snow facies (Nghiem et al., 2005;Culberg et al., 2021). Similar to ice slabs, the 55 formation and expansion of spatially coherent melt layers reduce the pore space within the upper snow and firn layers. They can also limit meltwater flow with minimum vertical percolation into the deeper firn layers, thereby potentially preconditioning the dry snow facies for the formation 60 of ice slabs and enhanced meltwater runoff from significantly higher elevations on accelerated timescales. If spatially coherent melt layers merge with ice slabs upslope of perennial firn aquifers areas, they may also simultaneously accelerate both meltwater runoff and meltwater-induced hydrofractur-65 ing during extreme melting seasons. The formation of spatially coherent melt layers overlying deeper perennial firn aquifers may result in the formation of shallow perched firn aquifers (Culberg et al., 2021) or may terminate gravitydriven meltwater drainage and seasonal recharging (Foun-70 tain and Walder, 1998), which may eventually completely refreeze stored meltwater into ice slabs or decameters-thick solid-ice layers overlying deeper glacial ice.

Summary and future work
In this study, for the first time, we have demonstrated the 75 novel use of the L-band microwave radiometer on NASA's SMAP satellite for mapping perennial firn aquifers and ice slabs together as a continuous englacial hydrological system over the percolation facies of the GrIS. We have adapted our previously developed empirical algorithm (Miller et al.,80 2020) by expanding our analysis of spatiotemporal differences in SMAP enhanced-resolution T B V imagery and temporal L-band signatures. We have used this analysis to derive a firn saturation parameter from a simple two-layer Lband brightness temperature model. And we have used the 85 firn saturation parameter to map the extent of the percolation facies. We have found that by correlating maximum and minimum T B V values, the firn saturation parameter, and the refreezing rate parameter with perennial firn aquifer and ice slab detections identified via the CReSIS AR and MCoRDS 90 instruments flown by NASA's OIB campaigns, we can calibrate our previously developed empirical algorithm (Miller et al., 2020) to map plausible extents.
We note that significant uncertainty exists in the mapped extents as a result of (1) correlating the SMAP-derived parameters with airborne ice-penetrating radar detections that are not coincident in time; (2) the lack of a distinct temporal L-band signature delineating the boundary between perennial firn aquifer areas, ice slabs areas, and adjacent percolation facies areas; and (3) the limited extent of the airborne ice-penetrating radar detections as compared to the rSIR grid cell extent and the effective resolution of the SMAP enhanced-resolution T B V imagery. 10 Miller et al. (2020) normalized SMAP enhancedresolution T B V time series and converted the exponential rate of T B V decrease over perennial firn aquifer areas to a binary parameter to map extent. In this study, we have converted the SMAP-derived parameters to binary parameters to map 15 the extent of both perennial firn aquifer and ice slab areas. Moreover, we have included additional analysis of the spatiotemporal differences in maximum and minimum T B V values, the firn saturation parameter, and the refreezing rate parameter. We have shown that spatiotemporal differences in 20 the SMAP-derived parameters are consistent with our assumption of spatiotemporal differences in the englacial hydrology and thermal characteristics of firn layers at depth. Future work will focus on simulating the temporal L-band signatures observed over perennial firn aquifer and ice slab 25 areas for a wide range of geophysical properties. Additionally, we will simulate the distinct temporal L-band signatures observed over spatially coherent melt layers and explore mapping the extent. Combining multi-layer depth-integrated L-band brightness temperature models (e.g., Jezek et al.,30 2015) that include embedded ice structure parametrizations (e.g., Jezek et al., 2018) with models of depth-dependent geophysical parameters can lead to an improved understanding of the extremely complex and poorly described physics controlling L-band emission over the percolation facies. The de- 35 velopment of more sophisticated empirical algorithms that incorporate multi-layer depth-integrated L-band brightness temperature models that are constrained by in situ measurements can help reduce the significant uncertainty in the current mapped extents and provide more accurate boundary de-40 lineation that can be used to further quantify interannual variability. https://doi.org/10.5067/B8X58MQBFUPA (last access: 1 April 2021; Howat, 2017), and the Digital Elevation Model, Version 1, is available at https://nsidc.org/data/nsidc-0645/versions/1 (last access: 1 April 2021; Howat et al., 2015). The coast-line data are available from GSHHG -A Global Self-55 consistent, Hierarchical, High-resolution Geography Database, https://doi.org/10.1029/96JB00104 (last access: 1 April 2021; Wessel and Smith, 1996). Ice surface temperature imagery (2015-2019) have been produced as part of the Multilayer Greenland Ice Surface Temperature, Surface Albedo, 60 and Water Vapor from MODIS V001 and are available at https://doi.org/10.5067/7THUWT9NMPDK (last access: 1 April 2021; Hall and DiGirolamo, 2019). OIB AR-and MCoRDSderived perennial firn aquifer detections (2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)   Author contributions. JZM initiated the study, adapted the empirical model, performed the analyses, and wrote the manuscript. RC processed and interpreted the airborne ice-penetrating radar surveys. All authors participated in discussions and reviewed manuscript drafts.

85
Competing interests. The contact author has declared that neither they nor their co-authors have any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

90
Acknowledgements. The authors thank Mike MacFerrin and the two anonymous reviewers for their constructive comments. We acknowledge the use of data from CReSIS generated with support from the University of Kansas, NASA Operation IceBridge (grant no. NNX16AH54G), NSF (grant nos. ACI-1443054, OPP-1739003, 95  Review statement. This paper was edited by Lars Kaleschke and reviewed by Michael MacFerrin and two anonymous referees.