Hydrologic Properties of a Highly Permeable Firn Aquifer in the Wilkins Ice Shelf, Antarctica

We present measurements of the density, hydraulic conductivity, and specific discharge of a widespread firn aquifer in Antarctica, within the Wilkins Ice Shelf. At the field site, the aquifer is 16.2 m thick, starting at 13.4 m from the snow surface and transitioning from water‐saturated firn to ice at 29.6 m. Hydraulic conductivity derived from slug tests show a geometric mean value of 1.4 ± 1.2 × 10−4 m s−1, equivalent to permeability of 2.6 ± 2.2 × 10−11 m2. A borehole dilution test indicates an average specific discharge value of 1.9 ± 2.8 × 10−6 m s−1. Ground‐penetrating radar profiles and a groundwater flow model show the aquifer is draining laterally into a large nearby rift. Our findings indicate that the firn aquifer in the vicinity of the field site is likely not in a steady state and its presence likely contributed to past ice shelf instability.


Introduction
Ice shelves, extensions of ice sheets and glaciers that have thinned sufficiently to become afloat on the ocean, are prevalent around the Antarctic Ice Sheet (AIS) and play a role in restraining ice-sheet discharge into the ocean (Siegert et al., 2019). More than 80% of Antarctica's ice discharge is released through ice shelf outflow and basal melting, making them an integral control of the mass balance of the AIS (Pritchard et al., 2012;Rignot et al., 2013). Observations show that there are extensive surface hydrologic systems and meltwater storage on ice shelves in Antarctica Kingslake et al., 2017;Lenaerts et al., 2017;Phillips, 1998) that can potentially accelerate their disintegration (Banwell et al., 2013;Scambos, 2004;Scambos et al., 2000). However, the meltwater volume, residence time, storage characteristics, and lateral/ vertical movement remain not well understood despite being critical to estimate impacts on mass balance, ice dynamics, and sea-level changes (e.g., Bell et al., 2018;Lenaerts et al., 2019;Smith et al., 2020).
Firn aquifers, well documented in mountain glaciers (e.g., Fountain & Walder, 1998) and more recently discovered in Greenland and Svalbard (Christianson et al., 2015;Forster et al., 2014), form when meltwater produced at the surface percolates into the firn and fills the available pore space above the firn-ice transition without refreezing during winter. Firn aquifers are located where there is sufficient pore space volume in the firn column for meltwater to be stored and high accumulation, which provides insulation that permits the saturated subsurface layer to remain at 0°C Kuipers Munneke et al., 2014). In Greenland, these conditions occur in the southeast, south, and northwest sectors where accumulation rates are~1-5 meters water equivalent (m w.e.) per year and melt rates are >650 mm WEyr −1 (Bell et al., 2018;Montgomery et al., 2020;Noël et al., 2018). For the AIS, areas with similar surface melt and snow accumulation signatures are rare, though recent modeling studies show widespread perennial firn aquifers on the Wilkins Ice Shelf (WIS) and elsewhere on the Antarctic Peninsula (van Wessem et al., 2016(van Wessem et al., , 2020. These aquifers can contribute to sea-level rise if connected to the ocean by slowly draining into crevasses (Koenig et al., 2014) and are especially important to understand on ice shelves where meltwater storage is likely a precursor to hydrofracture and ice shelf break-ups (Bell et al., 2018;Scambos et al., 2000).
In December 2018, we conducted fieldwork on the WIS using a combination of borehole drilling, hydrological tests, and ground-penetrating radar (GPR) profiles. We used hydrological and geophysical measurements combined with a groundwater flow model to quantify lateral water flow and assess the hydrologic balance of the aquifer. Further, these results can be used in future studies to examine the impact of firn aquifers on ice shelf stability.

Study Site
Our field site is located on the southwestern portion of the WIS on the Antarctic Peninsula ( Figure 1a). Measurements were taken~50 km from the edge of the ice shelf (−70.80S, −71.71W) from 3-12 December 2018. Average (2010-2017) annual accumulation and melt rates simulated from two regional climate models (MAR and RACMO2) at the field site are 590-800 and 250-425 mm WEyr −1 (Figure 1c and Table S1 in the supporting information) Datta et al., 2018Datta et al., , 2019Noël et al., 2018). A 17 km long rift in the ice shelf was observed~5 km from our field site ( Figure 1b).
Analysis of NASA Operation IceBridge 2014 radio-echo-sounding profiles collected over the WIS using the Multichannel Coherent Radar Depth Sounder (MCoRDS; CReSIS, 2020) indicated a bright reflector in the upper firn similar to high-amplitude reflections associated with firn aquifers for the Greenland Ice Sheet Miège et al., 2016). Based on the separation of the snow surface reflection and the bright subsurface return, we anticipated the top of the firn aquifer to be~13 m below the surface ( Figure S1; Studinger, 2014). MCoRDS processing steps to retrieve aquifer extent and depth to water are discussed in detail in Miège et al. (2016) and Brangers et al. (2020).

Borehole Drilling and Firn/Ice Cores
We used a custom-built lightweight electrothermal drill to drill three boreholes located~1 m apart, reaching depths of 14, 20, and 35 m to estimate density, stratigraphy, and use for other hydrological tests. Drill specifications are in the supplemental material of Miller et al. (2018). To determine gravimetric density, core sections (average diameter of 55 mm) were weighed and measured immediately after collection. The average uncertainties of these density measurements ranged from 9-11%, similar to previous studies (Text S5; Conger & McClung, 2009;Fornasini, 2008). We anticipate some meltwater to have drained during extraction and prior to weighing but cannot quantify the added uncertainty to the density measurements. The stratigraphy was recorded, and two main facies were identified: firn and ice lenses (Figures 2a and S2). No section of the retrieved core had temperatures substantially different from 0°C within the aquifer (accuracy~±0.5°C).

Hydraulic Conductivity Estimate
Slug tests are used to measure the hydraulic conductivity in the saturated zone of an aquifer. Before the slug tests, we placed a pressure transducer (HOBO Onset© U20-001-02) with an operational range between 0 and 400 kPa and a maximum error of 3 cm of water at the bottom of the borehole, which allowed us to record the water-level changes ( Figure S3a). To perform a slug test, we measured the time it took to displace a small volume of water by inserting/removing a solid cylinder (1,187.7 cm 3 ) into the water-filled part of the borehole and measuring the response time for the water level to return to its initial level ( Figure S3c). We use present measurements from six slug tests (with the least noise) from 9-11 December 2018 in the 20 m borehole.

Borehole Dilution Test
We used a saline dilution method in the 35 m borehole to estimate discharge through the aquifer and locate the base of the flow zone. This method is established for traditional groundwater studies (Pitrak et al., 2007) and was recently adapted for estimating water discharge within firn aquifers in the Greenland Ice Sheet (Miller et al., 2018). The method consists of injecting a saline solution of a known conductivity into the borehole and monitoring the conductivity variation over time throughout the borehole to estimate horizontal flow through permeable firn. A description of the method is in the supporting information (Text S1 and Figure S4). After injecting and homogenizing the salt solution in the borehole, vertical conductivity profiles were recorded with a vertical spacing of 30 cm in the water-filled part of the borehole: every 30 min for the first four profiles, followed by two profiles taken an hour apart, then two more profiles taken 2 hours apart ( Figure 2b). One additional profile was taken the following morning. The total duration of the experiment was~17 hours. This borehole dilution test also allows us to estimate specific discharge and linear velocity with depth (Figures 2b and 2c).

Geophysical Surveys and Water-Table Elevation Estimates
We used GPR to survey the top~50 m of the firn/ice column and identify the water-table depth around the drilling site and rift. We also surveyed~35 km of surface topography using a GPS receiver. Limited postprocessing was required because the water table was easily identified as a high-reflective high-amplitude signal. Postprocessing included shifting the traces vertically to align the surface with the first break position of the signal, geolocation, adjustment for surface elevation, and removing the mean trace for the profiles from each trace within the profiles to enhance layering details. We converted from two-way-travel time to depth using mean density from the cores above the water table and a relationship linking permittivity and density (Kovacs et al., 1995). Details of the GPR method are in the supporting information (Text S2).
We also conducted an 8-day continuous GPS survey (1-s epochs) at the drill site (3-11 December 2018) to determine the short-term ice flow vector and tidal range ( Figure S5). Ice flow speed was 92.8 m d −1 with a near-zero diurnal variation. The daily range of vertical ice motion, roughly equal to tidal motion, was 0.8 to 1.3 m.

Flow-Rate Modeling
We use SEEP2D (Jones, 1999), a 2-D finite element flow model within the Groundwater Modeling System package, to simulate water flow within the aquifer, based on the results of the slug tests, GPR, and firn core measurements to determine if the aquifer is in steady state and constrain recharge estimates. SEEP2D solves the steady-state (i.e., all recharge that reaches the aquifer is discharged into the rift) groundwater flow equation (Cherry & Freeze, 1979) that is based on mass balance and utilizes Darcy's Law calculated as where Q is groundwater flow (m 3 s −1 ), A is cross-sectional area (m 2 ), K is the hydraulic conductivity (m s −1 ), and ∂H/∂x is the slope of the water table. Details of the SEEP2D model setup and assumptions are in the supporting information (Text S3).
We ran SEEP2D simulations with the following recharge scenarios: Recharge scenarios were chosen to represent a wide range of possible climatic conditions on the ice shelf. We compare these steady-state output conditions (water table slope and aquifer thickness) to our observations to determine if these scenarios are plausible. One test of plausibility is to compare the amount of refrozen meltwater in the simulation to the column fraction of ice lenses in the firn core density profile.

Firn Core Characteristics
The deepest firn core was collected at 35 m, below the estimated firn-ice transition of~29.6 m (when the specific discharge reaches 0; Figure 2c). The average firn density above the aquifer was 650 kg m −3 and within the aquifer was 850 kg m −3 (Figure 2a). Two shallower firn cores were collected 1 m away from the deepest firn core and show small-scale spatial variability in density above the aquifer with average densities of 627 and 669 kg m −3 ( Figure S2), which agree within 9-11% uncertainty. The three firn cores also show variability in ice lenses, with an ice fraction above the aquifer varying from 17.9% to 20.3% and an average of 18.1%. We measured the depth to water table, which ranged from 13.39-13.46 m, averaging 13.43 m, proving a homogeneous water table depth at the site. We compare this measured depth with the GPR depth converted from two-way-travel time and find a good agreement (±20 cm).

Hydraulic Conductivity
We derived values of hydraulic conductivity from slug tests ranging from 1.0 × 10 −4 to 1.7 × 10 −4 with a geometric mean of 1.4 ± 1.2 × 10 −4 m s −1 . These correspond to permeabilities ranging from 1.8 × 10 −11 to 3.1 × 10 −11 m 2 and averaging 2.6 ± 2.2 × 10 −11 m 2 . Water level change values ranged from 10 to 15 cm ( Figure S3). The slug tests indicate a highly permeable aquifer (similar to unconsolidated sand) with similar hydraulic conductivity values to that of the firn aquifer found in Southeast Greenland (Miller et al., 2017).

Borehole Dilution Test
A time series of conductivity profiles from the 35 m borehole after the saline solution was injected is shown in Figure 2b. We note that this conductivity represents a vertical profile of salinity and not hydraulic conductivity. The background conductivity is below detection until~30 m, where it increased to 75 μS cm −1 . After the conductivity reached 200 μS cm −1 in the water-filled borehole (Text S1), we observed a gradual decrease in the conductivity above 20 m until background levels were reached. Below 20 m, the conductivity did not reach background levels during our experiment time. The decrease in conductivity over time in the profiles indicates lateral water flow, which dominates the freshening process in the borehole. We consider diffusion rates to be negligible in this process as they are~30 times smaller than the inferred advection of water. The decrease in dilution rate with depth (and therefore lateral flow) is primarily due to decreasing porosity. The profile reaches a pore close-off at~30 m, eliminating any dilution or lateral flow (Figure 2a).
The vertical specific discharge (Text S4) profile derived from the salt dilution indicates where and the rate at which water flows laterally in the borehole profile into connected pores in the surrounding aquifer ( Figure 2c). The bottom depth of the specific discharge profile where flow ceases (29.6 m) agrees with the bottom depth of the aquifer we found through coring (~29 m). The average specific discharge is 1.9 × 10 −6 m s −1 with a maximum value of 1.6 × 10 −5 m s −1 at the top of the aquifer due to high porosity 10.1029/2020GL089552
We also calculated the average linear velocity profiles using uniform porosities ranging from 0.1 to 0.3 (Figure 2d and Text S4; Koenig et al., 2014). The density measurements in the region just above the aquifer suggest that the porosity is near the high end of this range at the top of the aquifer (Figure 2a). The average linear water velocity ranged from 0.6-1.7 m d −1 depending on the porosity. Our maximum linear velocity value was 14.1 m d −1 assuming a porosity of 0.1. The calculated values of linear water velocity are substantially larger than the measured ice motion (0.25 m d −1 ); therefore, the water flow is faster than the surrounding ice flow.

GPR Survey of the Rift
The combined GPR-GPS radar profile from the field site across the southern terminus of the ice shelf rift (~5 km) indicates the surface of the ice shelf is initially~17 m above the WGS84 ellipsoid height and rises to~20 m within a few hundred meters of the rift location (3.8 km along the profile) before tipping downward toward the rift (Figures 1b and 3a). This profile shape is typical of rifts in ice shelves (e.g., Fricker et al., 2005). A survey of Landsat imagery (www.earthexplorer.usgs.gov) showed that the rift appeared in late 2009, shortly after a series of major disintegration events on the WIS (Humbert et al., 2010). Rifts (breaks through the entire ice shelf) generally occur when tensile stresses exceed the tensile strength of the ice plate (Braun et al., 2009;Scambos et al., 2009). The rift was filled with snow at the southern end.
The profile reveals firn and ice lenses parallel to the ice shelf surface (Figure 3a) and a continuous strong reflector at the same depth (±1 m) as the aquifer's upper surface at the core site. No layering is observed below the bright layer, consistent with the radar attenuation expected for an aquifer layer. As the profile approaches the rift, it shows the inferred aquifer surface deepening and crossing the firn ice layering. The inferred aquifer surface intersects the rift at an ellipsoidal height of about 0 m, close to the local sea level (Figure 3b). Decreasing aquifer height and consequently decreasing hydraulic head toward the rift strongly implies the lateral movement of water toward the rift and drainage through its sidewalls.

Groundwater Flow Modeling
To quantitatively evaluate that water is flowing into the rift (thereby removing some of the annual recharge to the aquifer) and constrain the parameters in which the aquifer system would be in steady state, we use the SEEP2D groundwater flow model. Table S2 shows the results of the SEEP2D modeling.
First, we present two extreme scenarios of a steady-state aquifer where all surface melt, ranging from 250 to 425 mmWE yr −1 depending on the RCM used, recharges the aquifer. The ∂H/∂x (slope) and A (thickness) parameters required for the high recharge scenario contribution would be physically impossible compared to our observations because the water table of the aquifer is not 26-46 m (Table S2) above the firn-ice transition (i.e., we see no ponding). Further, the presence of refrozen melt layers in the upper firn suggests that the aquifer is not being recharged by all surface meltwater.
We also examine a recharge scenario where 125 mm WEyr −1 (50%) of surface melt recharges the aquifer, using RACMO2 melt input. The other half of the melt input refreezes or densifies the snowpack, which could explain the ice layers in the stratigraphy above the aquifer. However, steady state is only reached if the hydraulic conductivity is twice the observed value, the ∂H/∂x has a 13 m gradient over 3,800 m distance, or if the aquifer is twice as thick. Major dynamical or structural changes would have to occur in a short distance for this scenario to be plausible.
Our final scenarios present low meltwater recharge values ranging from 42.5 to 62.5 mmWE yr −1 , consequently leading to 75-90% of the melt input refreezing in the firn above. With this small amount of melt recharging to the aquifer, our model output matches closely with field measurements (Figure 3b). However, the ice fraction above the aquifer measured at our study site is too small to represent this amount of refrozen meltwater per year (i.e., no thick ice layers).

Discussion and Conclusions
We provide the first measurements of hydraulic conductivity, density, and specific discharge of a firn aquifer on an ice shelf in Antarctica. We also show that the water in this firn aquifer at this field site likely is discharging into a nearby rift. These measurements were performed at only one location but allow us some insight into the meltwater storage in the WIS. Further field experiments spatially distributed (i.e., seismic refraction, magnetic resonance, hydrological tests, and firn stratigraphy and density profiles) would be necessary to fully quantify the overall structure and age of the aquifer, various drainage divides, and the magnitude of the total discharge into the rift.
The average hydraulic conductivity value for the WIS aquifer (from slug tests, 1.4 × 10 −4 ± 1.2 m s −1 ) is of the same magnitude but potentially below what was measured in Greenland (2.7 × 10 −4 ± 1.6 m s −1 ). In turn, the specific discharge from dilution tests (1.9 × 10 −6 ms −1 , σ = 2.8 × 10 −6 m s −1 ) is also less than the value (4.3 × 10 −6 m s −1 , σ = 2.5 × 10 −6 m s −1 ) from Greenland. This results from the smaller surface slope relative to the Helheim study area (~0°vs. 0.8°; Forster et al., 2014) and also suggests that there is less recharge on the ice shelf. The lower recharge rates can be explained by the lower melt rates compared to the Helheim Glacier region. In fact, both the modeled WIS accumulation and melt rates are less than half of the values in the southeastern Greenland region (590-800 and 250-425 [ Figure 1c] vs. 1,400-1,650 and 730 mm WEyr −1 ; Miège et al., 2016). Miller et al. (2017Miller et al. ( , 2018 found that the Greenland firn aquifer was highly permeable and had evidence of direct meltwater flow. This meltwater likely flowed into nearby crevasses and possibly down to the bed of the ice sheet into the ocean (Poinar et al., 2017). Similarly, we find that the WIS firn aquifer is highly permeable and is likely discharging into the nearby rift.
Our hydrological modeling results indicate that while plausible, it is unlikely that the WIS aquifer is in steady state with all meltwater discharged to a nearby rift. The high recharge scenario was deemed implausible because hydraulic gradients rose above the ice shelf surface. The medium recharge scenario provided a more plausible ratio of recharge to refreeze and approximately accounts for the ice fraction of 18% observed in the core above the aquifer. To achieve steady state, this scenario either requires a water table slope or hydraulic conductivity that is two times greater than our observed values. Since we had to extrapolate our single hydraulic conductivity measurement value throughout the entire system, the variability of K in the system could explain the discrepancy. It is also possible that refreezing at the base of the aquifer accounts for some of the unexplained water mass. The low recharge scenarios both match the observed water table (Figure 3b), though this implies that 75-90% of the meltwater input densifies or is refrozen in the firn column above the aquifer, which does not agree with the observed ice fraction above the aquifer (Figure 2a). However, meltwater is required to bring the firn to an isothermal state above the aquifer to allow recharge after the winter season, which varies annually based on seasonal snow thickness and temperature, and could explain some excess meltwater (Miller et al., 2020).
We assume that our study site on the WIS firn aquifer is slowly draining into a large rift that was formed in 2009. Earlier observations using airborne radar show a highly reflective high-amplitude internal reflector near the surface (1966 and1975 surveys;Operation IceBridge [OIB] surveys in 2014 shown in Figure 1) absent of any observed widespread surface flooding, even in the highest-melt years (Braun et al., 2009;Scambos et al., 2000;Vaughan et al., 1993). The presence and thickness of an aquifer are controlled primarily by the summer meltwater flux and the annual net snowfall and rainfall and to a lesser extent by freezing at the base of the aquifer column . Without surface flooding or complete refreeze seen from satellite imagery for several decades (Barrand et al., 2013;Braun et al., 2009;Scambos et al., 2000;Vaughan et al., 1993), the WIS has apparently maintained a near-balanced system with respect to water table height and therefore snowfall and recharge conditions. However, the specific dynamics of the past WIS aquifer are unknown, including if it was in steady state and has always had cracks in the ice shelf with which to evacuate water.
Firn aquifers offer sufficient water storage capacity to contribute to ice shelf disintegration events, as evidenced by the partial break-ups of the northern and northwestern WIS in 1993WIS in , 1998WIS in , and 2008WIS in -2009, where the bright reflector from airborne radar was present (Braun et al., 2009;Humbert & Braun, 2008;Scambos et al., 2000). Our study supports this conclusion, identifying the WIS aquifer as having high permeability and the drainage of the aquifer into the adjacent rift. However, the stability of an aquifer-bearing ice shelf will be dependent on the volume of the aquifer and relationship to the lateral flow of water within the aquifer. Discharge from firn aquifers has been modeled to be able to cause hydrofracturing to the bed of ice sheets with enough inflow of meltwater (Poinar et al., 2017). If a series of low melt years or high snowfall years decreased the relative height of the top of the aquifer column, the potential for hydrofracture is greatly reduced. The year-round availability of water at depth allows for enhanced fracturing whenever stresses change to favor tensile extensions or loss of compression, even in the winter (Scambos et al., 2009). With complete saturation of the vertical firn column, the hydrostatic head for hydrofracture is at a maximum value and has the potential to cause destabilization of ice shelves and iceberg calving through hydrofracturing Pollard et al., 2015;Scambos et al., 2003Scambos et al., , 2009). Continued and improved monitoring of firn aquifers would lead to a better understanding of the role these aquifers play in ice-shelf disintegration.

Data Availability Statement
All data used in this manuscript can be found at https://www.usap-dc.org/view/dataset/601390.