Mapping Firn Saturation Over Greenland Using NASA’s Soil Moisture Active Passive Satellite

—Mapping the spatial extent of recently identiﬁed englacial hydrological features (i.e., ice slabs and perennial ﬁrn aquifers) formed by meters-thick water-saturated ﬁrn layers over the percolation facies of the Greenland Ice Sheet using L -band microwave radiometry has recently been demonstrated. However, these initial maps are binary, and do not provide a parameter to estimate the spatial variability in the thickness and volumetric fraction of meltwater stored within the ﬁrn pore space. Here, we exploit enhanced-resolution vertical-polarization L -band brightness temperature ( T BV ) imagery (2015–2019) generated using observations collected over Greenland by NASA’s Soil Moisture Active Passive (SMAP) satellite and a simple two-layer L -band brightness temperature model. We map water-saturated ﬁrn layers via a “ﬁrn saturation” parameter, and interpret our results together with ice slab and perennial ﬁrn aquifer spatial extents, estimates of snow ac-cumulationsimulatedviatheRegionalAtmosphericClimateModel(RACMOp2.3),andairborneradarsurveyscollectedviaNASA’sOperationIceBridge(OIB)campaigns.Weﬁndthatvariableﬁrnsaturationparametervaluesaremappedinlowersnowaccumula-tioniceslabareasinwestern,northern,andnortheasternGreen-land,whereﬁrniscolderandwater-saturatedﬁrnlayersseasonallyrefreezeassolid-ice.Higherﬁrnsaturationparametervaluesaremappedinhighersnowaccumulationperennialﬁrnaquiferareasinsoutheastern,southern,andnorthwesternGreenland,whereﬁrnisnearthemeltingpoint,andmeters-thickwater-saturatedﬁrnlayersexist.Ourresultshaveimplicationsforidentifyingexpansiveenglacialreservoirsthatstoresigniﬁcantvolumesofmeltwaterinlocationsthatarevulnerabletomeltwater-inducedhydrofractur-ingandacceleratedoutletglacierﬂowyear-round.


I. INTRODUCTION
R ECENTLY launched low-frequency (1.4 GHz) satellite microwave radiometry missions, such as NASA's Aquarius Mission [1], ESA's Soil Moisture Ocean Salinity Mission (SMOS) mission [2], and NASA's L-band Soil Moisture Active Passive (SMAP) mission [3], enable mapping global sea surface salinity and terrestrial soil moisture from space. However, these missions also have potential for cryospheric applications over Earth's polar ice sheets as a result of large penetration depths that can probe as deep as the underlying bedrock [4], [5] or ocean. While higher frequency (6-37 GHz) satellite microwave radiometry missions have been extensively used to map seasonal surface and near-surface meltwater to depths of a few meters over ice sheets since early in the satellite era, e.g., [6]- [9], low-frequency satellite microwave radiometry missions are capable of mapping seasonal surface and near-surface meltwater, e.g., [10]- [14] as well as subsurface meltwater to depths of tens to hundreds of meters. This capability allows for mapping englacial hydrological features, such as ice slabs and perennial firn aquifers [15], [16], and subglacial hydrological features, such as large subglacial lakes [4], beneath the surface of the Greenland and Antarctic ice sheets.
The objective of this study is to further exploit enhancedresolution vertical-polarization L-band brightness temperature (T B V ) imagery (2015-2019) generated using observations collected over Greenland by the microwave radiometer on NASA's SMAP satellite [17] to map water-saturated firn layers from space. We calculate a "firn saturation" parameter derived from a simple two-layer L-band brightness temperature model [16] using both conventionally processed SMAP T B V imagery generated by a "drop-in-the-bucket" (GRD) algorithm, and enhancedresolution SMAP T B V imagery generated by the radiometer form of the Scatterometer Image Reconstruction (rSIR) algorithm [17]. We then interpret the rSIR mapping of the firn saturation parameter together with ice slab and perennial firn aquifer spatial extents mapped using SMAP rSIR T B V imagery and an empirical algorithm [16], estimates of snow accumulation simulated via the Regional Atmospheric Climate Model (RAC-MOp2.3) [18], and airborne radar surveys collected via NASA's Operation IceBridge (OIB) campaigns [19]. This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ The remainder of this article is organized as follows. A background on previous work is provided in Section II. A description of the methods used in this study is provided in Section III. Section IV provides an interpretation of our results. Finally, Section V concludes this article.

II. MAPPING ENGLACIAL HYDROLOGICAL FEATURES
The percolation facies of the Greenland Ice Sheet are characterized by heterogeneous vertical and lateral percolation of seasonal meltwater into porous snow and firn layers at depth during the melting season, [20]- [23]. Depending on the local climate, subsurface meltwater may refreeze as embedded ice structures [21], [23], or may be stored within englacial reservoirs. During the freezing season, vertically oriented ice pipes form when seasonal meltwater refreezes in percolation channels. Horizontally oriented ice layers and ice lenses form when seasonal meltwater refreezes on semi-impermeable layers deeper within the firn. Shallow water-saturated firn layers often fully refreeze, forming either a network of embedded ice structures consisting of discontinuous ice layers and ice lenses sparsely connected via ice pipes (i.e., spatially coherent melt layers) [24], or dense, low-permeability, solid-ice layers (i.e., ice slabs) [25], [26]. Seasonal meltwater may also be stored intermittently on top of existing ice slabs or other impermeable layers [27], or at the physical limit of vertical meltwater percolation in firn [28], forming perched firn aquifers. In contrast, deep water-saturated firn layers may be perennially stored and can extend from near the ice sheet surface to the firn-ice transition [29], [30]. Such meters-thick water-saturated firn layers that exist year-round and are known as perennial firn aquifers [31]- [33].
The percolation facies boundaries are interannually variable and delineated by the dry snow facies and the wet snow facies as well as the low-elevation (<∼1000 m a.s.l.) ablation facies [23]. The dry snow facies are located in the high-elevation (>∼2500 m a.s.l.) interior and represent the upper boundary of the percolation facies, which is defined by the dry snow line. This region experiences negligible surface melting, and is characterized by snow and firn layers extending from the ice sheet surface to the firn-ice transition. The wet snow facies are located in the lower elevations (<∼1500 m a.s.l.) and represent the lower boundary of the percolation facies. This region typically experiences extensive surface melting. The upper snow layer is completely saturated with seasonal meltwater that refreezes as a superimposed ice layer overlying deeper glacial ice.
Miller et al. [15] first demonstrated the novel use of the L-band microwave radiometer on NASA's SMAP satellite for mapping perennial firn aquifers in the percolation facies of the Greenland Ice Sheet from space [see Fig. 1]. In that study, exponentially decreasing temporal L-band signatures observed in SMAP rSIR T B V time series were correlated with perennial firn aquifer detections [32] identified via airborne ice-penetrating radar surveys collected by NASA's OIB campaigns [19]. An empirical algorithm to map spatial extent was developed by fitting exponentially decreasing temporal L-band signatures to a set of sigmoidal curves derived from the continuous logistic model. The relationship between the radiometric, and therefore the physical, temperature of perennial firn aquifer areas, as compared to other percolation facies areas, forms the basis of the empirical algorithm. That study hypothesized that the dominant control on the exponential rate of decrease in temporal L-band signatures over perennial firn aquifer areas is physical temperature versus depth. However, more recent work by Jezek [34] suggests the complex interplay between physical temperature, the volumetric fraction of meltwater stored within the water-saturated firn layers, latent heat, and volume scattering from embedded ice structures dominates depth-integrated L-band emission over perennial firn aquifer areas during the freezing season.
Miller et al. [16] adapted the empirical algorithm developed in [15] and demonstrated mapping ice slabs and perennial firn aquifers together as a continuous system over the percolation facies of the Greenland Ice Sheet from space [see Fig. 1]. The expanded relationship between the radiometric temperature of each of these englacial hydrological features forms the basis of the adapted empirical algorithm. Analyses of spatiotemporal differences in SMAP rSIR T B V imagery and temporal L-band signatures were used to derive a firn saturation parameter from a simple two-layer L-band brightness temperature model [9], [35]. The firn saturation parameter estimates the spatial variability in the thickness and volumetric fraction of meltwater stored within the firn pore space, and was integrated into the adapted empirical algorithm to map binary spatial extents. That study concluded that ice slabs and perennial firn aquifer areas represent distinct, but related, subfacies within the broader percolation facies that can be defined by differences in snow accumulation, which influences the englacial hydrology and thermal characteristics of the firn at depth.

A. Satellite L-Band Microwave Radiometry
NASA's SMAP satellite was launched on January 31, 2015 and flies at a 685 km altitude, 98.1 • inclination, polar orbit [3]. Operating at an L-band (1.41 GHz) frequency, the microwave radiometer on the SMAP satellite is currently collecting brightness temperature (T B ) observations of the horizontal, vertical, and 3rd and 4th Stokes parameter polarizations over Greenland at a nominal incidence angle of 40 • , and a total radiometric uncertainty of approximately 1.3 K [36]. The antenna spin rate is 14.6 r/min which, when coupled with the along-track motion of the satellite, produces a helical antenna scan pattern on the earth's surface with an along-track spacing of approximately 31 km between antenna rotations.
Conventionally processed SMAP GRD T B imagery generated using a drop-in-the-bucket algorithm provide low-noise levels, but also low spatial resolution. SMAP GRD T B imagery is generated by averaging all satellite observations within a discrete time interval whose center locations fall within the bounds of a particular grid cell. The rSIR algorithm was developed to reconstruct conventionally processed SMAP GRD T B imagery on a finer spatial resolution grid [17], [37]. It exploits the measurement response function (MRF) for each satellite observation, which is a smeared version of the helical antenna scan pattern. Using the overlapping MRFs, the rSIR algorithm reconstructs T B from the spatially filtered, low-resolution sampling provided by the satellite observations. In effect, it generates an MRF-deconvolved T B image. Combining multiple satellite orbital passes within a discrete time interval increases the sampling density, which improves the accuracy and resolution of the reconstruction [17]. Finer spatial resolution SMAP rSIR T B imagery generated using the rSIR algorithm have higher noise levels. Thus, the differing SMAP GRD and rSIR T B imagery provide a tradeoff between effective resolution and noise [17], [37]- [39].
Conventionally processed SMAP GRD and reconstructed rSIR T B imagery have been generated as part of the SMAP Radiometer Twice-Daily rSIR-Enhanced EASE-Grid 2.0 Brightness Temperatures, Version 1, Dataset [17]. The GRD and rSIR algorithms separately combine satellite orbital passes that occur between 8 A.M. and 4 P.M. local time-of-day over Greenland to generate imagery twice daily (i.e., morning and evening satellite orbital pass intervals, respectively). This minimizes fluctuations in the observed T B at high latitudes due to changes in the physical temperature from local daily temperature cycling, and provides improved temporal resolution, permitting resolution of diurnal variations [37]. SMAP GRD and rSIR T B imagery are projected on Northern Hemisphere Equal-Area Scalable Earth Grids (EASE-Grid 2.0) [40] at a 25 and 3.125 km grid cell spacing, respectively. The effective resolution for each grid cell is dependent on the number of satellite observations used in the gridding algorithm, and is coarser than the grid cell spacing. While the effective resolution of conventionally processed SMAP GRD T B imagery posted on a 25 km grid is approximately 30 km, the effective resolution of SMAP rSIR T B imagery posted on a 3.125 km grid is approximately 18 km, an improvement of up to 60%. The tradeoff is up to 70% higher noise [17].
As previously noted, for this study we use both GRD and rSIR T B V imagery over Greenland. As compared to the horizontalpolarization channel, the vertical-polarization channel exhibits decreased sensitivity to changes in the volumetric fraction of meltwater, which is attributed to reflection coefficient differences between channels [15], [16]. The Greenland Ice Mapping Project (GIMP) Land Ice and Ocean Classification Mask, Version 1, Dataset [43] is projected on a Northern Hemisphere EASE-Grid 2.0 at a 25 km GRD and 3.125 km rSIR grid cell spacing, respectively. The derived ice mask includes the Greenland Ice Sheet and its peripheral ice caps. We construct annual ice-masked SMAP GRD and rSIR T B V image time series between April 1, 2015 and March 31, 2019 that alternate morning and evening satellite observations.
Our analysis of temporal L-band signatures over the percolation facies suggests that deeper firn layers are saturated with seasonal meltwater quickly (i.e., time scales of weeks) following melt onset, and remain water-saturated throughout the melting season. We postulate that firn saturation results in a superimposed signal, where the rapidly varying signal induced by near-surface melting and refreezing that occurs as a result of daily temperature cycling is superimposed over the slowly varying signal from deeper water saturated firn layers.

B. Two-Layer L-Band Brightness Temperature Model
The firn saturation parameter [16] is derived from a two-layer L-band brightness temperature model [9], [35] that assumes a base layer underlying a water-saturated firn layer with a given depth and volumetric fraction of meltwater. Each of the layers is uniform. The ice sheet is discretely layered to calculate emission at an oblique incidence angle (1). Emission from the base layer is a function of both the macroscopic roughness and the dielectric properties. It occurs in conjunction with volume scattering at depth, and is locally dependent on englacial hydrological features, including embedded ice structures, spatially coherent melt layers, ice slabs, and perched and perennial firn aquifers. Buried supraglacial lakes may also influence base layer emission [16]. Reflectivity at depth (i.e., at the base layer-water-saturated firn layer interface), and at the ice sheet surface (i.e., at the watersaturated firn layer-air interface) is neglected. The contribution from each layer is individually calculated.
The two-layer brightness temperature model is represented analytically by where T B V,min is the minimum vertical-polarization brightness temperature emitted from the base layer. T B V,max is the maximum vertical-polarization brightness temperature at the ice sheet surface, and represents emissions from the maximum volumetric fraction of meltwater stored within the water-saturated firn layer. T is the physical temperature of the water-saturated firn layer, θ is the transmission angle, κ e is the extinction coefficient, and d is the depth.
We invert (1), and then solve for the firn saturation parameter (ξ) where ξ = κ e d. T B V,max 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. [4] 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) embedded 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 ) [42].
We calculate the firn saturation parameter from the SMAP GRD and rSIR T B V image time series. For each GRD and rSIR grid cell within the ice-masked spatial extent, we calculate T B V,min and T B V,max by smoothing T B V time series using a 14-observation (1-week) temporal averaging window, and then estimating the maximum and minimum values. We note that smoothing T B V time series may mask brief low-intensity surface melting events that occur in the high-elevation percolation facies, and slightly underestimate the mapped spatial extent [16]. T B V,min is chosen prior to T B V,max to avoid detecting changes in T B V values that may occur as a result of changes in ice sheet stratigraphy during or following the melting season (e.g., the formation of embedded ice structures, changes in the depth or volumetric fraction of meltwater stored within the firn pore space of perennial or perched firn aquifers). We set the physical temperature of the water-saturated firn layer to T = 273.15 K, and the transmission angle to θ = 40 • . We then calculate ξ. Following Miller et al. [16], we set the firn saturation parameter threshold to ξ T = 0.1, which corresponds to a uniform water-saturated firn layer with a depth of d = 1 m, and a volumetric fraction of meltwater of m v = 1%. If the calculated firn saturation parameter exceeds the firn saturation parameter threshold (ξ > ξ T ), the grid cell is mapped. The resulting total grid cell spatial extent represents the percolation facies.
We assume that thicker water-saturated firn layers with larger volumetric fractions of meltwater generate higher firn saturation parameter values. However, the observed thickness of the watersaturated firn layer is limited by the penetration depth. Miller et al. [16] estimated theoretical L-band penetration depths for a uniform layer of water-saturated firn using the empirically derived models for the complex dielectric constant described by Tiuri et al. [41], which are calculated as a function of m v . Theoretical L-band penetration depths calculated for a watersaturated firn layer range from between approximately 10 m for small volumetric fractions of meltwater (m v = 1%), and 1 cm for large volumetric fractions of meltwater (m v = 20%). Large volumetric fractions of meltwater stored within the firn pore space results in high reflectivity and attenuation at the watersaturated firn layer-air interface, and a radiometrically cold firn layer.

C. Empirical Algorithm
The empirical algorithm developed by Miller et al. [16] is derived from the two-layer L-band brightness temperature model and the continuous logistic model. The continuous logistic model is a differential equation that models the decrease in physical systems as a function of time using a set of sigmoidal curves. These curves begin at a maximum value with an initial interval of decrease that is approximately exponential. As the function approaches its minimum value, the decrease slows to approximately linear. When the function asymptotically reaches its minimum value, the decrease exponentially tails off and achieves stable values.
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) is known as the sigmoid function. The empirical algorithm uses SMAP rSIR T B V image time series and the sigmoid function. For each rSIR grid cell within the percolation facies spatial extent defined by the two-layer L-band brightness temperature model, we first normalize T B V time series where T B V,min is the minimum vertically polarized brightness temperature, and T B V,max is the maximum vertically polarized brightness temperature. Similar to the firn saturation parameter calculation procedure, we calculate T B V,min and T B V,max by smoothing T B V time series, and then estimating the minimum and maximum values. We then calculate the exponential rate of decrease by iteratively applying the sigmoid fit is the normalized vertical-polarization brightness temperature time series partitioned on the time interval t ∈ [t max , t min ], where t max is the time the function achieves its maximum value, and t min is the time the function achieves its minimum value. The initial normalized vertical-polarization brightness temperature T   Fig. 1(b)] [19]. The calibration parameter intervals used in the empirical algorithm are given in Table I. If the calculated calibration parameters are within the intervals, the rSIR grid cell is converted to a binary parameter to map the total spatial extent of ice slab and perennial firn aquifer areas. For a complete physical description of ice slabs and perennial firn aquifers within the Greenland Ice Sheet, and additional details on the empirical algorithm and the associated uncertainty, see Miller et al. [16].

D. Regional Atmospheric Climate Model
The Regional Atmospheric Climate Model (RACMO2) is a limited-area atmospheric circulation model from which a polar version (p) has been developed to simulate the climate over the polar ice sheets. It includes a multilayer firn model that integrates snow thermodynamics, and allows for surface melting, vertical percolation, retention and refreezing, and run-off [44]. Snow albedo is simulated as a function of snow grain size, cloud cover, solar zenith angle, and impurity concentration [45], and includes drifting, snow erosion, and sublimation [46]. Here, we use simulations from the most recent model version (RACMO2.3p2), which features updates in cloud physics, and an improved representation of impurity concentration, refreezing, albedo, and snowdrift [18]. RACMO2.3p2 simulations are currently generated at 5.5 km horizontal resolution between 1958 and 2019 on a domain covering the Greenland Ice Sheet. The model is forced at its lateral boundaries by six-hourly ERA-40 [47], ERA-Interim [48], and ERA5 [49] reanalyses.
Snow accumulation simulations are projected on a Northern Hemisphere EASE-Grid 2.0 at a 3.125 km grid cell spacing, and then ice sheet-masked using the GIMP Land Ice and Ocean Classification Mask, Version 1, Dataset [43]. We then calculate a climatological snow accumulation estimate between April 1, 2010 and March 31, 2019 for each rSIR grid cell within the percolation facies spatial extent.

E. NASA's Operation IceBridge Campaigns
The CReSIS Accumulation Radar was flown over Greenland on a P-3 aircraft in April and May between 2010 and 2017 as part of NASA's OIB campaigns [19]. The instrument operates at a center frequency of 750 MHz with a bandwidth of 300 MHz, resulting in a range resolution in firn of 0.5 m [50]. The collected observations have an along-track resolution of approximately 30 m with 15 m spacing between traces in the final processed radargrams. At a nominal flight altitude of 500 m above the ice sheet surface, the cross-track resolution varies between 20 m for a smooth surface and 54 m for a rough surface with no appreciable layover.
We construct radargram profiles over southwestern and southeastern Greenland. Each profile was collected prior to melt In the dry snow facies, the penetration depth is approximately 1 km in the cold snow and firn layers, and underlying glacial ice. However, the penetration depth is reduced to between meters and tens of meters in the warmer percolation facies by higher attenuation in the ice column, scattering from embedded ice structures and spatially coherent melt layers, and most significantly, by reflection at the upper surface of meltwater stored within perched and perennial firn aquifers, and buried supraglacial lakes.  1.80 × 10 6 km 2 for the rSIR mapping [see Fig. 2(b)]. Image reconstruction refines the ice sheet-land, ice sheet-ocean, and surrounding ice cap boundaries as well as the percolation faciesdry snow facies boundaries, which reduces Greenland's total initial ice-masked spatial extent by 50 000 km 2 (3%). The total spatial extent of the percolation facies is 5.84 × 10 5 km 2 for the GRD mapping, and 5.75 × 10 5 km 2 for the rSIR mapping. The percolation facies-dry snow facies and wet snow faciespercolation facies boundaries are also refined, which reduces the total percolation facies spatial extent by 9000 km 2 (2%).

A. GRD and rSIR Firn Saturation Parameter Mappings
The firn saturation parameter values exhibit spatial patterns that are consistent in both the GRD and the rSIR mappings [see Figs. 2 and 3]. Higher firn saturation parameter values are mapped over south eastern, southeastern, southern, and northwestern Greenland. Variable firn saturation parameter values are mapped over western, northern, and northeastern Greenland. However, western, northern, and northeastern the spatial patterns are consistent, the firn saturation parameter in the rSIR mapping achieves higher values as compared to the GRD mapping. This is a result of the refinement in the percolation facies boundaries, and the corresponding reduction in mixed emission within the percolation facies spatial extent. Maximum firn saturation parameter values of ξ max = 1.26 are exhibited in the GRD mapping [see Fig. 2(a)]. The mean value isx = 0.5, and the standard deviation is σ = ±0.21. Maximum firn saturation parameter values of ξ max = 2.8 are exhibited in the rSIR mapping [see Fig. 2(b)]. The mean value isx = 0.6, and the standard deviation is σ = ±0. 28. The rSIR mapping also exhibits tighter closed-contour intervals with higher firn saturation parameter values, especially in southeastern [see Fig. 3(c) and (d)] and southern Greenland, and steeper contour gradients. We infer that closed-contour intervals indicate thicker water-saturated firn layers with larger volumetric fractions of meltwater that are relatively uniform over expansive areas. We also infer that steeper contour gradients indicate spatial variability in the thickness and volumetric fraction of meltwater stored within the water-saturated firn layers over limited areas.
Our results indicate that the firn saturation parameter calculated from the SMAP rSIR T B V imagery effectively resolves boundaries and spatial patterns of water-saturated firn layers within the percolation facies that are not effectively resolved in SMAP GRD T B V imagery, and therefore, improves the accuracy of the mapping.  3. We note that exact snow accumulation ranges within the ice slab and perennial firn aquifer spatial extents are challenging to resolve as a result of the mismatch between the effective resolution of SMAP (∼18 km) and the horizontal resolution of RACMOp2.3 (5.5 km). Although image reconstruction improves SMAP s effective resolution, it still spans a spatial extent more than three times the horizontal resolution of RACMOp2.3. Snow accumulation distributions within each of the subfacies are difficult to interpret, especially at the percolation facies-dry snow facies and the wet snow facies-percolation facies boundaries. The spatial patterns are in general agreement. However, the firn saturation parameter is often shifted slightly upslope of higher snow accumulation areas throughout the percolation facies. Fig. 7 shows a scatterplot of the firn saturation parameter values versus snow accumulation. MacFerrin [26] previously V image time series [63] and the empirical algorithm [16]. (c) Estimates of snow accumulation (rainbow logarithmic colorscale) simulated via RACMOp2.3 [18].

B. rSIR Firn Saturation Parameter Mapping Comparisons
reported that snow accumulation is less than 572±32 mm w.e. yr-1 in ice slab areas. Forster et al. [31] reported that snow accumulation exceeds 800 mm w.e. yr-1 in perennial firn aquifer areas. Here, we find an unreasonably wide range of snow accumulation estimates in each of the subfacies. Fig. 7 illustrates the complexities of interpreting grid cell comparisons between satellite-derived mappings of geophysical parameters that often include significant uncertainty in the mapped spatial extents [16] and regional climate model-derived simulations.
The ice slab spatial extent extends over 76 000 km 2 (13%) of the percolation facies, and is primarily mapped over lower snow accumulation areas in central eastern, southern, western, and northern Greenland [see Figs. 4 and 5]. The perennial firn aquifer spatial extent extends over 64 000 km 2 (11%) of the percolation facies, and is primarily mapped over higher snow accumulation areas in southeastern, southern, and northwestern Greenland [see   Fig. 1). (b) 2015-2019 ice slab (cyan shading), perennial firn aquifer (blue shading), and percolation facies (purple shading) spatial extents mapped using the SMAP rSIR T B V image time series [63] and the empirical algorithm [16]. (c) Estimates of snow accumulation (rainbow logarithmic colorscale) simulated via RACMOp2.3 [18].
Bugt glaciers [see Fig. 8]. Steep contour gradients surround these closed-contour intervals, indicating rapid decreases in the firn saturation parameter values over limited areas. Higher firn saturation parameter values are mapped over higher snow accumulation areas in southern and northwestern Greenland. Closedcontour intervals indicate that the firn saturation parameter values range from between ξ = 1.0 and ξ = 1.6, and values generally decrease with increasing northern latitudes. Wider contour gradients surround these closed contour intervals, indicating more gradual decreases in the firn saturation parameter values over larger areas. The lowest firn saturation parameter values are mapped over the lowest snow accumulation areas in central eastern, central western, and northern Greenland. Contour gradients widen with decreasing firn saturation parameter values. Fig. 9 shows the firn saturation parameter values (2015-2019) together with the CReSIS Accumulation Radar southwestern (May 5, 2017) and southeastern (April 22, 2017) profiles collected as part of NASA's OIB campaigns. Surface elevation ranges from between ∼1700 and 2500 m a.s.l in the southwestern profile, and between ∼1500 and 2500 m a.s.l. in the southeastern profile. These radargram profiles display ice sheet stratigraphy in the dry snow facies and the percolation facies just prior to the melting season. Exceptionally bright surface-parallel reflectors that represent spatially coherent melt layers formed during recent extreme [52] and anomalous [51] melting seasons are exhibited throughout the upper firn layers of the radargram profiles [24]. The large dielectric contrast between the spatially coherent melt layers and the overlying, underlying, and interior firn layers results in high-reflectivity at the interfaces. However, electromagnetic energy still propagates through the high-reflectivity spatially coherent melt layers into the deeper firn layers.
The dry snow facies exhibit an alternating sequence of bright and dark surface-parallel reflectors that represent seasonal snow accumulation layers throughout the depth range displayed in both the southwestern and the southeastern profiles. Small density differences between summer and winter snow accumulation results in dielectric contrasts and interface reflectivity [53]. Upward dipping layers transitioning to relatively flat layers in the southwestern profile indicate decreases in snow accumulation toward the ice sheet periphery. Downward dipping layers in the southeastern profile indicate large increases in The higher-elevation percolation facies exhibit a faded alternating sequence of surface parallel reflectors that extend downslope ∼120 km (2450 m a.s.l.) in the southwestern profile, and ∼180 km (2100 m a.s.l.) in the southeastern profile. Firn layers warm further downslope where deep vertical percolation of seasonal meltwater occurs, and the volumetric fraction of embedded ice structures increases. Penetration depth decreases, and surface-parallel reflectors gradually disappear. The firn saturation parameter values increase at the dry snow line. The rate of increase over the southwestern profile is approximately one third the rate of increase over the southeastern profile.
The lower elevation southwestern profile exhibits thick dark surface-parallel regions of low reflectivity that represent ice slabs [see Fig. 9(b)]. The large dielectric contrast between ice slabs and the overlying and underlying firn layers results in high reflectivity at the interfaces. 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. Ice slabs extend downslope ∼50 km (1500-2100 m a.s.l.), range in depth from between ∼1 and 15 m beneath the ice sheet surface, and in thickness from between ∼5 and 20 m. The firn saturation parameter reaches the maximum value (ξ = 0.6) upslope ∼220 km (1700 m a.s.l.) of the ice slabs, and then decreases across the ice slabs toward the percolation facies boundary where meltwater runoff occurs. The lower elevation southeastern profile exhibits a bright lower reflector that

V. DISCUSSION
The simple two-layer L-band brightness temperature model does not adequately describe the physics controlling L-band emission over areas of the percolation facies where heterogeneous vertical and lateral percolation of seasonal meltwater into porous snow and firn layers at depth results in the presence of subsurface meltwater during the melting season. Field measurements collected in and southwestern [21] and southeastern [31] Greenland suggest seasonal meltwater is capable of vertical percolation to depths greater than 10 m beneath the ice sheet surface. When the volumetric fraction of meltwater stored within the water-saturated firn layer is m v = 1%, the theoretical L-band penetration depth is estimated to be approximately 10 m [16].
This estimate is consistent with the L-band penetration depth in damp temperate firn derived from airborne interferometric synthetic aperture radar and laser altimetry observations (10±4 m) [54]. To overcome capillary trapping and initiate vertical percolation, the required volumetric fraction of meltwater stored within the water-saturated firn layer is estimated to be at least m v = 7% [55]. At m v = 7%, the theoretical L-band penetration depth is reduced to less than 1 m [16]. Thus, as observed depth-integrated over SMAP's antenna footprint (i.e. ∼18 km effective resolution), the firn saturation parameter is likely underestimated once seasonal meltwater stored within the water-saturated firn layer exceeds the volumetric fraction required to initiate vertical percolation.
Perennial firn aquifers represent radiometrically cold subsurface meltwater reservoirs that exist at depths of between 1 and 40 m [32] beneath the ice sheet surface during the freezing season [15]. The mean depth is 22 m [32]. The thickness of the water-saturated firn layer ranges from between 4 and 25 m [30], [33]. The mean thickness is 14 m [29]. Field measurements collected in southeastern Greenland suggest that the volumetric fraction of meltwater stored within the pore space of perennial firn aquifers may be as high as 25% [29], which results in high permittivity. ( r ≈ 9 + 1i) that limits the transmission of electromagnetic radiation (10%) to the upper firn layers [15]. L-band emissions upwelling from below the upper surface of meltwater stored within the perennial firn aquifer are extinguished by reflection at the water-saturated layer interface. During the melting season, seasonal meltwater must percolate vertically from the surface to the subsurface to recharge the perennial firn aquifer [27]. The total thickness of the water-saturated firn layer includes the thickness of the seasonal water-saturated firn layer plus the thickness of the perennial firn aquifer. Plausible thicknesses range roughly between 5 and 65 m. The mean thickness is 36 m. Thus, the firn saturation parameter is likely further underestimated over the perennial firn aquifer spatial extent, even at m v = 1%. Similar to Jezek [34], we hypothesize the complex interplay between physical temperature, the volumetric fraction of meltwater stored within the water-saturated firn layers, and latent heat also dominates depth-integrated L-band emission over perennial firn aquifer areas during the melting season. 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 [15]. This subsurface heat source is not present in other areas of the percolation facies. Field measurements collected in southeastern Greenland indicate that firn layers overlaying the perennial firn aquifer are at, or near, the melting point year-round [29]. Our hypothesis is also consistent with temporal L-band signatures observed over the perennial firn aquifer spatial extent, which are radiometrically warm in the absence of seasonal meltwater during the freezing season [see Table I]. Field measurements collected in ice slab areas [25], [26] and in other areas of the percolation facies [21] in southwestern Greenland indicate that firn layers are colder as compared to perennial firn aquifer areas during the freezing season. Temporal L-band signatures observed over the ice slab spatial extent and over other areas of the percolation facies spatial extent are also radiometrically colder as compared to temporal L-band signatures over the perennial firn aquifer spatial extent [see Table I].
While the total thickness and volumetric fraction of meltwater stored within water-saturated firn layers is likely underestimated over much of the perennial firn aquifer spatial extent by the simple two-layer L-band brightness temperature model, the firn saturation parameter effectively identifies expansive englacial reservoirs that store significant volumes of meltwater year-round. Especially interesting are the locations mapped above the fast-flowing (>100 m yr −1 ) Helheim and Køge Bugt glaciers in southeastern Greenland. These large marine-terminating outlet glaciers were the second and fourth largest contributors of ice discharge from the Greenland Ice Sheet to the ocean between 1986 and 2020, respectively [56]. Recent modeling studies suggest that if perennial firn aquifers intersect with crevasse fields and contain a sufficient volume of meltwater, they are capable of initiating meltwater-induced hydrofracturing [57] through the full-thickness (∼1000 m) of the Greenland Ice Sheet [58]. Meltwater-induced hydrofracturing is capable of delivering meltwater stored in englacial reservoirs to the subglacial hydrological system, which may lead to the localized acceleration of outlet glaciers [59], [60], peripheral ice discharge [61], and mass loss to the ocean [62].

VI. CONCLUSION
In this study, we have demonstrated the novel use of the Lband microwave radiometer on NASA's SMAP satellite for mapping firn saturation over the percolation facies of the Greenland Ice Sheet from space. We have derived a firn saturation parameter from a simple two-layer L-band brightness temperature model that can be used to estimate the spatial variability in the thickness and volumetric fraction of meltwater stored within the firn pore space. We have demonstrated that as compared to conventionally processed SMAP GRD T B V imagery, enhanced-resolution SMAP rSIR T B V imagery improves the accuracy of the firn saturation parameter mapping. We have also demonstrated that spatial variability in the firn saturation parameter is correlated with the locations of SMAP-and NASA OIB-derived englacial hydrological features (i.e., ice slabs and perennial firn aquifers) as well as RACMOp2.3-simulated snow accumulation patterns. Further interpretation suggests that the firn saturation parameter likely underestimates the total thickness and volumetric fraction of meltwater stored within the water-saturated firn layers, however, effectively identifies expansive englacial reservoirs in locations that are potentially vulnerable to meltwater-induced hydrofracturing and accelerated outlet glacier flow year-round.
Future work will focus on expanding these results using individual years of enhanced-resolution SMAP rSIR T B V imagery to resolve spatiotemporal variability in the firn saturation parameter over the percolation facies. Given snow accumulation is a key climate parameter defining firn saturation (i.e., it seasonally controls firn pore space) as well as the formation, continued presence, and expansion of ice slabs and perennial firn aquifers, we will attempt to resolve the uncertainties in the estimates of snow accumulation simulated via RACMOp2.3 and derive more accurate estimates in the percolation facies and in each of the subfacies. Comparisons between yearly firn saturation parameter maps and satellite-derived ice flow may provide interesting new insights into the currently unknown influence of perennial firn aquifers on the overall mass balance and stability of the Greenland Ice Sheet.