Integrated Borehole, Radar, and Seismic Velocity Analysis Reveals Dynamic Spatial Variations Within a Firn Aquifer in Southeast Greenland

Perennial water storage in firn aquifers has been observed within the lower percolation zone of the southeast Greenland ice sheet. Spatially distributed seismic and radar observations, made ~50 km upstream of the Helheim Glacier terminus, reveal spatial variations of seismic velocity within a firn aquifer. From 1.65 to 1.8 km elevation, shear‐wave velocity (Vs) is 1,290 ± 180 m/s in the unsaturated firn, decreasing below the water table (~15 m depth) to 1,130 ± 250 m/s. Below 1.65 km elevation, Vs in the saturated firn is 1,270 ± 220 m/s. The compressional‐to‐shear velocity ratio decreases in the downstream saturated zone, from 2.30 ± 0.54 to 2.01 ± 0.46, closer to its value for pure ice (2.00). Consistent with colocated firn cores, these results imply an increasing concentration of ice in the downstream sites, reducing the porosity and storage potential of the firn likely caused by episodic melt and freeze during the evolution of the aquifer.


Introduction
Warming of polar ice sheets changes the hydrological structure of surface firn and ice. Distribution and routing of surface meltwater to the base of an ice sheet can alter the sub-and en-glacial hydrologic systems and affect ice flow (e.g., Bartholomew et al., 2010). In particular, englacial water storage in firn aquifers, defined here as a storage of liquid water above the firn-ice transition (Fountain & Walder, 1998), can have contrasting implications on the influence of surface melting on ice sheet dynamics. Firn aquifers could either temporarily buffer the delivery of water to the oceans Parry et al., 2007) or alternatively cause sudden release events when storage capacity is reached (Koenig et al., 2014). In Greenland, extensive firn aquifers have been observed in high-accumulation, high-melt (24-50 cm/yr; Miller et al., 2020) regions in the southeast in the last decade (e.g., Forster et al., 2014;Koenig et al., 2014;Miège et al., 2016). Meltwater stored within a firn aquifer is either drained from the system, though fracture/crevasse networks, or refreezes into ice. Refreezing densifies the firn column (Hubbard et al., 2016) and has been implicated as a precursor for the catastrophic collapse of Antarctic ice shelves . Densification also complicates the relationship between surface elevation change and changing ice mass, since an increase in surface density can cause a decrease in elevation (Reeh et al., 2005). ©2020. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

10.1029/2020GL089335
Key Points: • Characterizing the storage potential of firn aquifers is important for predicting water budgets in a warming world • Multiple geophysical methods are combined to evaluate the elastic properties and porosity of the firn aquifer • Seismic velocities detect an increase in ice content in the aquifer downslope, reducing the porosity and storage potential of the firn Supporting Information: • Supporting Information S1 • Table S1 • Table S2 • Table S3 • Table S4 Correspondence to: S. F. Killingbeck, eespr@leeds.ac.uk Therefore, knowledge of the geographic extent, storage capacity, water flow, and their temporal variations will enable an improved understanding of firn aquifers' impacts on glacial hydrology (Poinar et al., 2019), which, in turn, will be critical for accurately predicting the contribution of meltwater discharge and changing ice mass to global sea level rise.
Geophysical analysis is useful for quantifying firn hydrology without time-intensive drilling, providing the spatial distribution of firn properties with estimates of geophysical property uncertainty. Such quantitative Figure 1. (a) Study location,~50 km upstream of the Helheim Glacier terminus. Sites 1-12 (excluding 9) are seismic survey locations, with boreholes (FA16 4, 5, and 6) drilled in August 2016 at Sites 7, 8, and 12. The map hillshade and elevation contours are obtained from the ArcticDEM created by the Polar Geospatial Center from DigitalGlobe, Inc. imagery (Porter et al., 2018). (b) Elevation (above sea level) profile along the flow line and aquifer structure. (c) Seismic acquisition geometry used throughout. (d) Rayleigh wave dispersion curves for each site and (e) their implied depth sampling (Gazetas, 1982).
analysis can be improved, in terms of depth sensitivity and resolution, where independent geophysical observations are integrated into a joint analysis and combined with ground truth from cores or other point samples of aquifer thickness and hydraulic property measurements.
In 2011, a firn aquifer in the upstream region of Helheim Glacier (southeast Greenland) was studied with ice coring and ground-penetrating radar (GPR) methods Koenig et al., 2014). Airborne radar data have been used to measure the depth and extent of the aquifer (e.g., Chu et al., 2018;Miège et al., 2016), but the utility of radar sampling beneath the aquifer is limited because liquid water increases radar attenuation rates. Magnetic resonance sounding (MRS) methods have been used to estimate the volume of water stored within the aquifer; however, vertical resolution and the aquifer base were poorly resolved (e.g., Legchenko et al., 2018aLegchenko et al., , 2018b. Point samples, provided by boreholes, directly measure the density structure of the aquifer and its detailed hydrological properties (Miller et al., 2018), but these do not provide the spatial coverage of geophysical surveying.
Seismic refraction profiling has been used to obtain the compressional (P-) wave velocity (Vp) structure of the firn and ice to improve the estimate of the aquifer-ice transition depth (e.g., Montgomery et al., 2017). However, thus far, there is no constraint of the shear (S-) wave velocity (Vs) distribution in a firn aquifer.
Vs is related to the shear rigidity of a material, and the ratio Vp:Vs is a powerful indicator of consolidation and fluid content (Lee, 2002). Combined analysis of Vp and Vs could, therefore, improve discrimination between the saturated and unsaturated zones of the firn aquifer. The seismic data described in Montgomery et al. (2017) feature strong Rayleigh waves, a form of seismic energy often termed as "ground roll" that travels along the ground surface. While Rayleigh waves are typically considered noise in seismic refraction/reflection analysis, this wave component is useful to inform the Vs structure of the firn and ice (e.g., Diez et al., 2014;Young Kim et al., 2010).
In this paper, we present an integrated geophysical analysis of the firn aquifer on Helheim Glacier by deriving Vs from seismic Rayleigh waves, leveraging independent Vp, density, and thickness constraints using the "MuLTI" algorithm (Killingbeck et al., 2018) and estimating porosity of the aquifer using the derived Vp:Vs ratios and Biot-Gassmann theory (Biot, 1941;Gassmann, 1951). We use the results to infer the spatial variability of seismic velocity and porosity within the aquifer. Our approach is a novel complement to the conventional geophysical analysis of firn aquifers and provides valuable constraints for modeling their formation.

Field Site and Data Acquisition
GPR, seismic, and firn cores were acquired along a flow line of Helheim Glacier,~50 km upstream of its terminus, between July and August 2016 ( Figure 1a) over a known firn aquifer Miège et al., 2016). Active-source seismic data were acquired at 11 sites (numbered 1-12, excluding 9 where no seismic was acquired in 2016 and for consistency with Montgomery et al., 2017) along the flow line, and boreholes were drilled at Sites 7, 8, and 12. The water table depth was estimated as 10-20 m below the surface ( Figure 1b) from colocated GPR data (Miège et al., 2016). The seismic acquisition is detailed in Montgomery et al., (2017), and their seismic refraction analyses suggest the aquifer base (firn-ice transition) to be~28 m below the surface ( Figure 1b). Borehole dilution tests (Miller et al., 2018) show that water flows through the aquifer with an average specific discharge of 4.3 × 10 −6 m/s, decreasing to 0 m/s at the aquifer base. Miller et al. (2018) also show firn cores that record the density structure of the aquifer.

Rayleigh Wave Analysis
Multichannel analysis of surface waves (MASW; Park et al., 1999) is employed to estimate the Vs structure of the aquifer. Vs can be linked to the shear modulus (or stiffness), μ, and density, ρ, of a material by the following: Since the shear modulus of air and fluids is 0, the replacement of air by water in the saturated zone of the aquifer produces no change in stiffness; however, the same replacement increases the bulk density of the . The mean of the prior and its standard deviation are shown by the red and red-dashed lines, respectively; the mode solutions are shown by the black lines. The gray shaded area highlights data outside the reliable depth range estimated (Gazetas (1982).

Geophysical Research Letters
firn; hence, Vs would be expected to reduce if stiffness remained constant. Conversely, the presence of ice lenses in the aquifer increases its stiffness, therefore leading to an increase in Vs.
The zero-offset shot gathers, at both ends of Arrays 1 and 2 (Figure 1c), were transformed into the frequency-phase velocity domain, where the dispersive pattern of the Rayleigh waves was sampled (Park et al., 1999). The dispersion curves for each array are similar; hence, a homogenous layered subsurface was assumed for the underlying firn such that a single one-dimensional dispersion curve could be analyzed for each site. The dispersion curves suggest that shear-wave phase velocity increases downstream towards lower elevations (Figure 1d), while it decreases with depth through the aquifer (Figure 1e). The estimated depth range resolved by the Rayleigh wave dispersion curves is~3-40 m, based on Gazetas' (1982) approximation that the depth sampled is one-third of the wavelength used.

Inverting for Vs
Rayleigh wave dispersion curves can be inverted to provide vertical Vs profiles. Our analysis is based on a probabilistic Bayesian tool, "MuLTI" (Killingbeck et al., 2018), that jointly inverts Rayleigh wave dispersion curves with layer-depth information into a depth-dependent posterior distribution of Vs. Here we extend the algorithm into a variant "MuLTI III" that additionally allows the inversion to be constrained using depth-dependent density and Vp profiles, thus integrating all available geophysical observations. The data used for our analysis include Using a transdimensional Markov chain Monte Carlo method, we create a large ensemble drawn from the joint posterior distribution of the three parameters (Vs, Vp, ρ), with more likely models sampled with greater frequency. The posterior distribution is the product of the likelihood and the prior information. Here, the likelihood is based upon fitting theoretical dispersion curves to those observed, assuming a normal distribution with a mean given by the forward model calculated using the Geopsy algorithm (Wathelet et al., 2020) and standard deviation σ. The GPR and refraction data determine the vertical extent of three internal layers: unsaturated firn, aquifer, and ice. These data are also used to define the uniformly distributed prior distributions for Vs in each of the layers. The prior depth-dependent distributions for Vp and ρ are determined using seismic refraction (3) and density data (4), respectively, assuming normal distributions centered on their profile values with standard deviation σ Vp and σ ρ . The model is discretized using k piecewise constant functions based on Voronoi cells, where k is unknown but follows a distribution also sought by the method. Where acceptable in terms of likelihood, the method always chooses smooth (low k) models over rough models (high k). Figure 2 shows an example of the surface wave analysis and output products from MuLTI III for Site 7, colocated with borehole measurements from FA16-4 (see supporting information). We obtain complete marginal probability density functions (PDFs) of Vs, Vp, and ρ for which we show the mode (highest probability) value and uncertainty as one standard deviation for each depth (Figure 2c). The posterior distributions of Vp and ρ are comparable to their prior distributions, showing that the data used in the joint inversion are consistent with datas ets (3) and (4) individually. MuLTI III was applied to all sites; where boreholes were absent, the average density of the 3 internal layers was calculated from borehole measurements and applied with an uncertainty of 10%, typical of the uncertainty in the firn cores and permitting variability within the model space.

Profiles of Vs Distributions With Depth
The most likely Vs solutions are derived for all sites (Figure 3a). Average Vs values for the unsaturated firn above the water table (Layer 1), saturated aquifer (Layer 2), and underlying ice (Layer 3) are shown in Figure 3c. For Sites 1-7, located at higher elevations (>1.65 km) within the aquifer system, the average Vs   (1-7) and >1.65 km elevation (10-12). Layer 1 averages are calculated for unsaturated firn at depths >10 m; therefore, site 12 is not included in this average. All error bars represent one standard deviation.

10.1029/2020GL089335
Geophysical Research Letters lenses at Site 12 (Figure 4a), where the clear ice content increased by 9% and 12% compared to Site 7 and 8, respectively, calculated from 0 to 38 m depth. At Site 12, where the water table depth is only 10 m, Vs is 950 ± 120 m/s in the unsaturated firn and 1,280 ± 170 m/s in the underlying saturated layer.
Posterior distributions of the ratio of Vp:Vs were also output from MuLTI III. The most likely solutions are plotted for all sites (Figure 3b) and average Vp:Vs values derived for each layer (Figure 3d). Gagnon et al. (1988) show Vp:Vs ratios of a single ice crystal are 2.14 and 1.97, respectively, for measurements along and orthogonal to the c-axis. When observed from seismic data, Vp:Vs of 2 is typically recorded for standard ice in the upper ice column (Wittlinger & Farra, 2012); hence, this value is used as a reference ratio for ice in this study. Figures 3b and 3d show the average Vp:Vs of the saturated firn aquifer is 2.30 ± 0.54 in the upslope Sites 1-7 and 2.01 ± 0.46 in Sites 10-12, suggesting that porosity is likely lower in the downstream sites, with the Vp:Vs ratio closely approximating that of ice.  (Miller et al., 2018), showing specific discharge (SD), density, stratigraphy, depth to water table, aquifer base (Montgomery et al., 2017), Vp (Montgomery et al., 2017), and Vs (this study). The flow zone was defined in Miller et al. (2018) by positive specific discharge. (b) Average porosity estimation for all sites using the Biot-Gassmann theory, including water content estimates from Montgomery et al. (2017). Error bars are based on the upper and lower Vp:Vs estimates. Biot (1941) and Gassmann (1951) related the bulk Vp:Vs of a material, Vp:Vs bulk , to the velocity ratio of its matrix, Vp:Vs matrix , material total porosity (Ø), and the pore-filling interstitial fluid (here liquid water). The change in Vp:Vs bulk observed in a porous fluid-filled medium is

Porosity Estimation
As we obtain Vp:Vs bulk from our study and set Vp:Vs matrix ¼ 2.0 (Gagnon et al., 1988;Wittlinger & Farra, 2012), equation 2 can be rearranged to yield porosity (Figure 4b). In some cases, our derived Vp: Vs bulk < Vp:Vs matrix , resulting in negative porosities, an artifact of Vp:Vs bulk approaching the velocity of solid ice, and thus extending to lower Vp:Vs within uncertainty bounds. Furthermore, equation 2 does not account for differential pressure and compaction (as modified in Lee, 2002) and assumes that pore fluid does not affect the shear modulus of the sediment. Also, equation 2 does not account for velocity anisotropy effects, which may be present through the aquifer. Therefore, further work is required to address these limitations. Here, we report an underestimation for the downslope sites, 10-12, where the likely increase in ice content, particularly at the base (e.g., Site 12; Figure 4a), stiffens the material, hence decreasing Vp:Vs, thus biasing our Biot-Gassmann approach to lower porosities for these sites. Notably, Site 10 implies a negative porosity (−0.04 ± 0.22), indicating poorly constrained Vp:Vs data and potentially having an effective medium property that cannot be described by a Biot-Gassmann approach.
Nonetheless, there is a decrease in aquifer porosity downstream, below 1.65 km elevation, with average porosity 0.11 ± 0.22 in the upslope sites (1-7) and 0.01 ± 0.24 in the downslope sites (10-12), again consistent with the appearance of additional clear ice in cores (Figure 4a). Our porosity estimates generally agree, within error margins, with those of Montgomery et al. (2017), reported as upper-bound estimates of water content, which use the relative change in Vp in the aquifer verses the underlying ice.

Discussion and Conclusions
This study integrates multiple geophysical observations to evaluate the physical properties of a firn aquifer. Upstream, the infiltration of water to the firn matrix is characterized by a decrease in Vs and progressive increase in Vp:Vs; downstream, the increased clear ice content (9% increase from Site 7 to Site 12) increases Vs and decreases Vp:Vs, bringing it closer to that of pure ice. These additional ice lenses reduce the porosity and storage potential of the firn and potentially facilitate water flow consistent with Miller et al. (2018). From the flow measurements in Figure 4a, the ice lenses do not appear to impede lateral flow within the aquifer but could reduce vertical flow if the ice lenses cover a large area (e.g., Humphrey et al., 2012). These along-flow changes suggest that the aquifer transitions from a flooded firn layer with minor clear ice Figure 5. Schematic interpretation of the firn aquifer system. The enlarged firn images were adapted from Blunier and Schwander (2000).
lenses upstream to a more heterogeneous aquifer with many clear ice lenses downstream ( Figure 5). We interpret the presence of clear ice lenses as being diagnostic of episodic refreezing of the aquifer system, in particular at the aquifer base where water-saturated firn produces clear ice layers (Figure 4a; Miller et al., 2018). Recent studies have suggested that the firn aquifer formed earlier at the lower elevation sites (e.g., Miller et al., 2020), therefore undergoing more cycles of melt and freeze, having a longer history of advection and existence, compared to the upslope sites.
Our observations suggest a decrease in aquifer porosity below 1.65 km elevation: The porosity in the young aquifer may exceed 11% (average from Sites 1 to 7), whereas the mature aquifer may have porosity locally reduced to 1% (average from Sites 10 to 12), notwithstanding our error estimates. Compared to porosities derived from Vp alone (Montgomery et al., 2017), the downslope sites' (10-12) porosities are underestimated in this study, likely owing to the strong dependence Vs has on shear modulus, where the increase of ice in the matrix causes an increase in the shear modulus to that more like pure ice, where the porosity is 0.
This paper presents the first insights into the Vs distribution of a firn aquifer. It details how a complete seismic record can be used to extract elastic properties of the aquifer and inform porosity and stiffness. Our novel methodology enables the direct integration of multiple geophysical data sets, providing a comprehensive uncertainty analysis, circumventing ambiguous solutions that are often obtained from individual datasets. We provide an evidenced description of the development of the aquifer highlighting its complex spatial variability and the reduced water storage potential of older systems. This insight has important implications for predictive modeling and understanding existing and new firn aquifers being discovered globally. We recommend a combination of GPR, seismic, and borehole measurements for deriving accurate properties of an aquifer, combining the spatial coverage of geophysical data with the detailed ground-truth of intermittent borehole measurements.