**Research article**| 20 Oct 2022

# Exploring the capabilities of electrical resistivity tomography to study subsea permafrost

Mauricio Arboleda-Zapata Michael Angelopoulos Pier Paul Overduin Guido Grosse Benjamin M. Jones and Jens Tronicke

^{1},

^{2},

^{2},

^{1,2},

^{3},

^{1}

**Mauricio Arboleda-Zapata et al.**Mauricio Arboleda-Zapata Michael Angelopoulos Pier Paul Overduin Guido Grosse Benjamin M. Jones and Jens Tronicke

^{1},

^{2},

^{2},

^{1,2},

^{3},

^{1}

^{1}University of Potsdam, Institute of Geosciences, Potsdam, Germany^{2}Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research, Potsdam, Germany^{3}Institute of Northern Engineering, University of Alaska Fairbanks, Fairbanks, AK, USA

^{1}University of Potsdam, Institute of Geosciences, Potsdam, Germany^{2}Alfred Wegener Institute, Helmholtz Centre for Polar and Marine Research, Potsdam, Germany^{3}Institute of Northern Engineering, University of Alaska Fairbanks, Fairbanks, AK, USA

**Correspondence**: Mauricio Arboleda-Zapata (mauricio.arboleda@geo.uni-potsdam.de)

**Correspondence**: Mauricio Arboleda-Zapata (mauricio.arboleda@geo.uni-potsdam.de)

Received: 11 Mar 2022 – Discussion started: 20 Apr 2022 – Revised: 26 Aug 2022 – Accepted: 02 Sep 2022 – Published: 20 Oct 2022

Sea level rise and coastal erosion have inundated large areas of Arctic permafrost. Submergence by warm and saline waters increases the rate of inundated permafrost thaw compared to sub-aerial thawing on land. Studying the contact between the unfrozen and frozen sediments below the seabed, also known as the ice-bearing permafrost table (IBPT), provides valuable information to understand the evolution of sub-aquatic permafrost, which is key to improving and understanding coastal erosion prediction models and potential greenhouse gas emissions. In this study, we use data from 2D electrical resistivity tomography (ERT) collected in the nearshore coastal zone of two Arctic regions that differ in their environmental conditions (e.g., seawater depth and resistivity) to image and study the subsea permafrost. The inversion of 2D ERT data sets is commonly performed using deterministic approaches that favor smoothed solutions, which are typically interpreted using a user-specified resistivity threshold to identify the IBPT position. In contrast, to target the IBPT position directly during inversion, we use a layer-based model parameterization and a global optimization approach to invert our ERT data. This approach results in ensembles of layered 2D model solutions, which we use to identify the IBPT and estimate the resistivity of the unfrozen and frozen sediments, including estimates of uncertainties. Additionally, we globally invert 1D synthetic resistivity data and perform sensitivity analyses to study, in a simpler way, the correlations and influences of our model parameters. The set of methods provided in this study may help to further exploit ERT data collected in such permafrost environments as well as for the design of future field experiments.

In Arctic coastal regions, contemporary subsea permafrost thawing starts following the inundation caused by sea level rise and coastal erosion. Seawater is typically warmer than mean annual air temperatures, and the presence of saltwater (mostly through diffusive processes) lowers the freezing point of the seabed (Harrison and Osterkamp, 1978; Are, 2003). Additionally, groundwater flow of freshwater from inland areas might play an important role in thawing permafrost (Frederick and Buffett, 2015; Pedrazas et al., 2020), comparable to warm discharge from large rivers (Shakhova et al., 2017). Subsea permafrost is estimated to contain a large quantity of organic carbon (Sayedi et al., 2020), which can decompose microbially to generate carbon dioxide and/or methane after the permafrost thaws. Furthermore, gas hydrates are present in subsea permafrost and may act as an additional source of methane if they dissociate (Ruppel and Kessler, 2017). Understanding the development of permafrost degradation rates helps to fine-tune predictive models of greenhouse gas emissions that may represent a positive feedback for climate warming (Schuur et al., 2015). Furthermore, the correlation between permafrost degradation and coastal erosion proposed by Are (2003) and Overduin et al. (2012, 2016) might be used to refine coastal dynamics models.

Subsea permafrost is a perennially cryotic (<0 ^{∘}C) layer or body of sediments underneath a marine water column (Angelopoulos et al., 2020a). These sediments can be frozen or unfrozen. A layer of unfrozen ground in a permafrost area is known as talik, and in particular, the perennially cryotic unfrozen sediments forming part of the permafrost are known as cryopegs (Permafrost Subcommittee, 1988). Cryopegs can be isolated pockets or layers and are commonly found along Arctic coasts in saline marine sediments that are exposed following a marine regression, for example, due to isostatic uplift (O'Neill et al., 2020). Offshore, cryotic unfrozen sediment in between the water column and frozen ground is still generally referred to as a talik (Osterkamp, 2001). The subsea permafrost that contains ice is known as ice-bearing permafrost, and when the soil particles are cemented together by ice, it is termed ice-bonded permafrost (Permafrost Subcommittee, 1988). Because traditional geophysical methods such as electrical resistivity tomography (ERT) and seismic techniques can only distinguish between sediments with no or low ice content from those with high ice content (note that direct sampling is required for a quantitative estimation of ice content), we refer to them as unfrozen and frozen sediments in this study, and the interface that separates them is the ice-bearing permafrost table (IBPT). Imaging and determining the position of the IBPT (e.g., using geophysical or borehole data) is important for a better process understanding of subsea permafrost evolution and to infer degradation rates. For example, dividing the depth to the IBPT by the time since inundation results in the mean annual degradation rate (e.g., Are, 2003; Overduin et al., 2012, 2016).

Among the most used geophysical imaging techniques to study the subsea permafrost are different electromagnetic and seismic methods as well as ERT (Scott et al., 1990; Kneisel et al., 2008; Hubbard et al., 2013). Electromagnetic induction (EMI) methods are promising techniques to map both the top and bottom boundaries of the permafrost and might be used for a wide range of water depths (e.g., Sherman et al., 2017). EMI methods can properly work under low-resistivity seawater layers, while the use of ground-penetrating radar (GPR) is limited to freshwater bottom-fast ice environments characterized by high electrical resistivity values as found in delta areas (e.g., Stevens et al., 2009). Seismic methods have been employed widely in deep-water environments (e.g., Rekant et al., 2015; Brothers et al., 2016), and, more recently, researchers have used recordings of ambient seismic noise in shallow waters to map the IBPT (e.g., Overduin et al., 2015a). The ERT method is a suitable tool to investigate the resistivity distribution of the unfrozen sediments (e.g., talik and cryopeg) and for studying and delineating the IBPT position (e.g., Sellmann et al., 1989; Overduin et al., 2012, 2016; Angelopoulos et al., 2019, 2020b; Angelopoulos, 2022; Pedrazas et al., 2020).

In marine ERT surveying, floating electrodes are typically used to inject a current and measure potential differences that are used to calculate apparent electrical resistivity data. In the summer season of the Arctic, the measured values are often influenced by the resistivity and thickness of the water layer and by the unfrozen and frozen sediments. The ERT method can detect the IBPT but does not necessarily distinguish non-cryotic from cryotic taliks above the IBPT. The resistivity of seawater depends mainly on the amount of dissolved salts and temperature. The dissolved salts are mostly affected by water inflows from rivers and the cycles of sea ice melting, freezing, and brine release. The resistivity of the seawater is commonly in the range of 0.1 to 40 Ω m (e.g., Sellmann et al., 1989; Lantuit et al., 2011). On the other hand, the resistivity of the underlying sediments is influenced by porosity, pore size, grain size, water and ice content, porewater salinity, and temperature (Kneisel et al., 2008; Wu et al., 2017). For example, the resistivity of unfrozen sediments typically ranges from 1 to 25 Ω m (e.g., Sellmann et al., 1989; Overduin et al., 2012; Angelopoulos et al., 2019), while the resistivities of frozen sediments might vary from 10 Ω m up to more than 1000 Ω m (e.g., Overduin et al., 2012, 2016; Pedrazas et al., 2020; Rangel et al., 2021). The higher the ice content is, the more resistive the medium is (Pearson et al., 1986; Fortier et al., 1994; Kang and Lee, 2015). In cases where the resistivity of the frozen sediments is several orders of magnitude higher than the resistivity of the overlying unfrozen sediments, the electrical current injected through the electrodes is expected to be channeled through the less resistive layers (e.g., Spitzer, 1998), resulting in a limited penetration of the current system into the frozen sediment layers.

When analyzing ERT data collected in subsea permafrost environments, defining an appropriate inversion and model parameterization strategy is critical for deriving reliable resistivity models and interpreting these models in terms of the IBPT position. For example, when a priori information suggests that the nature of the contact between the unfrozen and frozen sediments is gradual, a grid-based model parameterization and a local inversion algorithm favoring vertical and/or horizontal smoothness in the final models might be an appropriate choice (e.g., Loke and Barker, 1996; Günther et al., 2006). Here, the experience of the interpreter might help to guess a specific resistivity threshold value to define the IBPT position (e.g., Overduin et al., 2016; Sherman et al., 2017; Angelopoulos et al., 2021). Additionally, one may also consider different gradient-based image-processing approaches to extract interfaces from the inverted resistivity model (e.g., Hsu et al., 2010; Chambers et al., 2012). Finally, when we have ample evidence of a sharp contact between the unfrozen and frozen sediments (e.g., Overduin et al., 2015b; Angelopoulos et al., 2020a), a layer-based model parameterization combined with a local and/or global inversion algorithms might be more suitable (e.g., Auken and Christiansen, 2004; Akça and Basokur, 2010; De Pasquale et al., 2019; Arboleda-Zapata et al., 2022).

In this study, we adapt the inversion and ensemble interpretation strategies as proposed by Arboleda-Zapata et al. (2022) to study submarine permafrost environments of the Arctic in terms of the resistivity distribution of the unfrozen and frozen sediments and the position of the IBPT, including estimates of uncertainties. We analyze and compare ERT data collected at two field sites in the Arctic characterized by different environmental conditions regarding seawater depth and resistivity, coastal erosion rate, and the sediment porewater chemistry. Additionally, we generate and interpret ensembles of globally inverted 1D electrical data to get a deeper understanding of the inverse problem for typical resistivity distributions in these kinds of environments. Finally, we also implement local and global sensitivity analyses to recognize the most influential parameters during 2D and 1D inversions.

Our field data were collected at two field sites: one offshore of the southern part of the Bykovsky Peninsula in the Siberian Laptev Sea (Fig. 1a) and the other one offshore of Drew Point in the Alaskan Beaufort Sea (Fig. 1d). To relate these data sets to the site-specific environmental settings, we summarize the main characteristics of each study area in a regional framework.

## 2.1 Regional setting of Bykovsky Peninsula

The Bykovsky Peninsula is located in northeastern Siberia in the vicinity of the Lena River Delta. The peninsula is mainly characterized by the presence of ice-rich sediments (volumetric ice content exceeding 80 %, also known as the Yedoma Ice Complex) that accumulated during the Late Pleistocene (Schirrmeister et al., 2002; Grosse et al., 2007). The Yedoma deposits extend to 15 m below sea level (Grigoriev, 2008). The sediments at or below sea level are composed of silt, sand, and gravel with variable grain size distributions (Grosse et al., 2007). During the Early to Middle Holocene, a general landscape transformation started, resulting in a thermokarst-dominated relief characterized by thermo-erosional valleys and thermokarst lakes (Schirrmeister et al., 2002; Grosse et al., 2007). The mean coastal erosion rates at different locations of the peninsula typically range between 0.4 and 1.5 m yr^{−1}, with maximum values of up to 10 m yr^{−1} mainly caused by storms and thermomechanical erosion of ice-rich sediments (Lantuit et al., 2011). The seawater around the peninsula is strongly influenced by freshwater and sediments originating from the Lena River (Lantuit et al., 2011). Additionally, the resistivity of the seawater is influenced by seasonal sea ice freezing and melting as shown by Lantuit et al. (2011), who report resistivity values for the eastern shore of the Bykovsky Peninsula of less than 1 Ω m in winter and above 10 Ω m in summer. Similar water resistivity values are also reported by Overduin et al. (2016) for the seawater near Muostakh Island. The depth of the seawater for the southern part of the Bykovsky Peninsula deepens 2 m within a distance of 100 m from the shoreline and increases to 5 m at about 2000 m from the coast (Lantuit et al., 2011; Fuchs et al., 2021).

## 2.2 Regional setting of Drew Point

Drew Point is located on the coast of the Alaskan Beaufort Sea. The local geology is characterized by glaciomarine, fine-grained, ice-rich sediments deposited in the Late Pleistocene (Black, 1964; Ping et al., 2011). The inland geomorphology is characterized by 3–5 m high coastal bluffs, thermokarst channels and lakes, and ice-wedge polygons on tundra plains with maximum elevations of ∼10 m (Barnhart et al., 2014; Jones et al., 2018). The average coastal erosion rate between 1979 and 2002 was around 9 m yr^{−1} (Jones et al., 2009) and increased for the period 2002 to 2016 up to approximately 20 m yr^{−1} (Jones et al., 2018). Lück (2020) reports brackish water resistivities observed during fieldwork in July 2018 of 0.4–0.5 Ω m, with weak stratification visible in the water column profiles. The depth of the seawater offshore of Drew Point deepens to 2 m within a distance of 500 m from the shoreline and increases to 3 m at distances of about 2000 m from the coast (Jones et al., 2018).

In marine ERT data acquisition, there is typically an excellent coupling between the floating electrodes and the seawater. This allows us to perform voltage measurements while the boat pulling the electrode streamer is in motion (preferably at constant speed) and, thus, to also efficiently measure profiles with a length on the order of kilometers. The sources of errors during data acquisition are mainly related to misalignments of the electrode streamer (e.g., due to water currents), the precision of electrode positioning (which is given relative to boat position), vertical oscillation of electrodes (e.g., due to wavy conditions), and surface area limitation of injection voltage. Furthermore, due to the large variety of environmental settings, one must tailor the survey parameters to each field site, which includes varying the electrode spacing, the transmitter voltage, the measurement duration, the boat speed, the digital resolution of the potential measurements, and the sampling frequency.

In Table 1, we compare the acquisition parameters for our ERT data from Bykovsky and Drew Point, which were collected during two fieldwork campaigns in July 2017 and July 2018, respectively. The two ERT data sets were collected using an Iris™ Syscal Pro Deep Marine system employing a streamer cable with 13 equally spaced floating electrodes. The resistivity measurements were acquired using a reciprocal Wenner–Schlumberger array configuration, where current was injected through the inner pair of electrodes and quasi-symmetric voltages were measured simultaneously with 10 channels using the outer pairs of electrodes (e.g., Overduin et al., 2012). The transmitter voltage was set at approximately 48 V at Bykovsky, while at Drew Point it was reduced to 24 V to avoid exceeding the electrode surface area limits. Additionally, different electrode spacings were used. The Bykovsky data were recorded using a 10 m spacing between electrodes, while at Drew Point 5 m spacing was chosen because the rapid coastal erosion rates suggested that the IBPT position at this field site should be shallower than at our Siberian field site for a given distance offshore. To collect the data along every profile, a cable was towed behind a small inflatable boat and voltages were measured as the boat moved at approximately constant speed of ∼1 m s^{−1} perpendicular to the shore. The Bykovsky soundings were collected at spacings of ∼7.7 m along a 418 m long profile, resulting in 540 measurements. At Drew Point, the soundings were collected at spacings of ∼4.7 m along a 854 m long profile, resulting in 1830 measurements. At both field sites, the electrode positions were estimated relative to the position of a GPS aboard the boat assuming a straight streamer cable. Complementary to the ERT measurements, we also recorded the water depth at each sounding location using a Garmin echo sounder attached to the boat (Fig. 1b, and e). Furthermore, we measured during each field campaign water resistivity and temperature at different depths close to our ERT profiles (the mean values are shown in the last two rows of Table 1) using a SonTek™ CastAway device also known as CTD.
In general, at our Drew Point field site, the seawater was shallower, less resistive, and slightly cooler than at our Bykovsky field site. Furthermore, the CTD measurements suggest low vertical and horizontal variations in the resistivity of the water layer at both field sites. For example, the largest variations are in the horizontal direction and are on the order of 1 and 0.04 Ω m for our Bykovsky and Drew Point data sets, respectively.

The measured apparent resistivities *ρ*_{a} are presented as pseudosections in Fig. 1c and f. Here the *x* coordinates represent the center position of each quadripole, and the vertical axes represent the relative penetration also known as levels; i.e., level =1 is the shortest quadripole, while level =10 is the quadripole with maximum electrode spacing. The range of *ρ*_{a} for Bykovsky is 5.9 to 43 Ω m and for Drew Point 0.9 to 7.6 Ω m. The lower *ρ*_{a} values at Drew Point are mainly due to the lower resistivity of the seawater at the Alaskan coast, which is less influenced by freshwater discharge from large rivers than at our Bykovsky field site. Additionally, we notice that levels larger than seven in our Bykovsky data set are characterized by higher variations due to noise or 3D subsurface structures. In contrast, the Drew Point data do not show obvious variations depending on the level number.

In this study, we follow the workflow of Arboleda-Zapata et al. (2022), who propose a layer-based model parameterization to globally invert 2D ERT data, which is used to generate an ensemble of representative model solutions. For completeness, we present a brief summary of this workflow in the following. For a more detailed analysis, we will also address complementary strategies such as 1D inversion tests as well as local and global sensitivity analyses.

## 4.1 Two-dimensional layer-based model parameterization

One of the most important steps in any geophysical inversion workflow is defining a model parameterization that can properly represent the studied geological environment. Because a priori information suggests a layered subsurface (i.e., unfrozen sediments overlying frozen sediments) at both of our field sites, we choose a layer-based model parameterization considering an unstructured mesh with local refinements along the interfaces separating individual layers. Additionally, because resistivity variations within each layer are negligible compared to the variations between different layers, we assume homogeneous layers; i.e., each layer is characterized by one resistivity value. For more complex geological settings, one might allow for lateral and/or vertical variations within the layers (e.g., Auken and Christiansen, 2004; Akça and Basokur, 2010). To parameterize the interface geometry that defines the contact between the individual layers, we may use different strategies, for example, based on spline interpolation (e.g., Koren et al., 1991), Fourier coefficients (e.g., Roy et al., 2021), or sums of arctangent functions (Gebrande, 1976).

Allowing for abrupt changes along the interfaces is considered to be a critical point in subsea permafrost environments where high structural variability is often found. In such environments, we expect sharp boundaries and variations along the interfaces due to inundated thermokarst structures (Angelopoulos et al., 2021), pingo-like features, bottom-fast ice versus floating ice regime transitions in winter, or changes in the ratio of coastal erosion vs. degradation rate; i.e., changing from a period of fast thawing and low coastal erosion to a period of fast coastal erosion and slow thawing can result in a heterogeneous structure of the IBPT (e.g., Overduin et al., 2016). Because we expect some of these processes and structures at our field sites, we adopt a strategy based on the sums of arctangent functions because it allows for abrupt and smooth changes along the interfaces (e.g., Roy et al., 2005; Rumpf and Tronicke, 2015). Following Arboleda-Zapata et al. (2022), the sums of arctangent functions for a single interface can be written as

where *z* is the depth, *n*_{nod} is the number of arctangent nodes, *z*_{0} is the average depth of the interface, *x*_{j} is the horizontal location of an arctangent node, and Δ*z*_{j} is the vertical throw attained asymptotically over a horizontal distance of *b*_{j}. Such sets of coefficients are used to obtain *z*(*x*) at horizontal distances *x*. Increasing the number of nodes allows us to resolve more complex interfaces. During preliminary experiments and parameter testing, we noticed that using three to seven nodes allows us to model rather complex interfaces. For both of our field studies, we fix the number of nodes to five, which results in 16 model parameters (1+3*n*_{nod}) per interface. Because we consider two interfaces, one for the seabed and the other for the IBPT separating three layers with homogeneous resistivities, our parameterization strategy results in a total of 35 model parameters.

## 4.2 Inversion strategy

During inversion, we search for a combination of model parameters (i.e., those describing the geometry of interfaces using Eq. 1 and the resistivities of the homogeneous layers) that minimizes the root mean squared logarithmic error (RMSLE). To reduce the space of possible solutions, we consider some constraints in our layer-based parameterization approach. For both case studies, we constrain the seabed position (Fig. 1b and e) allowing vertical variations of up to ±0.15 m, which is the approximate error level of our echo sounder data for water depths <5 m. Considering our CTD measurements, we allow the resistivity of the water to vary between 11 and 15 Ω m for our Bykovsky data and between 0.2 and 2 Ω m for our Drew Point data. Note that we consider additional freedom beyond the variabilities reported in Sect. 3 (1 Ω m for Bykovsky and 0.04 Ω m for Drew Point) to account for additional variations related to the different sensitivities and resolutions of our CTD and ERT data. Additionally, we set our search parameter range for the resistivity of the talik from 1 to 100 Ω m and that for the resistivity of the ice-bearing permafrost from 1 to 300 000 Ω m for both field studies.

Because we aim to find an inverse model independent of a reference or starting model, we use a global inversion strategy based on the particle swarm optimization (PSO) technique, which was originally introduced by Kennedy and Eberhart (1995). Over the last decade, the PSO algorithm has been widely used to invert different types of geophysical data sets because it has proven to be an effective tool for finding different local minima in objective functions with complicated topography (e.g., Tronicke et al., 2012; Fernández-Martínez et al., 2017).

In a first step, the PSO requires defining a set of particles where each particle represents a different model. The particles are initialized with random parameters bounded within realistic physical ranges. This defines our model space. The position of each particle is updated iteratively considering the best global position found so far by the entire swarm (i.e., the particle with the best fit performance in terms of the RMSLE), the best local position (i.e., the best fit performance in the history of each particle), and the inertia (i.e., the direction in which the particle is moving). These parameters are weighted and perturbed with random numbers drawn from a uniform distribution, which helps avoid getting trapped in a local minimum. For every particle and every iteration, we calculate the forward response using the Python library pyGIMLi (Rücker et al., 2017). Note that each particle contains model parameters that result in two different interface geometries: one representing the seabed and the other the IBPT. Thus, adding these interfaces to our finite-element mesh results in a different mesh geometry for each particle. To ensure good mesh quality, we constrain the minimum angle within each cell to 33.5^{∘} (Shewchuk, 1996). This parameterization strategy allows us to calculate the forward response with high precision and with a reasonable amount of time (Arboleda-Zapata et al., 2022). At the end, when the optimization reaches the maximum number of iterations or a minimum threshold in the objective function, we save the final best model. Using different seeds of the random number generator, we repeat this process until we obtain an ensemble *M*_{F0} consisting of several hundred independent models and an ensemble of corresponding residuals *δ*_{F0}. In this study, each residual vector is calculated as the difference between the observed and the corresponding modeled log-apparent resistivity values.

## 4.3 Ensemble interpretation

In a first step, to ease our ensemble analysis and interpretation in a pixel-wise fashion, all models in *M*_{F0} are interpolated using the nearest-neighbor algorithm on a densely discretized structured mesh (note that we use a unstructured mesh during inversion, Sect. 4.1). In a second step, we perform a cluster analysis using the *k*-means algorithm (MacQueen, 1967) and considering *M*_{F0} and *δ*_{F0} as input to group similar solutions from our ensembles. To find an optimal number of clusters *n*_{k}, we use the criterion proposed by Caliński and Harabasz (1974) supported by a visual inspection of the clustering results. Finally, we characterize in a pixel-wise fashion each found cluster *M*_{Fi} and *δ*_{Fi} (where $i=\mathrm{0},\mathrm{1}$, …, *n*_{k}, note *i*=0 correspond to the whole ensemble and *i*>0 to the clustered ensembles) by the median values ${\mathit{\mu}}_{\mathrm{1}/\mathrm{2}}\left({\mathit{M}}_{\mathrm{Fi}}\right)$ and ${\mathit{\mu}}_{\mathrm{1}/\mathrm{2}}\left({\mathit{\delta}}_{\mathrm{Fi}}\right)$ and the interquartile ranges IQR(*M*_{Fi}) and IQR(*δ*_{Fi}). Additionally, we describe *δ*_{Fi} in an overall fashion assessing the RMSLE(*δ*_{Fi}), the IQR(*δ*_{Fi}), and the quantile 90 % *q*_{90}(*δ*_{Fi}).

## 4.4 One-dimensional inversion

Often, we prefer 2D inversion algorithms in comparison to 1D strategies – especially for field data where the subsurface situation and its complexity are largely unknown. However, to investigate and understand, for example, the relationship between specific model parameters and the influence of a priori information and constraints, 1D approaches (also considering synthetic data examples) represent helpful interpretation tools (e.g., Sen and Stoffa, 1996; Malinverno, 2002).

In this study, we use 1D models consisting of five parameters, the depth of the seawater *z*_{w}, the depth of the contact between unfrozen and frozen sediments *z*_{pt} (i.e., IBPT), the water resistivity *ρ*_{w}, the resistivity of the unfrozen sediments *ρ*_{uf}, and the resistivity of the ice-bearing permafrost *ρ*_{p}. As for our 2D examples, we also consider PSO to invert our 1D synthetic data. Because for such 1D inversions the computational cost is significantly lower than for 2D problems, we can run several tests and create larger model ensembles. We use such a 1D approach to tackle some specific questions regarding the considered application. For example, we investigate how constraining the depth of the water layer and its resistivity affects the final ensemble of 1D model solutions. Additionally, the limited number of parameters in our 1D model parameterization strategy allows us to study in a simpler way the posterior correlation matrix as proposed by Sen and Stoffa (2013). Although in this study we do not consider cluster analysis to classify our 1D ensembles as implemented for our 2D analyses, this step may be adapted in future studies (e.g., investigating more complex model scenarios).

## 4.5 Sensitivity analysis

Sensitivity analysis is a powerful tool that can provide additional information to improve system or process understanding (Wainwright et al., 2014). In the context of subsea permafrost applications, several studies have shown the potential of the ERT method to image the IBPT position (e.g., Sellmann et al., 1989; Overduin et al., 2012). However, the sensitivity distribution of the ERT model parameters for such environments characterized by resistivity contrasts up to several orders of magnitude between unfrozen and frozen sediments is poorly understood. Adding sensitivity analysis to the interpretation workflow helps investigate the impact of our chosen model parameterization and the used constraints. Furthermore, such sensitivity studies might also help optimize ERT acquisition geometries and strategies before a field campaign.

In this study, we use 2D-local and 1D-global sensitivity analyses. To investigate which regions of the 2D discretized model have the greatest influence on our objective function, we consider the difference-based local sensitivity method of Günther et al. (2006), which is available within the pyGIMLi library (Rücker et al., 2017). For example, we assess the sensitivity of the shortest electrode configurations to understand if the corresponding measurements are influenced by both the water layer and the underlying unfrozen sediments. In turn, this helps to evaluate the reliability of imaging the uppermost water layer (e.g., for measurements where no CTD measurements are available). Furthermore, the longest electrode spreads (corresponding to the deepest levels in 2D pseudosections) and/or cumulative sensitivity distributions provide information on whether our ERT data are sensitive to the IBPT and/or the frozen sediments. For 1D model parameterizations and synthetic studies (considering *z*_{w}, *z*_{pt}, *ρ*_{w}, *ρ*_{uf}, and *ρ*_{p} as described in Sect. 4.4), we use the variance-based global sensitivity method of Sobol (Sobol, 2001; Saltelli et al., 2008) as implemented in the Python library SALib (Herman and Usher, 2017). Using this approach, we aim to understand how the total influence of the considered parameters might be affected by variations in *ρ*_{p} and *z*_{pt}.

In the following, we present the 2D inversion results for the Bykovsky and Drew Point data sets in two separate subsections. Each subsection is complemented with 1D inversion results of synthetic data simulated considering the site-specific environmental and electrode settings as well as with a 2D-local and a 1D-global sensitivity analysis.

## 5.1 Bykovsky

The geological and environmental settings of the Bykovsky area are described in Sect. 2.1, and a summary of the acquisition parameters and measured seawater properties is provided in Table 1. We invert the 540 apparent resistivity measurements recorded along a 418 m long profile (Fig. 1c) using a layer-based model parameterization as described in Sect. 4.1 and a PSO-based inversion strategy as outlined in Sect. 4.2. In the PSO, we use 60 particles and a maximum of 600 iterations as stopping criterion. To obtain a single inverted model (i.e., after one inversion run), we have to evaluate the forward response 36 000 times, which takes on average ∼40 h on a single core of a modern desktop computer. We repeat these inversion runs considering different initial seeds of the random number generator (note that this approach allows for a straightforward parallelization when multiple cores are available) until we obtain an ensemble *M*_{F0} consisting of 690 models.

### 5.1.1 Ensemble analysis

After the inversion, we interpolated all models to a refined structured mesh before performing any posterior statistical analyses (see Sect. 4.3). In Fig. 2a and b, we show the ${\mathit{\mu}}_{\mathrm{1}/\mathrm{2}}\left({\mathit{M}}_{\mathrm{F}\mathrm{0}}\right)$ and IQR(*M*_{F0}) models calculated from the Bykovsky model ensemble. The ${\mathit{\mu}}_{\mathrm{1}/\mathrm{2}}\left({\mathit{M}}_{\mathrm{F}\mathrm{0}}\right)$ model indicates that *ρ*_{uf} is ∼4 Ω m and *ρ*_{p} is ∼60 000 Ω m. However, when analyzing individual models, we note a bimodal distribution of *ρ*_{p} (some models with *ρ*_{p}<2000 Ω m and others with *ρ*_{p}>100 000 Ω m), which is also illustrated by increased IQR(*M*_{F0}) values for the lowermost layer. These observations already indicate different groups of models with distinct resistivity characteristics.

In the next step, we performed cluster analysis (Sect. 4.3) and found that our ensemble *M*_{F0} can be divided into two model families *M*_{F1} and *M*_{F2}. In Fig. 2b–c and e–f, we present the ${\mathit{\mu}}_{\mathrm{1}/\mathrm{2}}\left({\mathit{M}}_{\mathrm{Fi}}\right)$ and IQR(*M*_{Fi}) models (where $i=\mathrm{1},\mathrm{2}$). Comparing these models illustrates that *M*_{F1} and *M*_{F2} show a similar IBPT shape dipping toward the open sea (i.e., depth of the IBPT increases with increasing profile distances). However, the IBPT position in *M*_{F1} is shallower than in *M*_{F2}. We learn from this that for models favoring larger *ρ*_{p} values, the depth of the IBPT increases, highlighting the trade-off between these two parameters that cause model variations along the IBPT. According to the depth of the IBPT and its gradients in the profile direction, we laterally subdivide the models into three main parts. The first part is found at *x*<130 m and is characterized by a gentle dipping slope with minor convexities and concavities. The second part is found at $\mathrm{130}<x<\mathrm{280}$ m, where the IBPT is relatively flat with a minor change in depth at *x*∼200 m. Finally, the abrupt change at *x*=280 m marks the transition to the third part, where the IBPT reaches its deepest point and extends until the end of the profile at depths >20 m.

We assess the fit performance in a pixel-wise and in an overall fashion for the residuals associated with the ensemble containing all models *δ*_{F0}, as well as for the two clustered model families *δ*_{F1} and *δ*_{F2} (Fig. 3). Thus, we calculate ${\mathit{\mu}}_{\mathrm{1}/\mathrm{2}}\left({\mathit{\delta}}_{\mathrm{Fi}}\right)$ and IQR(*δ*_{Fi}) (where $i=\mathrm{0},\mathrm{1},\mathrm{2}$) in a pixel-wise fashion and present them as pseudosections in Fig. 3a–f. When comparing these pseudosections to each other, we notice that the ${\mathit{\mu}}_{\mathrm{1}/\mathrm{2}}\left({\mathit{\delta}}_{\mathrm{Fi}}\right)$ and IQR(*δ*_{Fi}) indicate similar fits of the data in terms of amplitudes and pseudosection patterns. The abrupt change from positive to negative residuals at *x*≃200 m coincides with the highest point of the bathymetric profile for *x*>150 m, which also corresponds to a general change in the gradients of the bathymetric profile (Fig. 1b). Therefore, a 3D subsurface structure (which cannot be explained by our 2D inversion strategy) and related 3D effects are a reasonable explanation of the discussed features in the residuals. For example, because the landscape was partly covered by lakes (which acted as a source of heat) prior to seawater submergence, lateral temperature gradients and heterogeneous sediment properties could affect subsurface resistivity and its 3D variations. The overall statistics RMSLE(*δ*_{Fi}), IQR(*δ*_{Fi}), and *q*_{90}(*δ*_{Fi}) (where $i=\mathrm{0},\mathrm{1},\mathrm{2}$) are presented as histograms in Fig. 3g–i. The histograms are characterized by bimodal distributions, especially evident in all shown RMSLE(*δ*_{Fi}) histograms. When comparing the histograms of *δ*_{F1} and *δ*_{F2}, we notice that they follow similar distributions (although in *δ*_{F1} there are less models). From these analyses of the residuals, we are not able to prefer one of the model families, and, thus, we perform some synthetic exercises to deepen our understanding of this inverse problem and the found model solutions.

### 5.1.2 One-dimensional inversion of synthetic data

To complement our understanding of the formulated inverse problem, we perform 1D inversions of a synthetic data set created considering a 1D subsurface model (see “Input model” in Table 2) as described in Sect. 4.4. The 1D subsurface model parameters were chosen by analyzing our 2D model solutions (e.g., Fig. 2b–c at *x*≃150 m). We calculate the forward response of 10 quadripoles considering the same electrode configurations as used for recording the Bykovsky field data (Table 1). We invert the simulated apparent resistivity data using two scenarios for constraining *z*_{w} and *ρ*_{w}, while the constraints for all other parameters remain unchanged (see Table 2). The resulting inverted models are shown in Fig. 4a and c. For all models, we have achieved RMSLE <0.028, which is equivalent to the noise level applied to the calculated synthetic data and comparable to the RMSLE achieved for the 2D inversion results of the Bykovsky field data. Comparing the results shown in Fig. 4a and c illustrates that constraining the water layer significantly decreases the non-uniqueness of the inverse problem. We also notice that the median model represents a good approximation to the input model except for *ρ*_{p}, which is overestimated as illustrated by the larger *ρ*_{p} of the quantile 25 % model compared to the input *ρ*_{p}. Additionally, from all models visualized in Fig. 4a and c, we calculate the corresponding posterior correlation matrices (Fig. 4b and d). For both cases, we see that [*ρ*_{uf}, *z*_{pt}] is the model parameter pair with the highest positive correlation, while the rest of the model parameter pairs are characterized by negative correlations with different amplitudes (except for pairs with *ρ*_{p} which show correlations approaching zero).

### 5.1.3 Sensitivity analysis

To understand the sensitivity distribution for our three-layer model (representing seawater and unfrozen sediments overlying frozen sediments), we calculate the cumulative sensitivity and the sensitivity for the shortest and widest quadripoles considering two model scenarios (Fig. 5). In the first scenario, we consider the same input model as for the 1D inversion exercise (Table 2). In the second scenario, we set *z*_{pt}=25 m while all other parameters remain unchanged. From the cumulative sensitivity plots (Fig. 5a and d), we learn that areas of sensitivities extend throughout the layer of unfrozen sediments for both scenarios. This suggests that we could interpret our inverted models even underneath the outer electrode positions; i.e., if the boat together with the electrode streamer is moving toward the right (i.e., increased *x* coordinates) to collect additional sounding curves, our interpretation of the inverted model should start at $x\sim -\mathrm{60}$ m. Note that a more conservative model interpretation might start at $x\sim -\mathrm{25}$ m, where we start having more significant cumulative sensitivities. When analyzing Fig. 5b and e, we see that the shortest quadripole is sensitive to both the water layer and the unfrozen sediments. For a wider electrode spacing and an IBPT located at a depth of 15 m (Fig. 5c), the sensitivities are focused around the inner electrodes but also with some minor contributions from the outer electrodes (note the reddish colors in the unfrozen sediments at $x<-\mathrm{60}$ m and at *x*>60 m), which may be critical when significant 2D or 3D resistivity variations are present. For a deeper IBPT (Fig. 5f), we notice that we are still sensitive at depths of ∼25 m; however, the lateral extensions of the sensitivity patterns within the unfrozen sediments appear to be reduced.

As noticed in our 2D sensitivity analysis, the high resistivity contrast between the unfrozen and frozen sediments seems to limit the penetration depth down to the IBPT. To complement and better understand our results of 2D sensitivity analysis, we investigate the global sensitivities (Sect. 4.5) of different 1D model parameterizations. Specifically, we use models where *z*_{w}=4.5 m, *ρ*_{w}=13.7 Ω m, and *ρ*_{uf}=4 Ω m are fixed, while *ρ*_{p} varies between 10 and 10 000 Ω m (eight values in total) and the IBPT is located at three different depths; i.e., *z*_{pt}=25, *z*_{pt}=15, and *z*_{pt}=5 m (Fig. 6a–c). Note that defining eight different values for *ρ*_{p} and three for *z*_{pt} results in 24 different 1D models. For the calculation of the total sensitivity for each of our five parameters in these 24 models, we set the parameters ranges to ${z}_{\mathrm{w}}=[\mathrm{4},\mathrm{6}]$ m, ${z}_{\mathrm{pt}}=[\mathrm{6.5},\mathrm{30}]$ m, ${\mathit{\rho}}_{\mathrm{w}}=[\mathrm{0.2},\mathrm{20}]\phantom{\rule{0.125em}{0ex}}\mathrm{\Omega}$ m, ${\mathit{\rho}}_{\mathrm{uf}}=[\mathrm{1},\mathrm{20}]\phantom{\rule{0.125em}{0ex}}\mathrm{\Omega}$ m, and ${\mathit{\rho}}_{\mathrm{p}}=[\mathrm{5},\mathrm{20}\phantom{\rule{0.125em}{0ex}}\mathrm{000}]\phantom{\rule{0.125em}{0ex}}\mathrm{\Omega}$ m. For these specific models and parameter ranges, our results (Fig. 6) suggest *ρ*_{w} is the most influential parameter followed by *ρ*_{uf}, which shows approximately half of the influence compared to *ρ*_{w}. The influence of *z*_{pt} is slightly larger than *z*_{w}, although *z*_{pt} is set up with a wider range than *z*_{w}. Furthermore, although we allow *ρ*_{p} to vary over 3 orders of magnitude, the result of this sensitivity analysis demonstrates that *ρ*_{p} is the parameter with the lowest influence, but it is not null as indicated by the results of our 2D sensitivity analyses (Fig. 5a and d). Such low sensitivity values help to explain the large variation of *ρ*_{p} in our 1D and 2D ensembles. Interestingly, we also notice in Fig. 6a–c that in general when increasing *ρ*_{p} (for *ρ*_{p}<100 Ω m) the total sensitivity index of the other parameters tends to decrease.

## 5.2 Drew Point

The geological and environmental settings of the Drew Point area are described in Sect. 2.2, and a summary of the acquisition parameters and measured seawater properties is provided in Table 1. We invert the 1830 apparent resistivity measurements recorded along an 854 m long profile (Fig. 1f) considering a layer-based model parameterization as described in Sect. 4.1 and a PSO-based inversion strategy as outlined in Sect. 4.2. In the PSO, because we notice that the inversion of the Drew Point data set converges much faster than in our Bykovsky example, we decide to lower the number of particles to 30 and the number of iterations to 400, thus allowing us to save some computational cost. Considering these settings, to obtain a single inverted model, we have to evaluate the forward response 12 000 times, which takes on average 57 h on a single core of a modern desktop computer. We repeat these inversion runs considering different initial seeds of the random number generator (using different processors in parallel) until we obtain an ensemble *M*_{F0} consisting of 416 models.

### 5.2.1 Ensemble analysis

After the inversion, we interpolate all models to a refined structured mesh before performing any posterior statistical analyses (Sect. 4.3). In Fig. 7a and b, we present the ${\mathit{\mu}}_{\mathrm{1}/\mathrm{2}}\left({\mathit{M}}_{\mathrm{F}\mathrm{0}}\right)$ and IQR(*M*_{F0}) models calculated from the Drew Point model ensemble. The irregular variations in the IQR(*M*_{F0}) model and the bimodal distribution of *ρ*_{p} (some models with *ρ*_{p}<500 Ω m and others with *ρ*_{p}>100 000 Ω m) already indicate different groups of models with distinct resistivity characteristics and IBPT positions.

In the next step, we performed cluster analysis (Sect. 4.3) and found that our ensemble *M*_{F0} can be divided into three model families (*M*_{F1}, *M*_{F2}, and *M*_{F3}). In Fig. 7b–d and f–h, we present the ${\mathit{\mu}}_{\mathrm{1}/\mathrm{2}}\left({\mathit{M}}_{\mathrm{Fi}}\right)$ and IQR(*M*_{Fi}) models (where $i=\mathrm{1},\mathrm{2},\mathrm{3}$). Comparing these models illustrates that *M*_{F1} and *M*_{F2} present a similar IBPT shape dipping toward the open sea. However, for *M*_{F3} the IBPT position is dipping toward the coast, which is not in agreement with our background knowledge of this field site. When comparing the *M*_{F1} and *M*_{F2} models in more detail, we note that the IBPT position in *M*_{F1} is shallower than in *M*_{F2}. Comparable to the Bykovsky example, models favoring high *ρ*_{p} values tend to show increased depths of the IBPT, resulting in thicker unfrozen sediments also near the coast. According to the depth of the IBPT and its gradients in the profile direction for *M*_{F1} and *M*_{F2}, we laterally subdivide the model into four main parts. The first part is found at *x*<100 m, and it is characterized by an intermediate convex slope. The second part is found at $\mathrm{100}<x<\mathrm{500}$ m, and the IBPT shows a gentle convex slope, whereas in the third part (at $\mathrm{500}<x<\mathrm{700}$ m) the IBPT is almost flat. Finally, the fourth part is found at *x*>750 m, where the IBPT may be located at depths ≥20 m.

We assess the fit performance for the residuals associated with the ensemble containing all models *δ*_{F0}, as well as for the three clustered model families *δ*_{F1}, *δ*_{F2}, and *δ*_{F3} (Fig. 8). We calculate ${\mathit{\mu}}_{\mathrm{1}/\mathrm{2}}\left({\mathit{\delta}}_{\mathrm{Fi}}\right)$ and IQR(*δ*_{Fi}) (where $i=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3}$) in a pixel-wise fashion and present them as pseudosections in Fig. 8a–h. When comparing these pseudosections to each other, we notice that ${\mathit{\mu}}_{\mathrm{1}/\mathrm{2}}\left({\mathit{\delta}}_{\mathrm{Fi}}\right)$ indicate similar fits of the data in terms of amplitudes and pseudosection patterns (although with slightly higher values for *δ*_{F3}). When comparing the IQR(*δ*_{Fi}) plots, we note that IQR(*δ*_{F0}) is characterized by several patches which are less prominent in the clustered residuals Fig. 8f–h. This indicates that our clustering results are properly grouping models with similar residuals. Furthermore, we associate the vertical feature at *x*=400 m in Fig. 8e–h to the variation in our models to locate the left edge of a bulge structure of the seabed (see Fig. 1e). This illustrates the applicability of exploring such residual statistics to identify possible drawbacks in our inversion results and, thus, allow us to re-evaluate our parameterization strategy. For example, we might consider to improve the inversion results by adding a node to our sums of arctangent functions around *x*=400 m. The overall statistics RMSLE(*δ*_{Fi}), IQR(*δ*_{Fi}), and *q*_{90}(*δ*_{Fi}) (where $i=\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{3}$) are presented as histograms in Fig. 8i–l. The histograms in Fig. 8i are characterized by bimodal distributions. Such bimodal distributions are less pronounced for the clustered families (Fig. 8j–l); however, small tails to the right are also evident for *δ*_{F1} and *δ*_{F2}. One may tend to reject the models falling in these tails, especially, when using the mean to estimate the central trend. However, because we consider robust statistical measures (e.g., median and IQR), we do not expect a significant impact from these models on our results and conclusions.

### 5.2.2 One-dimensional inversion of synthetic data

Following Sects. 4.4 and 5.1.2, we perform 1D inversions of a synthetic data set created considering a 1D subsurface model (see “Input model” in Table 3). The 1D model parameters were chosen by analyzing our 2D model solutions (e.g., Fig. 7b–c at *x*≈600 m). Note *ρ*_{p} is the same as in the 1D synthetic example from Sect. 5.1.2, which allows us to better compare the results of our 1D synthetic exercises. We calculate the forward response of 10 quadripoles considering the same electrode configurations as used for recording the Drew Point field data (Table 1). We invert the simulated apparent resistivity data considering two scenarios for constraining *z*_{w}, *ρ*_{w}, and *ρ*_{uf}, while the constraints for *z*_{pt} and *ρ*_{p} remain unchanged (see Table 3). The resulting inverted models are shown in Fig. 9a and c. For all models, we have achieved RMSLE <0.007, which is equivalent to the noise level applied to the calculated synthetic data and comparable to the RMSLE achieved for the 2D inversion results of the Drew Point field data. Comparing the results shown in Fig. 9a and c illustrates that the applied constraints improve the median model. However, we also observe an increase in the variability of the models around *z*_{pt} in Fig. 9c. Additionally, from all the models visualized in Fig. 9a and c, we calculate the corresponding posterior correlation matrix (Fig. 9b and d). For both cases, we see that the largest negative correlations are found for the model parameter pairs [*ρ*_{uf}, *ρ*_{w}] and [*ρ*_{w}, *z*_{pt}], while the most significant positive correlation is found for [*ρ*_{uf}, *z*_{pt}]. Note that the absolute correlations of these model parameter pairs are larger in Fig. 9d compared to Fig. 9c. Finally, we want to point out that the signs for the most significant parameter correlations are the same as the ones found for Bykovsky in Fig. 4d.

### 5.2.3 Sensitivity analysis

For the sensitivity analysis, we consider the two model scenarios indicated in Fig. 10. In the first scenario, we consider the same input model as for the 1D inversion exercise (Table 3). In the second scenario, we set *z*_{pt}=16 m while all other parameters remain unchanged. From analyzing the cumulative sensitivity plots (Fig. 10a and d), we infer that an interpretation of our inversion results should focus on the area around the inner electrodes; i.e., if the boat is moving toward the right to collect additional sounding curves, our interpretation of the inverted model should start at $x\sim -\mathrm{10}$ m. When analyzing Fig. 10b and e, we see that we are most sensitive to the water layer. Interestingly, when comparing Fig. 10c and f in detail, we realize that the sensitivity distribution in Fig. 10c reaches the IBPT interface while the sensitivity distribution in Fig. 10f is almost null for depths >12 m.

We perform the global sensitivity analyses (Sect. 4.5) considering 1D models described by five model parameters as used for the above-presented 1D inversions. We consider models where *z*_{w}=2 m, *ρ*_{w}=0.4 Ω m, and *ρ*_{uf}=5 Ω m are fixed, while *ρ*_{p} varies between 10 and 10 000 Ω m (eight values in total), and the IBPT is located at three different depths; i.e., *z*_{pt}=16, *z*_{pt}=10, and *z*_{pt}=4 m (Fig. 11a–c). For the calculation of the total sensitivity for each of our five parameters in the resulting 24 models, we set the parameter ranges to ${z}_{\mathrm{w}}=[\mathrm{0.5},\mathrm{3}]$ m, ${z}_{\mathrm{pt}}=[\mathrm{3.5},\mathrm{20}]$ m, ${\mathit{\rho}}_{\mathrm{w}}=[\mathrm{0.2},\mathrm{20}]$ Ω m, ${\mathit{\rho}}_{\mathrm{uf}}=[\mathrm{1},\mathrm{20}]$ Ω m, and ${\mathit{\rho}}_{\mathrm{p}}=[\mathrm{5},\mathrm{20}\phantom{\rule{0.125em}{0ex}}\mathrm{000}]$ Ω m. For these specific models and parameters ranges, our results (Fig. 11) suggest that *ρ*_{w} and *ρ*_{uf} are the most influential parameters, and the other three parameters (*z*_{pt}, *z*_{w} and *ρ*_{p}) are characterized in all cases by rather low total sensitivities. Furthermore, we also notice in Fig. 11a–c that *ρ*_{w} is the parameter showing the most significant changes when varying *ρ*_{p} and *z*_{pt}.

Knowledge of how fast permafrost thaws would improve predictive models of greenhouse gas release and coastal erosion, as well as coastal infrastructure design. The ERT method has been successfully used to image the unfrozen sediments overlying the permafrost layer in subsea permafrost environments, especially using smooth inversion approaches (e.g., Overduin et al., 2012; Pedrazas et al., 2020). In typical subsea permafrost environments, there might be a gradual transition zone consisting of a mixture of water and ice between fully unfrozen and frozen ice-bonded sediments. However, during ERT inversion, the nature of this transition can be either enlarged when using smooth inversion approaches or reduced to a single interface when using layer-based strategies. Whether we have a smooth or sharp transition between unfrozen and frozen sediments, there must be a threshold in the ice content that creates sufficient contrast in resistivity also influencing the penetration of the injected current and, thus, our apparent resistivity measurements (e.g., Kang and Lee, 2015). Because we wanted to target the interface defined by such a resistivity contrast (interpreted here as the IBPT), we considered a layer-based model parameterization to invert our ERT data. Additionally, we obtained estimates of uncertainties using an ensemble approach. For the sake of completeness, we provide the smooth inversion models for both of our field studies in Appendix A.

## 6.1 Insights from our parameterization and inversion strategies

We used a 2D layer-based model parameterization to globally invert marine ERT data and obtain different ensembles (e.g., after cluster analysis) of model solutions. We demonstrated with the two case studies that such ensembles allow us to reliably image the IBPT position with its associated uncertainties. The main advantage of using a layer-based model parameterization strategy is that we do not assume an arbitrary resistivity threshold or gradient to interpret the IBPT position, as we would need to do for our smooth inversion results (see Fig. A1a–b). This may be advantageous to compare ERT profiles collected at the same position in different years to track changes along the IBPT or in environments where the freezing point of the sediment porewater changes spatially. For example, offshore surveys that encounter submerged hypersaline lagoon deposits may show relatively low resistivity values for partially frozen sediments compared to colder ice-bonded permafrost with fresh porewater (Angelopoulos et al., 2021). Indeed, this interface may be related to a threshold in ice content. However, associating the IBPT with a certain ice content requires calibration by using borehole data or by additional geophysical information; e.g., by using the joint inversion approach of ERT and seismic refraction data of Wagner et al. (2019). Such thresholds may vary from site to site depending on properties of the sediments including temperature, grain size distribution, and the salinity of the porewater. Furthermore, we consider it convenient to use the sum of arctangent functions to parameterize the IBPT because there may be cases where the IBPT position varies steeply (as the ones we identified at the end of the median models in Figs. 2a–c and 7a–c) associated, for example, with submerged thermokarst structures or changes in the ratio of coastal erosion vs. degradation rate.

To find reliable and stable 2D model solutions, we performed experiments in which we ran the PSO several times to find the appropriate parameter settings for our layer-based ERT inversion. In a first stage, when performing our 2D inversions without considering any constraints, we found model solutions that were unrealistic according to our prior knowledge of our field sites. Therefore, we constrained our inversions considering our bathymetric and CTD measurements (see Sect. 4.2). Although the 2D inverted models without considering constraints are not shown in this study, we demonstrated with our 1D inversions how such constraints significantly improved the inversion results while reducing the number of possible solutions. Another exercise in the experimental phase consisted in using more than three layers in our parameterization strategy. However, we did not observe any significant improvement in the final median models when increasing the number of layers and, thus, restricted our inversion and analyses to three-layer scenarios.

One disadvantage of using a layer-based model parameterization relying on homogeneous layers is that it is not possible to resolve small-scale resistivity variations (e.g., horizontal heterogeneities at a spatial scale of meters). However, our workflow allowed us to inspect and evaluate model performance including the appropriateness of the model parameterization. For example, in the Bykovsky and Drew Point case studies, we observed in the residual pseudosections some regular lateral variations (see Figs. 3 and 8). This indicates that we were not completely explaining the data, because of either lateral subsurface resistivity variations or 3D effects. To tackle this problem, it could be beneficial to measure 3D bathymetric data around each ERT profile and collect additional parallel and perpendicular ERT profiles to better understand 3D resistivity variations at our field sites. Furthermore, direct measurement of the resistivity of water and unfrozen sediments (e.g., using additional water samples and drilling cores) might help to inform the model parameterization (e.g., account for lateral variations) and inversion strategies. We should notice that adding complexity to our model parameterization comes with the trade-off of increased computational cost to solve the inverse problems. An alternative to obtaining more complex resistivity models is to use our layer-based global inversion results as reference models to perform smoothness-constrained inversions (e.g., Günther et al., 2006).

The error level of ERT data is usually unknown – especially for marine data, where repeated or reciprocal measurements are not practical because the data are acquired while the boat is moving. This represents a challenge during the inversion when specifying an appropriate fit level. One alternative to get insights into the noise level is to perform repeated measurements in a static fashion (avoiding bending of the cable by wind or swells) for a certain section of the profile. For example, this can be achieved at the coast on a calm day where one end of the cable is secured to the beach and the other end is fastened to an anchored boat. However, such repeat measurements were not available for our field sites. Therefore, we set our stopping criterion by considering a fixed number of iterations rather than using a minimum threshold in our objective function. With this approach, we obtained model solutions characterized by different fit levels. For example, for the Bykovsky data, we found RMSLE values between 0.025 and 0.038 (Fig. 3g–i), while for our Drew Point data, we found RMSLE values between 0.007 and 0.016 (Fig. 8i–l). Although the RMSLE values for Drew Point are significantly lower than for Bykovsky, we found a family of models in the Drew Point study, which was considered as geologically unrealistic (Fig. 7d). This highlights the importance of estimating different ensembles of solutions with different fit levels and having an accurate estimate of data noise. Because the misfits for the model in Fig. 7d were higher, this family of models could potentially be discarded if they were found to exceed expected error levels without considering any prior knowledge of the environmental setting.

## 6.2 Parameter learning from 1D inversion

Our 2D inversion results showed large variations in the modeled resistivities of the permafrost, and we also noticed that, typically, the variabilities of IBPT position increase with depth. These observations indicate decreasing resolution capabilities of our ERT data with depth and limited penetration of the injected current in the frozen permafrost layers. To better understand these results in a more quantitative fashion, we reduced the number of parameters to five and performed selected 1D inversion experiments using synthetic data inspired by our 2D inversion results. Because such 1D inversions are significantly faster than 2D inversions, they represent an efficient way to explore the influence of constraining different parameters. For example, we noticed from our 1D inversion results that constraining the water layer significantly decreased the non-uniqueness of the inverse problem. This is essential for a reliable estimation of the IBPT position and for establishing petrophysical relations, for example, to estimate porewater salinity and ice content. Additionally, we noticed that the 1D inversion results for the Bykovsky data (Fig. 4c) provided similar uncertainties around the IBPT as the 2D inversion results at *x*=150 m (Fig. 2d–f). However, the 1D inversion results for the Drew Point data (Fig. 9c) showed uncertainties around the IBPT 3 times larger compared to the 2D inversion results at *x*≈600 m (Fig. 7f and g). This indicates that there is no general best way of using the results of such complementary synthetic 1D studies; the success and feasibility rather depends on the characteristics of the field site and analyzed data set. On the other hand, we can use our 1D inversion results to assess the posterior correlation matrix that, as we showed in our examples, can be helpful to identify interactions between the model parameters. Furthermore, comparing the changes across different posterior correlation matrices (e.g., associated with different model constraints) can help us detect changes in the parameters interactions and, thus, quantify the impact of our model constraints. Such straightforward but informative analysis provides a deeper understanding of the inversion process and the suitability of the entire inversion strategy.

Our 1D inversion results indicated some problems if we want to infer relative permafrost characteristics from ERT measurements. The 1D input models for our 1D synthetic examples (see Tables 2 and 3) assumed identical resistivities of the ice-bearing permafrost layer (*ρ*_{p}=4000 Ω m) and similar resistivities for the unfrozen sediments as found by our 2D inversion results. In contrast, the resistivity and depth of the seawater layer between both models were set according to field measurements at our field sites. Although the resistivities of the unfrozen and frozen layers were similar in both models, we noticed that for model scenarios derived from the Bykovsky site, the inverted *ρ*_{p} values were generally overestimated where already the *q*_{25} model indicates *ρ*_{p} values larger than the input *ρ*_{p} (Fig. 4c). On the contrary, for settings inspired by the Drew Point field site, the input *ρ*_{p} fell within the range defined by *q*_{25} and *q*_{75} models but showed more significant variabilities than our 1D Bykovsky experiment (Fig. 9c). These results demonstrated the influence of the depth and resistivity of the seawater layer in the inverted models which may be critical for subsequent interpretations. For example, assuming the same temperature and porewater salinity, the resistivity of the sediments increases with ice content (e.g., Pearson et al., 1986; Fortier et al., 1994; Kang and Lee, 2015). Thus, our results may lead us to conclude that the ice-bearing permafrost layer holds higher ice content at Bykovsky compared to Drew Point. Over- or underestimating the ice-bearing permafrost resistivity may lead to potentially erroneous interpretations, for example, related to the sediment's ice content, temperature, and composition. We would need complementary field information or further analyses like sensitivity assessments to avoid misleading interpretations.

## 6.3 System understanding with sensitivity analysis

We obtained an additional model understanding (e.g., in view of delineating confident and reliable model areas) by performing sensitivity analyses. From our examples, we learned that if the resistivity of the seawater were higher than the resistivity of the unfrozen sediments (as in the Bykovsky case study, Fig. 5), this would result in increased sensitivities inside the unfrozen sediments and, thus, to changes along the IBPT position. This type of situation may be prevalent in subsea permafrost areas affected by freshwater river discharge in summer. On the other hand, if the seawater were less resistive than the unfrozen sediments (e.g., as in the Drew Point case study, Fig. 10), we would be more sensitive to the water layer and, therefore, to bathymetric changes. This emphasizes the importance of accurate water depth measurements. We highlight the fact that although the local 2D sensitivities for the Drew Point data were rather small for the unfrozen sediments, the IQR of the models (Fig. 7f–g) showed equivalent variability around the depth of IBPT (1.5 to 2 m for depths ∼12 m) in comparison to the Bykovsky example (Fig. 2d–f), where the sensitivities showed a more pronounced influence within the unfrozen sediments.

This study used global sensitivity analysis considering only five parameters as needed for our 1D inversion examples. The Sobol approach proved to be a powerful method to distinguish the most influential parameters. After evaluating how the permafrost resistivity and the IBPT position may influence the rest of the parameters in our 1D three-layer examples, we noted some relevant differences. For example, in the Bykovsky example (Fig. 6), we noticed that for larger values of *ρ*_{p} and shallower *z*_{pt} the total influence of the rest of the parameters decreased. On the other hand, for the Drew Point example (Fig. 9), increasing *ρ*_{p} increased the total sensitivity of the rest of the parameters, while varying *z*_{pt} at shallower depths mainly increased the influence of *ρ*_{w}. We also want to highlight that *ρ*_{w}, *ρ*_{uf}, and *z*_{pt} were the parameters with the most significant total sensitivity in both examples and were also the parameters that formed model parameter pairs showing the largest correlation (see Figs. 4d and 9d). Encouragingly, *ρ*_{w} and *ρ*_{uf} can be informed from CTD casts and shallow sediment sampling, respectively. We must be aware that such a global sensitivity analysis is highly dependent on the predefined constraining parameter range and should be applied to address specific questions to allow, for example, parameter reduction or to guide our sampling strategies and experimental design.

## 6.4 Subsea permafrost features (Bykovsky vs. Drew Point)

The inverted ERT profiles yielded new insights into how subsea permafrost thaws because the Bykovsky Peninsula and Drew Point are characterized by distinct seawater properties and geological histories. The Bykovsky 2D inversion results at *x*=150 m, which corresponds to an inundation period of 357 years assuming an erosion rate of 0.42 m yr^{−1} (e.g., Lantuit et al., 2011), showed a median depth to the IBPT of ∼15 m (Fig. 2). This resulted in an average degradation rate of ∼0.04 m yr^{−1}. On the other hand, the Drew Point 2D inversion results at *x*≈600 m showed a median depth to the IBPT of ∼12 m. Note that this location coincides with the 1955 coastline position (see Fig. 1d–e), which corresponds to 63 years of inundation, yielding an average degradation rate of ∼0.19 m yr^{−1}. At Bykovsky, 63 years of inundation (again assuming an erosion rate of 0.42 m yr^{−1}) correspond to an offshore distance of ∼26 m, which corresponds to a median IBPT depth in the 2D inversion results at most 6 m (Fig. 2). Although the mean annual IBPT degradation rate slows with inundation time as the temperature gradient driving diffusive heat fluxes weakens (Angelopoulos et al., 2019), it is evident that the permafrost at Drew Point may thaw faster, presumably because Drew Point sediments are primed with salts in the pore space prior to inundation (Black, 1964; Sellmann, 1989).

Since salt diffusion is typically slower than heat diffusion (Harrison and Osterkamp, 1978), the IBPT degradation rate at Bykovsky should theoretically be faster than at Drew Point, provided that the permafrost sediments are similar. However, it appears that dissolved salts in the pore space of the sediments at Drew Point play an important role in lowering the permafrost freezing point and resulting in higher IBPT degradation rates than at Bykovsky. In fact, the top of onshore cryotic and saline unfrozen sediment layers (cryopegs) were observed near the Drew Point shoreline during coring (Bull et al., 2020; Bristol et al., 2021). This can lead us to interpret a faster IBPT degradation rate at Drew Point compared to Bykovsky in two ways: (1) a layer of submerged Drew Point sediments was already unfrozen upon inundation (e.g., *M*_{F2} in Fig. 7) and (2) the frozen layers at Drew Point contained less ice and had a lower freezing point. Jones et al. (2018) suggested that warming permafrost temperatures at Drew Point (3 to 4 ^{∘}C over the past several decades) have made saline permafrost more susceptible to erosion, potentially contributing to the enhanced coastal erosion rate (2.5 times that of the historical average) observed between 2007 and 2016. Warming by seawater submergence would presumably result in cryopeg spreading and IBPT degradation.

As shown in Fig. 1a and d, the coastal plains at our field sites consist of numerous thermokarst lakes and drained lake basins. When thermokarst lakes are breached by coastal erosion, the unfrozen sediments underneath the lake become integrated into the subsea permafrost environment, leading to bowl-shaped electrical resistivity structures. For example, Angelopoulos et al. (2021) showed steep IBPT gradients along ERT profiles parallel to the southern Bykovsky shoreline that traverse submerged thermokarst and undisturbed permafrost. These authors also suggested that drained lake basins, which have undergone thaw–refreeze cycles, are more susceptible to quicker thaw compared to undisturbed terrain. Comparing the first 400 m of our inverted median models for our field sites, we noticed that, in general, the IBPT at Drew Point is smoother than at Bykovsky. This might be the result of the higher erosion rates at Drew Point (>10 m yr^{−1}) than at Bykovsky (<1 m yr^{−1}) that expose coastal areas to inundation in a shorter time. Because of the longer inundation time at Bykovsky, we expect fluctuations in different environmental controls (e.g., water temperature, seawater salinity) that might result in step-like features as the one at *x*≈280 m. Furthermore, layered strata alternating between ice-rich and relatively ice-poor sediment may also contribute to step-like IBPT features. Similarly, in the Drew Point 2D inversion (Fig. 7a–c), there was a steep median IBPT decline observed at *x*≈750 m, where the IBPT deepened from ∼12 to ≥20 m. Although the resolution capabilities of our ERT data at these depths are limited, we suggest that thermokarst processes prior to seawater submergence may be responsible for the nature of this IBPT dip.

In this study, we illustrated how we could use ERT data to reliably estimate the IBPT position in shallow coastal areas of the Arctic. We found that using a layer-based model parameterization helps us target the IBPT position directly from the inversion of ERT data with the trade-off of omitting small-scale heterogeneities. To improve the inversion result, we noticed that constraining the water layer depth and resistivity reduces the non-uniqueness of the ERT inverse problem, improving the estimation of the resistivity of the unfrozen sediments (talik and/or cryopeg) and the IBPT position. However, even when constraining the water layer, we still found large variabilities in the resistivity of the frozen sediments. We suggest that constraining the resistivity of the unfrozen sediments (e.g., sediment sampling) during ERT inversion could improve resistivity estimates of the frozen layer and, thus, further permafrost's physical properties (e.g., ice content). Properly imaging the IBPT position may allow us to improve the estimation of the permafrost degradation rate, which might be used to better understand greenhouse gas emissions and coastal erosion processes. The workflow and methods presented in this study can guide future field campaigns and may be used as a reference for more detailed parameterizations and/or inversion strategies.

Although comparing different inversion strategies is beyond the scope of this study, for the sake of completeness, we show the smooth inversion results for our Bykovsky and Drew Point data sets in Fig. A1a–b. To invert these data sets, we use the inversion routine of the Python library pyGIMLi (Rücker et al., 2017). For both cases, we constrain the seawater layer by using the bathymetric profile data as collected by a Garmin echo sounder (see Sect. 3). The Bykovsky model (Fig. A1a) shows the highest resistivities to the left (near the coast), while resistivities drop in the offshore direction. For our Drew Point model (Fig. A1b), the highest resistivities are also present near the coast, but the resistivities decrease in the offshore direction more gradually than in the Bykovsky model. To derive the IBPT position from these models, we would need to assume or measure (e.g., through borehole data) the resistivity threshold that separates the talik from the ice-bearing permafrost layer (e.g., Overduin et al., 2016; Sherman et al., 2017; Angelopoulos et al., 2021). Because no ground-truth data are available and to avoid assuming a resistivity threshold, we decide to target such an interface using a layer-based model parameterization approach as explained in Sect. 4.1. Please refer to Angelopoulos et al. (2019), where a smooth inversion of the Bykovsky data set presented in this study is discussed in more detail, as well as Angelopoulos et al. (2021), where laterally constrained inversions of additional data sets offshore of Bykovsky are shown.

The data presented in this paper are freely available and can be downloaded from the PANGAEA website. Our Bykovsky and Drew Point data sets are published by Angelopoulos et al. (2018), https://doi.org/10.1594/PANGAEA.895887, and Angelopoulos et al. (2022), https://doi.org/10.1594/PANGAEA.948564, respectively.

This study was conceptualized by MAZ, MA, PPO, and JT. Data collection was performed by MA, GG, and BMJ. The calculations, visualization, and preparation of the original draft were performed by MAZ. The manuscript was reviewed and edited by all authors.

The contact author has declared that none of the authors has any competing interests.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

We thank Josefine Lenz and Juliane Wolter for their help during fieldwork at Drew Point and Mikhail N. Grigoriev, who helped to realize the entire Bykovsky study program. Finally, we sincerely thank the floatplane pilot Jim Webster for his years of service supporting fieldwork in Alaska.

Mauricio Arboleda-Zapata was supported by the Deutscher Akademischer Austauschdienst (grant no. 57395813) in the framework of the Research Training Group “Natural hazards and Risks in a Changing World” NatRiskChange 2043/2 funded by the Deutsche Forschungsgemeinschaft. Benjamin M. Jones was supported by the US National Science Foundation (grant no. OISE-1927553). Fieldwork was supported by a European Research Council project PETA-CARB grant (ERC-338335) to Guido Grosse and by Alfred Wegener Institute base funds.

This paper was edited by Adrian Flores Orozco and reviewed by two anonymous referees.

Akça, I. and Basokur, A. T.: Extraction of structure-based geoelectric models by hybrid genetic algorithms, Geophysics, 75, F15–F22, https://doi.org/10.1190/1.3273851, 2010. a, b

Angelopoulos, M.: Mapping subsea permafrost with electrical resistivity surveys, Nature Reviews Earth & Environment, 3, 6–6, https://doi.org/10.1038/s43017-021-00258-5, 2022. a

Angelopoulos, M., Westermann, S., Overduin, P. P., Faguet, A., Olenchenko, V., Grosse, G., and Grigoriev, M. N.: Conductivity, temperature and depth (CTD), snow and ice thickess and apparent resisitivity on the Bykovsky Peninsula, Lena Delta, in April and July 2017, PANGAEA [data set], https://doi.org/10.1594/PANGAEA.895887, 2018. a

Angelopoulos, M., Westermann, S., Overduin, P. P., Faguet, A., Olenchenko, V., Grosse, G., and Grigoriev, M. N.: Heat and Salt Flow in Subsea Permafrost Modeled with CryoGRID2, J. Geophys. Res.-Earth, 124, 920–937, https://doi.org/10.1029/2018JF004823, 2019. a, b, c, d

Angelopoulos, M., Overduin, P. P., Miesner, F., Grigoriev, M. N., and Vasiliev, A. A.: Recent advances in the study of Arctic submarine permafrost, Permafrost Periglac., 31, 442–453, https://doi.org/10.1002/ppp.2061, 2020a. a, b

Angelopoulos, M., Overduin, P. P., Westermann, S., Tronicke, J., Strauss, J., Schirrmeister, L., Biskaborn, B. K., Liebner, S., Maksimov, G., Grigoriev, M. N., and Grosse, G.: Thermokarst Lake to Lagoon Transitions in Eastern Siberia: Do Submerged Taliks Refreeze?, J. Geophys. Res.-Earth, 125, e2019JF00542, https://doi.org/10.1029/2019JF005424, 2020b. a

Angelopoulos, M., Overduin, P. P., Jenrich, M., Nitze, I., Günther, F., Strauss, J., Westermann, S., Schirrmeister, L., Kholodov, A., Krautblatter, M., Grigoriev, M. N., and Grosse, G.: Onshore Thermokarst Primes Subsea Permafrost Degradation, Geophys. Res. Lett., 48, e2021GL093881, https://doi.org/10.1029/2021gl093881, 2021. a, b, c, d, e, f

Angelopoulos, M., Arboleda-Zapata, M., Overduin, P. P., Jones, B. M., Tronicke, J., and Grosse, G.: Marine electrical resistivity data offshore of Drew Point, Alaska, PANGAEA [data set], https://doi.org/10.1594/PANGAEA.948564, 2022. a

Arboleda-Zapata, M., Guillemoteau, J., and Tronicke, J.: A comprehensive workflow to analyze ensembles of globally inverted 2D electrical resistivity models, J. Appl. Geophys., 196, 104512, https://doi.org/10.1016/j.jappgeo.2021.104512, 2022. a, b, c, d, e

Are, F.: Shoreface of the Arctic seas – A natural laboratory for subsea permafrost dynamics, in: Proceedings of the Eighth International Conference on Permafrost, Zürich, Switzerland, edited by: Philips, M., Springman, S. M., and Arenson, L. U., Swets & Zeitlinger, Lisse, 21–25 July 2003, 27–32, https://oceanrep.geomar.de/id/eprint/26899/ (last access: 20 June 2021), 2003. a, b, c

Auken, E. and Christiansen, A. V.: Layered and laterally constrained 2D inversion of resistivity data, Geophysics, 69, 752–761, https://doi.org/10.1190/1.1759461, 2004. a, b

Barnhart, K. R., Anderson, R. S., Overeem, I., Wobus, C., Clow, G. D., and Urban, F. E.: Modeling erosion of ice-rich permafrost bluffs along the Alaskan Beaufort Sea coast, J. Geophys. Res.-Earth, 119, 1155–1179, https://doi.org/10.1002/2013JF002845, 2014. a

Black, R. F.: Gubik formation of Quaternary age in northern Alaska, Tech. rep., US Government Printing Office, 1964. a, b

Bristol, E. M., Connolly, C. T., Lorenson, T. D., Richmond, B. M., Ilgen, A. G., Choens, R. C., Bull, D. L., Kanevskiy, M., Iwahana, G., Jones, B. M., and McClelland, J. W.: Geochemistry of Coastal Permafrost and Erosion-Driven Organic Matter Fluxes to the Beaufort Sea Near Drew Point, Alaska, Front. Earth Sci., 8, 598933, https://doi.org/10.3389/feart.2020.598933, 2021. a

Brothers, L. L., Herman, B. M., Hart, P. E., and Ruppel, C. D.: Subsea ice-bearing permafrost on the U.S. Beaufort Margin: 1. Minimum seaward extent defined from multichannel seismic reflection data, Geochem. Geophy. Geosy., 17, 4354–4365, https://doi.org/10.1002/2016GC006584, 2016. a

Bull, D. L., Bristol, E. M., Brown, E., Choens, R. C., Connolly, C. T., Flanary, C., Jones, B. F. J., Jones, C. A., Jones, M. W., McClelland, J., Mota, A., and Tezaur, I.: Arctic Coastal Erosion: Modeling and Experimentation (No. SAND2020-10223), Tech. rep., Sandia National Labs. (SNL-NM), Alburqueque, NM, USA, Sandia National Lab. (SNL-CA), Livermore, CA, USA, https://doi.org/10.2172/1670531, 2020. a

Caliński, T. and Harabasz, J.: A dendrite method for cluster analysis, Commun. Stat. Theory, 3, 1–27, https://doi.org/10.1080/03610927408827101, 1974. a

Chambers, J. E., Wilkinson, P. B., Wardrop, D., Hameed, A., Hill, I., Jeffrey, C., Loke, M. H., Meldrum, P. I., Kuras, O., Cave, M., and Gunn, D. A.: Bedrock detection beneath river terrace deposits using three-dimensional electrical resistivity tomography, Geomorphology, 177–178, 17–25, https://doi.org/10.1016/j.geomorph.2012.03.034, 2012. a

De Pasquale, G., Linde, N., Doetsch, J., and Holbrook, W. S.: Probabilistic inference of subsurface heterogeneity and interface geometry using geophysical data, Geophys. J. Int., 217, 816–831, https://doi.org/10.1093/gji/ggz055, 2019. a

Fernández-Martínez, J. L., Xu, S., Sirieix, C., Fernández-Muniz, Z., and Riss, J.: Uncertainty analysis and probabilistic segmentation of electrical resistivity images: the 2D inverse problem, Geophys. Prospect., 65, 112–130, https://doi.org/10.1111/1365-2478.12559, 2017. a

Fortier, R., Allard, M., and Seguin, M. K.: Effect of physical properties of frozen ground on electrical resistivity logging, Cold Reg. Sci. Technol., 22, 361–384, https://doi.org/10.1016/0165-232X(94)90021-3, 1994. a, b

Frederick, J. M. and Buffett, B. A.: Effects of submarine groundwater discharge on the present-day extent of relict submarine permafrost and gas hydrate stability on the Beaufort Sea continental shelf, J. Geophys. Res.-Earth, 120, 417–432, https://doi.org/10.1002/2014JF003349, 2015. a

Fuchs, M., Palmtag, J., Juhls, B., Overduin, P. P., Grosse, G., Abdelwahab, A., Bedington, M., Sanders, T., Ogneva, O., Fedorova, I. V., Zimov, N. S., Mann, P. J., and Strauss, J.: High-resolution bathymetry models for the Lena Delta and Kolyma Gulf coastal zones, Earth Syst. Sci. Data, 14, 2279–2301, https://doi.org/10.5194/essd-14-2279-2022, 2022. a

Gebrande, H.: A seismic-ray tracing method for two dimensional inhomogeneous media, in: Explosion seismology in Central Europe: data and results, 162–167, Springer-Verlag, New York, https://doi.org/10.1007/978-3-642-66403-8_23, 1976. a

Grigoriev, M. N.: Kriomorfogenez i litodinamika pribrezhno- shelfovoi zony morei Vostochnoi Sibiri (Cryomorphogenesis and lithodynamics of the East Siberian near-shore shelf zone), PhD thesis, Melnikov Permafrost Institute, Russian Academy of Sciences, Siberian Branch, Yakutsk, 2008. a

Grosse, G., Schirrmeister, L., Siegert, C., Kunitsky, V. V., Slagoda, E. A., Andreev, A. A., and Dereviagyn, A. Y.: Geological and geomorphological evolution of a sedimentary periglacial landscape in Northeast Siberia during the Late Quaternary, Geomorphology, 86, 25–51, https://doi.org/10.1016/j.geomorph.2006.08.005, 2007. a, b, c

Günther, T., Rücker, C., and Spitzer, K.: Three-dimensional modelling and inversion of DC resistivity data incorporating topography – II. Inversion, Geophys. J. Int., 166, 506–517, https://doi.org/10.1111/j.1365-246X.2006.03011.x, 2006. a, b, c

Harrison, W. D. and Osterkamp, T. E.: Heat and mass transport processes in subsea permafrost 1. An analysis of molecular diffusion and its consequences, J. Geophys. Res., 83, 4707–4712, https://doi.org/10.1029/JC083iC09p04707, 1978. a, b

Herman, J. and Usher, W.: SALib: an open-source Python library for sensitivity analisis, Journal of Open Source Software, 2, 97, https://doi.org/10.21105/joss.00097, 2017. a

Hsu, H. L., Yanites, B. J., chih Chen, C., Chen, Y. G., Chen, C. C., and Chen, Y. G.: Bedrock detection using 2D electrical resistivity imaging along the Peikang River, central Taiwan, Geomorphology, 114, 406–414, https://doi.org/10.1016/j.geomorph.2009.08.004, 2010. a

Hubbard, S. S., Gangodagamage, C., Dafflon, B., Wainwright, H., Peterson, J., Gusmeroli, A., Ulrich, C., Wu, Y., Wilson, C., Rowland, J., Tweedie, C., and Wullschleger, S. D.: Quantifying and relating land-surface and subsurface variability in permafrost environments using LiDAR and surface geophysical datasets, Hydrogeol. J., 21, 149–169, https://doi.org/10.1007/s10040-012-0939-y, 2013. a

Jones, B. M., Hinkel, K. M., Arp, C. D., and Eisner, W. R.: Modern erosion rates and loss of coastal features and sites, Beaufort Sea coastline, Alaska, Arctic, 61, 361–372, https://doi.org/10.14430/arctic44, 2008. a

Jones, B. M., Arp, C. D., Jorgenson, M. T., Hinkel, K. M., Schmutz, J. A., and Flint, P. L.: Increase in the rate and uniformity of coastline erosion in Arctic Alaska, Geophys. Res. Lett., 36, L03503, https://doi.org/10.1029/2008GL036205, 2009. a

Jones, B. M., Farquharson, L. M., Baughman, C. A., Buzard, R. M., Arp, C. D., Grosse, G., Bull, D. L., Günther, F., Nitze, I., Urban, F., Kasper, J. L., Frederick, J. M., Thomas, M., Jones, C., Mota, A., Dallimore, S., Tweedie, C., Maio, C., Mann, D. H., Richmond, B., Gibbs, A., Xiao, M., Sachs, T., Iwahana, G., Kanevskiy, M., and Romanovsky, V. E.: A decade of remotely sensed observations highlight complex processes linked to coastal permafrost bluff erosion in the Arctic, Environ. Res. Lett., 13, 115001, https://doi.org/10.1088/1748-9326/aae471, 2018. a, b, c, d

Kang, M. and Lee, J. S.: Evaluation of the freezing-thawing effect in sand-silt mixtures using elastic waves and electrical resistivity, Cold Reg. Sci. Technol., 113, 1–11, https://doi.org/10.1016/j.coldregions.2015.02.004, 2015. a, b, c

Kennedy, J. and Eberhart, R.: Particle swarm optimization, in: Proceedings of the IEEE International Joint Conference on Neural Networks, 27 November 1995–1 December 1995, Perth, WA, Australia, 1942–1948, https://doi.org/10.1109/ICNN.1995.488968, 1995. a

Kneisel, C., Hauck, C., Fortier, R., and Brian, M.: Advances in geophysical methods for permafrost investigations, Permafrost Periglac., 19, 157–178, https://doi.org/10.1002/ppp.616, 2008. a, b

Koren, Z., Mosegaard, K., Landa, E., Thore, P., and Tarantola, A.: Monte Carlo estimation and resolution analysis of seismic background velocities, J. Geophys. Res.-Sol. Ea., 96, 20289–20299, https://doi.org/10.1029/91JB02278, 1991. a

Lantuit, H., Atkinson, D., Overduin, P. P., Grigoriev, M., Rachold, V., Grosse, G., and Hubberten, H. W.: Coastal erosion dynamics on the permafrost-dominated Bykovsky Peninsula, north Siberia, 1951–2006, Polar Res., 30, 7341, https://doi.org/10.3402/polar.v30i0.7341, 2011. a, b, c, d, e, f

Loke, M. H. and Barker, R. D.: Rapid least-squares inversion of apparent resistivity pseudosections by a quasi-Newton method, Geophys. Prospect., 44, 131–152, https://doi.org/10.1111/j.1365-2478.1996.tb00142.x, 1996. a

Lück, A.: Subsea permafrost degradation and coastline retreat at Drew Point (Alaska), Master Thesis, University of Potsdam, Potsdam, Germany, 2020. a

MacQueen, J.: Some methods for classification and analysis of multivariate observations, in: Proceedings of the fifth Berkeley symposium on mathematical statistics and probability, 21 June–18 July 1965 and 27 December 1965–7 January 1966, University of California, Berkeley, 281–297, https://projecteuclid.org/euclid.bsmsp/1200512992 (last access: 22 October 2020), 1967. a

Malinverno, A.: Parsimonious Bayesian Markov chain Monte Carlo inversion in a nonlinear geophysical problem, Geophys. J. Int., 151, 675–688, https://doi.org/10.1046/j.1365-246X.2002.01847.x, 2002. a

Osterkamp, T.: Sub-sea Permafrost, in: Encyclopedia of Ocean Sciences, 2902–2912, Academic Press, New York, and London, https://doi.org/10.1006/rwos.2001.0008, 2001. a

Overduin, P. P., Westermann, S., Yoshikawa, K., Haberlau, T., Romanovsky, V., and Wetterich, S.: Geoelectric observations of the degradation of nearshore submarine permafrost at Barrow (Alaskan Beaufort Sea), J. Geophys. Res.-Earth, 117, F02004, https://doi.org/10.1029/2011JF002088, 2012. a, b, c, d, e, f, g, h

Overduin, P. P., Haberland, C., Ryberg, T., Kneier, F., Jacobi, T., Grigoriev, M. N., and Ohrnberger, M.: Submarine permafrost depth from ambient seismic noise, Geophys. Res. Lett., 42, 7581–7588, https://doi.org/10.1002/2015GL065409, 2015a. a

Overduin, P. P., Liebner, S., Knoblauch, C., Günther, F., Wetterich, S., Schirrmeister, L., Hubberten, H.-W., and Grigoriev, M. N.: Methane oxidation following submarine permafrost degradation: Measurements from a central Laptev Sea shelf borehole, J. Geophys. Res.-Biogeo., 120, 965–978, https://doi.org/10.1002/2014JG002862, 2015b. a

Overduin, P. P., Wetterich, S., Günther, F., Grigoriev, M. N., Grosse, G., Schirrmeister, L., Hubberten, H.-W., and Makarov, A.: Coastal dynamics and submarine permafrost in shallow water of the central Laptev Sea, East Siberia, The Cryosphere, 10, 1449–1462, https://doi.org/10.5194/tc-10-1449-2016, 2016. a, b, c, d, e, f, g, h

O'Neill, H. B., Roy-Leveillee, P., Lebedeva, L., Ling, F., O'Neill, H. B., Roy-Leveillee, P., Lebedeva, L., and Ling, F.: Recent advances (2010–2019) in the study of taliks, Permafrost Periglac., 31, 346–357, https://doi.org/10.1002/ppp.2050, 2020. a

Pearson, C., Murphy, J., and Hermes, R.: Acoustic and resistivity measurements on rock samples containing tetrahydrofuran hydrates: Laboratory analogues to natural gas hydrate deposits, J. Geophys. Res., 91, 132–14, https://doi.org/10.1029/jb091ib14p14132, 1986. a, b

Pedrazas, M. N., Bayani Cardenas, M., Demir, C., Watson, J. A., Connolly, C. T., and McClelland, J. W.: Absence of ice-bonded permafrost beneath an Arctic lagoon revealed by electrical geophysics, Sci. Adv., 6, eabb5083, https://doi.org/10.1126/sciadv.abb5083, 2020. a, b, c, d

Permafrost Subcommittee: Glossary of permafrost and related ground-ice terms, Associate Committee on Geotechnical Research, National Research Council of Canada, Ottawa, https://globalcryospherewatch.org/reference/glossary_docs/permafrost_and_ground_terms_canada.pdf (last access: 1 July 2021), 1988. a, b

Ping, C. L., Michaelson, G. J., Guo, L., Jorgenson, M. T., Kanevskiy, M., Shur, Y., Dou, F., and Liang, J.: Soil carbon and material fluxes across the eroding Alaska Beaufort Sea coastline, J. Geophys. Res.-Biogeo., 116, G02004, https://doi.org/10.1029/2010JG001588, 2011. a

Rangel, R. C., Parsekian, A. D., Farquharson, L. M., Jones, B. M., Ohara, N., Creighton, A. L., Gaglioti, B. V., Kanevskiy, M., Breen, A. L., Bergstedt, H., Romanovsky, V. E., and Hinkel, K. M.: Geophysical Observations of Taliks Below Drained Lake Basins on the Arctic Coastal Plain of Alaska, J. Geophys. Res.-Sol. Ea., 126, e2020JB020889, https://doi.org/10.1029/2020JB020889, 2021. a

Rekant, P., Bauch, H. A., Schwenk, T., Portnov, A., Gusev, E., Spiess, V., Cherkashov, G., and Kassens, H.: Evolution of subsea permafrost landscapes in Arctic Siberia since the Late Pleistocene: a synoptic insight from acoustic data of the Laptev Sea, Arktos, 1, 11, https://doi.org/10.1007/s41063-015-0011-y, 2015. a

Roy, A., Dubey, C. P., and Prasad, M.: Gravity inversion of basement relief using Particle Swarm Optimization by automated parameter selection of Fourier coefficients, Comput. Geosci., 156, 104875, https://doi.org/10.1016/j.cageo.2021.104875, 2021. a

Roy, L., Sen, M. K., McIntosh, K., Stoffa, P. L., and Nakamura, Y.: Joint inversion of first arrival seismic travel-time and gravity data, J. Geophys. Eng., 2, 277–289, https://doi.org/10.1088/1742-2132/2/3/011, 2005. a

Rücker, C., Günther, T., and Wagner, F. M.: pyGIMLi: An open-source library for modelling and inversion in geophysics, Comput. Geosci., 109, 106–123, https://doi.org/10.1016/j.cageo.2017.07.011, 2017. a, b, c

Rumpf, M. and Tronicke, J.: Assessing uncertainty in refraction seismic traveltime inversion using a global inversion strategy, Geophys. Prospect., 63, 1188–1197, https://doi.org/10.1111/1365-2478.12240, 2015. a

Ruppel, C. D. and Kessler, J. D.: The interaction of climate change and methane hydrates, Rev. Geophys., 55, 126–168, https://doi.org/10.1002/2016RG000534, 2017. a

Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., Saisana, M., and Tarantola, S.: Global sensitivity analysis: The primer, John Wiley and Sons, https://doi.org/10.1002/9780470725184, 2008. a

Sayedi, S. S., Abbott, B. W., Thornton, B. F., Frederick, J. M., Vonk, J. E., Overduin, P. P., Schädel, C., Schuur, E. A., Bourbonnais, A., Demidov, N., Gavrilov, A., He, S., Hugelius, G., Jakobsson, M., Jones, M. C., Joung, D. J., Kraev, G., Macdonald, R. W., David McGuire, A., Mu, C., O’Regan, M., Schreiner, K. M., Stranne, C., Pizhankova, E., Vasiliev, A., Westermann, S., Zarnetske, J. P., Zhang, T., Ghandehari, M., Baeumler, S., Brown, B. C., and Frei, R. J.: Subsea permafrost carbon stocks and climate change sensitivity estimated by expert assessment, Environ. Res. Lett., 15, 124075, https://doi.org/10.1088/1748-9326/abcc29, 2020. a

Schirrmeister, L., Siegert, C., Kunitzky, V. V., Grootes, P. M., and Erlenkeuser, H.: Late quaternary ice-rich permafrost sequences as a paleoenvironmental archive for the Laptve Sea Region in Northern Siberia, Int. J. Earth Sci., 91, 154–167, https://doi.org/10.1007/s005310100205, 2002. a, b

Schirrmeister, L., Grigoriev, M. N., Strauss, J., Grosse, G., Overduin, P. P., Kholodov, A., Guenther, F., and Hubberten, H.-W.: Sediment characteristics of a thermokarst lagoon in the northeastern Siberian Arctic (Ivashkina Lagoon, Bykovsky Peninsula), Arktos, 4, 13, https://doi.org/10.1007/s41063-018-0049-8, 2018. a

Schuur, E. A., McGuire, A. D., Schädel, C., Grosse, G., Harden, J. W., Hayes, D. J., Hugelius, G., Koven, C. D., Kuhry, P., Lawrence, D. M., Natali, S. M., Olefeldt, D., Romanovsky, V. E., Schaefer, K., Turetsky, M. R., Treat, C. C., and Vonk, J. E.: Climate change and the permafrost carbon feedback, Nature, 520, 171–179, https://doi.org/10.1038/nature14338, 2015. a

Scott, W. J., Sellmann, P. V., and Hunter, J. A.: Geophysics in the study of permafrost, in: Geotechnical and environmental geophysics, edited by: Ward, S. H., vol. I, Society of Exploration Geophysicists, https://doi.org/10.1190/1.9781560802785.ch13, 1990. a

Sellmann, P. V.: Seafloor temperature and conductivity data from coastal waters of the US Beaufort Sea, Tech. Rep. 1, US Army Corps of Engineers, Cold Regions Research and Engineering Laboratory, 1989. a

Sellmann, P. V., Delaney, A. J., and Arcone, S. A.: Coastal subsea permafrost and bedrock observations using dc resistivity, Tech. rep., US Army Corps of Engineers, Cold Regions Research and Engineering Laboratory, 1989. a, b, c, d

Sen, M. K. and Stoffa, P. L.: Bayesian inference, Gibbs' sampler and uncertainty estimation in geophysical inversion, Geophys. Prospect., 44, 313–350, https://doi.org/10.1111/j.1365-2478.1996.tb00152.x, 1996. a

Sen, M. K. and Stoffa, P. L.: Global optimization methods in geophysical inversion (2nd ed.), Cambridge University Press, 2013. a

Shakhova, N., Semiletov, I., Gustafsson, O., Sergienko, V., Lobkovsky, L., Dudarev, O., Tumskoy, V., Grigoriev, M., Mazurov, A. K., Salyuk, A., Ananiev, R., Koshurnikov, A., Kosmach, D., Charkin, A., Dmitrevsky, N., Karnaukh, V., Gunar, A., Meluzov, A., and Chernykh, D.: Current rates and mechanisms of subsea permafrost degradation in the East Siberian Arctic Shelf, Nat. Commun., 8, 15872, https://doi.org/10.1038/ncomms15872, 2017. a

Sherman, D., Kannberg, P., and Constable, S.: Surface towed electromagnetic system for mapping of subsea Arctic permafrost, Earth Planet. Sc. Lett., 460, 97–104, https://doi.org/10.1016/j.epsl.2016.12.002, 2017. a, b, c

Shewchuk, J. R.: Triangle: Engineering a 2D quality mesh generator and Delaunay triangulator, in: Applied Computational Geometry Towards Geometric Engineering, edited by: Lin, M. C. and Manocha, D., 1148, 203–222 Springer, https://doi.org/10.1007/BFb0014497, 1996. a

Sobol, I.: Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates, Math. Comput. Simulat., 55, 271–280, https://doi.org/10.1016/S0378-4754(00)00270-6, 2001. a

Spitzer, K.: The three-dimensional DC sensitivity for surface and subsurface sources, Geophys. J. Int., 134, 736–746, https://doi.org/10.1046/j.1365-246X.1998.00592.x, 1998. a

Stevens, C. W., Moorman, B. J., Solomon, S. M., and Hugenholtz, C. H.: Mapping subsurface conditions within the near-shore zone of an Arctic delta using ground penetrating radar, Cold Reg. Sci. Technol., 56, 30–38, https://doi.org/10.1016/j.coldregions.2008.09.005, 2009. a

Tronicke, J., Paasche, H., and Böniger, U.: Crosshole traveltime tomography using particle swarm optimization: A near-surface field example, Geophysics, 77, R19–R32, https://doi.org/10.1190/geo2010-0411.1, 2012. a

Wagner, F., Mollaret, C., Günther, T., Kemna, A., and Hauck, C.: Quantitative imaging of water, ice and air in permafrost systems through petrophysical joint inversion of seismic refraction and electrical resistivity data, Geophys. J. Int., 219, 1866–1875, 2019. a

Wainwright, H. M., Finsterle, S., Jung, Y., Zhou, Q., and Birkholzer, J. T.: Making sense of global sensitivity analyses, Comput. Geosci., 65, 84–94, https://doi.org/10.1016/j.cageo.2013.06.006, 2014. a

Wu, Y., Nakagawa, S., Kneafsey, T. J., Dafflon, B., and Hubbard, S.: Electrical and seismic response of saline permafrost soil during freeze – Thaw transition, J. Appl. Geophys., 146, 16–26, https://doi.org/10.1016/j.jappgeo.2017.08.008, 2017. a