1. Introduction
Significant reserves of methane (CH4) are held in the Arctic seabed, but the release of CH4 to the overlying ocean and, subsequently, to the atmosphere has been believed to be restricted by impermeable subsea permafrost, which has sealed the upper sediment layers for thousands of years [
1]. Subsea permafrost and shallow hydrate deposits are predicted to occupy almost 60% of the East Siberian Arctic Shelf (ESAS), which is the broadest and shallowest shelf of the world ocean [
2,
3,
4]. Until very recently, understanding of the current state and stability of the subsea permafrost system in the ESAS was primarily based on modeling results showing that degradation of subsea permafrost and destabilization of shelf Arctic hydrates was hypothesized to be negligible on a decadal century scale [
5]. Moreover, the northern thermokarst lakes were proposed as the major regional source of atmospheric CH4 [
6]. However, recent observational data collected from the ESAS demonstrated that this vast yet shallow region is a key area of the Arctic Ocean for atmospheric venting contributing annually two-three times more than the entire world ocean, and no less than terrestrial Arctic ecosystems, and [
7,
8], but unlike terrestrial ecosystems, the ESAS emits CH4 year-round due to its partial open during the winter when terrestrial ecosystems are dormant [
3,
9]. Observational data became a subject for theoretical and experimental work targeted at understanding the ESAS permafrost stability and its mapping.
In recent years, a number of surveys have been launched and are continuing in the ESAS region, with a primary focus on permafrost mapping, which is crucial to understand its current state and possible greenhouse effect-associated potential. The subsea permafrost (ice-bonded permafrost, IBP) now has temperatures 8–10 degrees higher compared to ground permafrost, which is close to the thawing conditions, and the latter has been shown to cause a massive release of bubble methane in mega-sips. The studies are mainly based on shallow drilling and logging, seismic and electrical imaging, temperature and geochemical observations, as well as heat transfer modeling [
10,
11].
Due to the significant effect of freezing/thawing transition on the electrical resistivity of fluid-saturated sedimentary rock, electrical and EM methods are particularly beneficial for IBP and gas hydrate imaging, and have found extensive application in permafrost-related studies [
12,
13,
14].
A large part of resistivity-based studies uses the direct current (DC) electrical resistivity tomography (ERT) method to image the seafloor rock structure in shallow water coastal Arctic environments [
15,
16,
17,
18]. These studies provide evidence for substantially inhomogeneous permafrost structure associated with intermittent thawed zones and resistive frozen regions and reveal the IBP top boundary down to 20 m depth. Being efficient in shallow water (1–3 m) coastal environments, ERT has essential limitations in terms of the depth of investigation (DOI) barely exceeding 30–40 m, although constraining particular inversion parameters, (seawater resistivity and thickness) may increase sensitivity to deeper structure and help reduce non-uniqueness when inverting for IBP resistivity and geometry [
19].
A deeper water setting requires the application of the technique providing higher DOI, and controlled source EM [
20,
21,
22,
23,
24] based on measuring transient EM field arising in response to injected square waveform current pulses due to induction (time-domain electromagnetics, or TDEM, sometimes called transient electromagnetics, or TEM), is one of such methods [
25,
26,
27]. Near-field TDEM sounding modification, which is based on the measurement of the late-stage transient EM process probing subsurface conductivity structure even at a zero source–receiver offset, offers the highest level of detail (resolution capability) among all electromagnetic survey techniques, except for ground penetrating radar [
28], being efficient in investigating structural features of the subsurface as well as in estimating reservoir (filtration-related) properties of specific formations. Marine TEM sounding modification is also an efficient sounding technique [
29,
30] that has been widely used in recent years in the studies of the sedimentary strata and permafrost rock in the shallow water of the Arctic shelf [
12,
31]. The basic limitations of this technique are associated with the need for accurately measuring the late-time portion of the transient response curve, where the signal is usually extremely low and biased due to rapid decay. On the other hand, TDEM data have a relatively high sensitivity to low resistivity regions associated with taliks, but are not always capable of providing sufficient depth resolution when determining the position of the upper and lower boundaries of resistive formation. The latter is crucial when the method is applied for IBP mapping, and this study is aimed at the simulation-based analysis of the TDEM capabilities in determining the parameters of the unfrozen sea bottom sediments, the IBP itself, and sub-permafrost structure strata.
This paper summarizes the synthetic response analyses and uncertainties of the time-domain electromagnetic (EM) technique to be used for the subsea permafrost mapping regarding use and estimates from numerical modeling data.
3. Results
Simulated apparent resistivities shown in
Figure 2d–f and
Figure 3d–f reveal the smoothed patterns with certain features of similarity to the true geometry of the resistivity models. However, these pseudo-sections are only the form of the response, where the vertical dimension is time, and these cannot be directly compared to the original resistivity images (panels a–c in
Figure 2 and
Figure 3). Nevertheless, the shape of the bluish region representing the response of the conductive portion of the model (water and unfrozen sediments), correlates with the thickness and conductivity of the upper 3 layers, including the IBP layer. The right-hand part of the pseudo-section images, corresponding to greater water depth and higher conductance of the seafloor sediments, have a larger portion of the early-time high-resistivity region (1–7 ms versus 1–3 ms time interval for a shallow-water setting), which is an indication of the fact that far-field behavior lasts longer in a less resistive environment [
26].
3.1. Synthetic Study of TDEM Response Sensitivity to Permafrost Parameters in ESAS Setting
To understand the effect of specific resistivity or geometry of the IBP formation, as well as the resistivity of unfrozen sediments above and below the IBP layer, we plotted relative response anomalies as time-distance pseudo-sections in
Figure 5. These were calculated from the apparent resistivity time-distance distributions shown in
Figure 2 and
Figure 3d–f as
where
x denotes distance and
t is the time passed since the current cut-off.
Despite the significant prevalence of seawater and bottom sediments’ conductances over those of the underlying IBP formation, all models show substantial anomalies relative to reference model A. Model B in which the IBP layer’s resistivity is ~3 times lower compared to model A, exhibits on average 10% negative anomaly (decreased apparent resistivity), with maximum absolute values seen at small distances (0–300 m), which is due to highest sensitivity to IBP present in this interval associated with the most shallow water layer. Model C, being essentially a no-IBP model, shows intense positive anomaly within an early time range (1–6 ms) followed by a strong negative anomaly in 7 ms−1 s intervals, with the amplitude decaying smoothly (to −5%) in the right-hand part of the plot. Inversely, Model D, having 2-fold increased IBP thickness, shows a weak −5% negative anomaly at early times, and then a clear 10–15% positive anomaly within the 10 ms−1 s time range, which gets slightly lower (around 5%) at distances greater than 1400 m. Finally, model E, being a perturbated version of model A, with resistivities varying smoothly in a lateral direction in all 3 upper layers, shows a mosaic pattern of the time-distance anomaly distribution, which is consistent with alternating series of increased and decreased resistivity regions. In terms of time distribution, the response shows coupled positive and negative anomalies which is the indication of change in the slope of the early-time apparent resistivity curve segment relative to the reference (model A).
Individual responses extracted at 6 specific locations are shown in
Figure 6. It can be seen that in a shallow-water setting, all 5 models differ significantly in terms of response curve shape in the time range of 10 to 100 ms (
Figure 6a–c). However, those differences become much smaller in the deep-water environment (
Figure 6d–f). Despite the upper layer being a conductor, the behavior displayed by the responses clearly indicates the prevalence of the far-field descending segment within the initial time range between 1 ms and 3–10 ms, depending on the local model, i.e., the conductance of the upper formation. Response minima, associated with the conductive layer(s), are typically observed at around 6–10 ms (x = 0 m), 10–20 ms (x = 350; 875 m) and 40–60 ms (x = 1300; 1575; 1975 m), which correlates with the increase in the formations’ thickness towards the greater distances. The minimum value varies between around 10 and 1.5 Ohm·m, consistently getting lower towards the right-hand part of the profile with a conductance increase. Locations x = 0, 350 and 875 m naturally exhibit significant differences between the responses corresponding to Models A–E (location x = 0 shows the largest anomalies, of about 30% within the intermediate time range), which is indicative of high sensitivity to subbottom structure in a shallow-water setting. Oppositely, a deeper-water setting (x = 1300, 1575 and 1975 m) implies much smaller misfits (not exceeding a few percent) between various curves, except for Model E at x = 1575 and 1975 m, which clearly deviates from other responses due to substantial differences in conductivity of water, seafloor sediments and IBP formation.
3.2. Inversion of the Synthetic Responses
In order to understand the practical capabilities and limitations of the galvanic excitation-based TDEM technology in the marine environments for IBP mapping, we applied 1D inversion code to synthetic apparent resistivity datasets calculated for models A and E. The resulting resistivity images in comparison with true models along with misfit pseudo-sections are presented in
Figure 7 and
Figure 8.
Both inverted models show reasonable consistency with the original (true) models, and major resistivity features have been successfully recovered. However, some model elements show erroneous configurations, which is mainly attributed to the underestimated thickness and resistivity of layer 2 (relatively conductive unfrozen seafloor sediments), causing the top of IBP to be mapped with depth error from a few meters in a shallow-water setting to 40 m in a deep-water environment, causing the reduction in sensitivity to the underlying structure. In
Figure 9, we evaluate the quality of the inverted models in terms of how they fit the true conductances calculated for individual layers 2 and 3, as well as the total conductance of layers 1 and 2. From the graphs, it can be easily seen that seawater’s (layer 1) conductance (20 to 200 S) is substantially higher than that of the seafloor unfrozen sediments (layer 2, 5–60 S) and largely prevails (by) over conductance of the IBP formation (layer 3, 0.2–2 S). Despite showing substantial errors in layer 2’s thickness and resistivity, the results for both models A and B exhibit a precise fit in terms of combined conductivity of seawater and unfrozen sea bottom sediments (green line indicates the true distribution, and the dotted one corresponds to inversion), which is in agreement with conductance equivalence known for inductive EM methods [
26].
The most prominent feature of the inverted models is accurate (to within a few percent) recovery of summated conductance of the seawater, layer 1, and unfrozen bottom sediments, layer 2 (green dotted and solid lines in
Figure 9 show good fit quality), while the recovered conductance of the layer 2 alone is slightly less accurate (blue line and dots) which is natural due to much lower conductance compared to seawater. Finally, target layer 3 shows the worst conductance fit quality, with errors reaching 5–10 fold of the true value, which generally increases from left to right as the conductance of the overlying formations increases from nearly 20 to over 200 S.
Obviously, this behavior is in agreement with the conductance equivalence principle stating that the recorded decay response is effectively controlled by the total conductance of the upper conductive strata, rather than by individual resistivities and thicknesses of its layers, assuming near-field approximation [
26]. Thus, even if the seawater resistivity and thickness are constrained by directly-measured salinity and bathymetry data, given the non-zero bias present in the response, the inversion targeted at mapping the IBP’s top boundary (i.e., the thickness of relatively conductive subbottom sediments) remains non-unique within the conductance-equivalent bounds. Nevertheless, it enables estimating its conductance with reasonable accuracy, which can be used further to constrain thickness if resistivity estimates are available from drilling or seafloor sampling. Conductance can also be used in coupled resistivity—heat transfer modeling through additional constraints.
In Models A and E, the true conductance of the IBP layer is 20 to 800 times lower than that of the overlying formation (including seawater); however, it still has an effect on calculated responses and is reliably and smoothly imaged by the inversion. In Model E inversion results (
Figure 8b), both talik zones with reduced resistivity located within distance intervals of 700–900 and 1500–2000 m, respectively, have been imaged by the inverted model.
3.3. Estimating Inversion Uncertainty for TDEM Stations Collected in Shallow-Water and Deep-Water Setting
Given the seawater depth range in ESAS varies from 10 to 100 m, and its typical resistivity value is around 0.3 Ohm·m), the conductance estimates fall between 30 and 300 S, which means that in order to be discriminated by the data, any features of interest in the underlying portion of the section must have a conductance of at least 20–30% of the conductance of overlying layers. Obviously, this implies a rapid decrease in subbottom sensitivity associated with an increase in seawater depth, a factor that mainly controls its conductance.
To estimate and illustrate the effect of sensitivity and inversion uncertainty (equivalent solutions) for the subbottom part of the resistivity structure under various noise settings, we have calculated the response dataset including over 100 k transient curves corresponding to a set of 1D 3-layer models with variable resistivities and thicknesses. The model includes the water layer (resistivity 0.25–0.6 Ohm·m, thickness 10–70 m), bottom sediments (resistivity 1–1000 Ohm·m, thickness 10–1000 m) and a deeper semi-infinite layer (resistivity 1–1000 Ohm·m); see
Table 2 for details. All possible combinations of parameter values varying within the above ranges produce a series of models, and a dataset of transient response curves, the latter is shown as a grey-colored cloud in
Figure 10. Among those, we have chosen two particular models/response curves, that imitate resistive permafrost (500 Ohm·m) overlaid by thawed bottom sediments and water, differing in water depth (15/60 m), water resistivity (0.4/0.3 Ohm·m), bottom sediment thickness (15/30 m) and resistivity (2/5 Ohm·m). Response for the shallower water model (15 m) is shown in
Figure 10 as a red dotted curve, while the one for deeper water (60 m) is indicated as blue points. Those two configurations simulate two different settings, characteristic of the most shallow part of the Laptev shelf (Bour Khaya Gulf), with less saline (and more resistive) water due to freshwater inflow from the Lena delta, and outer part of the shelf, with higher depth (ranging from 50 to 100 m and more) further to the north.
We calculated the misfits between each of these responses and each of all the responses forming the test dataset (grey cloud in
Figure 10, numbered by
j) as
where
is either “shallow water” or “deep water” transient response,
is
jth response from a complete simulated dataset, and
denote time samples, counting from 1 to
. Based on misfit values obtained, we categorized all the models into 4 groups, so that corresponding responses fit
curves to within 10, 5, 2 or 1 percent error
(averaged over the curve).
Figure 11 illustrates the margins of each category in model space, shown as red, blue, green and red regions, respectively, for 10%, 5%, 2% and 1% data error for each of two specific responses, “shallow-water” (
Figure 11a) and “deep-water” (
Figure 11b). One obvious conclusion is that a shallow-water setting results in a narrower resistivity range for layers 2 and 3, being of primary interest when resolving the structure in terms of lower resistivity (thawed sediments) or higher resistivity (frozen sediments). In most cases, 1%, 2%, 5% and 10% data errors are all acceptable to enable discriminating between the lower (less than 10 Ohm·m) and higher (more than 10 Ohm·m) resistivity in both layers 2 and 3; 1% data accuracy shows more certain inversion bounds (yellow).
In the case of a deep-water setting, the resistivity bounds for 10% and 5% data error are less definitive (
Figure 11b) compared to shallow water, which indicates a severe lack of sensitivity to deeper layers’ resistivities. Though 2% and 1% margins still exhibit some moderate sensitivity to the layer 3 resistivity, the measured transient time given is large enough. Note that intermediate layer 2, given its relatively small thickness, is resolved with the largest uncertainty. The above analysis shows the limitations of TDEM inversion precision, pointing to the necessity of high data quality, while in practice keeping the data error within a 2% threshold seems to be problematic, specifically in a late transient time range.
4. Discussion
Permafrost studies, unfolding in recent years on the Arctic shelves, require a reasonably accurate mapping tool, capable of imaging the sedimentary structure down to 200–300 m depth. Although TDEM has been employed for decades in a broad set of environments and applications [
25], and its capabilities and limitations are quite well understood [
20], this tool remains somewhat exotic in a marine setting. The marine environment implies a strong effect of the saline low-resistivity water layer and water-saturated seafloor sediments obscuring the resistivity structure beneath. Successful implementation of marine TDEM-based systems in deep-ocean environments has been achieved through the employment of submerged source–receiver arrays towed relatively close to the seabed [
42,
43], or receivers deployed on the seafloor [
34], or remotely operated vehicles carrying loop-based instruments [
30]. Another approach, promoted by Barsukov and Fainberg [
29], is based on vertical electric field measurements, carried out by a dipping receiver towed behind a horizontal dipole source. Locating the transmitter and receiver closer to the seafloor is obviously beneficial in terms of EM response sensitivity, however is costly and performance-limiting, while the present day large-area shallow-water surveys require fast and cheap solutions.
TDEM applications in permafrost studies remain quite limited compared to other fields, such as hydrocarbons and mineral exploration. Yet, since the 1980s there have been attempts to use TDEM for permafrost mapping, along with an investigation of its efficiency. The thesis by G.G. Walker [
44] extensively covers practical ground-based TDEM data collected in Alaska with a loop system. From apparent resistivity responses measured mainly within 10
−4 to 10
−2 s time intervals, by making use of a 1D inversion code, resistivity models have been obtained, showing a resistive permafrost layer with a thickness ranging up to 600 m. A synthetic study [
45] also simulating an onshore permafrost setting, found the ability of TDEM response to image lateral inhomogeneities both in permafrost geometry and resistivity of the conductive sediments beneath IBP.
Under typical onshore conditions, with a large-thickness high-resistivity top layer, TDEM data tend to display rapid decay immediately from the early transient time (100 mcs), making it hard to achieve measurable signal at late times (around 1 s), though remain sensitive to conductive formation beneath the resistive layer. On the opposite, in the case of subsea permafrost, the common setting involves low-resistivity layer(s) above the higher-resistivity target structure.
Despite the availability of some amount of onshore resistivity studies, to our knowledge, examples of marine TDEM permafrost studies are rare. Ice-based, yet marine, TDEM measurements taken in shallow (a few meters) waters near Muostakh Island (the coastal part of the Laptev Sea) [
36] yielded IBP thickness estimates of nearly 500 m; however, these were loop soundings along the relatively short line with significant spatial separations between the stations.
The recent TDEM onshore loop-system survey by Creighton et al. [
46] revealed a nearly 100 m thick low-resistivity (1–2 Ohm·m) talik beneath a thermokarst lake in Alaska, with an overall imaging depth of 200 m and no signs of IBP bottom boundary within this depth range (the IBP is likely much thicker in this area). This resistivity-thickness configuration of the conductive top layer (50 to 100 S in terms of conductance) is roughly equivalent to a marine setting with water depth ranging from 15 to 50 m and water resistivity varying between 0.3 and 0.5 Ohm·m. In the above study TDEM was capable to image resistive (tens Ohm·m) IBP beneath the talik, and, although a different array configuration was used, it is still fairly close to the one we use in our analysis in terms of depth of investigation and sensitivity.
Our dipole–dipole TDEM data from Laptev Sea collected in 2020 in the area with seawater depth varying between 7 and 15 m, and water resistivity between 0.47 and 0.55 Ohm·m, reveal relatively uniform permafrost layer having a thickness of 100–200 m, and distinctive sensitivity to sub-permafrost conductive sediments [
35], thus proving the technology efficacy in a shallow-water setting. Notably, the IBP thickness estimates obtained from inverted TDEM data turn out to be substantially lower compared to existing temperature-based models for this area [
10].
Thus, in this study we focus on a specific setting involving a marine towed dipole–dipole TDEM system, relatively shallow saline water covering subbottom sedimentary formations which may host a laterally non-homogenous permafrost layer, and apply the numerical study to evaluate the method’s sensitivity and inversion uncertainty to specific properties of the submarine IBP layer and overlying unfrozen sediments, assuming 1D approximation.
Synthetic data analysis included the inversion of simulated response data within a limited time range (1–100 ms, usually available with acceptable quality from the survey data) with constraints imposed on seawater parameters (these are known from bathymetry and salinity data) to evaluate the uncertainty of recovered resistivity images under conditions typical for ESAS environments with variable water depth and IBP characteristics.
Inversion of the response dataset calculated for the model having laterally-variable layer resistivity, including low-resistivity talik regions within the IBP formation (Model E), shows reasonable accuracy and proves the method’s ability to discriminate taliks from intact IBP, except for a particular zone (x = 1500–1600 m), where an overlying structure has conductance over 250 S.
Evaluation of model equivalence bounds based on compiling a large response set assuming the 3 layer structure (water, unfrozen sediments, permafrost) and calculating misfits with the true response, suggests the sufficient accuracy of the model recovery, given the data errors are within 1% threshold over the 1–1000 ms time range, allowing for discrimination between thawed (<10 Ohm·m) and frozen (>10 Ohm·m) states (yellow region in
Figure 11) of the basement (permafrost) layer. However, available ESAS survey data usually exhibit substantially higher errors [
47], while a 5% error floor is obviously not enough for such discrimination (blue region in
Figure 11), especially in deep-water settings.
One of the crucial problems is connected with response distortions [
48], mainly presenting as low-frequency bias that leads to a significant shortening of interpretable segments of measured responses, in turn causing the reduction of DOI.
Thus, TDEM acquisition/inversion technique yields large amounts of high-resolution data enabling in situ probing of permafrost thickness and lateral distribution aiming at a better understanding of the extent of the IBP formation in the study area. The peculiarities of transient EM field behavior observed from synthetic data determine the method’s capabilities depending on sea bottom depth, as well as resistivity features of the subbottom structure. Indeed, deeper water results in reduced DOI and sensitivity of TDEM response to the portion of section beneath the seafloor. Provided that transmitter and receiver are placed at the water surface (or slightly submerged), the high conductance of the seawater may play the dominant role, obscuring (or screening) the underlying formations so that the response curve becomes less sensitive to their parameters.
5. Conclusions
We carried out a modeling-based analysis of the towed dipole–dipole TDEM system capabilities with respect to subsea permafrost imaging. In compiling the synthetic resistivity models, we relied on known permafrost-related resistivity patterns [
12,
20,
40] as well as survey data on the bathymetry, water salinity and previous EM studies in the ESAS region [
31,
35,
36]. The study included forward modeling TDEM data (response) analysis for five synthetic models simulating different settings (presence/absence of ice-bonded permafrost layer, different positions of its top and bottom boundaries, different width and thickness of thawed bodies or taliks, variable sea-water depth and its resistivity), inversion of response for particular models and estimation of the inversion uncertainty for shallow-water and deep-water setting.
From the comparison of the modeled responses, it can be seen that the models can be essentially discriminated between each other despite the significant prevalence of the effect of seawater and unfrozen sediments, with all 5 models differing significantly in terms of response curve shape in time range 10 to 100 ms, except for the deep-water portion (water depth over 50 m), where this difference becomes smaller.
Inversion results are found to reproduce the main features of the true models, although the thickness of relatively conductive unfrozen seafloor sediments has some tendency to underestimation. Inverted images show accurate (to within a few percent) recovery of summated conductance of the seawater and unfrozen bottom sediments, while the recovered conductance of the latter alone is slightly less accurate. Finally, the target IBP shows lower conductance fit quality, which generally increases following the increase in conductance of the overlying formations. Nevertheless, talik zones with reduced resistivity have been imaged by the inversion, and IBP conductance still can be estimated with reasonable accuracy. Conductance estimates can be used further to constrain thickness if resistivity estimates are available from drilling or seafloor sampling, and also might be helpful in coupled resistivity—heat transfer modeling through additional constraints.
An uncertainty study involving multiple 1D 3-layer models suggests that the shallow-water (10–20 m) setting results in narrower resistivity ranges for target layers when resolving the structure in terms of lower resistivity (thawed sediments) or higher resistivity (frozen sediments). In most cases discrimination between lower (less than 10 Ohm·m) and higher (more than 10 Ohm·m) resistivity is possible. Deep-water setting (30–70 m) yields a substantial lack of sensitivity to deeper layers’ resistivities. However, the 2% response error margin still exhibits some moderate sensitivity to the lower layer resistivity, given the measured transient time is large enough.
The above analysis shows the limitations of TDEM inversion precision, pointing to the necessity of high data quality, namely, keeping the data error within the 2–3% threshold, especially in deep-water settings.
This study allowed us to understand the possible uncertainties in the geometry and resistivity of the reconstructed permafrost layer, depending on seawater depth and unfrozen layer thickness, and confirm the overall efficacy of TDEM technology for the subsea permafrost imaging, which is crucially important for understanding the current state of the subsea permafrost-hydrate system and possible future dynamics. Furthermore, the obtained modeling results can be used to modify and improve the technology for collecting permafrost mapping data in the ESAS.