Pulsed Electromagnetic Cross-Well Exploration for Monitoring Permafrost and Examining the Processes of Its Geocryological Changes

Abstract

This paper is dedicated to the topical problem of examining permafrost’s state and the processes of its geocryological changes by means of geophysical methods. To monitor the cryolithozone, we proposed and scientifically substantiated a new technique of pulsed electromagnetic cross-well sounding. Based on the vector finite-element method, we created a mathematical model of the cross-well sounding process with a pulsed source in a three-dimensional spatially heterogeneous medium. A high-performance parallel computing algorithm was developed and verified. Through realistic geoelectric models of permafrost with a talik under a highway, constructed following the results of electrotomography field data interpretation, we numerically simulated the pulsed sounding on the computing resources of the Siberian Supercomputer Center of SB RAS. The simulation results suggest the proposed system of pulsed electromagnetic cross-well monitoring to be characterized by a high sensitivity to the presence and dimensions of the talik. The devised approach can be oriented to addressing a wide range of issues related to monitoring permafrost rocks under civil and industrial facilities, buildings, and constructions.

1. Introduction

The study of the state of permafrost, which occupies at least 25% of the world’s land area [1], is now becoming an increasingly urgent scientific problem. This is owing to, first of all, climatic changes in recent decades, accompanied by an increase in the average annual air temperature on the Earth. Thawing permafrost, as a result, poses threats at a local, regional, and global scale.

In the Russian Federation, the territory of permafrost distribution is about 65% [2], which makes its study extremely relevant. The critical need is also determined by the latest events related to the environmental disaster at the oil storage facility in Norilsk on 29 May 2020 (Figure 1). There, a fuel storage tank was destroyed as a consequence of the settlement of its foundation, with about 20 thousand tons of diesel being released into the environment. The settlement is assumed to have occurred due to the permafrost degradation [3,4]. This necessitates the development and improvement of geophysical techniques for scrutinizing permafrost conditions.

Permafrost monitoring is one of the most effective and proven ways to deal with the changing state of permafrost. In this context, to solve the fundamental problems associated with climate change, various approaches have been proposed and applied: temperature borehole monitoring [6,7]; electrotomographic studies [8,9]; shallow seismic surveys [10,11]; ground-penetrating radar exploration [12]; satellite observations [13,14]; integration of methods [15], and others.

Permafrost monitoring has also found a wide application in relation to civil and industrial facilities, buildings, and constructions such as oil pipelines [16,17], oil and gas fields [18,19,20], buildings [21], highways [22,23,24,25], airports [26], railways, and power transmission lines [27].

The aim of the presented investigation was the scientific substantiation of a new technique for electromagnetic monitoring of permafrost under highways, which is a serious problem in the territory of Russia. The operation of highways built on frozen ground is characterized by a high capital intensity attributed to constantly occurring deformations of the roadbed, requiring reconstruction work and measures to stabilize the permafrost-affected environment [28]. As a rule, the cause of such deformations is the ground ice melting and the formation of thaw bowls at the road base, giving rise to thermokarst. The evolution of thermokarst goes on for several years, and sometimes for decades, which leads to significant financial costs to maintain the road in a condition that meets the necessary standards. For instance, 5 years after the start of operating the R-297 ‘Amur’ federal highway (Chita-Khabarovsk), in the Amur Region alone there were more than a hundred of the so-called irreversible deformations (subsidence) of the roadbed, with lengths of 10–40 m [29]. In the Trans-Baikal Territory, an example of the long-term subsidence is the road segment from the Amur highway to the Peschanka settlement on the outskirts of Chita (Figure 1), which has been experiencing deformations since the 1980s (Figure 2a).

In 2015, a large-scale expensive reconstruction of the site was carried out, with the replacement of ground at the road base. However, the misinterpretation of the engineering survey results without regard to the geophysical data led to the implemented measures being ineffective, with the roadbed deformations still taking place (Figure 2b).

Therefore, the research objectives were (1) to create a realistic geoelectric model of thawed permafrost by the example of that under the R-297 ‘Amur’ federal highway; (2) to construct a mathematical model of pulsed electromagnetic cross-well exploration by making use of the modern vector finite-element method, which has proven effective for complex geophysical problems; (3) to develop the respective high-performance parallel computing algorithm with an assessment of its performance; and (4) to simulate and analyze the electromagnetic signals.

2. Geoelectric Model and Method

The task of mapping a thaw bowl under the roadbed and determining its dimensions is efficiently solved by means of electrotomography [30]. Via 2D and 3D inversion of field survey data, we derived a realistic resistivity distribution in depth and in area. Figure 3 shows a 2D inversion geoelectric section through the deformation area indicated by an arrow in Figure 2, and a resistivity map at a depth of 15 m following the 3D inversion results (Figure 3b).

The thaw bowls at the road base have low resistivities of 15–30 ohm∙m against the 150–300 ohm∙m frozen background, and the lateral dimensions from about 20 × 25 to 20 × 60 m in plan, the vertical thickness being up to 10 m, and the upper edge depth being approximately 10 m. The clear delineation of the thaw bowl makes it possible to precisely plan engineering works and increase the efficiency of anti-deformation measures.

Note that frozen rocks, in general, have very diverse resistivity values: from a few dozens of ohm∙m to millions of ohm∙m, depending on the particle size (clay or sand), temperature, and pore moisture salinity [31]. In the present case, we dealt with high-temperature permafrost with temperatures close to 0 °C (−0.3 … −0.1 °C), typical of central and southern Transbaikalia. According to drilling data from well 318, to a depth of 8.5 m there are frozen loams (light, icy), with a fluid consistency during thawing. At depths 8.5–10 m, there is frozen coarse sand of massive cryotexture, with 2–3 cm interlayers of sandy loam and loam every 15–20 cm; the sand is moist when thawed. A large amount of clay particles and high temperatures in the frozen ground are the reason for its relatively low resistivities in the case at hand.

However, in addition to mapping thermokarst, it becomes necessary to monitor the processes of geocryological changes at the road base, for example, after taking measures to regulate the permafrost of concern. The processes of permafrost thawing or freezing will be accompanied by a change in the electrical conductivity, which enables deploying resistivity exploration methods, including electromagnetic ones, for monitoring.

We propose a technique for monitoring permafrost under highways, based on the pulsed electromagnetic cross-well exploration of the geological space. With account taken for the above-described resistivity inversion results, the following geoelectric model and cross-well system parameters are considered in sectional view (Figure 4). At the top, there is a 2 m high trapezoidal roadbed with embankment (gray). Its upper base is 15 m long, and the length in the plan is 30 m (not shown in the Figure 4). The base angles of the trapezoid are 30°. The electrical resistivity of the road element under examination is equal to 300 ohm∙m. Below, there is a seasonally thawing rectangular layer (blue) with a height of 2 m and a resistivity of 100 ohm∙m. Permafrost (yellow) is characterized by a resistivity of 250 ohm∙m. The indicated geoelectric parameters are in a good agreement with the results of electrotomography field data interpretation.

Under the road and seasonally thawing layer, a thaw-affected area (talik) with a trapezoidal cross-section has formed (red). We investigated the cases of two different taliks both having a resistivity of 15 ohm∙m: a large one (5 m upper base length, 10 m distance from it to the embankment base) and small (2 m upper base length, 4 m distance).

For the cross-well electromagnetic survey of the medium under consideration, two vertical wells 15 m deep from the embankment base were drilled on both sides of the road, at a 25 m distance from each other (green). The first well contains the signal source (T), whereas the second one—the receiver (R). The source is a loop with a radius of 0.02 m, and the receiver is a coil of the same radius. A direct current with an amplitude of 1 A flows through the loop. At the time t = 0, the current is abruptly switched off, after which the temporal variation of the electromotive force (EMF) induced in the receiver coil is recorded. During the cross-well exploration, the source and receiver are located at the same depth level and are synchronously moved along the wellbores.

3. Mathematical Model

3.1. Boundary-Value Problem

To obtain a mathematical model describing the sounding process with a pulsed source of electromagnetic field excitation in a three-dimensional region with a complex physical and geometric structure, we use the Maxwell system of equations. As the boundary of the computational domain Ω, we will apply the surface ∂Ω distant from the source to an extent that the values of the sought fields can be assumed to be zero at ∂Ω. Since the electromagnetic field source will be excited by a current pulse, the initial values of the fields are taken to be zero. By employing the state equations, we exclude the electric and magnetic induction vectors from the Maxwell system of equations. Consequently, we get the following mathematical model (1)–(8):

$$\mathrm{rot}\overrightarrow{\mathrm{E}}=-\mathsf{\mu }\frac{\partial \overrightarrow{\mathrm{H}}}{\partial \mathrm{t}},$$

(1)

$$\mathrm{rot}\overrightarrow{\mathrm{H}}=\mathsf{\epsilon }\frac{\partial \overrightarrow{\mathrm{E}}}{\partial \mathrm{t}}{+\mathsf{\sigma }\overrightarrow{\mathrm{E}}+\overrightarrow{\mathrm{J}}}_0,$$

(2)

$$\mathrm{div}\mathsf{\epsilon }\overrightarrow{\mathrm{E}}=\rho ,$$

(3)

$$\mathrm{div}\mathsf{\mu }\overrightarrow{\mathrm{H}}=0,$$

(4)

$${\overrightarrow{\mathrm{E}}|}_{{\mathrm{t}=\mathrm{t}}_0}=0,$$

(5)

$${\overrightarrow{\mathrm{H}}|}_{{\mathrm{t}=\mathrm{t}}_0}=0,$$

(6)

$${\overrightarrow{\mathrm{E}}|}_{\partial \mathsf{\Omega }}=0,$$

(7)

$${\overrightarrow{\mathrm{H}}|}_{\partial \mathsf{\Omega }}=0,$$

(8)

where $\overrightarrow{\mathrm{E}}$ is the electric field strength, $\overrightarrow{\mathrm{H}}$ is the magnetic field strength, ${\overrightarrow{\mathrm{J}}}_0$ is the external electric current density, $\mathsf{\sigma }$ is the specific electrical conductivity, $\rho$ is the electric charge density, $\mathsf{\epsilon }$ is the dielectric permittivity, and $\mathsf{\mu }$ is the magnetic permeability.

Using the forward time-domain Fourier transform:

$$\mathrm{f}(\mathsf{\omega })=\underset{-\infty }{\overset{\infty }{{\displaystyle \int }}}{\mathrm{f}(\mathrm{t})\mathrm{e}}^{-\mathrm{i}\mathsf{\omega }\mathrm{t}}\mathrm{dt},$$

we transfer from the time domain to the frequency domain. The mathematical model (1)–(8) will take the following form:

$$\mathrm{rot}\overrightarrow{\mathrm{E}}=\mathrm{i}\mathsf{\omega }\mathsf{\mu }\overrightarrow{\mathrm{H}},$$

(9)

$${\mathrm{rot}\overrightarrow{\mathrm{H}}=-\mathrm{i}\mathsf{\omega }\mathsf{\epsilon }\overrightarrow{\mathrm{E}}+\mathsf{\sigma }\overrightarrow{\mathrm{E}}+\overrightarrow{\mathrm{J}}}_0,$$

(10)

$$\mathsf{\epsilon }\overrightarrow{\mathrm{E}}=\rho ,$$

(11)

$$\mathsf{\mu }\overrightarrow{\mathrm{H}}=0,$$

(12)

$${\overrightarrow{\mathrm{E}}|}_{\partial \mathsf{\Omega }}=0,$$

(13)

$${\overrightarrow{\mathrm{H}}|}_{\partial \mathsf{\Omega }}=0.$$

(14)

Since it is necessary to know the electromagnetic field values at several measurement points ${\mathrm{X}}_{\mathrm{i}}$ (i = 1,…, n), for each measurement point we will construct an interpolant ${\mathrm{I}}_{\mathrm{i}}(\mathsf{\omega })$ based on a set of solutions to the problem (9)–(14) at different angular frequencies ${\mathsf{\omega }=\mathsf{\omega }}_{\mathrm{j}}$ (j = 1,…, k). Taking the inverse frequency-domain Fourier transform:

$$\mathrm{f}(\mathrm{t})=\frac{1}{\mathsf{\pi }}\mathrm{Re}(\underset{0}{\overset{\infty }{{\displaystyle \int }}}{\mathrm{f}(\mathsf{\omega })\mathrm{e}}^{\mathrm{i}\mathsf{\omega }\mathrm{t}}\mathrm{d}\mathsf{\omega })$$

to the interpolant ${\mathrm{I}}_{\mathrm{i}}(\mathsf{\omega })$, we obtain the required electromagnetic field components at the time $\mathrm{t}$ and measurement point ${\mathrm{X}}_{\mathrm{i}}$ [32].

Eliminating the magnetic field strength from the Equations (9)–(14), we have the boundary-value problem:

$$\mathrm{rot}\frac{1}{\mathsf{\mu }}{\mathrm{rot}\overrightarrow{\mathrm{E}}+\mathrm{k}}^2{\overrightarrow{\mathrm{E}}=-\mathrm{i}\mathsf{\omega }\overrightarrow{\mathrm{J}}}_0,$$

(15)

$$\overrightarrow{\mathrm{E}}{|}_{\partial \mathsf{\Omega }}=0,$$

(16)

where ${\mathrm{k}}^2{=\mathrm{i}\mathsf{\omega }\mathsf{\sigma }-\mathsf{\omega }}^2\mathsf{\epsilon }$.

A consequence of (10) is the charge conservation law:

$$\mathrm{div}((\mathsf{\sigma }-\mathrm{i}\mathsf{\omega }\mathsf{\epsilon })\overrightarrow{\mathrm{E}}+{\overrightarrow{\mathrm{J}}}_0)=0.$$

(17)

3.2. Vector Variational Formulation

Let $\mathsf{\Omega }$ be a three-dimensional, possibly inhomogeneous in physical properties domain with a Lipschitz-continuous boundary $\partial \mathsf{\Omega }$. We introduce the function spaces [33,34]:

$$\begin{array}{c}\mathrm{H}(\mathrm{rot};\mathsf{\Omega })=\left\{\overrightarrow{\mathrm{v}}\in {[\mathrm{L}}^2(\mathsf{\Omega }){]}^3:\mathrm{rot}\overrightarrow{\mathrm{v}}\in {[\mathrm{L}}^2(\mathsf{\Omega }){]}^3\right\},\\ {\mathrm{H}}_0(\mathrm{rot};\mathsf{\Omega })=\{\overrightarrow{\mathrm{v}}\in {\mathrm{H}(\mathrm{rot};\mathsf{\Omega }):\overrightarrow{\mathrm{v}}\times \overrightarrow{\mathrm{n}}|}_{\partial \mathsf{\Omega }}=0\text{}},\end{array}$$

with the norm:

$$\parallel \overrightarrow{\mathrm{u}}{\parallel }_{\mathrm{rot},\mathsf{\Omega }}^2=\underset{\mathsf{\Omega }}{{\displaystyle \int }}\overrightarrow{\mathrm{u}}\cdot {\overrightarrow{\mathrm{u}}}^{\ast }\mathrm{d}\mathsf{\Omega }+\underset{\mathsf{\Omega }}{\overset{}{{\displaystyle \int }}}\mathrm{rot}\overrightarrow{\mathrm{u}}\cdot {\mathrm{rot}\overrightarrow{\mathrm{u}}}^{\ast }\mathrm{d}\mathsf{\Omega }.$$

Let us define the scalar product:

$$(\overrightarrow{\mathrm{u}},\overrightarrow{\mathrm{v}})=\underset{\mathsf{\Omega }}{\overset{}{{\displaystyle \int }}}\overrightarrow{\mathrm{u}}\cdot {\overrightarrow{\mathrm{v}}}^{\ast }\mathrm{d}\mathsf{\Omega }.$$

For the problem (15)–(16), we formulate the following variational statement [35]:

For ${\overrightarrow{\mathrm{J}}}_0\in {\mathrm{L}}^2(\mathsf{\Omega })$ find $\overrightarrow{\mathrm{E}}\in {\mathrm{H}}_0(\mathrm{rot};\mathsf{\Omega })$ such that $\forall \overrightarrow{\mathrm{v}}\in {\mathrm{H}}_0(\mathrm{rot};\mathsf{\Omega })$:

$$(\frac{1}{\mathsf{\mu }}{\mathrm{rot}\overrightarrow{\mathrm{E}},\mathrm{rot}\overrightarrow{\mathrm{v}})+(\mathrm{k}}^2{\overrightarrow{\mathrm{E}},\overrightarrow{\mathrm{v}})=-\mathrm{i}(\mathsf{\omega }\overrightarrow{\mathrm{J}}}_0,\overrightarrow{\mathrm{v}}).$$

(18)

For the introduced space, the further embedding property holds:

$$\mathrm{grad}\mathsf{\varphi }\in \mathrm{H}(\mathrm{rot};\mathsf{\Omega }),\forall \mathsf{\varphi }\in {\mathrm{H}}^1(\mathsf{\Omega }).$$

(19)

The variational formulation (18) is fulfilled for all $\overrightarrow{\mathrm{v}}\in \mathrm{H}(\mathrm{rot};\mathsf{\Omega })$; according to (19), we take $\overrightarrow{\mathrm{v}}=\mathrm{grad}\mathsf{\varphi },\mathsf{\varphi }\in {\mathrm{H}}^1(\mathsf{\Omega })$, in which case (18) has the form:

$$(\frac{1}{\mathsf{\mu }}\mathrm{rot}\overrightarrow{\mathrm{E}},\mathrm{rotgrad}\mathsf{\varphi }{)+(\mathrm{k}}^2\overrightarrow{\mathrm{E}},\mathrm{grad}\mathsf{\varphi })=-\mathrm{i}(\mathsf{\omega }{\overrightarrow{\mathrm{J}}}_0,\mathrm{grad}\mathsf{\varphi });\forall \mathsf{\varphi }\in {\mathrm{H}}^1(\mathsf{\Omega }).$$

Taking into account the property $\mathrm{rotgrad}\mathsf{\varphi }=0$, we get:

$${((\mathrm{i}\mathsf{\omega }\mathsf{\sigma }-\mathsf{\omega }}^2\mathsf{\epsilon })\overrightarrow{\mathrm{E}}+\mathrm{i}\mathsf{\omega }{\overrightarrow{\mathrm{J}}}_0,\mathrm{grad}\mathsf{\varphi })=0;\forall \mathsf{\varphi }\in {\mathrm{H}}^1(\mathsf{\Omega }),$$

(20)

$${((\mathrm{i}\mathsf{\omega }\mathsf{\sigma }-\mathsf{\omega }}^2\mathsf{\epsilon })\overrightarrow{\mathrm{E}}+\mathrm{i}\mathsf{\omega }{\overrightarrow{\mathrm{J}}}_0,\mathrm{grad}\mathsf{\varphi })=\underset{\mathsf{\Omega }}{\overset{}{{\displaystyle \int }}}{((\mathrm{i}\mathsf{\omega }\mathsf{\sigma }-\mathsf{\omega }}^2\mathsf{\epsilon })\overrightarrow{\mathrm{E}}+\mathrm{i}\mathsf{\omega }{\overrightarrow{\mathrm{J}}}_0)\cdot \mathrm{grad}\mathsf{\varphi }\mathrm{d}\mathsf{\Omega },$$

$$\underset{\mathsf{\Omega }}{\overset{}{{\displaystyle \int }}}{((\mathrm{i}\mathsf{\omega }\mathsf{\sigma }-\mathsf{\omega }}^2\mathsf{\epsilon })\overrightarrow{\mathrm{E}}+\mathrm{i}\mathsf{\omega }{\overrightarrow{\mathrm{J}}}_0)\cdot \mathrm{grad}\mathsf{\varphi }\mathrm{d}\mathsf{\Omega }=\underset{\mathsf{\Omega }}{\overset{}{{\displaystyle \int }}}{\mathrm{div}[(\mathrm{i}\mathsf{\omega }\mathsf{\sigma }-\mathsf{\omega }}^2\mathsf{\epsilon })\overrightarrow{\mathrm{E}}+\mathrm{i}\mathsf{\omega }{\overrightarrow{\mathrm{J}}}_0]\mathsf{\varphi }\mathrm{d}\mathsf{\Omega }.$$

(21)

It follows from (21) that Equation (20) is a variational analog of the conservation law (17). Consequently, the solution to the variational problem (18) satisfies the charge conservation law (17) in the weak sense.

3.3. Discrete Analogue of the Variational Problem

Traditionally, basis functions forming the space ${\mathrm{H}}^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega })$are divided into two types (in what follows, the index h in the name of a space denotes its discrete subspace). The first type was proposed in [33], and the second one in [34]. Following the embedding property (19), part of the space ${\mathrm{H}}^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega })$ consists of functions that are gradients of scalar functions. If scalar functions of an order not exceeding p are used to determine a basis of order p, then the resulting basis is a basis of the first type; if scalar functions are of order p + 1 at most, this is a basis of the second type. Thus, the second-type basis can be obtained from the first-type basis, supplementing it with basis functions that are gradients of scalar functions of order p + 1.

To create a discrete analogue of the variational problem, we will approximate the elements of the space $\mathrm{H}(\mathrm{rot};\mathsf{\Omega })$ by the elements of the discrete subspace ${\mathrm{H}}^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega })$. As its basis functions, we will take the second-type third-order vector elements on a tetrahedral grid [36].

To improve the spectral properties of the matrices obtained after discretizing the original problem, the basis functions can be orthogonalized. Their full orthogonalization would lead to a sharp increase in the number of non-zero matrix elements. In [36], partial orthogonalization is proposed, that is to split the basis functions into many groups and then perform orthogonalization within each group.

The high-order vector basis functions can be associated with the edge, face, or tetrahedron itself. It depends on determining the degree of freedom for the specific basis function: an integral along the edge, over the face, or over the entire geometric element, respectively. Since for the high-order bases one edge, face, or element is associated with a number of functions, this peculiarity will be utilized as a delimiter into groups. Such determination of orthogonalization groups does not add to the number of non-zero matrix elements, with its portrait being unaffected.

The article [36] addresses the standard scalar product for orthogonalization: The basis functions within one group are orthogonalized relative to the bilinear form applied to derive the variational formulation.

Now, we introduce a discrete analogue of the variational problem:

For $\overrightarrow{\mathrm{J}}\in {\mathrm{L}}^2(\mathsf{\Omega })$ find ${\overrightarrow{\mathrm{E}}}^{\mathrm{h}}\in {\mathrm{H}}_0^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega })$ such that $\forall {\overrightarrow{\mathrm{v}}}^{\mathrm{h}}\in {\mathrm{H}}_0^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega })$:

$$(\frac{1}{\mathsf{\mu }}{\mathrm{rot}\overrightarrow{\mathrm{E}}}^{\mathrm{h}}{,\mathrm{rot}\overrightarrow{\mathrm{v}}}^{\mathrm{h}}{)+(\mathrm{k}}^2{\overrightarrow{\mathrm{E}}}^{\mathrm{h}}{,\overrightarrow{\mathrm{v}}}^{\mathrm{h}}{)=-\mathrm{i}(\mathsf{\omega }\overrightarrow{\mathrm{J}},\overrightarrow{\mathrm{v}}}^{\mathrm{h}}).$$

As far as the constructed discrete subspaces are characterized by the inclusion property:

$${\mathrm{u}}^{\mathrm{h}}\in {\mathrm{H}}^{\mathrm{h}}(\mathrm{grad};\mathsf{\Omega })\to {\mathrm{gradu}}^{\mathrm{h}}\in {\mathrm{H}}^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega }),$$

the approximation of the electric field strength ${\overrightarrow{\mathrm{E}}}^{\mathrm{h}}$ will satisfy the charge conservation law in the weak form:

$${((\mathrm{i}\mathsf{\omega }\mathsf{\sigma }-\mathsf{\omega }}^2{\mathsf{\epsilon })\overrightarrow{\mathrm{E}}}^{\mathrm{h}}{,\mathrm{gradu}}^{\mathrm{h}})=0,\forall {\mathrm{u}}^{\mathrm{h}}\in {\mathrm{H}}_0^{\mathrm{h}}(\mathrm{grad};\mathsf{\Omega }).$$

Via expanding the electric field strength vectors ${\overrightarrow{\mathrm{E}}}^{\mathrm{h}}$ into the basis of the discrete subspace ${\mathrm{H}}_0^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega })$ and choosing the basis functions ${\overrightarrow{\mathrm{w}}}_{\mathrm{m}}^{\mathrm{h}}$ as the function $\overrightarrow{\mathrm{u}}$, we turn from the discrete variational problem to the equivalent system of linear equations:

$$[(\mathrm{A}+\mathrm{B})+\mathrm{iC}]\overrightarrow{\mathrm{E}}=\overrightarrow{\mathrm{f}},$$

(22)

where $\overrightarrow{\mathrm{E}}$ is the weight vector in the expansion into the basis ${\overrightarrow{\mathrm{E}}}^{\mathrm{h}}$, while the elements of the matrices A, B, C and those of the right-hand side vector $\overrightarrow{\mathrm{f}}$ are given by the relations:

$$\begin{array}{c}{[\mathrm{A}]}_{\mathrm{i},\mathrm{j}}=(\frac{1}{\mathsf{\mu }}\mathrm{rot}{\overrightarrow{\mathrm{w}}}_{\mathrm{i}}^{\mathrm{h}},\mathrm{rot}{\overrightarrow{\mathrm{w}}}_{\mathrm{j}}^{\mathrm{h}}{)}_{\mathsf{\Omega }},\\ {[\mathrm{B}]}_{\mathrm{i},\mathrm{j}}{=-(\mathsf{\epsilon }\mathsf{\omega }}^2{\overrightarrow{\mathrm{w}}}_{\mathrm{i}}^{\mathrm{h}}{,\overrightarrow{\mathrm{w}}}_{\mathrm{j}}^{\mathrm{h}}{)}_{\mathsf{\Omega }},\\ {[\mathrm{C}]}_{\mathrm{i},\mathrm{j}}{=(\mathsf{\sigma }\mathsf{\omega }\overrightarrow{\mathrm{w}}}_{\mathrm{i}}^{\mathrm{h}}{,\overrightarrow{\mathrm{w}}}_{\mathrm{j}}^{\mathrm{h}}{)}_{\mathsf{\Omega }},\\ {[\mathrm{f}]}_{\mathrm{i}}{=(\mathsf{\omega }\overrightarrow{\mathrm{J}},\overrightarrow{\mathrm{w}}}_{\mathrm{i}}^{\mathrm{h}}{)}_{\mathsf{\Omega }},\end{array}$$

${\overrightarrow{\mathrm{w}}}_{\mathrm{i}}^{\mathrm{h}}$ being the basis functions of the space ${\mathrm{H}}^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega })$.

3.4. Algorithm for Solving the System of Linear Equations

Let us introduce the following notation:

• ${\mathrm{H}}^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega };2)$ is the discrete subspace formed by the second-type basis.
• ${\mathrm{H}}^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega };1)$ is the discrete subspace formed by the first-type basis.
• ${\mathrm{N}}^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega };2)=\{\overrightarrow{\mathrm{u}}\in {\mathrm{H}}^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega };2);\nabla \times \overrightarrow{\mathrm{u}}=0\text{}}$ is the kernel of the curl operator in the space ${\mathrm{H}}^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega };2)$.
• ${\mathrm{P}}_1$ is the matrix whose columns form the coordinates of the basis functions in the subspace ${\mathrm{H}}^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega };1)$ in ${\mathrm{H}}^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega };2)$.
• ${\mathrm{P}}_2$ is the matrix whose columns form the coordinates of the basis functions in the subspace ${\mathrm{N}}^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega };2)$ in ${\mathrm{H}}^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega };2)$.

It should be noted that the kernel of the curl operator consists of functions that are gradients of scalar functions. Consistent with the method for constructing the spaces ${\mathrm{H}}^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega };2)$ and ${\mathrm{H}}^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega };1)$ [33,34], the following relationships take place:

$$\begin{array}{c}{\mathrm{H}}^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega };1)\cup {\mathrm{N}}^{\mathrm{h}}{(\mathrm{rot};\mathsf{\Omega };2)=\mathrm{H}}^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega };2),\\ {\mathrm{H}}^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega };1)\cap {\mathrm{N}}^{\mathrm{h}}{(\mathrm{rot};\mathsf{\Omega };2)=\mathrm{N}}^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega };1).\end{array}$$

It follows from assigning the subspaces that the matrices ${(\mathrm{P}}_2^{\mathrm{T}}{\mathrm{AP}}_2)$ and ${(\mathrm{P}}_1^{\mathrm{T}}{\mathrm{AP}}_1)$ will correspond to discrete variational problems formulated with respect to the vector Helmholtz equation in the spaces ${\mathrm{N}}^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega };2)$ and ${\mathrm{H}}^{\mathrm{h}}(\mathrm{rot};\mathsf{\Omega };1)$. In view of this, the finite-element technology can be used to obtain these matrices. For an approximate solution to the linear equation system (22), we formulate an iterative algorithm, which is an extended version of the multiplicative algorithm proposed in [37]. One iteration of the new algorithm is written as:

$$\begin{array}{l}{\mathrm{g}}_1{=\mathrm{P}}_1^{\mathrm{T}}{\mathrm{r}}_{\mathrm{i}-1},\\ {\mathrm{z}}_1{=\mathrm{S}}_1{((\mathrm{P}}_1^{\mathrm{T}}{\mathrm{AP}}_1{),\mathrm{g}}_1{,\mathsf{\gamma }}_1),\\ {\mathrm{x}}_{\mathrm{i}-1/2}{=\mathrm{x}}_{\mathrm{i}-1}{+\mathrm{P}}_1{\mathrm{z}}_1,\\ {\mathrm{r}}_{\mathrm{i}-1/2}{=\mathrm{b}-\mathrm{Ax}}_{\mathrm{i}-1/2},\\ {\mathrm{g}}_2{=\mathrm{P}}_2^{\mathrm{T}}{\mathrm{r}}_{\mathrm{i}-1/2},\\ {\mathrm{z}}_2{=\mathrm{S}}_2{((\mathrm{P}}_2^{\mathrm{T}}{\mathrm{AP}}_2{),\mathrm{g}}_2{,\mathsf{\gamma }}_2),\\ {\mathrm{x}}_{\mathrm{i}}^{\ast }{=\mathrm{x}}_{\mathrm{i}-1/2}{+\mathrm{P}}_2{\mathrm{z}}_2,\\ {\mathrm{r}}_{\mathrm{i}}^{\ast }{=\mathrm{b}-\mathrm{Ax}}_{\mathrm{i}}^{\ast },\\ {\mathrm{z}}^{\ast }{=\mathrm{S}(\mathrm{A},\mathrm{r}}_{\mathrm{i}}^{\ast },\mathsf{\gamma }),\\ {\mathrm{x}}_{\mathrm{i}}{=\mathrm{x}}_{\mathrm{i}}^{\ast }{+\mathrm{z}}^{\ast },\\ {\mathrm{r}}_{\mathrm{i}}{=\mathrm{b}-\mathrm{Ax}}_{\mathrm{i}},\end{array}$$

where ${\mathrm{x}}_{\mathrm{i}}=\mathrm{S}(\mathrm{A},\mathrm{b},\mathsf{\gamma })$ is the iterative solver for the linear equation system $\mathrm{Ax}=\mathrm{b}$, and $\mathsf{\gamma }$ determines how many times the residual norm of the approximate solution ${\mathrm{x}}_{\mathrm{i}}$ is less than the norm of the right-hand side (the initial approximation being zero). As the solver $\mathrm{S}$, we will take the symmetric QMR [38]. It is a variation of the QMR method, focused on solving linear equation systems with symmetric matrices. Further, we assume that $\mathsf{\gamma }=0.5$, ${\mathsf{\gamma }}_1=0.5$, ${\mathsf{\gamma }}_2=0.01$.

4. Parallel Implementation and High-Performance Computing

For a three-dimensional simulation of the sounding process with a pulsed current source, we deploy the Fourier time transform. Initially, it is necessary to solve many independent and laborious three-dimensional problems, whose solutions are taken into account together only when performing the inverse Fourier transform. Such an approach provides an opportunity to carry out the parallelization in a natural way, solving each frequency problem in a separate computational thread.

The approach proposed shows a significant advantage over the parallelization of matrix-vector multiplication (the most time-consuming operation in the solution process), since in this case the connection between different computational threads has to be accomplished at least once per iteration in the algorithm for solving the system of linear equations. Moreover, in this connection, the time for transferring large amounts of data between various computing units becomes significant. However, when the solution of one frequency problem requires a fairly large amount of random-access memory, depending on the characteristics of the computer technology exploited, a situation may arise when all the available memory is already used, but there are unused computing cores. Under such conditions, for a more complete use of the available computing resources, it becomes relevant to apply parallel algorithms to the matrix-vector multiplication procedure. The traditional way to parallelize the procedure of multiplying a large sparse matrix by a vector is selecting a block structure in the matrix and performing a set of multiplications of the matrix blocks by the vector ones in parallel, and then sequentially assembling the final product from the vector blocks.

In the presented research, the division of the original matrix into blocks is based on the structure of the high-order vector basis functions employed to construct a discrete analogue of the differential equation being solved. These functions are defined on a tetrahedron and can be associated with the edge, face, or tetrahedron itself, which depends on setting the degree of freedom for the specific basis function: an integral along the edge, over the face, or over the entire geometric element, respectively. As the functions under discussion (one edge, face or element) are related to several functions, we deploy this property as a divider into groups.

Each of these groups can be divided into two more subgroups: a subgroup of the basis functions that are gradients of scalar functions, and a subgroup of the basis functions that cannot be expressed in such terms. Thus, we have 6 groups of the basis functions, which allow us to split the matrix into 36 matrix blocks. Such a division always exists and does not depend on the computational grid in use. It makes the parallelization effect more stable and predictable in comparison with the formation of the matrix blocks through algorithms for decomposing the computational grid into subdomains. The procedure for parallelizing matrix multiplication is auxiliary with respect to frequency parallelization, the total number of involved computational threads being equal to the product of the number of frequency problems solved simultaneously by the number of threads used for matrix-vector multiplication.

The testing of the parallel algorithm for simulating the sounding process with a pulsed current source was carried out on the NKS-1P supercomputer at the Siberian Supercomputer Center of SB RAS in consequence of the resource intensity of the computations. The supercomputer comprises a set of computational nodes of two types, divided by task queues: Intel KNL and Broadwell (

Table 1. Parameters of the Intel KNL and Broadwell computational nodes.

Characteristic	KNL	Broadwell
Central processor	1 × Intel Xeon Phi 7290 (1.5 GHz, 72 cores × 4 threads)	2 × Intel Xeon E5-2697v4 (2.6 GHz, 16 cores × 2 threads)
Number of hardware threads	288	32
Random-access memory	96 GB	128 GB
Peak performance	3456 Gflop/s	1331.2 Gflop/s

).

We analyzed the algorithm performance for two types of the cluster nodes, while varying the number of frequencies at which the signals of the sounding system were simulated in parallel, as well as the number of computational threads on which the problem was solved for a specific frequency [39] (Figure 5).

Due to the limited amount of random-access memory and the large number of computational threads, it is impossible to utilize all the threads on the KNL nodes when parallelizing the algorithm by frequency tasks. It becomes expedient to parallelize both frequency tasks and matrix multiplication operations within individual tasks. The best performance when using the KNL nodes is achieved when all the 288 hardware threads are involved, in which case the problem is solved in parallel for 48 frequencies, and matrix multiplication operations are performed on 8 threads. The total time for solving the problem is about 6 h.

Fewer computational threads are available when applying the Broadwell nodes, but with higher performance. The large amount of random-access memory allows for solving each frequency problem in a separate computational thread, thereby achieving the maximum performance [40]. In this case, the minimum time for solving the problem is approximately 4 h, with all the 32 hardware threads being used for frequency parallelization, and matrix-vector multiplication being performed on one computational thread.

Furthermore, the algorithm performance analysis results indicate that for both node types it is optimal, if possible, to make use of 8 threads when parallelizing matrix-vector multiplication operations. Apparently, with an increase in the number of threads, the time required to synchronize data between the threads goes up; to solve one problem, such synchronization has to be performed many times. In addition, it should be borne in mind that the blocks into which the matrix is divided have different sizes, which also slows down the synchronization. In the case of frequency parallelization, such resource-intensive synchronization does not need to be carried out; therefore, the time for solving the problem monotonically decreases with an increase in the number of the threads being deployed.

5. Numerical Simulation Results

The simulated temporal variations of the EMF induced in the receiver coil are presented on a logarithmic scale for four depths: 1, 5, 10, and 15 m (Figure 6). The zero depth corresponds to the embankment base in Figure 4.

The EMF in the reference model has the simplest form: Its transition through zero is evidenced only at a depth of 15 m at late times. As for both talik-containing models, the corresponding EMFs differ from the mentioned one significantly: Zero crossing is manifested on all the graphs. At the same time, the crucial difference between the latter two lies in the locations of the extrema when the depth being explored changes. More specifically, if the talik is large, when increasing the depth, the extremum time does not seem to change and the EMF value becomes slightly smaller. For the small talik, a rise in the depth leads to both a shift in the extremum of the EMF graph towards short times and an increase in the measured signal. The largest difference between the EMFs regarding the zero crossings and absolute values is observed in the middle range of depths (for instance, at 5 m). At a depth of 15 m, the size of the heterogeneity (talik) ceases to be discernible, since the EMF dependences become almost identical. The main conclusion from the analysis is the capability of the proposed cross-well electromagnetic system to reliably indicate the presence of a talik under a highway on permafrost, and estimate its size as well.

Let us consider the relative differences between the EMF in the reference model and those with a talik (Figure 7; 1, 3—small talik and 2, 4—large talik) as a function of depth along the well. The first two (1 and 2) are obtained at the short time t₁ = 2∙10⁻⁷ s, while the second pair (3 and 4) is attributed to the late time t₂ = 10⁻⁶ s. A careful examination reveals the pairwise differences between the graphs 1, 2 and between 3, 4 to be significant, especially for the time t₂. The exception is the maximum depths, where the measuring system sensitivity to the taliks is practically absent subsequent to their remoteness.

The suggested pulsed electromagnetic cross-well system is sensitive to the presence and size of a talik under a road embankment in permafrost areas. The simulation results are compliant with [23], where the same federal highway was explored by surface electromagnetic induction surveys. The principal regional-scale result obtained there was that the induction responses to thawed permafrost decay with depth at rates correlated with temperatures, which is why the induction-decay rates are suggested to be targeted for monitoring permafrost temperature changes.

6. Conclusions

The results of our investigation were aimed at the scientific substantiation of a geophysical technique for monitoring the processes of geocryological changes at a highway base. To study the processes of permafrost thawing or freezing, which are accompanied by a change in the electrical conductivity, we proposed pulsed electromagnetic sounding. Numerical simulation was conducted in order to investigate the capability of employing a pulsed electromagnetic cross-well exploration of the geological environment.

We developed a mathematical model that describes the process of cross-well sounding with a pulsed source of electromagnetic field excitation in a three-dimensional region complex in physical parameters and geometric structure. For the discretization, the Fourier time transform was employed, which makes it possible to substantially simplify the analysis of the current pulse shape. This is due to the application of the fundamental solution in time and its convolution with the excitation pulses of different shapes, with no need for solving a great number of three-dimensional problems.

Through the vector finite-element method, we obtained a discrete analogue of the original problem in spatial coordinates. The discrete variational problem has an important property, namely, its solution satisfies the charge conservation law in the weak sense. To construct the matrix and the right-hand side of the linear equation system, we utilized hierarchical vector basis functions of the third order, defined on a tetrahedral grid.

To improve the speed of solving the equation system, we deployed a partial orthogonalization of the basis functions and a specially developed extended multiplicative algorithm based on the properties of the discrete function space used to construct the discrete variational problem. The testing and computational experiments show the high efficiency of the devised computational scheme.

For the fullest usage of the computational resources of the NKS-1P supercomputer at the Siberian Supercomputer Center of SB RAS, the authors developed a two-level approach to the parallelization of the problem being solved. We divided the initial problem into frequencies to be computed on individual cores of the cluster nodes with a sufficient amount of random access memory. In its turn, to make full use of the computational resources of the nodes with a limited amount of memory, in addition to the previous technique we parallelized matrix-vector multiplication. The computational experiments exhibited a strong performance.

We simulated temporal variations of the electromotive force in the receiver coil at different depths when the source and receiver were moved along the wellbores opposite each other. In addition, we calculated the relative differences between the electromotive force in the reference model and those with a talik at short and late times as a function of depth. It follows from the simulation results that the proposed pulsed electromagnetic cross-well exploration system has a high sensitivity to the presence and size of a talik under a highway over permafrost.

The proposed pulsed electromagnetic cross-well monitoring system has an advantage over the conventional surface electrotomography in cases where there is no possibility to ground the electrodes. Conversely, its own practical application can be limited by the local topography, bogginess, and vandalism.

The devised algorithm and its software implementation are one of the main computing tools for the scientific substantiation and optimal design of the configurations of pulsed cross-well sounding systems applied to solving a wide range of problems of monitoring the state of permafrost under civil and industrial facilities, buildings, and constructions.

In the future, we will improve the geoelectric model on the basis of interpreting more field geophysical data, and the current computational experiments open new prospects for monitoring in permafrost regions.