Geophysical constraints on the properties of a subglacial lake in northwest Greenland

. In this study, we report the results of an active-source seismology and ground-penetrating radar survey performed in northwestern Greenland at a site where the presence of a subglacial lake beneath the accumulation area has previously been proposed. Both seismic and radar results show a ﬂat reﬂector approximately 830–845 m below the surface, with a seismic reﬂection coefﬁcient of − 0 . 43 ± 0.17, which is consistent with the impedance an intermittent bottom A strong and reﬂections a of sediments although its thickness and material properties are un-certain. Finally, we use these results to conduct a ﬁrst-order assessment of the lake origins using a one-dimensional thermal model and hydropotential modeling based on published surface and bed topography. Using these analyses, we narrow the lake origin hypotheses to either anomalously high geothermal ﬂux or hypersalinity due to local ancient evaporite. Because the origins are still unclear, this site provides an intriguing opportunity for the ﬁrst in situ sampling of a subglacial lake in of relative to hydrology.

While the presence and nature of subglacial lakes underlying the Antarctic ice sheet has been studied for more than 50 years, the existence of subglacial lakes below the Greenland ice sheet (GrIS) is a relatively recent discovery and comparatively little 35 is known about their properties and origin. Detection of subglacial lakes has relied on a variety of methods, including airborne radio-echo sounding (RES) (Robin et al., 1970;Siegert et al., 1996;Langley et al., 2011;Palmer et al., 2013;Young et al., 2016;Bowling et al., 2019) satellite altimetry measurements (Fricker et al., 2007;Palmer et al., 2015;Siegfried & Fricker, 2018;Willis et al., 2015), and active source seismic experiments (Horgan et al., 2012;Peters et al., 2008). Using these techniques, approximately 400 subglacial lakes have been detected in Antarctica (Wright & Siegert, 2012), of which 124 are 40 considered "active" (Smith et al., 2009). In Greenland, subglacial lakes were first detected in RES data by Palmer et al. (2013), who identified two small (roughly 10 km 2 ) flat regions of anomalously high basal reflectivity below the northwestern GrIS.
These features, named "L1" and "L2", were discovered below 757 m and 809 m of ice respectively. Recently, Bowling et al. (2019) greatly expanded the inventory of subglacial lakes in Greenland to approximately 54 candidates based on a combination of RES and satellite altimetry data. The new inventory shows that, in contrast to subglacial lakes in Antarctica which tend to 45 form under thick (> 4 km) warm-based ice in the continental interior, the majority of subglacial lakes in Greenland are found under relatively thin (1 -2 km) ice near the margins of the GrIS. Bowling et al. (2019) find that most subglacial lakes in Greenland appear to be stable features, showing temporally consistent RES signatures and an absence of vertical surface deformation over the decadal time scales of observation. Of the 54 candidate lakes, only 2 showed signs of vertical surface deformation indicative of active draining or recharge. While evidence of active subglacial lake dynamics in Greenland has 50 been sparse (Palmer et al., 2015;Willis et al., 2015), approximately 40% of Antarctic subglacial lakes appear to be active features (Bowling et al., 2019).
The formation and location of the detected subglacial lake features in Greenland remains elusive because many are located in regions where observations and modeling suggest that the base of the ice is frozen to its bed (MacGregor et al., 2016). 55 Complicating our understanding of the nature of subglacial lakes is the fact that uniquely identifying lakes in airborne radar data is challenging since basal reflectivity is sensitive to both the physical properties and the roughness of the material underlying the ice (e.g., Jordan et al., 2017). Amplitude anomalies of radar echoes in the range of +10 to +20 dB are often interpreted as subglacial lakes, although flat regions of saturated sediment may produce similar anomalies. Furthermore, the total volume of water stored in subglacial lake systems is unknown since airborne and space based remote sensing observations 60 are incapable of measuring lake depth (i.e., water column thickness). Seismic investigations provide an independent means of confirming the presence of subglacial lakes and are capable of measuring lake depth and underlying geological structures which can provide valuable clues into their formation and total volume.
In this study, we report results of an active source seismology and ground penetrating radar (GPR) survey performed in data. Our survey was conducted above subglacial lake candidate L2 which is within the accumulation area of the ice sheet.
Both seismic and GPR results show a flat reflector approximately 830 -845 m below the surface, with a seismic reflection coefficient of -0.42 +/-0.15, which is consistent with the acoustic impedance contrast between a layer of water below glacial ice. Additionally, in the seismic data we observe an intermittent lake bottom reflection arriving between 14 -20 ms after the 70 lake top reflection, corresponding to a lake depth of approximately 10 -15 m. Strong coda following the lake top and lake bottom reflections is consistent with a package of lake bottom sediments although its thickness and material properties are uncertain. Finally, we use these results to conduct a first-order assessment of the lake origins using a one-dimensional thermal model and hydropotential modeling based on published surface and bed topography (Morlighem et al., 2017). Because the origins are still unclear, this site provides an intriguing opportunity for the first in situ sampling of a subglacial lake in 75 Greenland, which could better constrain mechanisms of subglacial lake formation, evolution, and relative importance to glacial hydrology.

Field experiment 80
In June 2018, we conducted a geophysical survey in northwestern Greenland above the candidate subglacial lake feature L2 (Palmer et al., 2013). This feature sits within a 980 km 2 drainage basin and is roughly adjacent (< 10 km) to the nearest drainage divide ( Fig. 1A and 1B). The mean annual snow accumulation rate at the field site is ~30 cm/yr ice equivalent (determined from RACMO2 model results, Noël et al., 2018). In order to confirm the presence of the subglacial lake and investigate its physical properties, we collected data using both active source seismology and GPR. The active source seismic experiment 85 consisted of a moving line of 24 40 Hz vertical component geophones spaced 5 m apart. For each line, we collected data at 4 shot locations using an 8 kg sledgehammer impacted against a 1.5 cm thick steel plate. At each source location at least 5 hammer shots were stacked into a single shot gather in order to increase the signal to noise ratio. The first shot location of each line was offset 115 m from the first geophone, and subsequent shot locations were moved 115 m along the line. After data was collected for each of the 4 shot locations, the line was moved 230 m east along the traverse and data collection was repeated. 90 The seismic line was moved a total of 10 times, totalling 40 separate shot locations. Using this geometry, we obtained reflection points at the ice bottom spaced every 2.5 m along a traverse totalling 2400 m ( Fig. 1C and 1D). The GPR survey was conducted across a ~5.5 km transect roughly parallel to the seismic survey (Fig. 1C), using a 10 MHz monopulse radar system.

Seismic and GPR imaging 95
We created a longitudinal seismic reflection image by bandpass filtering data between 100 -200 Hz and applying a normal moveout (NMO) correction with a velocity of 3700 m/s, which was found to be the average velocity of the ice column from https://doi.org/10.5194/tc-2020-321 Preprint. Discussion started: 13 November 2020 c Author(s) 2020. CC BY 4.0 License. NMO analysis of the primary bed reflection. High frequency spatial noise with wavenumber greater than 0.05 m -1 was removed with f-k filtering. Shot gathers with offsets of −115 m and 230 m from the first geophone contained an air wave arrival that was muted by zeroing a 10 ms window with a moveout of 315 m/s. The seismic reflection profile ( Fig. 2A) shows a clear ice 100 bottom reflection across the entire transect arriving with a two-way travel time between 400 -460 ms. The ice bottom multiple (i.e., a wave that has travelled twice between the surface and ice bottom) is also visible between 800 -920 ms. From hereon, we refer to the ice bottom reflection as R1 and its multiple as R2. At transect distances between 0 -1700 m, the R1 reflection is flat and relatively uniform in character, which we interpret to be the signal of the top of the subglacial lake. In this region, R1 arrives at 457 ms, which corresponds to a depth of 845 m, assuming an average VP of 3700 m/s within the ice. At larger 105 transect distances, the reflections arrive earlier with increasing distance, which likely reflects the bed topography adjacent to the subglacial lake. An additional reflection with opposite polarity of R1 is observed arriving between 14 -20 ms after R1, which is consistent with a lake bottom reflection (Fig. 2B). This signal is intermittently observed but is most continuous at transect distances between 660 -1200 m. The travel time differential between the lake top and lake bottom reflection is used to measure the thickness of the water column across the seismic section. Assuming VP in the lake of 1498 m/s (Table 1), the 110 lake is between 10 -15 m deep (Fig. 2C). An uncertainty of +/-50 m/s on lake velocity would correspond to a lake depth uncertainty of +/-0.5 m. A strong coda following the lake bottom reflection is apparent which is likely caused by a thin (~ 10 m) sediment package underlying the lake (see Discussion).
To verify our seismic results, we also created a GPR image by converting the radar data to depth using a radar velocity of 172 115 m/µs (see Supporting Information). The subglacial lake is apparent as a flat reflector at an elevation of ~510 m across the majority of the transect (Fig. 3A). The surface topography slopes gently to the west across the transect, hence the lake is slightly deeper towards the east (Fig. 3B). The lake is beneath 840 m of ice at transect distances between 2 to 4.5 km, which roughly corresponds to the location of the seismic survey. The transition from the lake top to the adjacent bed is observed at approximately 4100 m along the transect. 120

Basal reflectivity
In order to confirm that the flat region of R1 is consistent with a reflection from liquid water below ice, we calculate the reflection coefficient at the base of the ice by analysing the amplitudes of the first (R1) and second (R2) reflections. When both R1 and R2 are visible, the basal reflection coefficient cR can be determined as a function of incidence angle using Eq. 125 (1) where AR1 and AR2 are the amplitude of the first and second ice bottom reflections, respectively, a is the absorption coefficient, and L is the raypath length of the R1 reflection (e.g., Peters et al., 2008).
At a given geophone, two factors control the amplitude ratio between R1 and R2. First, R1 and R2 reflect off of the lake with slightly different angles, which changes the relative amount of energy partitioned into each reflection. Second, since R2 travels farther than R1, its amplitude is diminished due to geometrical spreading and attenuation. However, at incidence angles in this 130 study, the difference in reflection coefficients between R1 and R2 is negligible. Additionally, the path lengths of R1 and R2 vary by < 5% between their shortest and farthest offsets. Therefore, to calculate the reflection coefficient cR we use the normal incidence approximation and compare amplitude ratios AR2/AR1 on individual seismograms. In order to minimize the influence of the air wave on AR2/AR1 ratio we exclude data from geophones with offsets between 135 -155 m, where there is potential interference between R1 and the air wave. Measurements of AR1 and AR2 are made prior to f-k filtering. 135 The relationship between the absorption coefficient a and the seismic quality factor Q is given by Eq. (2), where c is the seismic velocity, and f is frequency (Bentley & Kohnen, 1976). While in principle, the spectral ratio of the R1 and R2 reflections can be used to determine the attenuation (Q -1 ) of the glacial ice (Dasgupta & Clark, 1998;Peters et al., 2012) the low signal to noise ratio of the R2 reflection prevents us from making a robust measurement. Here, we estimate the absorption coefficient a 140 based on the study of Peters et al. (2012), who reported Q = 355 +\-75 in the upper 1 km of ice in Jakobshavn Isbrae, western Greenland. Using Eq. (2) with c = 3.7 km/s, and assuming a frequency of 100 Hz (the predominant frequency observed in the reflections), this corresponds to an absorption factor a = 0.23 +/-0.06.
Assuming an absorption factor of a = 0.23, the average reflection coefficient of the lake bottom across the transect is -0.42 +/-0.15 (Fig. 4A). In Fig. 4B, we plot cR calculated for each shot gather above the lake as a function distance along the transect. 145 For comparison we show the expected reflection coefficients of several different geologic materials underlying glacial ice. The reflection coefficients were modeled using the two-term approximation of the Zoeppritz equations (e.g., Aki & Richards, 2002;Booth et al., 2015) with the material properties shown in Table 1. In contrast to other likely geological materials at the base of the ice, liquid water is expected to have a strongly negative reflection coefficient. The reflection coefficient modeled for lithified sediments underlying ice is similar in amplitude to liquid water but opposite in sign, thus, without polarity information 150 sedimentary rock strata could be mistaken for a lake signature. Here, we observe an opposite polarity of R1 and the seismic source which is consistent with the negative reflection coefficient expected for a subglacial lake. Thus, liquid water is the only material that is consistent with both the amplitude and sign of the observed reflections.

Lake geometry and volume 155
If our interpretation of the observed reflections as signals from the lake top and bottom is correct, it implies that L2 could hold a significant volume of water. Assuming the imaged lake depth of approximately 15 m is representative of average lake depth throughout the roughly 10 km 2 surface area determined by RES, we estimate the total volume of liquid water to be 0.15 km 3 (0.15 Gt of water). While this is only a small fraction of the 217 +/-32 Gt of ice that Greenland is estimated to lose each year to glacier discharge and surface melting (IMBIE Team Report, 2019), the net storage capacity of all of Greenland's subglacial 160 lakes could be appreciable. To verify our interpretation of the lake top and bottom reflections, we modeled synthetic seismic waveforms of shot gather 12, which contained some of the clearest reflections. This shot gather corresponds to transect distances between 660 -720 m in the seismic reflection image. Synthetic seismograms were computed using Specfem2D (Tromp et al., 2008) for two simple layered models of a 12 m thick lake underlying 850 m of glacial ice. In the first model the lake is underlain by a thick layer of sediments that extends to the bottom of the model domain. In the second model there is 10 165 m of sediments overlying a discontinuity with the bedrock below. The seismic velocity profiles for the two cases are shown in the insets in Fig. 5B and 5C. The source used in the simulations was a Ricker wavelet with a dominant frequency of 100 Hz. shot gather contains a coda following the lake bottom reflection that is absent in the synthetics that do not include a discontinuity at the base of the sediment package (Fig. 5B). When a discontinuity between the sediment and underlying bedrock is included a strong sediment bottom reflection is introduced which more closely matches the observations (Fig. 5C). In the observed data it is difficult to clearly identify a sediment bottom reflection since the complex coda could be caused by reverberations within a thin sediment sequence, or many superposed reflections from individual discontinuities. However, if 175 the first positive peak following the lake bottom reflection represents the base of the sediment, we can estimate a sediment thickness of 8.5 m assuming a sediment VP of 1700 m/s (Table 1).

Lake origin
While our results confirm that L2 is indeed a subglacial lake, its presence is perplexing given its location with a mean annual surface temperature of -22 o C, and its position beneath a relatively thin column of glacial ice. In contrast to many well studied 180 subglacial lakes below Antarctica such as Lake Vostok that lie below ~ 4 km of ice, the basal temperature at our field site is expected to be well below the pressure dependent melting point of ice. Distinguishing between the different hypotheses of subglacial lake formation has implications for the stability and dynamics of the GrIS since they predict different basal thermal and hydrological conditions. Thus, constraining the temperature of L2 is an important goal. We determine the range of possible basal temperatures using a 1D steady state advection-diffusion heat transfer model solved using the control volume method (Patankar, 1980). The modeling assumes an ice density r = 920 kg/m 3 , a heat capacity * = 2000 J/kg/K, and a thermal conductivity = 2.3 W/m/K. The basal geothermal heat flux is varied between 50 -60 mW/m 2 , which is consistent with estimates derived from magnetic data (Martos et al., 2018) and thermal isostacy modeling (Artemieva, 2019). Fig. 6 shows results for surface temperatures Ts of -20 o C and -22 o C and accumulation rate ranging from 0 to 0.3 m/yr ice equivalent.
When advection is ignored (i.e., no ice accumulation), most scenarios predict frozen bed conditions with the exception of the 190 relatively warm surface condition (TS = −20 o C) and high heat flow ( = 60 mW/m 2 ) scenario (Fig. 6A). When ice accumulation is considered, all scenarios predict frozen bed conditions (Fig. 6B). For an ice accumulation rate of 0.3 m/yr, which most closely matches the conditions of the field site, the basal temperature is expected to be between approximately -12 o C and -14 o C.

195
There are several possible explanations for the existence of liquid water in L2, including hypersalinity, recharge by surface meltwater, high geothermal flux, and latent heat from freezing.
(1) Hypersalinity: If the lake is hypersaline it would depress the freezing temperature. Because the ice surrounding the lake would be frozen to the bed, it would form a closed hydrologic system that could remain isolated on geologic timescales. In 200 this scenario, the lake could represent a body of ancient marine water that was trapped as glacial ice advanced over the area and potentially further enriched in salt through cryogenic concentration processes (Lyons et al., 2005(Lyons et al., , 2019. Similar hypersaline lakes with salt concentrations several times higher than sea water are known to exist below the McMurdo Dry Valleys in Antarctica (Hubbard et al., 2004;Lyons et al., 2005Lyons et al., , 2019Mikucki et al., 2009) and in the Devon Ice Cap, Canada (Rutishauser et al., 2018). However, because the current elevation of the lake is more than 500 m above sea level, it is unlikely 205 to be trapped sea water as in the McMurdo Dry Valleys; however, we cannot rule out an ancient evaporite deposit as is proposed for the Devon Ice Cap (Rutishauser et al., 2018). However, the geologic map of Greenland does not indicate likely evaporites in this area (Dawes, 2004).
(2) Surface meltwater: The lake may be part of an open hydrological system that is continually recharged by surface meltwater. 210 This lake, however, is in the accumulation area of the ice sheet, near the ice divide (Fig. 1B) and there are no obvious sources for significant surface recharge. To determine possible pathways for surface recharge, we estimate the local hydraulic head based on surface and bed elevations (Fig. S4) and find no pathways given the present resolution of bed and surface topography.
It is possible that a subglacial pathway exists that is smaller than the resolution of BedMachine (Morlighem et al., 2017). If the hydrological system is connected and the rate of recharge matches or exceeds the rate of freezing, a lake could persist 215 despite sub-freezing temperatures in the lower part of the ice. At other locations in Greenland, observations of vertical surface deformation and collapse features have suggested that surface meltwater plays a prominent role in subglacial lake formation and dynamics (Palmer et al., 2015;Willis et al., 2015). Thus, subglacial lakes in Greenland could become more widespread as the GrIS responds to a warming polar climate. https://doi.org/10.5194/tc-2020-321 Preprint. Discussion started: 13 November 2020 c Author(s) 2020. CC BY 4.0 License.
(3) High geothermal flux. Anomalously high basal heat flux may promote melting of the GrIS from below (Fahnestock et al. 2001;Rogozhina et al., 2016). If this is the case, the local geothermal heat flux must greatly exceed regional estimates of the geothermal heat flux beneath the northwestern GrIS, which are typically in the range of 50 -60 mW/m 2 (Artemieva, 2019;Martos et al., 2018;Rogozhina et al., 2016). Based on the one-dimensional model shown in Fig. 6, a geothermal flux on the order of 100 mW/m 2 would be necessary to sustain the lake. While high heat flux in this region would be unexpected based 225 on the cratonic bedrock geology and lack of recent volcanism, a local region of high heat flux could be promoted by the presence of upper crustal granitoids rich in radiogenic heat producing elements or hydrothermal fluid migration through preexisting fault systems (e.g., Jordan et al., 2018).
(4) Latent heat from freezing. For either the isolated lake of actively freezing brine (as in Hypothesis 1) or the hydrologically 230 connected continuous flow (Hypothesis 2), or if the lake is a relic of a larger freshwater body that is slowly freezing, the thermal profile of the ice would show a curvature change at depth due to a latent heat source at the bottom boundary. Given a latent heat of freezing of 334 J/g, freezing a layer 1 m thick to the bottom of the ice over one year is roughly equivalent to increasing the geothermal flux by 10 mW/m 2 . Sustaining a freezing rate of several m/yr to generate the latent heat necessary to maintain warm basal ice is less likely than locally elevated geothermal anomaly. We, therefore, rule out cryoconcentration 235 and narrow the lake origin hypotheses to either anomalously high geothermal flux or hypersalinity due to local ancient evaporite. Measuring the thermal profile above this lake would provide important information to assess these hypotheses as the basal ice temperature would be near 0 o C for the freshwater lake created by high geothermal flux hypothesis and substantially below zero for the hypersaline lake created by evaporite.

240
A freshwater lake and a hypersaline lake have different physical properties and thus may have different signatures that could be detected in geophysical surveys. Radar reflections from an ice/brine boundary undergoing freezing and cryoconcentration of the brine is known to cause scattering and decrease the reflectivity (Badgeley et al., 2017) which we do not see in our data; this provides a second justification to rule out cryoconcentration. Further, because the seismic velocity and density of water depends on temperature and salinity, we would expect that lakes formed by different mechanisms would have slightly different 245 basal reflection coefficients cR. To estimate the difference in cR for different lake conditions, we calculate the basal reflectivity as described in Section 2.2 for low and high salinity end member scenarios. The VP in water is calculated as a function of pressure, temperature, and salinity using Equation (7) from Coppens (1981). At 0 o C, the difference between cR of a freshwater lake and a saline lake with 100 ppt (approximately three times the salinity of ocean water) is 0.06 at normal incidence. While a hypersaline lake could remain liquid at temperatures below 0 o C, reducing temperature and increasing salinity have the 250 opposite effect on seismic velocities thus reducing temperature below 0 o C would decrease the contrast in cR. Because the difference in reflection coefficient calculated for end member scenarios of 0.06 is smaller than the uncertainty on our measurements (approximately 0.15) we are unable to distinguish between different models of lake salinity and temperature https://doi.org/10.5194/tc-2020-321 Preprint. Discussion started: 13 November 2020 c Author(s) 2020. CC BY 4.0 License. from our basal reflectivity analysis. On the other hand, because the electrical resistivity of water is strongly dependent on salinity, magnetic sounding could provide useful constraints on lake composition. Future active source seismic experiments 255 focused on measuring Q within the overlying glacial ice could help us constrain the thermal profile above the lake (Peters et al., 2012). Additionally, repeated seismic reflection or GPR surveys calculated along the same transect could provide clues into whether or not lake levels are changing over time (e.g., Church et al., 2020). Direct sampling with drilling would provide the best measurements on subglacial lake properties and could also yield useful biological and paleoenvironmental information.

Conclusions
We conducted an active source seismic reflection and GPR survey in northwestern Greenland above a site that was previously identified as a possible subglacial lake. We observed a horizontal reflector at the base of the ice across the majority survey with a seismic reflection coefficient of -0.42 +/-0.15, confirming the presence of a subglacial lake. Additionally, we observed a lake bottom reflection near the center of our seismic profile consistent with a lake depth of approximately 15 m. From 265 previous observations of the lateral extent of the lake based on RES (Peters et al. 2013), we estimate the subglacial lake holds a total of 0.15 Gt of water. Strong coda arriving after the lake-bottom reflection suggests that the lake is underlain by a sedimentary package but its thickness and material properties are uncertain. To the authors knowledge, this is the first time a ground based geophysical survey has confirmed the existence of a subglacial lake in Greenland and provided constraints on its depth. Understanding the nature and origins of recently detected subglacial lakes in Greenland is important since wet basal 270 conditions enable glacial ice to flow more easily which can further promote ice loss. Our analysis of the seismic, radar, as well as thermal and hydropotential modeling narrow the lake origins to either locally high geothermal flux or an ancient evaporite deposit. Future work, such as additional geophysical investigations or drilling expeditions, should focus on constraining the temperature and salinity of the lake which will provide clues to its origin. Additionally, monitoring the vertical motion of the ice surface above the subglacial lake system with satellite altimetry could provide further insight into the subglacial hydrology 275 (e.g., Siegfried & Fricker, 2018).