Groundwater Flow Through Continuous Permafrost Along Geological Boundary Revealed by Electrical Resistivity Tomography

In continuous permafrost regions, pathways for transport of sub‐permafrost groundwater to the surface sometimes perforate the frozen ground and result in the formation of a pingo. Explanations offered for the locations of such pathways have so far included hydraulically conductive geological units and faults. On Svalbard, several pingos locate at valley flanks where these controls are apparently lacking. Intrigued by this observation, we elucidated the geological setting around such a pingo with electrical resistivity tomography. The inverted resistivity models showed a considerable contrast between the uphill and valley‐sides of the pingo. We conclude that this contrast reflects a geological boundary between low‐permeable marine sediments and consolidated strata. Groundwater presumably flows toward the pingo spring through glacially induced fractures in the strata immediately below the marine sediments. Our finding suggests that flanks of uplifted Arctic valleys deserve attention as possible discharge locations for deep groundwater and greenhouse gases to the surface.

the thaw-protecting active layer (Mackay, 1998). By definition, this pingo will be of the open-system type because it is fed by groundwater not enclosed by permafrost (Liestøl, 1996). Pingos persist for as long as permafrost conditions remain, and even after the through-talik has potentially frozen over and the spring discharge has ceased. Consequently, open-system pingos indicate current or previous presence of active through-taliks (Yoshikawa, 2013).
Both active through-taliks and open-system pingos require artesian pressure in the sub-permafrost groundwater system (French, 2017). In areas of continuous permafrost, such pressures may be produced by recharge from meltwater infiltrating the ground below glaciers with basal melting (e.g., Liestøl, 1977;Scheidegger & Bense, 2014) or, where permafrost is relatively young, by freezing expansion associated with basal permafrost aggradation . While artesian pressure is a prerequisite for the transport of deep groundwater toward the surface, a sufficiently hydraulically conductive pathway is also needed. Permeable geological units (e.g., Haldorsen et al., 1996) and faults (e.g., Rossi et al., 2018;Scheidegger et al., 2012;Scholz & Baumann, 1997;Wu et al., 2005) comprise current examples of such migration pathways.
In Svalbard, many pingos are found along valley flanks (Humlum et al., 2003), and several of these occur where no links to hydraulically conductive geological units or faults are known (Ballantyne, 2018). We propose that a combination of low-permeability Holocene marine sediments and underlying bedrock fractures resulting from pre-Holocene glacial loading and unloading may constitute a previously overseen explanation for springs located at valley margins. Figure 1 illustrates our conceptual model for spring formation at valley margins with cross-sections of the side of a typical glacially cut valley on Svalbard ranging from glaciation to present day conditions. During the various glaciation cycles, glacial loading and unloading has caused ground compression and decompression along with fracturing (Figures 1a and 1b; e.g., Neuzil, 2012). Glacial fracturing of the subsurface is likely to be most abundant within the valleys, because of the greater pressures generated there (Leith et al., 2014a(Leith et al., , 2014b. Following deglaciation, low-permeability marine and deltaic sediments are deposited on top of the fracture zone (Gilbert et al., 2018), confining groundwater flow (Figure 1c). Given the right conditions, a spring forms at the end of the hill slope (Figure 1d and Fitts, 2002) when the sea retreats. In Late Holocene, temperatures drop to form continuous permafrost (Humlum, 2005;Mangerud & Svendsen, 2017), but the ground stays unfrozen below the spring site due to hydrological advective heat transfer ( Figure 1e). As permafrost thickness increases, the active through-talik forms along the fractured zone, because it comprises the most hydraulically conductive pathway (Figure 1f  , respectively, cause compression and decompression of the ground that results in fracturing (Leith et al., 2014a(Leith et al., , 2014b. The fractures produced this way are more abundant below valley bottoms. (c) Low-permeable marine and deltaic sediments are deposited in the fjord valley (Gilbert et al., 2018) constituting a low-permeability cover on top of the conductive fracture zone. (d) After relative sea-level fall, a spring forms at the end of hill slope. (e) Continuous permafrost forms, but the ground stays unfrozen below the spring site due to advective heat transfer and latent heat release during freezing of spring water within the pingo. (f) Comprising the most hydraulically conductive pathway, groundwater is transported toward the spring along the fractured layer.
Surface-based electrical methods have been widely used to map and characterize frozen and unfrozen ground in permafrost environments . In most locations, frozen ground can be expected to have a significantly higher electrical resistivity (>1,000 Ωm,  than unfrozen (<500 Ωm, Palacky, 1988). However, clay-rich and saline permafrost environments may possess significantly lower resistivities. Frozen clay and other fine-grained sediments can host microfilms of unfrozen water even at temperatures below −5°C (Scott et al., 1990) and show electrical resistivities below 100 Ωm (Harada & Yoshikawa, 1996;Keating et al., 2018;Minsley et al., 2012). Upon ground freezing, groundwater brine exclusion takes place as solutes are expelled to the residual water (Cochand et al., 2019). Saline, unfrozen, and electrically conductive groundwater may occur as microfilms within frozen ground (Keating et al., 2018) or as larger inclusions (i.e., cryopegs; Gilbert et al., 2019;Gilichinsky et al., 2003).
We investigate the above conceptual model by elucidating the geological and hydrogeological context at the margins of a valley-flank, active open-system pingo by measuring the electrical resistivity in the ground.

Study Site
The study site was Førstehytte Pingo (FHP), one of five open-system pingos in Lower Adventdalen, found in central Spitsbergen, the biggest island in the Svalbard archipelago ( Figure 2a). As for the rest of Svalbard, continuous permafrost dominates Adventdalen due to a cold and dry climate. Permafrost thicknesses range from <200 m in the valley floor to >450 m in the adjacent mountains (Christiansen et al., 2005;Humlum et al., 2003;Liestøl, 1977). With one exception, all five pingos are active and perennially discharge brackish (the spring at FHP has a Cl − concentration of ∼1,100 mg L − ), methane-rich waters in orders of 10 −1 L s −1 Hornum et al., 2020). The two most up-valley pingos, Innerhytte and River pingos, have formed in fractured shale and their positions are likely explained by an underlying fault (Figures 2a and Rossi et al., 2018). Moving westwards into the lowest part of Adventdalen, FHP is the first of three pingos (the other two being Longyear (LYRP), and Lagoon (LP) pingos) that all have formed in Holocene marine muds (Yoshikawa & Harada, 1995). All three pingos align with the Northeastern flank of Adventdalen and thus locate close to the boundary of well-consolidated sedimentary rock that outcrops on the mountainside ( Figure 2b). The elongated shapes of LP and FHP are both parallel with the aforementioned alignment.
Below the valley floor of Adventdalen, a succession of Late Weichselian to Holocene glacio-marine and deltaic sediments up to 60 m thick overlies well-consolidated rocks of Cretaceous age or older that dip in a southwestern direction ( Figure 2, Gilbert et al., 2018). Together, all units comprise a groundwater system with a very low permeability, and most fluid flow is restricted to fractures in the consolidated strata. Such fractures are pronounced in the Festningen Sandstone member (Figures 2b and Olaussen et al., 2020) and in the consolidated sedimentary rock immediately below glacio-marine succession (Figures 2b and Gilbert et al., 2018). The latter fracture zone is likely of glaciogenic origin (e.g., Neuzil, 2012).

ERT -Data Collection and Processing
Direct current resistivity measurements were conducted near FHP. Four 2D electrical resistivity tomography (ERT) surveys were implemented using the Wenner-α configuration (cf., Reynolds, 2011) during threeweeks in September 2017 (Figure 3). At this time of year, the thawed active layer allowed for easy installation of electrodes and good electrical connectivity with the ground. The ERT surveys were performed with an ABEM-SAS-1000 Terrameter coupled with an ABEM-ES10-64 Electrode Selector. The layout for a single survey consisted of four 100-m-long cables in a roll-along configuration, each with 21 electrode takeouts. Only uneven electrode take-outs were used and the last takeout on a cable was aligned with the first takeout on the subsequent one so that the combined cable was 400 m long and connected to 41 stainless steel electrodes with 10 m spacing. All possible four-electrode Wenner-α configurations were measured in both normal and reciprocal mode to assess measurement error (Binley & Kemna, 2005;Kim et al., 2016). This resulted in 260 unique electrode configurations and a maximum of 520 measurements for each line. Fewer measurements were available when the instrument could not connect to all electrodes and therefore failed to measure some electrode configurations. The current injected to the ground varied between 200-1,000 mA. Thirteen surveys were carried out along the four transects covering most of the pingo margin (Figure 3). At each transect, two to four surveys were undertaken and provided 300 m overlap between consecutive surveys. This resulted in total transect lengths of 500, 600 or 700 m, respectively comprising 375, 490 and 605 unique electrode configurations.
Electrode positions were mapped with a handheld GPS device (Garmin GPSMAP® 76C). When measuring the coordinate position within a limited time (<1 h), this device showed to have a relatively high precision (<0.1 m) but low accuracy (<2 m). We adjusted for the low accuracy by noting particular electrodes, whose locations could be accurately pinpointed on the orthomap ( Figure 3) and translated the coordinates accordingly. Because of the relatively poor vertical precision of handheld GPS measurements, we inferred the topography along the ERT lines by projecting the electrode positions on a 5-m-resolution DEM of the field area (not shown, Norwegian Polar Institute, 2020).
To ensure good quality of the resistivity data used for the inversion, we first performed statistical data cleaning. The final product of this pre-processing was four files, one for each transect, containing up to one measurement for each unique electrode configuration. Details of the data pre-processing can be found in the supporting information (Text S1).
HORNUM ET AL. Overview of Lower Adventdalen that shows the location of pingos, the study site (red square), and the Holocene marine limit. Topographic data used to create the map by courtesy of Norwegian Polar Institute (2020). (b) Geological cross-section across Adventdalen and the study site. Fractures in the sandstone unit and below the succession of glacio-marine and deltaic sediments interrupt the dominant low-permeability of the groundwater system. Crosssection modified from Hodson et al., 2020. 2D inversion of the measured apparent resistivities were carried out using the graphical user interface of ResIPy 3.0.1, an open-source software for inversion and modeling of geoelectrical data (Blanchy et al., 2020) that builds on the R2 code (version 4.02, Binley, 2019) for the inversion of DC resistivities. We employed a triangular mesh for the inversion. The mesh was composed of a fine mesh that defined the region of the final resistivity model encompassed by a coarse mesh. The transect length defined the lateral extent of the fine mesh and the coarse mesh extended five times to both sides. The fine-to-coarse mesh boundary was at 50 m below ground level (m b.g.l.) and the coarse mesh extended to a depth of 30% the total lateral mesh extent. The resolution of the fine mesh was defined by a characteristic length of 4.38 and a growth factor of 4. This resulted in fine meshes with 1,705, 1,490, 1,582 and 1,741 triangles for transects A, B, C and D, respectively. We used the inversion type "normal regularization with linear filtering" and the convergence criterion was defined by a root-mean-square error of <1.2%. The certainty of the electrical resistivities predicted by the inversion was quantified by the depth of investigation (DOI) method (Oldenburg & Li, 1999). This technique uses the difference of inverted resistivity models with varied starting resistivities to calculate a DOI-index map. We used a two-sided difference scheme with starting resistivities 0.1 and 10 times the logarithmic mean of the apparent resistivities. DOI values above 0.2 indicate that the model is weakly constrained by the measurements (Oldenburg & Li, 1999). Similarly, lower values indicate that inverted resistivities are well constrained by data and allows for greater faith in the resistivity model.
All measured and inverted electrical resistivity data resulting from this research is public available from the Zenodo repository (Hornum, 2021 Figure 4 shows the electrical resistivity models produced by the inversion of the measured values and the DOI-index predicted for these models. To facilitate further spatial understanding, we also produced a 3D animation, which is available as supporting information (Movie S1). In addition to the resistivity models produced from our own survey, the animation also shows a resistivity model from FHP presented by Ross et al. (2007).

Results
The electrical resistivity models predicted significantly varying values and patterns at different sides of FHP. Based on the differences of the predicted resistivity values, we divided the transects into three segments (I, II, and III), which are summarized in Table 1 and described in detail below.
Segment I covers transects A, B, and the eastern part of Transect C and is situated between FHP and the mountainside. The resistivity model shows that the subsurface here is generally characterized by high resistivity values that range from 1,000 to 5,000 Ωm. Relatively large and elongated zones up to ∼200 × 60 m (width × height) of very high resistivities (5,000-50,000 Ωm) are also common, but these do not extend to depths shallower than ∼10 m b.g.l.
Segment II possesses the most complex resistivity pattern of this survey. This segment locates south of the southeastern end of FHP and covers the western part of Transect C and the southeastern part of Transect D.
In approximately the deepest 15-25 m, Segment II is characterized by low resistivity values that range from 20 to 100 Ωm. A relatively sharp boundary (<5 m) marks the transition to a lateral zone of moderate to high resistivity values (500-5,000 Ωm). This resistivity range generally dominates the shallowest 25-35 m of the subsurface, but not at the boundary to Segment II (Transect C), where low resistivity values extend to near the surface. The moderate to high resistivities are distributed in a heterogeneous way, and vary between the extreme ends of the range at several points along the extent of Segment II.
Moving on to Segment III and the southwestern flank of FHP, a further decrease in the ground resistivity can be observed. Segment III covers the northwestern part of the Transect D and locates between FHP and the valley center. Low resistivity values of up to 50 Ωm characterize the lower part of Segment III and gradually decrease upwards to very low resistivity values of down to 1.8 Ωm. A sharp boundary can be observed in the shallow part of Segment III toward Segment II in the form of the contrasting resistivities, but at greater depths the resistivity values are close to identical in both segments.
Showing mostly values below 0.2, the DOI-index maps on Figure 4 indicate that the majority of the inverted resistivity values are relatively well constrained by the measurements. Segment III forms an exception to this pattern by mostly showing higher DOI values except for in the lower part. The predicted resistivities in Segment III were thus generally not well constrained by the measurements.

Discussion -Implications of Resistivity Models
Indicating a robust inversion, DOI values below 0.2 dominated the majority of the resistivity models (Figure 4) and suggested that most of the predicted values represent true ground conditions. However, the higher DOI values dominating Segment III indicated that the resistivity values predicted here were not well constrained by the measurements. DOI values below 0.2 are indeed also present in the lower part of Segment III, but as these are situated below consistently high DOI values, the entire segment should be interpreted with greater caution. When low resistivity values dominate shallow ground conditions, the depth of current flow is reduced and measurements are less sensitive to deeper layers of the subsurface (Binley, 2015). For Segment III, this implied that the predicted low resistivities may conceal zones of higher resistivities. To quantify the potential concealment, we conducted a series of forward modeling experiments with ResIPy, which are described in detail in the supporting information (Text S2). From these experiments, we conclude that low resistivities dominate at least the shallowest 15 m b.g.l. and likely extent to more than 25 m b.g.l.
The relatively strong differences observed on the resistivity models surrounding FHP (Figure 4) indicate varying conditions in the subsurface. In the following, we consider salinity, lithology and phase of state as possible explanations for these differences.
HORNUM ET AL.   In most other settings, very low electrical resistivities as those of Segment III would be inconsistent with frozen ground (e.g., . However, permafrost drill cores from Adventdalen unambiguously show that the majority of the ground is frozen (Gilbert et al., 2018(Gilbert et al., , 2019. Completely unfrozen ground could therefore not explain the low resistivities in Segment III (Figure 4), because the low resistivities completely dominate the ground, rather than appearing as zones within higher resistivity values. Instead, we attributed the low resistivities of segment III to the Holocene marine sediments of which FHP is also composed (Yoshikawa & Harada, 1995). Although such low resistivities (1.8-50 Ωm, Figure 4) would not be expected for most permafrost environments (e.g., Draebing & Eichel, 2017;Lewkowicz et al., 2011;Sjöberg et al., 2015), they are consistent with previous salinity measurements (1-71 ppt., Gilbert et al., 2019) and electrical measurements of marine sediments in Adventdalen (Harada & Yoshikawa, 1996;Keating et al., 2018;Ross et al., 2007). We infer that an unfrozen saline water content of <5% documented in other parts of Adventdalen (Gilbert et al., 2019;Keating et al., 2018) likely also explains the low resistivities of Segment III.
The high and very high resistivities measured on the other side of FHP (Segment I, Figure 4, Table 1) did not comply with the above explanation. Instead, the modeled values (1,000-50,000 Ωm) pointed to permafrost with a limited unfrozen water content  and as such would be difficult to explain if the ground consisted of the aforementioned marine sediments. We instead interpret the high resistivities to reflect a different lithology, which, given the geological context, is likely to be shale or mudstone (Figure 2b). The significant resistivity range may have resulted from differences in fracture abundance, lithology or ground ice concentration, but borehole calibration or other investigations are needed before an unequivocal interpretation can be made.
Constituting the transition between Segments I and III, Segment II presumably spans a geological boundary. At the same time, this segment passes closely to recent spring locations that may affect subsurface thermal regimes and influence subsurface resistivities. Despite the presence of two different lithologies, the resistivity pattern of Segment II is unlikely to be explained solely by lithogy since, we would expect a less complex pattern with an upwards increase in resistivity. This is because the older lithology is the more resistive. Instead, the moderate to high resistivities distributed heterogeneously throughout the segment are best explained as zones with high ice concentrations related to the accumulation of frozen spring water (as is typical near a pingo) within the marine sediments. Similar discontinuous ice bodies were exposed at the summit of the pingo following a large mass movement before the study was undertaken. This explanation is consistent with ERT surveys by Ross et al. (2007) that documented similar complex resistivity patterns for the internal structures of LP and FHP.
10.1029/2021GL092757 8 of 11 The insert at the bottom shows the location of the transects (see also Figure 3) and the transect from Ross et al. (2007) is also included. The number of iterations and final root-mean-square error are written in the lower left corner of each transect. Based on the observed resistivities, we divided the transects into three segments. A description of these are summarized in Table 1. Very low resistivities in the top gradually increasing to low resistivities at the base.

Table 1
Summary of Electrical Resistivity Patterns Observed on the Resistivity Models (Figure 4) Following the above interpretation, FHP locates exactly at the boundary between the consolidated strata and the marine valley infill. Assuming that this is not a coincidence, the geological boundary must be considered as an explanation for the location of FHP. The conjunction of the pingo and the geological boundary might be explained by groundwater recharge in the highlands discharging at the foothill. However, such an explanation would not be consistent with the high electrical ground resistivities found toward the mountainside, and there would be no obvious reason why the pingos in lower Adventdalen (LP, LYRP and FHP in Figure 2a) do not locate on both sides of the valley. Instead, the conjunction of FHP and the geological boundary is in agreement with the conceptual model presented in Figure 1, wherein glacially induced fractures in the sedimentary strata comprise a hydrological pathway for deep groundwater to reach the surface. Considering the regional southwestern dip of the sedimentary strata (Figure 2b), this view readily explains why the pingos locate on the northeastern margin. The conceptual model is further supported by the resistivities of Segment III that indicated accumulation of frozen spring water at the inferred geological transition, and by the geochemistry of pingo spring waters in Adventdalen, which point to a deep groundwater origin Hornum et al., 2020).
To our knowledge, no other investigation at any of the open-system pingos in Svalbard that are found along valley flanks (e.g., Humlum et al., 2003) have mapped the geological context in detail. Still, we hypothesize that a similar mechanism may also contribute to the formation of some of the open-system pingos found in the Arctic in similar geological contexts. In Adventdalen, Svalbard, this would readily explain, for example, the elongated shapes of LP and FHP and their alignment with the valley flank. Further research into the hydrogeological setting of other spring sites in high Arctic valleys would confirm the appropriateness and representativeness of our conceptual model.

Conclusions
This study is the first to show a direct relationship between a geological boundary and an open-system pingo. The strong electrical resistivity contrast observed between the uphill and valley sides of FHP likely reflects a lithological difference: The high resistivities observed toward the mountainside are consistent with frozen sedimentary rocks with a limited groundwater content, while permafrost with a low but saline content of groundwater explains the low resistivities on the valley-side. Groundwater presumably flows to the pingo springs through fractures in the sedimentary strata induced during glacial loading and unloading. This view is supported by a heterogeneous resistivity pattern suggesting accumulation of frozen spring water at the inferred geological transition; by spring water geochemistry that indicates a deep groundwater origin; and by the consistently high electrical ground resistivities toward the mountainside of FHP, which does not favor groundwater recharge from the mountains above. The numerous pingos on Svalbard that also locate along valley margins are possibly associated with this boundary as well, and if so, these are explained by groundwater flow enhanced by glacial fractures. Our findings indicate that shallow fractures in the Late Weichselian landscape relief may constitute a previously overlooked groundwater pathway. The fracture zone may link deep groundwater systems to the surface, where low-permeable sediments cover this surface. On a circumpolar scale, flanks of uplifted valleys deserve particular attention as possible pathways for subsurface fluid migration.

Data Availability Statement
The measured and inverted electrical resistivity data supporting this research is public available from the Zenodo repository via https://doi.org/10.5281/zenodo.4479529 (Hornum, 2021).