Exploring the capabilities of electrical resistivity tomography to study subsea permafrost

. Sea level rise and coastal erosion have inundated large areas of Arctic permafrost. Submergence by warmer 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 5 (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 10 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 as well as for the design of future field

. The higher the ice content, the less conductive is the medium (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 more conductive layers (e.g., Spitzer, 1998) resulting in a limited penetration of the current system into the frozen sediment layers. 70 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 75 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 80 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) 85 to study submarine permafrost environments of the Arctic in terms of the resistivity distribution of the unfrozen and frozen sediments including estimates of uncertainties, also around the IBPT. We analyze and compare ERT data collected at two field sites in the Arctic characterized by different environmental conditions (e.g., regarding seawater depth and resistivity, coastal erosion rate, and the sediments porewater chemistry), which have been recorded using different acquisition settings (e.g., different electrode spacing and spatial sampling). Additionally, we generate and interpret ensembles of globally inverted 90 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 analysis to recognize the most influential parameters during 2D and 1D inversions.
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 efficiently measure also profiles with a length in the order of km. 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 are given relative to boat position), vertical oscillation of electrodes (e.g., due to wavy conditions), and 130 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 135 during two fieldwork campaigns in July 2017 and July 2018, respectively. The two ERT data sets were collected using an Iris TM 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 pair of electrodes (e.g., Overduin et al., 2012). The transmitter voltage was set at approximately 48 V at Bykovsky while, at Drew 140 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 (∼ 1.1 m/s) perpendicular to the 145 shore. The Bykovsky soundings were collected at a spacing of five to seven meters (540 measurements in total, 418 m long) and the Drew Point soundings at one to two meters (1, 830 measurements in total, 854 m long). 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 f). Furthermore, we measured water conductivity and temperature at different depths (the mean values are shown 150 in the last two rows of Table 1) using a Sontek TM CastAway device (also known as CTD). In general, at our Drew Point field site, the seawater was shallower, more conductive, and slightly cooler than at our Bykovsky field site.
The measured apparent resistivities ρ a are presented as pseudosections in Fig. 1c and g. Here the x coordinates represent the center position of each quadripole, and the vertical axis (levels) represent the relative penetration; i.e., level = 0 is the shortest 155 quadripole, while level = 9 is the quadripole with maximum electrode spacing. The range of ρ a for Bykovsky is 5.9 Ωm to 45 Ωm and for Drew Point 0.9 Ωm to 6.3 Ωm. The lower ρ a values at Drew Point are mainly due to the higher conductivity of the seawater at the Alaskan coast, which is less influenced by freshwater discharge from large rivers than at our Bykovsky field site. Table 1. Acquisition parameters for our ERT data sets and further site-specific information for our two field sites.

Bykovsky, Siberia Drew Point, Alaska
Number of electrodes 13 13 Electrode spacing (m) 10 5 Transmitter voltage (V) 48 24 Sounding separation (m) 5 to 7 1 to 2 Length of profile (m) 418 854 Number of data points 540 1830 Water resistivity (Ωm) 13.7 0.5 Water temperature ( • C) 7 5.5 Additionally, we found it useful to plot the data as 1D sounding curves ( Fig. 1d and h) because this allows for additional visual inspections of data quality. For example, we noticed for our Bykovsky data that levels seven and eight are characterized by 160 higher noise levels than all other levels. In contrast, the Drew Point data not show obvious variations in data quality depending on the level number (i.e., along the 1D sounding curves).

Methodology
In this study, we follow the workflow of Arboleda-Zapata et al. (2022) who propose a layer-based model parameterization to 165 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.

2D layer-based model parameterization
One of the most important steps in any geophysical inversion workflow is defining a model parameterization that can prop-170 erly 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 may be characterized by one resistivity value. For more complex geological settings, one might allow for lateral 175 and/or vertical variation 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 functions based, for example, on spline interpolation (e.g., Koren et al., 1991), Fourier coefficients (e.g., Roy et al., 2021), or sums of arctangent functions  (Gebrande, 1976). For this study, we adopt a strategy based on the sum of arctangent functions because it allows for abrupt and smooth changes along the interfaces (e.g., Roy et al., 2005;Rumpf and Tronicke, 2015). Allowing for abrupt changes is 180 considered to be important in permafrost environments where high structural variability is often found; for example, 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).

Inversion strategy
During inversion, we search for a combination of model parameters (i.e., those describing the geometry of interfaces and the resistivities of the homogeneous layers) that minimizes the root mean squared logarithmic error (RMSLE). 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 190 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).

195
In a first step, the PSO requires defining a set of particles (where each particle represents a different model) that are initialized with random parameters (bounded within realistic physical ranges). This defines our model space. In the following iterations, the position of each particle is updated 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 200 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 in a different finite-element model mesh using the python library pyGIMLi (Rücker et al., 2017). This allows for local refinements around the interfaces and, thus, 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 205 best model. Using different seeds of the random number generator, we repeat this process until we obtain an ensemble M F 0 consisting of several hundred independent models and an ensemble of corresponding residuals δ F 0 . In this study, each residual vector is calculated as the difference between the observed and the corresponding modeled log-apparent resistivity values.

Ensemble interpretation
In a first step, to ease our ensemble analysis and interpretation in a pixel-wise fashion, all models in M F 0 are interpolated using 210 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 F 0 and δ F 0 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 F i and δ F i (where i = 0, 1, ..., n k , note i = 0 correspond to the whole 215 ensemble and i > 0 to the clustered ensembles) by the median values µ 1/2 (M F i ) and µ 1/2 (δ F i ) and the interquartile ranges IQR(M F i ) and IQR(δ F i ). Additionally, we describe δ F i in an overall fashion assessing the RMSLE(δ F i ), the IQR(δ F i ), and the quantile 90 % q 90 (δ F i ).

1D inversion
Often, we prefer 2D inversion algorithms in comparison to 1D strategies; especially for field data where the subsurface situ-220 ation 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 225 frozen and unfrozen 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 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 230 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 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).

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 high resistivity contrasts (up to several orders of magnitude) 240 between frozen and unfrozen 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.
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 250 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 .

Results
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.

260
The geological settings of the Bykovsky area are described in Sect. 2.1 and a summary of the acquisition parameters are 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. We parameterize our model with two interfaces (considering five nodes for each interface), one for the seabed and the other for the IBPT. In total, this results in 36 model parameters, which describe the structure of the interfaces and the layer resistivities. In 265 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 hours 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 F 0 consisting of 690 models.

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  models, we note a bimodal distribution of ρ p (some models with ρ p < 2, 000 Ωm and others with ρ p > 100, 000 Ωm) which is also illustrated by increased IQR(M F 0 ) values for the lowermost layer. These observations already indicates 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 F 0 can be divided into two model 280 families (M F 1 and M F 2 ). In Fig. 2b-c and e-f, we present the µ 1/2 (M F i ) and IQR(M F i ) models (where i = 1, 2). Comparing these models illustrates that M F 1 and M F 2 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 F 1 is shallower than in M F 2 . We learn from this that if ρ p increases, the depth of the IBPT increases (resulting in thicker unfrozen sediments), also near the coast.
According to the depth of the IBPT and its gradients in profile direction, we laterally subdivide the models into three main 285 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 130 < x < 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, which is characterized by a rather deep IBPT (with depths > 20 m) and extends until the end of the profile. 290 We assess the fit performance in a pixel-wise and in an overall fashion for the residuals associated to the ensemble containing all models δ F 0 , as well as for the two clustered model families δ F 1 and δ F 2 (Fig. 3). Thus, we calculate µ 1/2 (δ F i ) and IQR(δ F i ) (where i = 0, 1, 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 µ 1/2 (δ F i ) and IQR(δ F i ) 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 295 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 (a heat source) prior to seawater submergence, lateral temperature gradients and heterogeneous sediment properties could affect subsurface resistivity and its 3D variations. The overall statistics RMSLE(δ F i ), IQR(δ F i ), 300 and q 90 (δ F i ) (where i = 0, 1, 2) are presented as histograms in Fig. 3g-i. The histograms are characterized by bimodal distributions, especially evident in all shown RMSLE(δ F i ) histograms. When comparing the histograms of δ F 1 and δ F 2 , we notice that they follow similar distributions (although in δ F 1 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.

1D 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  rameters were chosen by analyzing our 2D model solutions (e.g., Fig. 2b-c at x ≃ 150 m). We calculate the forward response 310 of 10 quadripoles (which is equivalent to one sounding curve in Fig. 1d) 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 315 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 because the quantile 25 % and 75 % models (q 25 , q 75 ) indicate resistivities larger than 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 [z pt , ρ uf ] is the model parameter pair with the highest positive correlation 320 while the rest of the model parameter pairs are characterized by negative correlation with different amplitudes (except for pairs with ρ p which show correlations approaching zero).

Sensitivity analysis
To understand the sensitivity distribution for our three-layer model (representing seawater and unfrozen sediments overlying 325 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 can interpret our inverted models even underneath the outer electrode positions; i.e., if the boat together 330 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 ∼ −60 m (a more conservative model interpretation might start at x ∼ −25 m). 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 335 unfrozen sediments at x < −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 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 higher range than z w . Furthermore, although we allow ρ p to vary over three orders of magnitude the result of this sensitivity analysis demonstrates that ρ p is the parameter with the 350 lowest influence, but it is not null as indicated by the results of our 2D sensitivity analyses ( Fig. 5a and d). 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 tend to decrease.

Drew point
The geological settings of the Drew Point area are described in Sect. 2.2 and a summary of the acquisition parameters is 355 provided in Table 1. We invert the 1830 apparent resistivity measurements recorded along an 854 m long profile (Fig. 1g) considering a layer-based model parameterization consisting of two interfaces and five nodes for each interface (see Sect. 4.1).
In total, this results in 36 model parameters, which describe the structure of the interfaces and the layer resistivities. In the PSO, we use 30 particles and a maximum of 400 iterations as stopping criterion. To obtain a single inverted model, we have to evaluate the forward response 12, 000 times, which takes on average 57 hours 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 obtaining an ensemble M F 0 consisting of 416 models.

Ensemble analysis
After the inversion, we interpolate all models to a refined structured mesh before performing any posterior statistical analy-365 ses (Sect. 4.3). In Fig. 7a and b, we present the µ 1/2 (M F 0 ) and IQR(M F 0 ) models calculated from the Drew Point model ensemble. The irregular variations in the IQR(M F 0 ) 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.

370
In the next step, we performed cluster analysis (Sect. 4.3) and found that our ensemble M F 0 can be divided into three model families (M F 1 , M F 2 , and M F 3 ). In Fig. 7b-d and f-h, we present the µ 1/2 (M F i ) and IQR(M F i ) models (where i = 1, 2, 3).
Comparing these models illustrates that M F 1 and M F 2 present a similar IBPT shape dipping toward the open sea. However, for M F 3 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 F 1 and M F 2 models in more detail, we note that the IBPT position in M F 1 is shallower than 375 in M F 2 . 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 profile direction (for M F 1 and M F 2 ), 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 100 < x < 500 m and the IBPT shows a gentle convex slope whereas in the third part (at 500 < x < 700 m) the IBPT is almost flat. Finally, the fourth part is found at x > 750 380 m, where the IBPT may be located at depths ≥ 20 m.
We assess the fit performance for the residuals associated to the ensemble containing all models δ F 0 , as well as for the three clustered model families δ F 1 , δ F 2 , and δ F 3 (Fig. 8). We calculate µ 1/2 (δ F i ) and IQR(δ F i ) (where i = 0, 1, 2, 3) in a pixel-wise fashion and present them as pseudosections in Fig. 8a-h. When comparing these pseudosections to each other, we 385 notice that µ 1/2 (δ F i ) indicate similar fits of the data in terms of amplitudes and pseudosection patterns (although with slightly higher values for δ F 3 ). When comparing the IQR(δ F i ) plots, we note that IQR(δ F 0 ) 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. 1f). This illustrates the applicability of exploring such residual 390 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(δ F i ), IQR(δ F i ), and q 90 (δ F i ) (where i = 0, 1, 2, 3) are presented as histograms in for the clustered families ( Fig. 8j-l), however, small tails to the right are also evident for δ F 1 and δ F 2 . One may tend to reject 395 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 not expect a significant impact from these models on our results and conclusions.

400
Following Sect. 4.4 and Sect. 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  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 405 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 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 [ρ w , ρ uf ] and [z pt , ρ w ] while the most significant positive correlation is found for [z pt , ρ uf ] (note that the absolute correlations of these model parameter pairs are larger in Fig. 9d). 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. 415

Sensitivity analysis
For 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 inver-420 sion 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 ∼ −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. 425 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 Ωm and 10, 000 Ωm (eight values in total), and the IBPT is located at three different depths; i.e., z pt = 16 m, z pt = 10 m, 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 w = [0. ρ p = [5, 20, 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 .

435
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   using layer-based strategies. Whether we have a smooth or sharp transition between frozen and unfrozen 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). In this study, we considered a layer-based model parameterization to target (during ERT inversion) the interface defined by such resistivity contrast (interpreted here as the IBPT) including estimates of uncertainties.

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 450 strategy is that we do not assume an arbitrary resistivity threshold or gradient to interpret the IBPT position as needed for the interpretation of smooth inversion results. 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 455 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 with direct sampling. 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. 2 and 7 ) associated, for example, with submerged 460 thermokarst structures or changes in the ratio of coastal erosion vs. degradation rate.
To find reliable and stable model solutions, we performed several experiments where we ran the PSO several times to find the appropriate parameter settings for the inversion. We noticed that constraining the seabed interface using the water depth data as collected with an echo sounder (allowing for depth variations in the range ±0.15 m, which is the approximate error 465 level of the echo sounder for water depths < 5 m) and the resistivity data as recorded with a CTD device (allowing variation between the minimum and maximum values of the measured data) significantly improved the inversion results and reduced the non-uniqueness of the inverse problem (as also demonstrated by the results of 1D synthetic inversion experiments). Although in this study we only analyzed the results considering a three-layer model parameterization (water and unfrozen sediments overlying an ice-bearing permafrost layer), in the experimental phase, we also explored the possibility of using more than three 470 layers. 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 work-475 flow 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, either because of 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 480 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 global inversion results (considering homogeneous resistivity layers) as reference models to perform local inversions relying 485 on a grid-based model parameterization.
The error level of ERT data is usually unknown; especially, for marine data, where repeated or reciprocal measurements are not possible 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 490 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 RM-

495
SLE 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 smaller 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. 500

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 505 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, which is essential, for example, for a reliable estimation of the IBPT position and for establishing petrophysical relations (e.g., to estimate porewater salinity and ice content). Additionally, we noticed that the 1D inversion results for the 510 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 three 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 515 that, as we showed in our examples, can be helpful to identify interactions between the models 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.

520
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 Table 2 and Table 3) assumed identical resistivities of the ice-bearing permafrost layer (ρ p = 4, 000 Ω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 field measurements at our field sites. Although the resistivities of the unfrozen and frozen layers were similar in both models, we 525 noticed that for model scenarios derived from the Bykovsky site, the inverted ρ p values were overestimated because q 25 and q 75 models indicate ρ 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 resistivity of the sediments 530 increases with ice-content (e.g., Pearson et al., 1986;Fortier et al., 1994;Kang and Lee, 2015), such results may lead us to conclude that the ice-bearing permafrost layer holds higher ice content at Bykovsky compared to Drew Point (assuming the sediments have the same temperature and porewater salinity). Over or under-estimating the resistivity may lead to potentially erroneous interpretations, for example, related to the sediments ice content, temperature, and composition. We would need complementary field information or further analyses like sensitivity assessments to avoid misleading interpretations.

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 540 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, we noticed that if the seawater were less resistive than the unfrozen sediments (e.g., as in the Drew Point case study, Fig. 10), we were 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 545 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 frozen 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 550 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 (which were the parameters with the largest total sensitivity 555 in both examples) were also the parameters that formed model parameter pairs showing the largest correlation (see Fig. 4d and Fig. 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 pre-defined 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. 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 575 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 to interpreting a faster IBPT degradation rate at Drew Point compared to Bykovsky in two ways: 580 1) a layer of submerged Drew Point sediments was already unfrozen upon inundation (e.g., M F 2 in Fig. 7); 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 even larger cryopeg spreading and IBPT degradation.

585
As shown in Fig. 1a and e, 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 U-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 590 thermokarst and undisturbed permafrost. These authors also suggested that drained lake basins, which have undergone thawrefreeze 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 per year) than in Bykovsky (< 1 m per year) that expose coastal areas to inundation in a shorter time. Because of the longer inundation time at Bykovsky, we expect fluctuations 595 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 deepens from ∼ 12 m 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 600 this IBPT dip.

Conclusions
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 605 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 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 610 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.