Surface nuclear magnetic resonance observations of permafrost thaw below floating, bedfast, and transitional ice lakes Open Access

The Arctic climate is changing with implications for weather, ecosystems, biogeochemical cycles, and human habitability. Along with warmer air temperatures, comes the potential to warm permafrost, resulting in mobilization of subsurface carbon that may escape to the atmosphere in gas form and contribute to the permafrost carbon feedback (Schuur et al., 2015). Arctic lakes have been identified as an important pathway for carbon to move between the subsurface and the atmosphere (Walter et al., 2007) and as a harbinger for change as permafrost thaw below lakes is highly sensitive to annual maximum ice thickness and mean annual lake bottom temperature (Arp et al., 2016).

Lakes that are initiated by permafrost thaw and subsidence are known as “thermokarst lakes,” and a zone of thaw expanding below a thermokarst lake is known as a “talik.” Taliks may also form below nonthermokarst water bodies, or independent of lake systems (i.e., terrestrial settings). Seasonal ice conditions on any given Arctic lake fall into three regimes: floating ice, bedfast ice, or transitional ice (Engram et al., 2018). Floating ice lakes are those that retain a water column below the ice throughout the entire winter and promote the thaw of permafrost below the lake. Bedfast ice lakes are shallower than the maximum potential ice thickness; therefore, the lake ice will meet the lake bottom over $>95%$ of the lake area at some point during the cold seasons — bedfast ice lakes preserve the permafrost below the water column even though the lake acts to warm the permafrost. If lake ice freezes to $>60%$ of the maximum ice thickness on average over multiple years, taliks will most likely not develop (Mackay, 1992; Arp et al., 2016). Transitional lakes fall between these two end member categories over a given duration of observation, and they may be transitioning from bedfast to floating ice, or vice versa, or they may be intermittent and not have an obvious trend over a period of observation. In general, we hypothesize that under a warmer and snowier climate, lake ice will be thinner, resulting in more transitional ice lakes than bedfast ice lakes and promoting widespread subaqueous permafrost thaw.

Lake ice regime transitions have direct implications for ecology (Prowse and Brown, 2010), hydrology (Walvoord and Kurylyk, 2016), energy balance (Jeffries et al., 1996), permafrost (Arp et al., 2016), and water resources (Jones et al., 2009). The ice regime of a lake may be determined by direct measurement with ice drilling equipment (Arp et al., 2012) or through analysis of radar remote sensing data (Surdu et al., 2014; Atwood et al., 2015). Measurement in one year can only detect floating ice or bedfast ice conditions. It is necessary to observe a trend of bedfast and floating lake ice conditions over time to categorize lakes as transitional. Recent results have shown high interannual variability of lake ice regime pointing to a complex set of interacting factors that drive the long-term transition to floating ice lakes (Engram et al., 2018). Over a 25-year time horizon, only localized transitions to increased floating ice were observed, whereas most of the Alaskan Arctic showed inconclusive patterns of changing ice regime. We are motivated by the hydrologic function and role that lakes are thought to play as carbon sources and sinks. If the permafrost below is frozen, the carbon remains sequestered. However, once thawed, it becomes mobile and may contribute to the permafrost carbon feedback. Understanding how a lake talik responds to different lake ice conditions may help to predict the risk of large carbon stores becoming released.

We also want to understand how a potential increase in floating lake ice regimes will result in a larger number of lakes initiating taliks. Thermal modeling under conditions similar to the modern North Slope of Alaska suggest that a talik could expand vertically into the permafrost 10–25 m below a lake over 500 years under a steady climate depending on the lake bed temperature and substrate (Ling and Zhang, 2003). Conceptually, lakes that have a mean annual lake bottom temperature $<0°C$ consistently over time should have no thawed permafrost below them, whereas substrates below floating ice lakes should be thawed from the bottom of the water column to the base of the talik, or to the regional permafrost thickness. Transitional lakes are expected to have talik development related to the direction of ice regime change, potentially even including layers of variably frozen and thawed material related to periods of either floating or bedfast ice conditions. All of these assumptions may be readily modeled; however, they are difficult to measure for two main reasons. First, very few Arctic lakes have ever been drilled to directly measure the talik properties (Brewer, 1958; Mackay, 1997). Second, although permafrost is defined by temperature, there are conditions in which a lake talik may be below 0°C, but substantial pore water may remain unfrozen due to depressed freezing points (Kleinberg and Griffin, 2005). In these cases, liquid pore water $<0°C$ is present, even though the material would be classified as permafrost, and biogeochemical reactions with carbon cycling implications may still occur.

Active-source geophysical measurements have been used for a variety of investigations of thermokarst lake talik properties and geometry (Shwamborn et al., 2002; Yoshikawa and Hinzman, 2003). Waterborne seismic measurements have been able to image the base of a talik in deep, continuous permafrost of Siberia (Shwamborn et al., 2000), though lake depth and substrate can be physical limitations to the success of this type of observation. Airborne electromagnetic imaging has been successful at observing taliks, particularly related to discontinuous permafrost (Minsley et al., 2012); however, it is ambiguous to interpret resistivity values as either contrasts in the substrate or freeze/thaw state. Other methods such as ground-penetrating radar have proven valuable for very shallow investigations related to lake ice and lake sediments (Shwamborn et al., 2002; Jones et al., 2013), but they do not have the depth sensitivity to fully explore the talik geometry and physical properties. Surface nuclear magnetic resonance (NMR) has been demonstrated on thermokarst lakes in the past (Parsekian et al., 2013) with sensitivity to talik thickness and unfrozen water content in the sublake sediments. The use of surface NMR in permafrost and Arctic environments has been relatively limited, although successful recent experiments have measured internal glacial water in Greenland (Legchenko et al., 2018) and unfrozen water in permafrost where the freezing point is suppressed in coastal Svalbard (Keating et al., 2018).

In this study, we use surface NMR soundings because of the direct sensitivity of this measurement to unfrozen water, therefore reducing some of the interpretation ambiguity. We focus on lakes with three major lake ice regime types (floating, bedfast, and transitional) to understand the related distribution of unfrozen water in sublake sediments. As part of our effort to understand how lakes control permafrost thaw state, the specific question we seek to address in this study is: What relationship exists between lake ice regime and permafrost thaw below lakes?

We have investigated floating, bedfast, and transitional ice lakes, as well as terrestrial permafrost sites, in two regions on the North Slope of Alaska near Inigok (INI) and Teshekpuk (TES) Lake (Figure ¹). Terrestrial sites with no known thaw history were used for control/comparison with on-lake measurements. The region around INI is located on the inner coastal plain with rolling hills topography and a typical relief of 20 m. Lakes near INI appear to be often controlled by undulations in a Pleistocene-aged dune field, and as a result, they are typified by shallow margins and deep floating ice center basins that likely support talik formation. In contrast, the region north of TES Lake is on the outer coastal plain with nearly flat topography having a relief of $<5m$⁠. This landscape is occupied by numerous drained lake basins and second (or later)-generation thermokarst lakes that are generally oriented northwest. Most of our lake study sites indicated in Figure ¹ are named following a numbering scheme developed during previous research (Arp et al., 2011), although some named sites use letter abbreviations: Peatball Lake (PBL), Peatball Bluff (PBB), and Independent Fox (INDFOX).

Based on temperature measurements in a network of boreholes (Figure ²), we know that the thickness of permafrost along the coast of the Beaufort Sea is 320–400 m, whereas at approximately 75 km inland on the inner coastal plain, the permafrost is 286 m (Clow, 2014). A borehole was drilled and logged at INI, but temperature logs were not measured. The closest permafrost thickness observation available is the Northern Inigok borehole, approximately 25 km northeast of INI. Drilling logs (Figure ²) indicate a thin mantle of recent unconsolidated deposits at all sites, underlain by sandstone and siltstone formations (Gryc, 1988). The resistivity measured in the borehole was between 1 and 100 ohm-m, and the resistivity generally decreased with depth (Figure ²).

Gradients in earth’s magnetic field may influence the surface NMR measurement, and consequently, we have overlain interpolated airborne magnetic observations in Figure ¹ (Connard et al., 1999). Though the magnetic field gradient data were measured years before our NMR observations, the International Geomagnetic Reference Field (IGRF) model determined by the National Oceanic and Atmospheric Administration (NOAA) magnetic field calculator indicates a slow total field rate of change of $−3.8nTyear−1$ for TES Lake and $−5.8nTyear−1$ for INI. There are likely some differences in the magnetic field properties between the historic observations (Figure ¹) and the current day; however, it is reasonable to assume that the overall pattern of low spatial gradients is consistent.

The abundance and distribution of lakes and drained lake basins in both study areas is a testament to the importance of lakes impacting the regional permafrost characteristics on the Arctic Coastal Plain (Hinkel et al., 2005; Jorgenson and Shur, 2007). Based on an existing synthetic aperture radar (SAR) classification of lake ice regime (Engram et al., 2018), the TES region has approximately 51% by area of stable floating ice lakes in comparison with approximately 94% stable floating ice lakes in the INI region. There are approximately 38% by area of stable bedfast ice lakes in the TES region, and there are only approximately 4% stable bedfast ice lakes in the INI region. The remainder of 11% and 2% by area for TES and INI, respectively, are either transitional or intermittent.

The surface NMR method uses electromagnetic pulses generated using a wire loop on the earth’s surface to excite subsurface hydrogen protons in unfrozen water aligned with earth’s magnetic field at equilibrium. After excitation, the protons relax back to alignment with earth’s field, inducing a measurable time-varying field observed by switching the transmitter loop to the receiving mode. The depth of investigation is partially controlled by the loop diameter, with the general rule that the maximum resolvable depth is approximately equal to the loop diameter. The frequency of the excitation pulse is tuned to the Larmor frequency $fL$ of hydrogen, so the measurement is directly sensitive to unfrozen water. The Larmor frequency is determined locally for each measurement as the product of the total field strength and the gyromagnetic ratio of water. For the following sections related to surface NMR measurement, all references to “water” imply liquid water, given that the relaxation time of ice is too short to measure with surface NMR. The duration of the excitation pulse impacts the location in the subsurface where water achieves its maximum excitation (i.e., 90° from equilibrium), and longer pulse lengths can lead to a deeper depth of investigation. The most common surface NMR pulse sequence is the free-induction decay (FID) in which a single excitation pulse is used, due to engineering limitations that restrict the number of useful subsequent refocusing pulses due to available power. The pulse moment is the amplitude and duration of the transmitted excitation waveform, reported in A s, and larger pulse moments enable observation of water at greater depth. The exponential relaxation of hydrogen protons, measured after the transmitted pulse is shut off, is then inverted to retrieve a vertical distribution of volumetric water content (VWC) and $T2*$⁠. The relaxation parameter $T2*$ obtained using FID pulse sequence is closely related to the transverse relaxation parameter $T2$ more commonly associated with NMR, and it contains information about the mean pore size (Grunewald and Knight, 2011). Large-scale gradients in earth’s magnetic field or contrasts in magnetic susceptibility between pore water and grains can contribute to enhanced relaxation (i.e., shorter relaxation times). Please refer to Weichman et al. (2002), Walsh (2008), and Behroozmand et al. (2015) for additional details and comprehensive reviews of the surface NMR principles.

We made surface NMR measurements using a GMR (Vista Clara, Mukilteo, WA) with either 75 or 90 m circle loops located approximately at the center of each lake of interest. We used similar measurement parameters for two additional control sites on terrestrial upland permafrost that showed no evidence of ever previously being a lake. Due to the low cultural noise conditions in the region, noise canceling loops were not required. When making NMR measurements near developed areas, noise canceling loops are commonly used to measure ambient noise so that it can be later subtracted from data measured on the transmitter/receiver loop. We used an FID excitation pulse duration (⁠$tP$⁠) of 0.02 s for all except one data set when $tP=0.04s$ was used to increase the depth of investigation. Earth’s magnetic field inclination (⁠$I$⁠) values were retrieved from the NOAA IGRF model based on the geographic location of each surface NMR sounding. Immediately before each NMR measurements, we used a Geometrics G816 magnetometer to measure the magnitude of earth’s background magnetic field $B0$ three times each at five locations within each loop to measure for calculation of $fL$ (Table 1). Each surface NMR data set required an average of 20 min of measurement time to acquire, and there was no evidence of time-varying $B0$ magnitude (i.e., changing Larmor frequency in the NMR data) during the measurements. Based on the power discharged with each pulse, the instrument automatically determines the number of pulse moments $q$ and provides a logarithmic sampling. Due to this automated procedure, the selected $q$ values are slightly different for each site.

We used MRSMatlab to invert the surface NMR data (Müller-Petke et al., 2016). After stacking, the data were all preprocessed with a 200 Hz band-pass filter and trimmed to an 0.8 s record length. The sensitivity kernel for each inversion was built based on local magnetic field properties, measurement parameters, and the average subsurface resistivity $ρ$ distribution (Braun, 2007; Braun and Yaramanci, 2008) measured by time-domain electromagnetics (Parsekian et al., 2018) that were similar to values measured in nearby boreholes (Figure ²). The QT data cube refers to the time series $T$ for each $q$ stored as a matrix for inversion in a single step (Müller-Petke et al., 2011). All smooth, single-relaxation time inversions were run using a rotated complex data space where the imaginary component of the data is combined with the real part leaving a single QT data cube containing all signal information. This is preferable over inverting only the real part because it accounts for an imaginary signal, and it is more stable than inverting the real and imaginary data spaces separately (Müller-Petke et al., 2011). The smooth inversion was chosen in this case over the blocky minimum layer inversion used on thermokarst lakes by Parsekian et al. (2013) because (1) based on time-domain electromagnetic observations, the talik thickness would be well beyond the depth of sensitivity of the surface NMR and (2) at transitional ice lakes that may have talik with complex thaw geometry, gradual transitions of liquid water content would be plausible and therefore invalidate the assumptions of the blocky inversion. After inversion, we calculated resolution matrices following (Müller-Petke and Yaramanci, 2008) that take into account the regularization as well as the measurement configuration (loop diameter and $q$’s).

To analyze the data measured on floating ice lakes, we carry out the typical inversion on the entire unmodified data set first, and then we do a second inversion to suppress the lake water signal following Falzone and Keating (2016). We use the forward modeling module in MRSMatlab to calculate the NMR signal for the water column only based on field observations (Table 1) of ice thickness $dICE$ and lake water depth (⁠$dTOTAL−dICE$⁠). Then, we subtract the forward-modeled signal from the field measured signals to leave a new data set with only the sublake-sediment NMR signal remaining. The subtracted signal is then inverted as previously described, constrained for a maximum $T2*$ value of 0.5 s and a maximum water content of $0.5m3m−3$⁠. The forward modeled water signal for PBL was $1m3m−3$ water content and a 0.3 s $T2*$⁠, whereas for INDFOX, we used $1m3m−3$ water content and 0.6 s $T2*$⁠. We based the $T2*$ relaxation times on the optimization of the model relative to the observed data finding the largest $T2*$ value that still removed all long signals from the field data.

At the center of each lake that was floating ice at the time of field observation, we measured lake $dICE$ in three or more locations within an approximately 15 m radius using a powered ice drill (Kovacks Enterprise, Roseburg, OR). The ice thickness and lake depth were measured with a Kovacs ice thickness gauge, which consists of a spring-loaded brass arm that pulls tight against the bottom ice surface of a 5 cm diameter auger hole. For each lake that was bedfast at the time of field observation, we measured $dICE$ in at least one location.

Lake ice classifications based on SAR remote sensing data are taken from the database produced by Engram et al. (2018). To develop that database, annual satellite SAR data from the time of maximum ice thickness were analyzed to quantify backscatter intensity from lakes, where high intensity is associated with the high dielectric contrast at the ice/water interface and low intensity is related to the low dielectric contrast at the ice/substrate interface (Jeffries et al., 1994; Morris et al., 1995). Each data set resulted in a population of categorized pixels for each lake ranging in size from 10 to 75 m (depending on the available satellite sensor). By using a time series of observations, each lake can then be categorized as bedfast ice, floating ice, or transitional (Figure ³).

There was no observable NMR signal at any of the terrestrial control sites or bedfast ice lakes as shown in Figure ⁴, plotted as the absolute value of the complex signal. When inverted (Figure ⁴), all of the data sets produced $<0.01m3m−3$ water content consistently throughout the depth profile. Due to the uniformly small inverted water content $<0.01m3m−3$ for these sites, we have omitted a figure showing these results. The slightly higher noise observed in the PBB data set (Figure ⁴) compared with the INI and TES data sets (Figure ⁴, ⁴, and ⁴) can be attributed to the larger loop used in the measurement at PBB. The residual plots are shown in Figure ⁴, and they are calculated as the difference between the $dobs$ and $dest$ immediately above. If the data fit the model exactly (i.e., a perfect inversion result), the residual plot would be uniformly dark-blue colored, whereas any part of the modeled signal that is not explained by the data will be a warm color. Visual inspection of the residual plot is used in conjunction with the calculated residual value (Table 1) for each data set to judge each inverse result as reliable. The resolution matrices, shown in Figure ⁴, can be used to determine the maximum depth of investigation by looking at where the weighting factor (warm colors) declines with depth at the largest $q$ value (Müller-Petke and Yaramanci, 2008).

Measured QT data cubes (Figure ⁵) show the high-amplitude, long-relaxation-time signal associated with lake water, in comparison to the data sets after the lake water is subtracted. From the data, it is clear that a measurable NMR signal is present at all pulse moments at INDFOX and only through 1 A s at PBB. This side-by-side comparison also illustrates the effectiveness of the water subtraction approach. The modeled (inverse result) QT data cubes in Figure ⁵ are visually similar to the measured data, and this is corroborated quantitatively through the residuals (Figure ⁵; Table 1). We do observe some subtle patterns remaining in the residuals, particularly in the unmodified data inversions (Figure ⁵ and ⁵), suggesting that the inversion was not able to find a model that entirely describes the data. Residuals in the subtracted inversions (Figure ⁵ and ⁵) are more randomly distributed, suggesting that once the water signal is removed, the inversion is able to find a model that more closely fits the data. The resolution matrices show sensitivity to 50–55 m for both soundings (Figure ⁵), and we note that the same kernel was used for unmodified and subtracted data sets at each lake.

The floating ice lake water content profiles based on the unmodified data sets (Figure ⁶ and ⁶) show low unfrozen water content layers at the top, followed by zones of approximately 100% water content corresponding to the lake water column. For comparison, the direct-measured total lake depth (⁠$dTOTAL$⁠) for INDFOX was 5.65 m, whereas the value retrieved from surface NMR was 5.90 m. Similarly, the $dTOTAL$ for PBL was 2.29, and the value retrieved from surface NMR was 2.16 m. Below the lake, the INDFOX water content profile is variable (0.05–$0.80m3m−3$⁠) and variable with depth. The high-variability and thin high-water-content layers (1.6 m thick at $0.8m3m−3$⁠) are interpreted to be inversion artifacts. For PBL, the water content below the lake is primarily a single thin high-water-content layer that is 3.3 m thick and has $0.84m3m−3$ water content. Because this unmodified data set is constrained to $1.1m3m−3$ water content due to the requirement of fitting the known lake water, this allows sublake water content to also be high; there is no way to have two constraints on the same inversion. Below the thin, high-water-content layer, the water content is approximately constant at $0.035m3m−3$ at depths $>11.7m$⁠. The $T2*$ profiles (Figure ⁶ and ⁶) are characterized by long relaxation times associated with lake water of 0.95 and 0.19 s for INDFOX and PBL, respectively. Below the lake, INDFOX $T2*$ ranges from 0.06 to 0.78 s until 40 m depth, when the $T2*$ increases substantially. This increase of $T2*$ at depth is too long to be physically realistic in a porous media, and it is therefore interpreted as an inversion artifact related to the challenges of having restricted single exponential $T2*$ and low water content values. It can be considered unreliable to fit multimodal relaxation curves that would be best suited for fine-textured sediments given the low signal-to-noise ratio of the surface NMR data (Müller-Petke et al., 2016). At PBL, the $T2*$ values range from 0.06 to 0.13 s, and they are variable across the entire depth range.

After the modeled lake water signal is subtracted from each data set, the layer of previously high water content is nearly entirely suppressed in both data sets (Figure ⁶ and ⁶), although there are still remnant small spikes in the zone where signal was subtracted that are interpreted to be associated with uncertainties in measurement and modeling of the lake water layer. After subtraction, the INDFOX profile is less variable than the unmodified result and only contains physically plausible water content values ranging from 0.045 to $0.43m3m−3$⁠. The subtracted PBL result shows a $0.48m3m−3$ water content layer from 7.6 to 15 m depth, and water content approximately $0.02m3m−3$ above and below it. Although this layer does have physically plausible water content, it is surprising that the sediment immediately below the water column — assumed to be unconsolidated and unfrozen — would have such a low water content. We hypothesize that fine-grained lake sediments may result in relaxation that is too fast to measure given the dead time of the surface NMR instrument (0.0055 s), or that vertical magnetic field gradients influence the relaxation time. Below 15 m depth, we interpret a change in lithology to sandstone, again with small pores and fast relaxation. The $T2*$ profiles below the lake water after subtraction range from 0.21 to 0.128 s at INDFOX and 0.13 to 0.470 s at PBL.

The difference in the water content profiles below these two floating ice lakes is interpreted to be due to the contrast in substrate between the “Sand Sea” (i.e., the paleo dune fields surrounding INI) and the marine-depositional coastal plain (Carter, 1981). The INI area has thicker accumulations of unconsolidated sediments than the TES region (Figure ²); in addition, the borehole at INI is at a topographic low point (40 m asl), whereas INDFOX is elevated (47 m asl) suggesting that unconsolidated sediments $>13m$ thick would be interpreted from the well log (Figure ²). The bottom of the high water content zone in the subtracted PBL sounding at 15 m depth likely corresponds with the 12 m thick unconsolidated thickness observed in the Drew Point borehole.

As seen in the L31 data set (Figure ⁷), there is no NMR signal at small pulse moments, consistent with the lack of measurable water in the top 15 m. However, a very clear NMR signal is present between 0.4 and 4 A s pulse moments giving a clear indication of unfrozen water in this range. In contrast, the TIL3 data (Figure ⁷) show clear, though short, NMR signals through the 1.855 A s pulse moment, consistent with unfrozen water throughout the near surface. The data set for ITIL16 (Figure ⁷) shows a clear NMR water signal until approximately the 0.63 A s pulse moment, and little signal below that consistent with permafrost or $T2*$ times that are too short to measure. At ITIL11, the data show clear NMR signals at all pulse moments consistent with thaw across the entire sensitive depth range. The inverted model results (Figure ⁷) are visually similar to the data, confirmed quantitatively by the residuals (Figure ⁷; Table 1). The residuals for TIL3 (Figure ⁷) suggest that the inversion result may slightly overestimate the water content. The resolution matrices show the depth of investigation sensitivity to between 45 and 60 m for the transitional lake ice data sets (Figure ⁷).

The transitional lakes were entirely frozen at the time of the measurement, and the depth of the lake is approximately 1.6 m (Table 1). The unfrozen water content profiles, shown in Figure ⁸, for the transitional lakes in the TES region (TIL3 and L31) are very different from each other. At lake TIL3, there is evidence of thawed sediments averaging $0.12m3m−3$ water content from the surface to at least a 12.8 m depth and unfrozen water content $>0.02m3m−3$ below. Similar to PBL, this may be interpreted as unfrozen, unconsolidated sediments near the surface, and a transition to sandstone at a depth that has NMR relaxation too fast for surface NMR to measure. In contrast, L31 has just $0.018m3m−3$ water content from the surface to a depth of approximately 15 m, then a 13 m thick layer with maximum water content of $0.16m3m−3$⁠, below which the measured water content is $<0.02m3m−3$⁠. This could be interpreted as an upper layer of permafrost and an isolated talik between 15 and 28 m deep. Other possible explanations would be that there are fine-textured lake sediments below L31 that have NMR relaxation times faster than can be measured with surface NMR, or that the SAR classification of the lake ice regime is erroneous for this site due to geomorphic processes causing partial lake drainage in the mid-2000 s (Jones et al., 2009).

Water content profiles from transitional lakes in the INI region show similar low water content (i.e., ice thickness) to approximately 1.8 m, and both have a zone of elevated water content in the near surface: $0.33m3m−3$ between 1.8 and 5.4 m below the surface for ITIL16 and $0.23m3m−3$ between 1.8 and 7.0 m below the surface for ITIL11. ITIL16 had a 0.05 m water layer below the ice, and ITIL11 was bedfast at the time of measurement, so it is not surprising that no bulk water signal is observed. Below 5.4 m at ITIL16, the water content is generally $<0.01m3m−3$⁠; however, at ITIL11, the sublake water content is $>0.05m3m−3$⁠, including a 13 m thick layer with a maximum $0.28m3m−3$ water content.

The $T2*$ relaxation profiles in the TES region lakes (Figure ⁸) are smooth and below 0.1 s. For TIL3, the zone associated with thaw in the top 12.8 m has $T2*$ of approximately 0.045 s. The slight increase in $T2*$ below 12.8 m is associated with a drop to very low water content; therefore, these relaxation times may not be reliable. The elevated water content zone 15–28 m below L31 has $T2*$ of 0.09 s, and due to the low water content above and below this zone, the $T2*$ values are likely unreliable. The ITIL16 $T2*$ relaxation time is 0.094 s in the zone of elevated water content, and it is likely unreliable outside this zone due to the low water content. At ITIL11, there is a water signal throughout the entire depth range and the $T2*$ values range from 0.085 to 0.122 s.

Our results show that (1) bedfast ice lakes had no detectable unfrozen water content, (2) floating ice lakes had highly detectable water content from the lake and the sediments below the lake and it was possible to isolate the signal from the unfrozen water within the talik, and (3) transitional ice lakes were variable, but in all cases, unfrozen water was detected in the sediment below the lake indicating the presence of a talik. In a changing climate, we anticipate warmer air temperatures and increased snowfall that lead to thinner lake ice; as a consequence, some bedfast ice lakes may transition to floating ice lakes (Arp et al., 2012; Surdu et al., 2016; Engram et al., 2018). Changing lake ice regimes may have environmental consequences to hydrologic cycles, energy balances, ecological habitats, and permafrost thaw. Our geophysical interpretations of liquid water content below lakes can be used to infer the state of thaw; therefore, we can begin to understand the implications of lake ice regime change on permafrost thaw. Our observations support the overarching hypothesis that bedfast ice lakes will be underlain by permafrost and floating ice lakes will have taliks of varying depths. Furthermore, we hypothesize that transitional ice lakes will have some complex freeze/thaw patterns as a function of depth determined by their lake ice regime history and whether they are changing from bedfast to floating or vice versa. Recent analyses of thousands of lakes on the North Slope using SAR satellite records reveal lake ice regime classification during the period from 1992 to 2016 (Engram et al., 2018), an unprecedented view into recent trends of lake ice regime change. However, permafrost thaw may be a slow process and we cannot know the system behavior before 1992, so our interpretations must be conservative and acknowledge the data limitations. Furthermore, our geophysical data set is considerably smaller than a satellite remote sensing data set due to the level of effort required for overland Arctic travel and the time needed to make surface NMR measurements. Therefore, even though this set of surface NMR measurements spanning a variety of lake ice regime types across two regions of Arctic Alaska is unique, we have only begun to capture the characteristics of system variability. Perhaps not surprisingly, we can clearly see that there is no evidence of permafrost thaw below bedfast ice lakes — they have the same absence of NMR signal as terrestrial permafrost of Pleistoscene age (Lenz et al., 2016). At the other extreme, we clearly see evidence of unfrozen water in the sediments below floating ice lakes as expected. The observations of water content distribution below transitional lakes are more challenging to interpret, though they do suggest more complexity than is predicted by the current generation of talik growth models (Ling and Zhang, 2003; West and Plug, 2008; Rowland et al., 2011). Direct observation or geophysical imaging of talik properties is limited in the North Slope region, and available information focuses primarily on floating ice lakes where steady increases of temperature with depth are observed (Brewer, 1958) that do not reveal information about variability in water content. Isolated taliks in discontinuous terrestrial permafrost have been observed using surface NMR (Parsekian et al., 2013).

The talik structure in the TES region is more difficult to interpret due to the likelihood that surface NMR is unable to detect some liquid water due to the short relaxation times based on interpretation of the PBL data set (Figure ⁶). The lake sediments to 6 m below the water column show no signal, and these should certainly be saturated and liquid below a floating ice lake, and based on the ¹⁴C minimum lake age estimates of 1400 cal years before present (Lenz et al., 2016) and probable thaw rate below lakes (Ling and Zhang, 2003), the talik should be $>15m$ (Creighton et al., 2018). Due to the interpretation challenges, we cannot confidently explain the water content structure of transitional lake L31 (Figure ⁸) by a partially refrozen talik. Based on historical aerial photograph analysis of L31 since 1955 (Jones et al., 2009), this lake has been expanding over time and may have partially drained into an adjacent lake between 2002 and 2007, explaining the occurrence of bedfast ice conditions in recent years (Figure ³). The L31 has been bedfast for only three of the past 11 measured years since partial drainage (a change from 2.2 to 1.4 m lake depth), and this seems unlikely to result in 15 m of permafrost aggradation from the top down. Deeper than 28 m in L31, we expect that a lithologic transition is present with small pores that prevent surface NMR measurement. The observed 26 m depth to a transition between sandstone and shale would be consistent with the observed change in liquid water content in the L31 profile.

Based on estimates of PBL age (Lenz et al., 2016) and transient electromagnetic measurements (Creighton et al., 2018), we expect the talik thickness to be in excess of 60 m. If we assume that the talik below PBL is $>60m$ thick with unconsolidated and sandstone substrate (Figure ²), this would suggest that geology plays a role in what water we can see because surface NMR would likely be unable to measure water with very short relaxation times (Walsh and Ferre, 2011) in underlying geologic formations. Observed $T2*$ relaxation times in other sandstone units are in a range from 0.064 (Hein et al., 2017) to 0.080 s (Legchenko et al., 2002) indicating that the sandstone lithology should not preclude surface NMR measurements. This corresponds with our measured $T2*$ of 0.026 s 7.6–14.7 m below the surface at PBL. Porosity deformation has been observed in North Slope sandstone units (Smosna, 1989) that could result in smaller pore sizes and reduced $T2*$ values; however, deformation processes are normally associated with deeper rock. The DRP borehole (19 km north–northeast) reports that shale may be present starting 26 m below the surface. Laboratory NMR relaxation measurements on water-saturated shale have found $T2$ relaxation times on the order of 0.01 s (Borysenko et al., 2009). Because $T2*≤T2$ given a surface NMR dead time of approximately 0.005 s, it is highly likely that no signal could be detected in shale formations even if liquid water is present. The discrepancy between the 15 m depth maximum vertical extent of unfrozen water (Figure ⁶) and the 26 m depth to shale in the DRP log may be explained by variation in structure in the 19 km distance between the lake and the borehole. Shallow cores at PBL describe the mineral component of lake sediment texture as silty (Lenz et al., 2016), suggesting that saturated sediments are present with relaxation times too short to observe with surface NMR. The subtracted result for PBL reveals relaxation times in the top 5 m of 0.012 s — although these are dubious given the small observable water content, they nonetheless suggest that short relaxation times are present that cause the measurement to miss a substantial volume of unfrozen water. Given the close proximity of L31 (12 km east-northeast from PBL), this explanation is likely the same for the upper 15 m of that water content profile (Figure ⁸).

The surface NMR FID measurement may be influenced by heterogeneous regional or pore-scale magnetic fields that can lead to incorrect estimation of the water content and relaxation time (Grunewald and Knight, 2011). We examine the definitive geomagnetic reference field-corrected magnetic field anomalies in our two study regions (Figure ¹), and we observe that there are no obvious strong gradients that would affect our data as would be expected for a region underlain by thick sedimentary packages. Drillers’ logs for DRP report the presence of “ferromagnesian minerals” (i.e., iron- and magnesium-rich) at depths $>46m$ (Brockway, 1983) that could induce pore-scale magnetic field gradients that would lead to nonexponential relaxation and overestimation of water content. False aquifers have been interpreted when surface NMR excitation frequencies are transmitted away from the true subsurface water Larmor frequency (Legchenko, 2005). These off-resonance effects are anticipated for 0.04 s excitation pulses if the frequency offset is $>5Hz$ (Walbrecker et al., 2011) or in some cases, less than 5 Hz (Grombacher and Knight, 2015). For shorter excitation pulses, the off-resonance effect is reduced (Grombacher and Knight, 2015). For each of our data sets with reliable liquid water NMR signals (i.e., floating ice and transitional ice lakes only), we observe frequency offsets to be $<3.3Hz$ (Table 1). Given the relatively short excitation pulse of 0.02 s used and the small frequency offsets, we anticipate that the off-resonance transmission has not introduced artifacts into our results.

We observe that there is a $<0.07m3m−3$ water content layer immediately below the water column on floating ice lakes when inverting the unaltered data set. For INDFOX, this layer is approximately 2 m thick (Figure ⁶), and for PBL, it is approximately 6 m thick (Figure ⁶), although due to the smooth inversion, these values are estimates. Previous observations on thermokarst lakes (Parsekian et al., 2013) used a “blocky” minimum-layer inversion;therefore, they did not resolve the low-water-content layer immediately below the lake. Three explanations for the observed low-water-content layer are (1) very fine particle sizes that cause relaxation times so short that the surface NMR dead time prevents accurate resolution of total water content, (2) magnetic properties of the sediments that increase the NMR relaxation rate (Keating and Knight, 2007), or (3) vertical magnetic field gradients that we cannot measure.

Although we know that large-scale gradients in earth’s magnetic field are absent at our field sites (Connard et al., 1999; Figure ¹), the pore-scale magnetic properties of lake sediments may still influence NMR relaxation (Grunewald and Knight, 2011). Lakes in interior Alaska have reported magnetic susceptibility values ranging from approximately 25 to $250×10−5SI$ (Begét et al., 1990) up to 16 m below the lake bottom, with most values falling in the range of $25−75×10−5SI$⁠. This range would be considered intermediate magnetic susceptibility, and in large pores, there would be little effect on the surface NMR signal (Grunewald and Knight, 2011). This leads to the possibility that the low apparent water content of the layer immediately below the lake is caused by the intermediate magnetic susceptibility and small particle/pore sizes. We cannot test this directly with the available data; however, sediment cores between the shore and the center of the lakes could illuminate the spatial distribution of the particle texture and the measurement of susceptibility. The Carr Purcell Mieboom Gill or spin echo pulse sequence would be effective for eliminating any possible magnetic contribution to relaxation (Grunewald and Walsh, 2013), and although we did not collect multipulse data during this campaign, it would be interesting to explore in future work to image lake sediments.

Lakes and ponds have been the target of past surface NMR investigations for validation of inverted parameters (Müller-Petke et al., 2011), thermokarst lake talik measurements (Parsekian et al., 2013), artificial infiltration (Walsh et al., 2014; Falzone and Keating, 2016), and mountain hydrology (Liefert et al., 2018). Although the NMR signal between saturated sediments and water column is a contrast of approximately $3×$ (e.g., sediment porosity of approximately $0.3m3m−3$⁠; $bulk water=1.0m3m−3$⁠), the influence of long NMR relaxation signals from bulk water to mask shorter relaxation signals of sediments has been demonstrated (Walsh et al., 2014; Falzone and Keating, 2016). Although the lake water signal is desirable for some validation purposes (Müller-Petke et al., 2011), many cases require only measurement of the sublake environment, and they therefore attempt to reduce the influence of bulk water in the inversion. For example, Parsekian et al. (2013) constrain for the water column in a blocky inversion, whereas Falzone and Keating (2016) present several approaches including subtraction of forward modeled data. Liefert et al. (2018) work on shallow lakes and were therefore able to make all observations when the ice conditions were bedfast.

As with all inverted surface NMR data, the water content and $T2*$ profiles are based on nonunique model fits to the data (Behroozmand et al., 2015). Changing the regularization parameter will result in smoother or more complex models, and although we derived our choice of regularization from L-curve tests, it is not possible to perfectly know what the correct value is. The effect of this is generally either an underestimate or overestimate of the water content and the possibility of subsurface layer structure to be more “complex” than reality. By showing uncertainty estimates (Parsekian and Grombacher, 2015), it is possible to interpret a likely range of water content, as well as zones where the best model may indicate overcomplexity.

One unusual observation related to the lake water in INDFOX (0.8 s) and PBL (0.3 s) is the $T2*$ relaxation time for lake water obtained by optimizing the lake water subtraction. The relaxation time for PBL (Figure ⁵), somewhat shorter than the expected range of 0.6–1.5 s (Schirov et al., 1991), is similar to the field observation of 1.0 s (Müller-Petke et al., 2011). The low $T2*$ for the lake water columns cannot be explained as an inversion artifact because the relatively short signal is clearly seen in the raw data (Figure ⁵). Of the three algorithms proposed by Falzone and Keating (2016) for removal of the bulk water signal, two assume that the lake water has a much longer relaxation than the sediments. Because this is apparently not the case for PBL where the sediments have a relaxation that is only slightly shorter than the lake water in some pulse moments, we are restricted to using only the a priori approach to bulk water removal. Fortunately, this is highly suitable to our case because the lakes are known to have relatively flat bathymetry, and the thickness of the ice and water layers is directly observed at the time of NMR measurement. After suppressing the water signal, the new “raw” data show that the sediment signal for INDFOX and PBL has increasing the initial amplitude and relaxation time between the lake bottom and the pulse moment 0.8 A s for both data sets, after which point the relaxation time and amplitude decrease at larger pulse moments. For INDFOX, the inverted water content result shows a smoother profile with depth that has less erratic behavior than the unsuppressed result (Figure ⁶). The water content profile for PBL shows a smoother, thicker layer of approximately $0.5m3m−3$ between 7.6 and 14.7 m compared with the $0.83m3m−3$ spike at 10.2 m in the unsuppressed result. Because the suppressed inversion no longer has to allow for a maximum water content constraint of bulk water, the optimization is limited to only the range of plausible water content values of sediments. Because there is no surface relaxation component in the lake water, and gradients in the earth’s magnetic field are low (Figure ¹), there is not obvious direct evidence to support an explanation for the observed short signal at PBL. We leave the possibility open that there may be a vertical magnetic field gradient due to the difference in susceptibility between water and the substrate that we cannot directly measure. Based on observations of iron- and magnesium-rich minerals at depths >46 m in drill cuttings (Brockway, 1983), it is reasonable to assume that the substrate has a different susceptibility to the lake water, even though the subsurface is consolidated or unconsolidated sedimentary in nature to $>400m$ depth (Figure ²).

Although there is a growing interest in using surface NMR to study cryosphere targets in the Arctic, it is important to highlight specific limitations of the measurement. The inclination of earth’s magnetic field at a measurement location controls the sensitivity kernel and, therefore, the maximum depth of investigation. Above the Arctic circle (66°34′N latitude), when the magnetic field inclination is 76°–90°, the depth of investigation is reduced. For example, consider the 7 A s pulse moment using a 75 m loop at $B0$ inclination of 80° (i.e., the parameters from this study), the maximum sensitivity depth is approximately 11% smaller in comparison to the same measurement at a location where $B0$ inclination is 40° based on kernel calculations using Müller-Petke et al. (2016). Additionally, the low electrical resistivity of the subsurface on the Arctic Coastal Plain (Figure ²) reduces the depth of investigation of surface NMR. Although the resistivity values are somewhat lower than might be expected for a frozen material, our observations (Table 1) are corroborated by the borehole data (Figure ²).

Besides the lake taliks observed in this study, other zones of unfrozen sediment do exist that are excellent surface NMR targets, such as when saline pore waters depress the freezing point in coastal environments (Keating et al., 2018), or englacial firn aquifers (Legchenko et al., 2018). However, in many Arctic regions, liquid water in sediments and geologic formations may be small, particularly in the top 100 m where surface NMR is useful and permafrost can be present. When porous media are frozen, the freezing point of water is depressed. This results in thin films of liquid water remaining along the grain surfaces due to adsorption forces (Watanabe and Mizoguchi, 2002) and soluble salt concentrations (Tice et al., 1985). In cores from Alaska’s North Slope, Kleinberg and Griffin (2005) demonstrate that formations with total porosity of approximately 0.35 may have liquid-filled porosity of 0.03–0.10 at $−14°C$⁠, corresponding to $T2$ relaxation times of $<0.005s$⁠. Surface NMR is expected to be able to resolve water with $T2*$ relaxation times as short as 0.006 s (Walsh and Ferre, 2011) — given that $T2*≤T2$ for any given formation, the short relaxation times effectively preclude measurement of unfrozen water in porous permafrost materials, consistent with our observations at terrestrial control sites (Figure ⁴). Even sites with pingos (groundwater-driven permafrost features) and ideal noise conditions have faced challenges when measuring surface NMR (Yoshikawa et al., 2002), though the early surface NMR instruments had even longer dead times, and it may be worthwhile to revisit attempts to image pingo groundwater sources with the most modern low-dead-time NMR instruments. Additionally, isolated taliks related to drained thermokarst lakes (Mackay, 1997), rivers (Minsley et al., 2012), or disturbance (Yoshikawa et al., 2002; Nossov et al., 2013) are likely to be measurable using surface NMR and future research related to these features is encouraged. In these scenarios, where liquid water content is low and relaxation times are short, logging NMR methods are more appropriate than surface NMR soundings. This has recently been demonstrated in shallow boreholes focused on active-layer and near-surface permafrost (Minsley et al., 2016; Kass et al., 2017), but it would also be plausible to install deeper borings with nonmetallic casing in permafrost, or even temporary borings in lake sediments.

We found that floating ice lakes had detectable liquid pore water below them, whereas bedfast ice lakes showed no signal similar to the terrestrial permafrost control sites. The transitional ice lakes had variable vertical thaw profiles that we attribute to periods of bedfast ice and floating ice that cause partial thaw. Although we identified a talik below floating ice lakes, we did not interpret a talik thickness due to the limited depth of investigation of the surface NMR measurements given the relatively conductive substrate and steep magnetic field inclination. These results support the hypothesis that lake ice regimes exert controls on thaw dynamics of sublake permafrost, although we also highlight the need for development of a larger geophysical data set to fully elucidate this relationship. Our results indicate that when characterizing sublake permafrost thaw using remotely sensed lake ice regime (i.e., SAR classification), the appearance of a lake as bedfast does not necessarily mean that permafrost is present immediately beneath.

This work was supported by NSF award #1417345 with additional funding from the U.S. Geological Survey, the U.S. Fish and Wildlife Service, and the Bureau of Land Management. We thank the associate editor and three anonymous reviewers for their helpful comments that substantially improved the manuscript. We thank E. Babcock, M. Engram, A. Bondurant, R. Daanen, and A. Gadake for help with the field work.

The data used in this analysis are publicly available at the NSF Arctic Data Center, http://doi.org/10.18739/A20K35 and https://doi.org/10.18739/A21J9771X.