Electrical Resistivity Tomography Investigation of Permafrost Conditions in a Thermokarst Site in Fairbanks, Alaska

1.   Introduction

The presence of permafrost in cold regions plays an essential role in shaping the Arctic landscape through its impact on small and large-scale topography changes. Permafrost shapes the terrain through a complex interplay of freeze-thaw processes and ground stability. From subtle, gradual shifts to dramatic, large-scale transformations, permafrost is an indispensable factor driving alterations in the topography of cold regions ^[1]. Permafrost provides physical support to the ground surface and regulates moisture content, soil temperature, and subsurface hydrology ^[2]. Because of these unique functions, its degradation strongly influences the Arctic and sub-Arctic communities.

Interior Alaska is characterized by discontinuous permafrost with ground temperatures close to 0 ℃ ^[3]. Multiple studies ^[4,5,6,7] suggest that climate changes during the 20th and 21st centuries have resulted in significant degradation of permafrost across extensive regions within the Arctic. This thawing process releases greenhouse gases and causes infrastructural challenges and socio-economic consequences to communities in cold regions. ^[8] hinted that one of the major issues confronting the Arctic region due to climatic warming and thawing of ice-rich permafrost is the creation of thermokarst terrain. Thermokarst is formed where the edges of isolated permafrost are on the verge of thawing. In such locations which consist predominantly of ice (greater than or equal to 50% ice), the thawing of permafrost changes the ice to water, which forms a mud slurry. The result of this process is the degradation of the environment in the form of ground subsidence, the formation of thermokarst topography, and landslides since the mud slurry cannot provide enough support to sustain the weight of the overlying soil ^[8]. Also, we have other factors, such as groundwater movement ^[9], among the major drivers of permafrost instability. The presence of thermokarst significantly changes the ground surface and modifies surface and subsurface hydrological regime. Thermokarst is noticeable in interior Alaska and can significantly affect infrastructure such as roads and highways.

A better understanding of the processes and dynamics of thermokarst development is crucial to determine the appropriate remedial measures to help reduce or prevent the impact in communities located in cold regions. Understanding this requires a detailed characterization of permafrost's physical and thermal conditions and imaging of the subsurface morpho-structures beneath the thermokarst. Moreover, the active layer characterization is also crucial in permafrost studies since most subsurface ecological, biological, and pedogenic activities occur within this shallow zone ^[10]. The use of geospatial and remote sensing techniques in conjunction with computational models and resources has aided scientists in extrapolating point-scale measurements of active layer thickness using empirical and statistical approaches ^[10,11]. Research by ^[12], however, has shown that the remote sensing approach cannot characterize permafrost's fine-scale patterns and transient features. Furthermore, to gain a comprehensive understanding of frozen ground, researchers use field boreholes equipped with instruments for measurements. Drilling in frozen terrain can be costly and challenging, particularly in remote areas ^[13]. Although boreholes offer valuable subsurface information at specific locations, their invasive nature and limited spatial distribution are limitations compared to geophysical methods. Therefore, noninvasive geophysical investigation methods provide non-discrete spatial data that are ideal in permafrost investigation. These methods are cost-effective and allow for relatively quick data collection compared to boreholes.

Geophysical methods in permafrost studies and subsurface characterization can detail the spatial variations in permafrost conditions and properties vital for understanding the seasonal variations and other transient features that reflect permafrost degradation and thermokarst formation ^[13]. Researchers have made significant advancements applying electrical resistivity tomography (ERT) in mapping the distribution of permafrost layers, identifying ice-rich zones and ice-wedges, detecting variations in thaw depth, and delineating the boundaries of permafrost features ^{[14,15,16,17,18]}.

While ERT has demonstrated its effectiveness in permafrost investigations, its application to studying thermokarst dynamics in cold regions remains relatively limited. In this study, we addressed this knowledge gap by utilizing 2D and 3D ERT grids to characterize the permafrost conditions and subsurface morpho-structures of a thermokarst site (Figure 1) in a discontinuous permafrost region in Fairbanks, Alaska, and determine optimum ERT parameter settings for effective permafrost investigation. Analysis of the ERT data showed prominent resistivity anomalies that correspond to subsurface morpho-structures related to the formation and development of the thermokarst. Among the major resistivity anomalies beneath the thermokarst are very-low resistivity anomalies (lower than 50 Ωm), presumably taliks, and very-high resistivity anomalies (>700 Ωm) that are apparently ice-wedges and frozen ice-rich soils.

Figure 1.  Location map of the study site. Coordinates are in UTM Zone 6 (Base map source: Google Earth).

DownLoad: Full-Size Img PowerPoint

2.   Site description

2.1.   Location and site description

The study area for this research is a thermokarst site located on the University of Alaska Fairbanks (UAF) campus in Fairbanks, Alaska (Figure 1). The campus is situated at approximately 64.8556 N and 147.8341 W, offering a unique opportunity to investigate permafrost dynamics and related phenomena in an accessible environment. UAF has encountered challenges with recurrent thermokarst affecting its roads and parking lots ^[19]. By conducting this study on the UAF campus, researchers can offer valuable insights into the effects of permafrost thaw on the surrounding environment and infrastructure, potentially leading to improved strategies for mitigating or adapting to permafrost-related challenges.

The thermokarst site is located behind UAF's virology lab (VL) and extends into the black spruce forest from the parking lot, crossing an asphaltic surfaced road. This location is susceptible to thermokarst depression formation, with continuous reports of such depressions occurring at the site for over ten years ^[19]. The thermokarst depression at this site (Figure 2) measures approximately 1 m in diameter and about 1.5 m in depth, displaying a distinct steep slope towards the southern margin. The transect along this site exhibits various anthropogenic disturbances, including ski trails, the installation of fences, and drainage pipes.

Figure 2.  Thermokarst on UAF's upper campus site.

DownLoad: Full-Size Img PowerPoint

2.2.   Geological setting

The Fairbanks area exhibits a diverse geological makeup shaped by tectonic and erosional activities ^[20,21]. A combination of metamorphic and igneous rocks characterizes the geology of the Fairbanks area. The basement rocks are comprised of ancient metamorphic rocks such as schist and gneiss, indicating a history of intense heat and pressure. Intrusions of granite and other plutonic rocks are also present, suggesting magmatic activity during the region's geological past ^[22]. Péwé and Prindle et al. ^[20,21,22] have conducted extensive studies, offering detailed cryostructural and stratigraphic data that greatly contribute to our understanding of the geology in the Fairbanks area.

The predominant soil type in the Fairbanks area is redeposited loess and Pleistocene aged silt, which have a thickness ranging from 14–45 m ^[23]. These deposits contain a significant amount of ice, approximately 80–92% by volume ^[23]. Below these deposits are layers of creek gravel that lie on a highly fractured Precambrian Birch Creek schist basement ^[20]. The depth of the bedrock in the Fairbanks area can vary considerably across the region. Historical records indicate that bedrock was encountered at a depth of 52 m during drilling operations for the Arctic Biology building in 1969 ^[24].

2.3.   Climate and vegetation

The climate of the Fairbanks area is characterized as continental, exhibiting distinct seasonal variations and extreme temperature ranges ^[25]. According to historical climate data, the mean annual air temperature in Fairbanks is approximately −2.4 ℃. The coldest months are December and January, with average temperatures ranging from −21.9 ℃ to −19.8 ℃. Conversely, July is the warmest month, with average temperatures reaching around 20 ℃ (https://www.ncdc.noaa.gov/cdo-web/). Temperature extremes in Fairbanks can be particularly severe, with recorded lows as low as −51 ℃ and highs as high as 38 ℃ ^[26]. According to the National Center for Environmental Information, the climatic records as of January 2023 indicate that the lowest recorded temperature occurred on January 3, 2022, with an average lowest temperature of −42 ℃. Conversely, the highest average daily temperature, reaching 23 ℃ was observed on July 3, 2022.

Vegetation found in the study area represents the characteristic sub-Arctic taiga forest of the Alaskan interior. This forest is characterized by white or black spruce trees, which grow above a thin layer of moss, aspens, and birches. The southern section of the site lacks vegetation (parking lot), while a thin moss layer and black spruce trees characterize the northern part. Evidence of degrading ice-wedge polygons along the northern half was noted based on observation of the drunken forest.

3.   Materials and methods

In this study, 2D and 3D-ERT grid images integrated with borehole data were used to characterize the subsurface condition of the thermokarst site. Borehole data on selected locations and literature input on the resistivity of permafrost were used to constrain the interpretation of the ERT images.

3.1.   ERT survey

In recent decades, the application of geophysics has become increasingly valuable in studying permafrost. Advancements in geochronology, physics, and the understanding of physical properties related to permafrost have contributed to a growing body of knowledge in this field ^[27]. Geophysical methods offer advantages in characterizing the subsurface continuously over large areas at relatively low costs ^[28], often with depths of investigation reaching a few tens of meters ^[29]. In permafrost environments, geophysical techniques focus on detecting differences in physical properties between the frozen ground containing ice or unfrozen water and the non-frozen subsurface materials ^[13]. These properties vary based on pore size, ice structure, water or ice content, pore water chemistry, ground temperature ^[30], and overburden pressure.

Geophysical methods like ground-penetrating radar (GPR) and ERT have been extensively used to study permafrost features and the active layer ^[31]. The ERT method has been widely utilized to investigate permafrost features like the active layer, ground ice content, and water content ^{[14,17,32,33]}. ERT is a geophysical method used to image and characterize subsurface structures and properties based on their electrical resistivity. The principle behind ERT is based on the fact that different materials have different electrical resistivity ^[34]. The method involves introducing electrical currents into the ground through pairs of electrodes and measuring the resulting voltage differences between other electrode pairs ^[34]. By analyzing the distribution of electrical potential in the subsurface, it is possible to infer the resistivity distribution and delineate subsurface features. In this study, we used 2D-ERT and 3D-ERT grids to investigate the thermokarst site.

3.1.1.   2D-ERT survey geometry and electrode array

In May 2022, a 2D-ERT survey was conducted to acquire nine ERT profiles around the margins of the thermokarst (Figure 2) using Syscal pro units with 48 and 72 electrodes. The measurements took place during the early stage of the thawing season when substantial runoff was observed. The main objective of the study was to investigate the internal structure and morphology of the thermokarst. Specifically, the focus was on the initial phase of the thawing season, which is known to coincide with significant runoff. This phase is particularly important as it marks the typical development of the thermokarst in the study area. By examining the internal characteristics of the thermokarst during this specific period, the study aimed to gain insights into its structure and morphology, contributing to a better understanding of its formation and behavior.

The profiles were divided into two groups. The first group consisted of seven profiles that extended from south to north (SN), while the second group included two profiles spanning from west to east (WE) along the margins of the thermokarst as shown in Figure 1. In the SN transect, there were seven profiles, with five of them measuring 142 meters in length. These profiles were equipped with 72 electrodes spaced at 2-meter intervals. Profiles 2 and 3 covered a length of 94 meters with 48 electrodes spaced at 2-meter intervals. Also, the WE profiles (Profiles 8 and 9) extended 94 m. Furthermore, the survey transects experienced significant topographic variations, starting from the parking lot, crossing a roadway, and transitioning gradually to a gently sloping topography (Figure 1).

The dipole-dipole (DD) and the Schlumberger (SC) array were utilized in this study because the DD array is more sensitive to resistivity changes in the horizontal direction than in the vertical direction ^[34]. Therefore, it is well-suited for mapping vertical structures, such as the expected contacts between the talik and ice-rich cores of permafrost, but it is not suitable for delineating horizontal structures like the permafrost table and base ^[16]. The SC array combines features of both the Wenner and the DD arrays, offering a moderate level of sensitivity to both horizontal and vertical geological structures ^[34]. Given that permafrost studies may involve both types of structures, the SC array could be a suitable middle compromise between the DD array, which is sensitive to vertical structures, and the Wenner array, which is sensitive to horizontal structures.

3.1.2.   3D visualization of ERT transects

In this study, we utilized the series of seven parallel ERT transects, with varying spacing from 2–7 m. These transects served as the basis for generating 3D ERT images using Surfer Golden software. From this 3D ERT dataset, we then extracted 2D ERT slices at different depths for further analysis. It is important to note that the recommended maximum spacing between parallel and adjacent profiles for generating a pseudo-3D survey for environmental and engineering investigations at large depths is four times the electrode spacing, according to ^[35], and a study by ^[36] showed a 4:1 ratio of the spacing between adjacent profiles to electrode spacing in a profile is adequate to generate an effective pseudo-3D ERT. To visualize and create the 3D ERT data in Surfer, we merged the seven 2D ERT data files into a single ERT file in Microsoft excel. This composite dataset was subsequently inverted using the Surfer program ^[37]. The resulting output comprised horizontal depth slices representing the resistivity profiles, allowing for a comprehensive understanding of subsurface conditions.

3.1.3.   Inversion

The RES2DINV software ^[38] was used to perform the inversion process. This software utilizes an inversion algorithm ^[39] to convert measured resistivity profiles into a resistivity model capable of representing various geological features in both their lateral and vertical distribution. The main objective is to generate a resistivity model that closely matches the observed apparent resistivity values. To ensure the accuracy of the model, the root mean square error (RMSE) ^[40] method, was applied.

The RES2DINV software offers two inversion regularization techniques: L₂-norm optimization and L₁-norm optimization. L₂-norm optimization produces smoother results and is suitable for regions with gradual resistivity changes. On the other hand, L₁-norm optimization is more effective in areas characterized by abrupt resistivity variations, noisy data, and sharp transitions in subsurface materials ^[41]. Given the specific conditions of our study sites, the L₁-norm optimization approach was employed to handle sharp transitions and prevent overfitting of the resistivity model. Multiple damping factors were tested to find an ideal balance between smoothing the data and preserving observed data consistency. The inversion process was accomplished with maximum and minimum damping factors that ranged from 0.01–0.4 ^[42], where the maximum damping factor was five times the minimum damping factor in each test. Given that the resolution of the resistivity survey decreased rapidly with increasing depth, the damping factors were raised for each subsequent deeper layer using the software's built-in automatic calculation feature. Additionally, a mesh refinement to half the electrode spacing was incorporated to improve the resolution of the anticipated resistivity differences between frozen and unfrozen materials. The best inversion setting for each resistivity survey was obtained by testing four different inversion parameter options and assessing the RMSE presented by the final inversion; The four tested options are:

➢ Standard setting: No changes from the program's default settings, including standard least-squares constraints (L₂-norm) and the standard mesh size.

➢ Robust setting: This included using robust constraints for sharp boundary inversion (L₁-norm) at standard cut-off factors while retaining all other settings from the standard setting.

➢ Robust setting with finest mesh: This involved setting the mesh type for the forward modeling sub-routine to the finest setting while retaining all other settings from the robust setting. A finer mesh increased accuracy but also increased computation time.

➢ Robust setting with finest mesh and no side block effect: This involved removing the impact of side blocks on the inversion process, which removed the presence of high or low resistivities on the model's edge. All other settings were kept the same as in the robust setting with finest mesh.

The inversion process terminated when one of the following conditions was fulfilled: the desired misfit between data was attained, or the maximum iteration count was achieved. Using the RES2DMOD software ^[43], the expected apparent resistivities for the simulated data set were calculated using the same electrode configuration used for collecting the actual data. The confidence level of the resistivity models from each simulation was assessed using the depth of investigation (DOI) test ^[44].

3.1.4.   DOI

To assess the accuracy and reliability of the observed features at depth, it is crucial to determine the profile's DOI. Traditional techniques involve computing the maximum signal at depth ^[45] or determining the median depth of investigation ^[44,46] for a homogeneous medium. ^[47], however, have shown that these methods are inefficient for quantifying the investigation depth in a subsurface with heterogeneous high-resistivity features. Instead, we used DOI index calculations to ascertain the depth below which the physical properties of the subsurface no longer affect the data ^[44].

The DOI test is a technique for evaluating the depth of investigation by conducting at least two inversions of the data using different constraints. The concept of DOI generally refers to the depth to which resistivity data are sensitive. This implies that inverted structures are no longer tied to the subsurface data below this depth and must be interpreted cautiously. The DOI method has been applied in various studies, such as investigating mountain permafrost ^[48] and studies of rock glaciers and permafrost in the Bernese Alps, Switzerland ^[49]. In the calculation of the DOI, each model is given a reference resistivity, computed as the average of the logarithm of the apparent resistivity values. The result of this inversion is the reference resistivity model (m₁). The reference resistivity is then altered by assigning each cell a value one-tenth or ten times the original reference resistivity.

In some cases, even more drastic changes of one hundredth or one hundred times the reference resistivity are necessary. Making the reference resistivity more conductive or resistive depends on the properties of the subsurface features of interest. If the targets are resistive, the reference resistivity should be more conductive, and if they are conductive, the reference resistivity should be more resistive. The output of the inversion process is referred to as the altered resistivity model (m₂) and comparing it to the reference model constitutes a one-sided test. In instances where the subsurface contains resistive and conductive targets (such as alternating ice-rich and ice-poor zones or ice wedges surrounded by ice-poor permafrost), two altered resistivity models (one made more conductive and the other more resistive) are compared, which is referred to as a two-sided test. The two-sided test was conducted for all transects in this study. Equation 1 is an example of a DOI calculation for a one-sided test:

where ${m}_{1}\left(x, z\right)$ and ${m}_{2}\left(x, z\right)$are the model cell resistivities obtained from the first and second inversions, respectively, and ${m}_{r1}$and ${m}_{r2}$ refer to the resistivity values assigned as the references. When the altered resistivity models match the reference resistivity models, the DOI value tends to approach zero, meaning that the resistivity of the cell constrains the data well. On the other hand, if the model cells contain limited information about cell resistivity, the DOI value tends to approach 1, indicating that these regions are unreliable. ^[44] proposed a cut-off value of 0.1 as a cautious threshold. In the inversion process, a reference model was used along with a self-damping factor (α s) of 0.1 to control its effect on the inversion process. The horizontal and vertical damping factors (α x and α z), respectively, were set to 1.0, ensuring that the inversion process gave equal weight to both horizontal and vertical features. The mesh for the finite-element calculation of the DOI matrix must extend to the sides of the survey area and reach a depth of about 3–5 times the median DOI for the array type utilized in the data set. The pseudo-depths for each layer are based on the median DOI for various n-factors and array lengths ^[46]. Similar to ^[46], the top layer was set to half the median depth of investigation for the shortest electrode spacing, and each subsequent layer was made 10% thicker to account for the decreasing resolution with depth ^[50]. The mesh's geometry and self-damping factor were determined to be accurate since the DOI values at the bottom of the models exceeded 0.2. To verify this, the scaled DOI test proposed by ^[44] was applied to all simulations. The scaled DOI is calculated using Equation 2, where DOI_max is the maximum DOI value obtained from Equation 5. After computing the scaled DOI values, the results were interpolated using the kriging technique available in Surfer and superimposed on the unaltered resistivity model.

3.2.   Boreholes

Following the completion of the ERT data collection, a subsequent phase of the investigation involved drilling four boreholes along Profile 4 of the VL. The objective of the borehole drilling was to obtain direct subsurface samples for further analysis and characterization. Two drilling techniques, hand-held drilling and core drilling (Figure 3), were employed during this phase. Hand-held drilling was performed using an auger, which effectively crushed the drilled material into fines. Core drilling utilized a diamond-tipped drill bit, enabling intact core samples to be extracted. Both drilling techniques were implemented at various locations of interest along Profile 4. The drilling operations extended to a depth of approximately 5.5 m. At intervals of 0.5 m during the drilling process, the extracted material was carefully collected, analyzed, and logged for detailed examination. The analysis of the drilled material encompassed various parameters, including lithology, color, texture, and moisture properties. These borehole data were used to constrain the interpretation of the ERT profiles.

Figure 3.  Borehole drilling along Profile 4: a) hand-held drill (January 15, 2023) and b) core drilling (March 9, 2023).

DownLoad: Full-Size Img PowerPoint

4.   Results

4.1.   Optimum ERT parameter settings

The optimum parameters for inversion modeling were determined for all resistivity profiles following the procedure described in Section 3.1.3. The summarized results obtained using the robust with finest mesh setting and no side block effect method for the DD and SC arrays are presented in Figures 4 and 5, respectively.

Figure 4.  Robust setting with the finest mesh and no side block effect for the DD array for Profiles 4 and 7.

DownLoad: Full-Size Img PowerPoint

Figure 5.  Robust setting with the finest mesh and no side block effect for the SC array for profiles 4 and 7.

DownLoad: Full-Size Img PowerPoint

Among the four tested parameter combinations, the robust parameter setting with the finest mesh and no side block effect provided the best fit for the observed profiles and for further interpretation. While the other robust setting models produced satisfactory results, the standard least-squares-constrained model did not perform satisfactorily. The presence of anomalies within the area of interest can be distinguished better using these settings. The selected mesh setting (finest mesh) and the use of the robust model and data regularization resulted in a much lower RMSE value (Tables 1 and 2) when compared to the other models generated. For this study, the robust setting was the preferred choice due to the sharp boundary interfaces in the measured target.

Table 1.  RMSE for Profile 4.

Inversion parameter	RMSE DD (%)	RMSE SC (%)
Standard setting	12.5	6.6
Robust setting	6.6	4.4
Robust setting with finest mesh	6.6	4.4
Robust setting with finest mesh and no side block effect	6	4.1

---

 | Show Table

DownLoad: CSV

Table 2.  RMSE for Profile 7.

Inversion parameter	RMSE DD (%)	RMSE SC (%)
Standard setting	3.7	1.7
Robust setting	1.6	0.66
Robust setting with finest mesh	1.6	0.65
Robust setting with finest mesh and no side block effect	1.4	0.34

---

 | Show Table

DownLoad: CSV

4.2.   2D-ERT profiles

Figures 6 and 7 summarize the inversion results and DOI contour for the SC array, respectively. The intensity of the dark grey colors in Figure 7 was adjusted based on DOI values > 0.1. Dark grey colors indicate a high DOI index (DOI value > 0.1), which implies less reliable inversion results. All resistivity models generated contain areas where the DOI values exceed 0.1, either located at the right and left ends of the profiles or in proximity to conductive and resistive structures in the subsurface. Notably, there is a distinct area with high sensitivity to DOI in the central margin of all the profiles (anomalies T1 to T7), particularly beneath the conductive structure. This presents a challenge in obtaining reliable information about the materials in those locations. In Profile 1, this low-resistive anomaly (T1) may not extend to an elevation of approximately 170 m, and the data at depth does not provide a clear indication of the vertical extent of this structure. Similar observations can be made for Profiles 6 and 7, where the resistivity of the structure beneath this conductive body cannot be confirmed for both arrays. Towards the northern margin (80–130 m horizontal distance) of Profiles 1, 4, and 5 (anomalies I1, I4, I4-1, and I5), isolated structures of high resistive bodies (>5000 Ωm) are present. Data cannot confirm this anomaly's vertical extent since the bottom falls within the 0.1 contours.

Figure 6.  Inverted resistivity plots for the SC array along the NS transects.

DownLoad: Full-Size Img PowerPoint

Figure 7.  DOI contours for the SC array along the NS transects.

DownLoad: Full-Size Img PowerPoint

The inversion results for the WE transect (Profiles 8 and 9) are shown in Figure 8 and DOI in Figure 9. The ERT results revealed interesting findings regarding the subsurface characteristics. Specifically, in the section spanning from 35–91 m on Profile 8, the top 2 m exhibited high resistivity values. In contrast, Profile 9 demonstrates low resistivity values within the upper 2 m. Below this zone, a uniform resistivity pattern was observed across the two profiles. The resistivity values in this region remained consistently low, measuring approximately 90 Ωm. Furthermore, in Profile 9, a high resistive anomaly was identified towards the western margin, specifically between 20 m and 30 m mark horizontal distance. Small patches of low resistive anomalies were present across the cross section for the two profiles.

Figure 8.  Resistivity profiles for the WE transect.

DownLoad: Full-Size Img PowerPoint

Figure 9.  DOI contours for the WE transect profiles.

DownLoad: Full-Size Img PowerPoint

4.3.   3D ERT grid and 2D-ERT depth slices

The pseudo-3D ERT grid generated from the seven parallel ERT profiles provides valuable insight into our study site's subsurface features and patterns. Several patterns and features can be identified by examining the grid's horizontal depth slices in Figure 10. The ERT results reveal a heterogeneous pattern of high resistivity features at the shallowest depth slice (0.5–2 m). Moving deeper into the subsurface, the mid-depth slices (4–7.5 m) show a gradual decrease in resistivity and thin patterns of low resistive anomalies within the subsurface. Also, the area of the thermokarst displays a similar resistivity pattern. Specifically, a significant area featuring low resistivity values (10–50 Ωm) is noticeable around the thermokarst (7.5–11 m slice) margin. The north area of the thermokarst (slice 10 m and 11 m) is characterized by a high resistivity plunging feature (>1000 Ωm). Resistivity decreases along the entire margin of this depth slice. Another notable feature observed in the Pseudo-3D ERT grid is a continuous layer of high resistivity anomaly extending horizontally across the study area. This layer is present in multiple horizontal depth slices (0.5–4 m) and may correspond to a feature identified in the 2D ERT grid underneath the roadway.

Figure 10.  A grid of 2D resistivity slices in the horizontal plane at depths (X). The red circled area is the location of the existing thermokarst.

DownLoad: Full-Size Img PowerPoint

Overall, the 3D ERT grid provided valuable information about the spatial distribution of resistivity and highlights the presence of distinct features within the subsurface that can assist in further interpretation of the profiles.

4.4.   Boreholes

The stratigraphic sequence observed within the four boreholes on Profile 4 is similar. The sequence begins with a layer of organic material approximately 0.5 m thick. Below the organic layer is a dark-brown, fine-grained, silty soil containing a significant amount of mica. This sediment layer extends to at least 5.5 m. In two boreholes, BH-2 and BH-4, ice-rich material and ice wedges were encountered at 3.8 and 3.9 m, respectively. These ice features indicate the presence of frozen conditions within the subsurface. BH-3 intercepted saturated soil at a depth of approximately 2 m. This indicates the presence of a layer with high water content. The stratigraphic column for BH-1 generally showed dry, silty soil below the organic layer. Figure 11 summarizes the stratigraphic columns of the borehole logs, outlining the layers encountered at different depths in each borehole. These borehole results were used to constrain the interpretation of the ERT profiles.

Figure 11.  Stratigraphic cross-section of the boreholes.

DownLoad: Full-Size Img PowerPoint

5.   Discussion

The ERT surveys yielded compelling evidence that ERT has the potential to map variations in permafrost properties in the subsurface; however, the efficiency of differentiating between various frozen ground features and properties is contingent on the selection of the most suitable inversion parameters. Each simulation carried out a DOI test to measure the data quality. The areas with a DOI value greater than 0.1 contours were considered unreliable. Based on the results from this investigation, we limited the interpretation of the inversion results to areas of the model that were deemed reliable (with DOI values < 0.1).

The Fairbanks region typically has an active layer ranging from 0.6 to 2 m below the surface ^[15]. To interpret the ERT sections, we relied on approximate resistivity scale values from previous studies in the area and borehole data. Resistivity values above 450 Ωm were considered as an indication of permafrost, based on the findings of ^[15]. Resistivity values below 450 Ωm were evaluated on a case-by-case basis due to the changes in resistivity as permafrost approaches 0 ℃. In the data interpretation process, anomalies with resistivity values above 450 Ωm that extended beyond 2 m below the surface were considered as evidence of permafrost. Anomalies with resistivity lower than this threshold were excluded from the permafrost layer and attributed to seasonal frost within the active layer.

The borehole data and ERT data from Profile 4 were compared to constrain the interpretation of the ERT profiles, and the values are shown in Table 3 and Figure 12. The borehole data provided detailed information on the stratigraphic sequence encountered at various depths, while the electrical resistivity data provided insights into the subsurface resistivity properties and variations. The resistivity profile obtained for Profile 4 exhibited significant differences in resistivity values, indicating contrasting subsurface conditions; however, the active layer was not well resolved across the profile from the parking lot to the thermokarst location due to the presence of base materials of the parking lot and the road in that area.

Table 3.  Summary of materials properties and resistivity.

Soil/feature	Resistivity value on Profile 4 (SC) (Ωm)	Depth of top of the feature on Profile 4 (m)	Depth of top of the soil/feature in boreholes (m)
			BH-1	BH-2	BH-3	BH-4
Saturated silt	< 100	3	–	–	2	–
Dry silt	300–400	0.5	0.5	0.45	–	0.42
Ice wedge	> 700	3.9	–	–	–	3.9

---

 | Show Table

DownLoad: CSV

Figure 12.  Comparison between boreholes and stratigraphic cross-section of Profile 4.

DownLoad: Full-Size Img PowerPoint

Borehole logs indicated relatively uniform soil types, except for the presence of ice wedge and ice-rich materials along the profile. In particular, BH-4 encountered an ice wedge at a depth of 3.9 m, which correlated well with the resistivity value greater than 700 Ωm (Figure 12 and Table 3). Hence, ice wedges in the study area are characterized by a resistivity higher than 700 Ωm. The vertical extent of the ice wedge could not be fully determined as the borehole was terminated at 5.7 m and the ERT data at the base of the anomaly is not reliable for interpretation as the DOI was greater than the 0.1 contour. BH-2 identified ice-rich soils at a depth of 3.8 m, which was not clearly delineated in the resistivity profiles. BH-3 provided confirmation of saturated material at a depth of 2 m, and the resistivity profiles showed a zone with a resistivity value which, based on literature, corresponds to a saturated material, at a depth of about 3 m. Also, the stratigraphic column of BH-1 showed the presence of unsaturated silt across the cross section from 0.45–5.5 m depth, while the resistivity profile showed the presence of saturated material at a depth of 3.6 m. Minor variations were noted in the depth correlation of some features and materials in the borehole data and the corresponding ERT anomalies, which can be attributed to the difference in timing between drilling (winter) and electrical resistivity measurements (summer), leading to slight variations in subsurface moisture conditions. Based on literature and the borehole data, the electrical resistivity measurements in the study indicated the likely presence of permafrost at the observation site, as significantly high resistivity values (>450 Ωm) were encountered at various depths along the profile.

During the interpretation of the ERT profiles in Figure 6, several significant findings were observed. One significant finding was the decrease in resistivity from 200 Ωm to approximately 100 Ωm between elevations of 170–190 m, which is depicted by the light green color. This decrease in resistivity can be attributed to the transition from slightly moist silt to saturated silt. The smooth reduction in resistivity values suggests a change in the subsurface soil moisture content. The presence of low-resistive bodies (anomalies T1–T7) across all transects, located at a horizontal distance mark of 50–90 m, suggests an area of saturated soil in the subsurface. The resistivity values in these bodies are below 50 Ωm, indicating a higher moisture content in the subsurface. Furthermore, localized high resistivity values (anomalies I1, I4, I4-1, and I5) along the northern margin of the transects (Profiles 1, 4, and 5) from 80–110 m suggests the presence of ice wedges. Ice wedges can significantly influence the resistivity measurements within their vicinity. The identification of these ice wedges is valuable as they contribute to the overall understanding of the permafrost dynamics in the study area. Patches of high resistivity values to the north of the transects within the active layer indicate the presence of frozen conditions within the seasonally frozen layer.

The inversion results (Figure 8) for the WE transect (Profiles 8 and 9) exhibited a uniform resistivity distribution below the active layer. This uniformity suggests the presence of relatively homogeneous soil moisture, which is attributed to the fact that these profiles are aligned parallel to the road and are along the topographic low (more surface water infiltration) bordering the road.

The analysis of the 3D ERT grid provided valuable insights into the resistivity distribution within the thermokarst area. Lower resistivity values indicated the presence of thawed ground, while higher resistivity values corresponded to frozen ground areas. Notably, the grid also identified distinct features, such as potential ice-rich zone or an ice-wedge, which appeared as high-resistivity anomalies in the 10 m and 11 m depth slices (Figure 10). A noteworthy observation between the pseudo-3D ERT grid and vegetation at the study site can be noted within the top 0.5 m and 2 m depth slices. In areas with dense vegetation cover (50–92 m), we observed a consistent pattern of high resistivity values in the horizontal depth slices of the pseudo-3D ERT grid. We hypothesize that this high resistivity pattern results from the protective role of vegetation cover on the active layer, reducing temperature fluctuations. Our study aligns with previous research underscoring the importance of vegetation cover in protecting the active layer and permafrost ^[1].

Conversely, areas with no vegetation cover were characterized by a slight decrease in resistivity values in the pseudo-3D ERT grid. We attribute this decrease in resistivity pattern to the lack of insulation, which leads to greater temperature fluctuations within the active layer. The discussion on electrode spacing highlights its critical role in the effectiveness of ERT surveys. The choice of a 2 m electrode spacing in this investigation was influenced by specific subsurface targets at significant depths, including ice wedges and taliks in permafrost regions. This strategic use aligns with the aim of capturing subsurface complexities at great depth, suitable for characterizing geological structures and groundwater resources. The review by ^[40] emphasizes the effective use of a 2 m electrode spacing for mapping bedrock depth and delineating subsurface lithological changes. On the other hand, high-resolution ERT, characterized by 0.5 m electrode spacing, emerges as a specialized tool. This configuration excels in near-surface investigations, showcasing its utility in mapping shallow geological features. Specifically, it is particularly effective for studies focusing on the active layer in permafrost regions, where detailed insights into near-surface conditions are crucial.

The 2D-ERT and 3D-ERT data presented in this study uncover a complex form of subsurface degradation resulting from thawing ice-rich soils and ice wedges at the vertical and lateral margins of a thermokarst-prone site. While earlier studies such as ^[25] have acknowledged the existence of a thin permafrost shelf underneath deteriorating margins in peatlands, our ERT measurements offer a more comprehensive and conclusive account of this lateral degradation in the subsurface.

A quantitative analysis comparing 2D ERT and 3D ERT grids reveals the distinctive characteristics and performance metrics. Traditionally, 2D ERT provides a cost-effective solution for subsurface imaging, excelling in scenarios where structures are predominantly vertical. However, limitations arise in capturing detailed lateral variations, as it may oversimplify certain subsurface geometries. On the other hand, the 3D ERT grid significantly improves lateral resolution, offering a more detailed representation of subsurface structures. By incorporating information from different depths, the 3D ERT grid minimizes vertical exaggeration seen in 2D profiles, providing a realistic portrayal of complex subsurface geological features. This comparative assessment highlights the nuanced strengths of each technique, emphasizing their application based on specific project goals and geological needs.

Our study indicates three potential factors contributing to the extensive permafrost degradation and development of thermokarst. These factors include the infiltration of surface water into the subsurface, topography, and the presence of near-surface ground ice. The results of the ERT survey revealed several low resistivity anomalies that were prominent on the ERT profiles. These anomalies have been labeled on Profiles 1 and 4 (see Figure 13) and were noticed extending from the surface to the bottom of the survey area. When compared to the resistivity of the surrounding soil within these zones, the absolute values are much lower, which may be attributed to the variation in water content within the soil. As water content increases, bulk resistivity decreases ^[18]. Of particular interest were the plunging nature of these anomalies, indicating active permafrost degradation due to infiltration of near-surface water. The anomaly labeled A1 on Profile 1 plunges from the surface to the bottom of the survey line; other anomalies labeled on Profile 4 (B1, B2, B3, and B4) follow a similar pattern. The width of these plunging features increases with distance from the surface, which suggests that the infiltration of near-surface water is the primary driver of permafrost degradation resulting in subaerial talik formation. Due to the significant topographic changes in this study site, the migration of snowmelt in combination with rainfall from up slope of the VL are transported downhill into the black spruce forest. Field observation shows the entrapment of water in thaw pits along the survey line (Figure 14a). In sloping areas, the development of thermokarst is often more complex and dynamic, as the gradient and orientation of the slope can influence the movement of water and soil. In the spring, as temperatures rise and the snow melts, the water infiltrates the soil and can initiate the thawing of the underlying permafrost.

Figure 13.  Possible surface water infiltration zones along Profiles 1 and 4 shown as purple lines.

DownLoad: Full-Size Img PowerPoint

Figure 14.  Surface expression facilitating thermokarst development along the study site: a) ponding of water along the profile line and b) steep depression and drunken trees along the profile line.

DownLoad: Full-Size Img PowerPoint

On sloping terrain, the melted snow can accumulate and flow downhill, concentrating in depressions or channels and further accelerating the thawing process. Snow melt can strongly influence the development of thermokarst in permafrost areas by affecting the distribution and characteristics of permafrost, the water balance, and the rate and direction of water movement ^[51].

Our site is also characterized by the formation of troughs and pits in polygonal shapes resulting from the thawing of ice-wedge polygons. The trees situated within and around the thermokarst pit are tipped and damaged in a polygonal manner (Figure 14b). When the active layer deepens, the ice wedges are exposed to warmer temperatures, causing them to thaw and ultimately degrade. As the ice wedges degrade, the surrounding soil can collapse, leading to the formation of depressions and the development of thermokarst ^[52]. The deepening of the active layer can also increase the amount of water that infiltrates the soil and permafrost, further exacerbating the degradation of ice wedges and the formation of thermokarst. The infiltration of water can cause the ice wedge to melt more rapidly, leading to increased subsidence and erosion.

By understanding these relationships, scientists and landowners can better predict the impacts of climate change and land use on the formation and expansion of thermokarst and develop strategies to mitigate their effects on ecosystems and human communities. As Alaska is projected to experience climate warming in the coming decades ^[26], the non-uniform subsurface composition of discontinuous permafrost is expected to react in unexpected ways. Understanding the morphology of these thawed-frozen margins presented in this study can aid in thermal and hydrologic modeling to predict how permafrost and subsurface groundwater flows might respond to a warming climate.

These characteristic degradation features may cause talik formation, which could create positive feedbacks and facilitate additional permafrost thaw if warm surface or shallow subsurface water interacts with permafrost ^[53]. Overall, ERT is a highly effective geophysical technique for imaging the subsurface characteristics and morphologies of thermokarst. Using ERT, researchers can gain valuable insights into the distribution of permafrost features, such as the presence of ice-rich permafrost and the structure of sediment layers within the thermokarst area. These findings can help inform our understanding of the formation and evolution of thermokarst and can be used to inform management and mitigation strategies for these unique landforms in the face of climate change. As geophysical technology continues to evolve and become more accessible, we can expect to see further advancements in our understanding of thermokarst dynamics and the implications for ecosystems and communities in permafrost regions.

6.   Conclusions

This study used a combination of inversion parameters including L₁ norm, the DOI, and two array types (DD and SC) to analyze ERT data. The use of additional methods such as the borehole data also provided valuable information for the characterization of the permafrost conditions at the study site. The study found that ERT is an effective tool in detecting the presence and distribution of near-surface anomalies, which were interpreted to be ice wedges, and taliks. The DOI analysis was useful in evaluating the reliability of the 2D-ERT resistivity models by highlighting areas that were not tied to resistivity data at depth. Borehole data were essential for interpreting and validating the geophysical results and providing insight into the subsurface conditions of the study area. The results from the study also showed that the combination of these techniques can effectively capture the spatial variability of thermokarst and provide a more detailed understanding of the mechanisms causing thermokarst formation, which is attributed to ice wedge degradation of the study site. Among the tested configurations, the SC array with the robust parameter setting, finest mesh, and no side block effect demonstrated the best fit for imaging subsurface features.

In conclusion, the distinctive resistivity anomalies revealed in the ERT profiles highlighted the presence of various permafrost-related phenomena, including ice wedges, seasonally frozen layers, and variations in moisture content. Notably, the identification of surface water infiltration zones and their role in thermokarst initiation and development provides critical insights into the complex interplay of topography, near-surface ground ice, and water infiltration in the context of permafrost degradation. The observed plunging nature of resistivity anomalies along specific profiles signifies active permafrost thawing attributed to the infiltration of near-surface water, leading to the formation of subaerial taliks. These findings contribute significantly to our understanding of the mechanisms driving thermokarst formation in the study site.

One valuable lesson learned from ERT measurements, along with borehole data, is the importance of considering the temporal aspect and potential variations between different measurement times. Minor differences observed between the individual models generated from ERT and borehole data highlight the dynamic nature of subsurface conditions and the need for cautious interpretation. By comparing measurements taken at different times, it became apparent that subsurface properties and moisture distribution can vary due to seasonal changes. These variations can introduce discrepancies between the ERT models and borehole data results. Therefore, it is crucial data acquisition is carried out within the same season in order to have a comprehensive understanding of the subsurface from the various methods incorporated in this study. This approach leads to more reliable and accurate interpretations of the subsurface characteristics from the integrated investigation method. It also highlights the importance of planning and coordinating data collection efforts to synchronize measurements and reduce the likelihood of temporal discrepancies.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgment

We sincerely thank Dr. Margaret Darrow, Luke Woosley from UAF, and Mr. Abraham Emond from the US Alaska Division of Geological & Geophysical Surveys for generously sharing their invaluable insights during this research endeavor. We also thank the UAF Permafrost Laboratory for providing equipment during this study's data collection phase.

Conflict of interest

The authors declare no conflict of interest.