The article considers the results of mathematical modeling of transient electromagnetic cross-borehole monitoring data for civil and industrial cryolithozone facilities containing thaw zones (taliks) in their vicinity. A solution to the direct problem is presented based on the Sumudu integral transform and the vector finite-element method for two types of borehole sources: induction coils and a less common electric current line, taking into account the frequency dispersion of the electrical conductivity of permafrost. Three-dimensional numerical modeling of the transient signals is performed in realistic geoelectric models of a gas-producing borehole and a residential building on piles. Based on the modeling results, we have revealed the main features and differences of transient cross-borehole monitoring with coils or a line as the source.
INTRODUCTION
In world practice, increasing attention is being paid to climate change issues and permafrost degradation in particular, with the development of unified approaches [Boike et al., 2022]. Concerning Russia, a significant part of the country’s territory is located in the area of the most intense climate change – in the Arctic latitudes [Mokhov and Parfenova, 2021]. Warming in the Arctic is occurring several times faster than the rate of increase in global surface temperature, owing to the rapid reduction in the sea-ice area of the Arctic Ocean [Parfenova et al., 2022].
Subsequent to the thawing of permafrost, northern cities, towns, settlements, infrastructure objects, and oil-and-gas production facilities are being destroyed [Nefedova and Solovyev, 2020; Brushkov et al., 2023]. When permafrost thaws, large amounts of greenhouse gases rise into the atmosphere [Xu et al., 2023], which aggravates the situation with global warming [Mokhov and Smirnov, 2022]. An immediate response of tundra vegetation to climate change in the Arctic is noted [Heijmans et al., 2022]. In addition, during the development of the Arctic territories, huge amounts of hydrocarbon fuel and other toxic elements penetrate permafrost, causing additional harm to the vulnerable environment [Lupachev et al., 2022]. However, there are also a few advantages associated with global warming: for instance, a significant increase in the duration of the navigation period on the Northern Sea Route [Parfenova et al., 2022].
At the end of 2023, the Climate Doctrine of the Russian Federation was approved. Among its most important points are the necessity to develop a state monitoring network and the urgency of providing scientific information for proactive adaptation to climate change and timely decision-making.
Both in Russia and abroad, the most common approach to monitoring objects of the Earth’s cryosphere is thermometry, both in surface and borehole configurations [Sysolyatin et al., 2022; Frolov et al., 2023]. Many practically important results for permafrost monitoring are obtained from ground-penetrating radar (GPR) survey data [Westermann et al., 2010; Edemsky et al., 2024] as well as from electrical exploration with direct (electrical tomography) and alternating (frequency sounding) current [Pavoni et al., 2021]. The integration of thermometry and electrical tomography data is successfully employed [You et al., 2016; Tourei et al., 2024], including cross-borehole exploration [Phillips et al., 2023], with the elaboration of joint inversion algorithms [Tomaškovičová and Ingeman-Nielsen, 2024]. Similarly, GPR measurements are coupled with temperature [Ding et al., 2023] and electrical resistivity tomography [You et al., 2017] data.
There are works on the application of seismic techniques: active and passive seismic measurements [Albaric et al., 2021]. Changes in the active layer of permafrost are studied through the seismoacoustic emission method [James et al., 2019; Köhler and Weidle, 2019]; the depth of the top of Quaternary frozen deposits is specified [Neradovskii, 2021]. To enhance the reliability of the obtained results, electrical resistivity data are involved along with seismic data [Wu et al., 2017; Tourei et al., 2024], with an assessment of their sensitivity to the content of water and ice [Mewes et al., 2017].
The direction of establishing interconnections between permafrost distribution and terrestrial landscapes is evolving [Zotova, 2021]. Stresses and deformations are measured [Wen et al., 2016], and relationships between geomechanical properties and the temperature of frozen rocks are determined [Wang et al., 2023].
There are increasingly more investigations on the use of unmanned aerial vehicles for monitoring the cryolithozone [Ding et al., 2023], including the modification of thermal photogrammetry [Ponti et al., 2024]. Remote sensing of the Earth is actively applied to mapping the stability of permafrost [Zhao et al., 2023].
In addition to purely fundamental research, the properties of permafrost under the foundations of civil and industrial facilities, buildings, and structures are examined: power transmission line supports [You et al., 2016], bridges and railways [You et al., 2017], road embankments [Tian et al., 2021], thawing zones around production boreholes [Filimonov and Vaganova, 2021], buildings on piles [Shang et al., 2022], oil-and-gas trunk pipelines [Cao et al., 2023], and gas-field sites with their geocryological zoning [Brushkov et al., 2023].
Recently, the direction of transient electromagnetic (TEM) sounding of the cryolithozone has been developing [Lorentzen et al., 2024], including its cross-borehole version [Glinskikh et al., 2021, 2024]. For an induction coil as the source, algorithms for mathematical modeling of TEM signals have been devised based on the Sumudu integral transform [Epov et al., 2023] and artificial neural networks [Epov et al., 2024a] as well as rapid modeling tools via a numerical–analytical solution [Nikitenko et al., 2023]. Also, algorithms have been created for transforming recorded transient signals into values of apparent electrical resistivity [Nikitenko et al., 2024]. Field experiments with a prototype of the equipment have been designed [Bukhtiyarov and Glinskikh, 2022; Glinskikh et al., 2023]. We have performed numerical modeling of TEM cross-borehole exploration signals for the monitoring of permafrost rocks under the foundations of industrial facilities [Mikhaylov et al., 2023]. Moreover, approaches to the TEM monitoring data inversion have been elaborated [Nechaev et al., 2024; Epov et al., 2024b].
The presented research is aimed at developing the theory of modeling and inversion of the TEM signals when considering an electric current line as the source. This source has so far been addressed only in individual theoretical investigations, for example [Karinsky, 2004; Karinsky and Daev, 2011; Karinskiy et al., 2023].
A higher signal level and a simpler engineering solution compared to an induction coil should be expected. Transient electromagnetic cross-borehole exploration is regarded further for the two different types of sources: induction coils and a current line.
MATHEMATICAL MODEL OF ELECTROMAGNETIC MONITORING
The mathematical model of electromagnetic monitoring is described through the following boundary-value problem:
where E(t) is the electric field strength vector; J0 is the current density vector in the source; σ(t) is the frequency-dispersive specific electrical conductivity; * is the convolution operation; ε0 is the dielectric permittivity, and μ0 is the magnetic permeability; ∂Ω is the boundary of the computational domain Ω, so far from the electric field source that the electric field strength on it can be neglected. Then, we assume the frequency dispersion of the electrical conductivity of permafrost in the frequency domain (Fourier transform) to be described by the Cole–Cole formula [Olhoeft, 1979; Lee, 1981]:
where σ0 is the direct-current electrical conductivity; m is the polarizability; c = 1 is the exponent; τ is the relaxation time. In environments other than permafrost, the conductivity does not depend on frequency (or time).
Let us apply the Sumudu integral transform to the boundary-value problem:
The Sumudu transform was proposed in [Watugala, 1993] as an alternative to the Laplace transform. Its important property is preservation of scale and dimension of a function [Belgacem, 2006; Belgacem and Karaballi, 2006]. The Sumudu image of a real function is a real function. In this regard, unlike the cases with the Laplace or Fourier transform, there is no need to go to the complex plane. This reduces the resource intensity of the computational algorithm when finding the Sumudu image. Unlike the Laplace transform, there is no explicit formula for performing the corresponding inverse transformation for Sumudu; it requires solving the Fredholm integral equation of the first kind [Epov et al., 2023, 2024a].
By means of the Sumudu transform over time, we convert the mathematical model of electromagnetic monitoring to the following form:
The Sumudu image of the Cole–Cole formula is represented as follows:
To solve the boundary-value problem in partial derivatives with respect to the Sumudu image of the electric field vector, we utilize the vector finite-element method.
VECTOR VARIATIONAL FORMULATION
Let Ω be a three-dimensional domain, possibly inhomogeneous in physical properties, with a Lipschitz-continuous boundary ∂Ω. We introduce functional spaces [Nedelec, 1980; Nédélec, 1986]:
with the norm
The scalar product is defined as
For problem (1), (2) a variational formulation is expressed [Hiptmair, 2002]: For J0 ∈ L2 (Ω) find (curl; Ω), such that for ∀v ∈ H0 (curl; Ω)
The introduced space has the embedding property:
Owing to this property, the electric field satisfying the vector variational formulation also satisfies the charge conservation law in weak sense [Nechaev et al., 2008]. Thereby, when solving the monitoring problem, the source type – coil or current line – is taken into account naturally and does not require special procedures. Only the different source geometry is preset.
Conventionally, the basis functions that form the discrete subspace Hh (curl; Ω) are divided into two types. The first type was introduced in [Nedelec, 1980], while the second in [Nédélec, 1986]. According to the embedding property, part of the space Hh (curl; Ω) consists of functions that are gradients of scalar functions. If scalar functions of order no greater than p are used to define a basis of order p, then the resulting basis is a basis of the first type; if scalar functions have order no greater than p + 1, this is a basis of the second type. Thus, the secondtype basis can be obtained from the first-type basis by supplementing it with basis functions that are gradients of scalar functions of order p + 1.
To construct a discrete analog of the variational problem, the elements of the space H (curl; Ω) are approximated by the elements of the discrete subspace Hh (curl; Ω). Second-type third-order vector elements on a tetrahedral grid are taken as the basis functions [Webb, 1999]. To improve the spectral properties of the matrices acquired after discretization of the original problem, the basis functions can be orthogonalized; complete orthogonalization would lead to a sharp increase in the number of nonzero elements of the matrix. In [Webb, 1999], it was proposed to perform partial orthogonalization – to divide the basis functions into many groups and then perform orthogonalization within each group. High-order vector basis functions can be associated with an edge, a face, or the tetrahedron itself. This depends on how the degree of freedom of a particular basis function is defined: by an integral along the edge, over the face, or over the entire geometric element, respectively. Since one edge, face, or element for high-order bases are associated with several functions, this property is employed as a separator into groups. Defining orthogonalization groups in this way does not add to the number of nonzero elements of the matrix, nor does it change its portrait. In [Webb, 1999], the standard scalar product is used for orthogonalization. In the present work, the basis functions within one group are orthogonalized with respect to a bilinear form that is utilized to construct the variational formulation.
We introduce a discrete analog of the variational problem:
For J0 ∈ L2 (Ω) find , such that :
By expanding the electric field strength vectors Eh over the basis of the discrete subspace and choosing the basis functions as the function vh, a transition is made from the discrete variational problem to an equivalent system of linear algebraic equations:
where X is the vector of weights in the expansion over the basis Eh and the elements of matrices A and B and of the right-hand side vector f are determined by the relations
where are the basis functions of the space Hh (curl; Ω). A modified multiplicative algorithm is applied to solve the resulting system of linear equations [Glinskikh et al., 2021].
REALISTIC GEOELECTRIC MODELS OF CRYOLITHOZONE OBJECTS IN RUSSIA
Following the analysis of specialized scientific publications, we have selected typical geocryologic parameters of permafrost rocks for a number of regions of Russia – with civil and industrial facilities, buildings, and structures that are exposed to the processes of permafrost thawing/freezing.
This study examines sections of the Yamalo-Nenets Autonomous District – on the Yamal Peninsula and in the city of Salekhard [Misurkeeva et al., 2021; Buddo et al., 2022]. Strong lateral variability of the permafrost section is noted, including its lower boundary. The permafrost layer thickness varies from about 70 to 180 m; the resistivity of the frozen rocks ranges from 1500 Ohm·m in the upper part to about 50 Ohm·m in the lower part. The seasonally thawed layer has a thickness of about 2 m. The model may contain local conductive taliks with a typical resistivity of the first dozens of Ohm·m.
At the next step, to get a detailed understanding of the formation and structure features of talik zones in the vicinity of civil and industrial facilities, buildings, and structures in the cryolithozone of Russia, a representative number of open publications have been analyzed. The studied objects include buildings [Zaplavnova et al., 2022; Kosyakina et al., 2023] and structures [Ermakov et al., 2021] on pile foundations, oil-and-gas fields [Sergeev et al., 2015; Cherepanov, 2018; Marinenko et al., 2019], large fuel and energy facilities [Vasiliev et al., 2021], and multiscale structures in permafrost [Hjort et al., 2022].
In reliance on the publication analysis, as well as on materials of field electrical exploration studies, realistic geoelectric models have been formulated. They comprise elements of buildings and structures and thawing zones with talik formation. Polarization parameters in frozen rocks are taken into account.
In what follows, we consider two models:
A production borehole at a gas field on the Yamal Peninsula.
A five-story panel residential building on piles in the city of Salekhard.
Geoelectric model of a production borehole at a gas field on the Yamal Peninsula. The geoelectric model of a production borehole at a gas field on the Yamal Peninsula (Fig. 1) includes permafrost with a resistivity of ρ0 = 1500 Ohm·m (conductivity σ0 = 1/1500 S/m) and air with a resistivity of 106 Ohm·m (conductivity = 1/106 S/m) above the Earth’s surface. The production borehole contains natural gas with a resistivity of 200 Ohm·m (conductivity = 1/200 S/m).
The outer diameter of the first (outer) casing string is 0.426 m (to simplify the model, several actually used casing strings with cement stone between them are combined into one, its wall thickness being 0.138 m). The height of the column from the Earth’s surface is 50 m. The outer diameter of the second (inner – production) casing string of infinite length is 0.15 m; the wall thickness is 0.01 m. The resistivity of the casing string metal equals 0.01 Ohm·m (conductivity 100 S/m, chosen as a compromise in terms of calculation time and the obtained modeling results). At a depth of 5 m, owing to the thawing of permafrost during the extraction of warm natural gas, a talik with 20 Ohm·m resistivity (conductivity 1/20 S/m) forms around the upper casing string. Next, three situations are analyzed: a small (1×5 m), medium (2×10 m), and large (3×15 m) talik. The Earth’s surface is shown as a green dotted line (Fig. 1). All the objects in the model are characterized by axial symmetry. The following parameter values in the Cole–Cole formula are typical of frozen rocks of the Yamal Peninsula: m = 0.7; c = 1; τ = 2.5·10–4 s.
Below are the results of numerical experiments on calculating TEM monitoring signals for a sequential assessment of talik influence:
– in the reference permafrost environment (1/1500 S/m) with two casing strings (100 S/m), natural gas in the borehole (1/200 S/m), air (1/106 S/m), without a talik;
– with the addition of the small talik (20 Ohm·m, 1 × 5 m);
– with the medium talik (20 Ohm·m, 2×10 m);
– with the large talik (20 Ohm·m, 3×15 m).
For signal excitation, two types of sources are used separately – induction coils and an electric current line. The pulse type is rectangular “switch-off.”
The following measuring configuration for crossborehole exploration is utilized. The distance between two monitoring boreholes is 30 m, and that between the monitoring and gas-producing borehole is 15 m. The monitoring borehole depth is 20 m (Fig. 1).
In the case of induction coils, the transmitters and receivers are located in different monitoring boreholes at the same depths from 1 to 20 m, with a step of 0.5 m. The magnetic moment of the transmitter coil is 1 A·m2; the z-component of the magnetic field strength Hz in the center of the measuring coil is analyzed.
For a current line, the cross-borehole configuration is as follows. The line is located vertically in the first monitoring borehole at depths from 1 to 20 m. In the measuring borehole, at depths of 1–20 m every 0.5 m, the z-component of the electric field strength Ez on the borehole axis is under analysis. The current in the line is 1 A.
In both cases, the imaginary line connecting the monitoring boreholes passes through the center of the gas borehole (Fig. 1).
Geoelectric model of a five-story panel residential building on piles in Salekhard. The geoelectric model contains a five-story panel residential building on piles, permafrost with resistivity ρ0 = 1000 Ohm·m (conductivity σ0 = 1/1000 S/m), and air with a resistivity of 106 Ohm·m (conductivity = 1/106 S/m) above the Earth’s surface (Fig. 2). The model includes a 2-m-thick seasonally thawed layer with a resistivity of 50 Ohm·m (conductivity = 1/50 S/m). The five-story panel building is 15 m high, 11.5 m wide, and 34.2 m long and has a resistivity of 100 Ohm·m (conductivity = 1/100 S/m). The building stands on 39 reinforced-concrete piles 9 m long (three rows of 13 piles), 1 m of which is above the ground, creating a vented underfloor space 1 m high with a resistivity of 106 Ohm·m (conductivity = 1/106 S/m). The distances between the piles along the long side of the building are 2.6 or 3.2 m and 5.75 m along the short side. The cross section of the reinforced-concrete pile is a square with a side of 0.4 m and 100 Ohm·m resistivity (conductivity = 1/100 S/m). A cubic talik with a resistivity of 10 Ohm·m (conductivity 1/10 S/m) forms close to the corner pile, flush with its end. Taliks of three sizes are considered: a small talik with a side of 1 m, a medium one with a side of 2 m, and a large one with a side of 3 m. The following parameter values in the Cole–Cole formula are typical of permafrost rocks in the city of Salekhard: m = 0.6; c = 1; τ = 2.5·10–4 s.
For the rectangular “switch-off” pulse, TEM signals are simulated for a sequential assessment of the talik effect:
– in the reference permafrost environment (1/1000 S/m) with the seasonally thawed layer (1/50 S/m), air (1/106 S/m), vented underfloor space (1/106 S/m), five-story panel building, and reinforced-concrete piles (1/100 S/m);
– with the addition of the small cubic talik (10 Ohm·m, 1×1×1 m);
– with the medium talik (10 Ohm·m, 2 × 2 × 2 m);
– with the large talik (10 Ohm·m, 3 × 3 × 3 m).
The distance between the monitoring boreholes is 13.5 m; the depth of the boreholes is 10 m. The distance from the monitoring borehole to the edge of the building is 1 m. The induction coils are positioned at the same depths in two boreholes from 1 to 10 m with a step of 0.5 m. The electric current line is located vertically in the first monitoring borehole at depths from 1 to 10 m. Measurements in the second borehole are made at 1–10 m depths, every 0.5 m. The imaginary line connecting the monitoring boreholes passes through the lateral edge of the short side of the building (Fig. 2c).
NUMERICAL SIMULATION OF SIGNALS FROM INDUCTION COILS AND CURRENT LINE
Production borehole at a gas field on the Yamal Peninsula. The absolute values of the vertical component of the magnetic field |Hz| are scrutinized during the TEM monitoring with induction coils (Fig. 3). The key difficulty of monitoring in this case is to identify talik zones against the background of highly conductive metal of the casing strings.
The range of recording times is from 5·10–8 to 5·10–6 s; the main range of |Hz| is 10–7–10–6 A/m. One may note a differentiation of the signals depending on time (in the horizontal direction). In the diagram without a talik (Fig. 3a), a transition of Hz through zero is observed in the vicinity of a time of 3·10–7 s in the lower part of the borehole up to approximately 2·10–7 s in its upper part. With the appearance of a talik and a growth in its size, the area of the transition of Hz through zero shifts toward increasing time: up to 4·10–7 s at the bottomhole and up to 3·10–7 s near the Earth’s surface. In this case, the area of maximum signal values expands to the left of the transition through zero and narrows to the right of it. The characteristic bend of the zero-crossing at depths of up to 5 m corresponds to the area above the talik and is due to the influence of the air. In the range of large times (2·10–6–5·10–6 s), the changes in the monitoring diagrams are negligibly small. In addition, from a comparison of Fig. 3a and Fig. 3b it seems problematic to identify the 1×5 m talik by means of the monitoring diagrams.
The overall pattern of the TEM monitoring changes drastically when we apply an electric current line as the source (Fig. 4). The range of recording times gets wider: 5·10–8–10–4 s; the main range of absolute values of the electric-field vertical component |Ez| is 10–2–10–1 V/m.
The analysis of the monitoring diagrams indicates, first, the signal differentiation to become more pronounced in the vertical direction than in the horizontal direction. In the absence of a talik, at early times (up to 4·10–7 s), an oval area of maximum signal values is evidenced. With an increment in the recording time, a concentric area with signal values of up to 4·10–2 V/m may be delineated. With the formation of a talik and growth in its size, an expansion of the area with the lowest signal values (about 10–2 V/m) takes place from top to bottom. In addition, on successive monitoring diagrams there occurs a shift of the oval and concentric areas to the bottomhole zone, with a simultaneous decrease in the signals to 2.5·10–2 V/m at medium and large recording times. It is worth noting that the smallest talik considered – 1×5 m – is distinguished according to the diagrams for the electric line somewhat more clearly than according to those for the induction coils. For the line, in comparison with the coils, the taliks of a larger size are more pronounced; signal transition through zero is absent. Moreover, the signal of the electric line attenuates over time much more slowly.
Five-story panel residential building on piles in the city of Salekhard. Let us consider similar monitoring diagrams for the residential building in Salekhard: with induction coils (Fig. 5) and with an electric line (Fig. 6). The main monitoring issue consists in the ability to identify a forming talik against the background of the dense pile field.
In the case of induction coils (Fig. 5), the range of recording times is 5·10–8–5·10–7 s, the |Hz| signal varying from 10–6 to 10–5 A/m. The key features of the monitoring diagrams are similar to those in the case of the production borehole: A transition of Hz through zero is noted, but in the vicinity of a shorter time of 6·10–8 s. As a talik forms and increases in size, the transition shifts to a time of 8·10–8 s. In general, as the talik grows, the region with the largest |Hz| values to the left of the Hz zero-crossing expands and to the right of it narrows, especially at medium depths and recording times. In addition, at times of about 4·10–7–5·10–7 s, the analyzed signal drops slightly.
When an electric line is used as the source (Fig. 6), the range of recording times is 5·10–8–10–4 s, the signal |Ez| changes from 10–2 to 10–1 V/m. A characteristic feature of all the monitoring diagrams is a lenticular shape of the anomaly in the whole time range, at depths of approximately 2 to 8 m, which is caused by a significant three-dimensional heterogeneity of the geoelectric model. In the absence of a talik, the highest signal values of about 10–1 V/m are observed at depths of 5–8 m. As the talik size increases, this area becomes progressively less pronounced. At the same time, in the lower part of the borehole (8–10 m), the signal noticeably decreases to 3–4·10–2 V/m, and in the upper part there is an even greater decrease in the signal to 10–2 V/m. In contrast to the horizontal variability of |Hz| for the induction coils, the expansion of the talik is generally characterized by the vertical differentiation of |Ez|.
Thus, the results of the three-dimensional numerical modeling show, first, that the formation of talik zones manifests itself on the diagrams of TEM signals in an essentially different way when induction coils or an electric line are used as the source. Second, the signals from a line attenuate over time several orders of magnitude more slowly, which might be a key factor at significant distances between monitoring boreholes. Third, despite the fact that a current line is an integral source compared to a set of induction coils, this source achieves no lower resolution in relation to the target objects of TEM monitoring of the cryolithozone in its cross-borehole configuration.
CONCLUSIONS
We have devised algorithmic and software tools for numerical modeling of transient cross-borehole signals with the excitation of an electromagnetic field by induction coils or a current line in a borehole. The algorithms are based on a combination of the Sumudu integral transform with the vector finite-element method. Drawing on the analysis of both open publications and materials of field electrical exploration studies, two realistic geoelectric models of the Russian cryolithozone have been created – a production borehole at a gas field on the Yamal Peninsula and a five-story panel residential building on piles in the city of Salekhard. Inherently three-dimensional heterogeneous models describe in detail structural elements as well as zones of thawing and talik formation, taking into account polarization parameters in permafrost. We have carried out numerical modeling of crossborehole electromagnetic monitoring signals in the specified geoelectric models. The diagrams of talik zone monitoring with induction coils and an electric line are shown to differ significantly. More specifically, coils are characterized by signal transition through zero, with the main changes occurring along the recording time scale. A line does not have a transition through zero, and the greatest dynamics are traced along the borehole depth. Transient electromagnetic cross-borehole monitoring of civil and industrial objects of the cryolithozone is shown to be possible both with the use of induction coils in two boreholes and with an electric line as the source. The appearance of a talik and change in its size can be visually traced by time-lapse transient signal diagrams plotted as a function of registration time and depth. The application of a current line may be preferable when we identify small taliks against the background of highly conductive monitoring objects (for instance, the casing of a gas-producing borehole). The signal from a line attenuates much more slowly over time, which gives an advantage in electromagnetic monitoring with widely spaced boreholes. In addition, a simpler engineering and more economically advantageous implementation of a borehole current line is expected as compared with a large set of induction coils.
ACKNOWLEDGMENTS AND FUNDING
The study was supported by grant No. 22-17-00181 (“Pulsed electromagnetic sounding of permafrost: theoretical and experimental development of a high-resolution geophysical method, scientific substantiation and creation of an innovative technology for monitoring of the cryolithozone”) from the Russian Science Foundation (https://rscf.ru/en/project/22-17-00181/).