Detecting Permafrost in Plateau and Mountainous Areas by Airborne Transient Electromagnetic Sensing

: Transportation has become a key bottleneck which restricts economic development in Western China. However, during the construction of the western railway, the permafrost problem has plagued railway construction on the Qinghai–Tibet Plateau, and has not yet been resolved. Accurately identifying permafrost by geophysical method is the most e ﬀ ective means to solve this problem. However, the mountainous and plateau terrain in Western China impose huge challenges in collecting geophysical data. To address this issue, this paper proposes an airborne transient electromagnetic method to collect geophysical electromagnetic data to identify permafrost in the mountains and plateaus of Western China. Based on Maxwell’s equations, the forward model of the airborne electromagnetic was derived, and the ﬁnite element method was used to calculate the two-dimensional (2D) space electromagnetic responses of di ﬀ erent permafrost geo-electrical models. Furthermore, a coupling function was constructed to estimate the distribution of the resistivity of the permafrost by the least-squares ﬁtting algorithm. Comparison between inversion resistivity distribution and the geo-electrical model showed that the proposed airborne transient electromagnetic method was valid for exploring the permafrost in the mountains and the Qinghai–Tibet Plateau in Western China.


Introduction
In order to reduce the economic gap between Eastern and Western China, the National Council established a group for the Development of the Western Region in 2000, and the Vice Premier Jiabao Wen served as the Deputy Leader [1]. Although 20 years have passed, it has not yet achieved the desired results. The main reason is that most of the larger industrial bases are built in Eastern China, whereas only energy bases are built in Western China. The energy bases cannot drive much local economic development. Recently, a new campaign, the "Master Plan for the New Land, Sea, and Western Corridor", was promulgated by Chinese government in August 2019, with the aim of promoting the deep integration of transportation to serve the economic development of Western China [2]. As is well known, railway transportation is the most important tool in Western China. Therefore, railway construction plays a pivotal role to develop economics in Western China. However, a large amount of permafrost with high thermal sensitivity in these areas would bring great potential dangers to railway transportation [3][4][5][6][7][8][9]. Geophysical data can provide important information to detect the permafrost [10][11][12]. Hence, identifying the geological conditions of the railway subgrade and locating the permafrost are crucial for the railway construction in Western China. Firstly, research on The purpose of the development of the airborne electromagnetic method was mainly to overcome difficult terrain factors and enhance work efficiency. The US [19] and Australia [20][21][22] are the most successful in applying airborne electromagnetic methods to explore metal, minerals, and resources in remote rural areas with little population. This was a good reference for the investigation of the frozen soil in the Qinghai-Tibet Plateau in China, and it was a good choice to collect geophysical data by airborne electromagnetic method. Although there are few related works in the literature, we tried to do theory research to verify its feasibility to provide the theoretical basis for its application in the future.
Due to the steep terrain and high altitude, it is quite difficult to carry out ground surface geophysical exploration along the survey line. Compared with these contact geophysical methods, the non-contact electromagnetic methods were better solution. The resistivity difference between permafrost and surrounding rock was notable. As a result, one could analyze the resistivity of the permafrost area using the geophysical electromagnetic method to detect the permafrost [23]. The purpose of the development of the airborne electromagnetic method was mainly to overcome difficult terrain factors and enhance work efficiency. The US [19] and Australia [20][21][22] are the most successful in applying airborne electromagnetic methods to explore metal, minerals, and resources in remote rural areas with little population. This was a good reference for the investigation of the frozen soil in the Qinghai-Tibet Plateau in China, and it was a good choice to collect geophysical data by airborne electromagnetic method. Although there are few related works in the literature, we tried to do theory research to verify its feasibility to provide the theoretical basis for its application in the future.

Qinghai-Tibet Plateau
Due to the steep terrain and high altitude, it is quite difficult to carry out ground surface geophysical exploration along the survey line. Compared with these contact geophysical methods, the non-contact electromagnetic methods were better solution. The resistivity difference between permafrost and surrounding rock was notable. As a result, one could analyze the resistivity of the permafrost area using the geophysical electromagnetic method to detect the permafrost [23]. According to field test results [24,25] between using the ground surface electromagnetic method and the airborne electromagnetic method, the airborne electromagnetic method was superior to the former because the noise caused by the violent terrain for ground surface electromagnetic was much more serious than for airborne electromagnetic methods. Therefore, it will broaden the application to conduct airborne electromagnetic permafrost detection in mountainous and plateau regions. However, the detection accuracy of airborne electromagnetic method was still an unresolved problem [26].
Recently, the transient electromagnetic method (TEM) was found to be able to accurately detect the depth of permafrost and effectively locate the area of permafrost [27,28]. For example, Yi et al. [29] simulated the freeze-thaw process of thermal karst ponds and lakes in northern Siberia using TEM and found that TEM could perfectly describe the process of freeze-thaw. However, very little work has explored permafrost detection using TEM and other airborne electromagnetic methods, especially for research of permafrost regions in Western China. To address this issue, this paper attempted to study, by using the airborne transient electromagnetic method to explore, the shallow geological problems of permafrost detection in the western plateau cold region by numerical simulation and hoped to provide reliable data for the evaluation of permafrost disasters.

Forward Modeling of Airborne Transient Electromagnetic Method
The transient electromagnetic method (TEM) is a kind of time domain method which employs coils with currents as the transmitters to cause the conductive geological body and generate the inductive electromagnetic field. The geological targets can be analyzed by measuring the inductive electromagnetic field. The equipment carried by an aircraft to explore the earth is called the airborne transient electromagnetic method.

Research on Forward Modeling Theory
The electromagnetic field excited by a magnetic source satisfied the following Maxwell equation [30,31] where, ε, µ, ω, and σ are the permittivity, magnetic permeability, frequency, and conductivity of the medium, respectively; H and E are the electric and magnetic fields, respectively; and J m is the magnetic current density. The electromagnetic field was divided into a primary field (H p , E p ) and a scattering field (H s , E s ), and their relationship could be described as The primary field satisfies the following equations: Combining Equations (1)-(6), the following expressions can be obtained where ε s = ε − ε p , σ s = σ − σ p , and µ s = µ − µ p indicate the difference of the permittivity, resistivity, and permeability between the background field and the anomalous geological bodies. If the earth is isotropic media, Equations (7) and (8) could be simplified as where J s = σ s E p andŷ = σ + iεω. Based on Equations (9) and (10), one could note that the source was governed by the primary electromagnetic field and the conductivity difference in the area of the anomalous geological bodies, which avoided the source singularity. To simplify the expressions in Equations (9) and (10), let H and E instead of H s and E s . The curl operator was introduced into Equations (9) and (10) to obtain the following expressions Based on Equations (11) and (12), a 2D geo-electric model, as shown in Figure 2, was adopted to perform numerical simulations. Let us assume that the electrical property did not vary in the y-direction. In this case, the electrical field, magnetic field, and source should be processed by the Fourier transform whereF is the Fourier transform of H or E. Two coupled governing differential equations forÊ y and H y can be defined as The other components were obtained from the space-derivatives of E y and H y as follows Electronics 2020, 9, 1229 5 of 12 Electronics 2020, 9, x FOR PEER REVIEW 5 of 12 The Galerkin method [32] was employed to solve the above equations. Firstly, the isoparametric elements method was used to establish the finite-element equations of the 2D model in Figure 1, because, in the space discretization, the isoparametric elements were flexible to approximate the complex structure in the subsurface formation. In addition, the higher-order interpolation function was able to accurately calculate the space-derivative of ̑ and ̑ in the finite-element equations. Therefore, a quadratic element with eight nodes, as shown in Figure 3, was adopted for calculating the finite-element equations. In each element, ̑ and ̑ were expressed by the quadratic interpolation function, and the space-derivative could be described with a linear interpolation function. In each element, , , were constant, and the interpolation function was represented with the local coordinates ( , ) [33]. The global coordinates and each component F in the calculation domain could be written as where e j N is the quadratic interpolation function of the i-th node in the e-th element. The Galerkin method [32] was employed to solve the above equations. Firstly, the isoparametric elements method was used to establish the finite-element equations of the 2D model in Figure 1, because, in the space discretization, the isoparametric elements were flexible to approximate the complex structure in the subsurface formation. In addition, the higher-order interpolation function was able to accurately calculate the space-derivative of E y and H y in the finite-element equations. Therefore, a quadratic element with eight nodes, as shown in Figure 3, was adopted for calculating the finite-element equations. In each element, E y and H y were expressed by the quadratic interpolation function, and the space-derivative could be described with a linear interpolation function. In each element, σ, µ, ε were constant, and the interpolation function was represented with the local coordinates (ξ, η) [33]. The global coordinates and each component F in the calculation domain could be written as where N e j is the quadratic interpolation function of the i-th node in the e-th element.
Electronics 2020, 9, x FOR PEER REVIEW 6 of 12 Applying the Galerkin method to Equations (13) and (14), we obtained the following linear equations: Applying the Galerkin method to Equations (13) and (14), we obtained the following linear equations: where the second integrals on the right hand sides are line integrals around the boundary ∂D e of the element regions; D e , E l , and H l are the tangential components for ∂D e , respectively; and N e is the total number of the finite elements. Finally, the matrix equation was obtained as where X = E y (1), E y (2) · · · E y (N), H y (1), H y (2) · · · H y (N) T is the electrical and magnetic fields located in the discrete nodes in the model; N is the number of discrete nodes; and B can be regarded as the source item, which is the product of the primary field and the conductivity difference between the background conductivity and the anomalous geological bodies' conductivity. In fact, σ s is not equal to zero only in the areas of the anomalous geological bodies; hence, E p can be solved only in the nodes of the anomalous geological bodies.
After calculating the electromagnetic field in the k y -domain, the results should be transformed from k y -domain to real space by the inverse Fourier transformation [34]. The electromagnetic fields in the y-domain in the real space were evaluated by F(x, y = y 0 , z, ω) = 1 2π k y −k yF x, k y , z, ω e ik y y dk y , k y → ∞ The integration in Equation (25) was carried out by the cubic spline interpolation of the solutions [35] in the logarithmic ware-number domain k y .
The electromagnetic field obtained from the above calculations was the frequency domain, and to obtain the electromagnetic field in the time domain we should transform the electromagnetic field in the frequency domain into the electromagnetic field in the time domain. Compared with other transform methods, we could obtain much more accurate results in the later period in the domain by sine and cosine transformation [36], as shown in Equation (26).
Electronics 2020, 9, 1229 7 of 12 The time derivative of H(s) could be described as Equation (27). In fact, measurement data by the receiving coil was ∂t , which is called induced voltage.

Numerical Simulation One
A large amount of permafrost is distributed in the mountains and on the Qinghai-Tibet Plateau in Western China. The permafrost has high-temperature characteristics [37], which means that it has thermal sensitivity, high ice content, and an easy change from the solid state of ice to the flow state of water. These characteristics would make railway foundations face great danger in spring and autumn seasons [38] especially; the danger of thawing is that it is easy to occur during the period of the spring season due to ice melting. It is crucial to detect the permafrost for the Qinghai-Tibet railway to reduce potential danger during the railway construction. Considering that the permafrost overlays are very common in Western China plateaus, and a low resistivity layer is formed due to ice melting, a typical geo-electrical model, shown in Figure 4, was designed to simulate the responses of the airborne transient electromagnetic methods. In Figure 4, there are two permafrost distributions with different volumes, i.e., 50 m × 100 m and 50 m × 50 m, respectively, and the thickness of the cover layer is 5 m. According to the electrical characteristics of the permafrost, the resistivity in the melting zone of the permafrost was set to 10 Ohm-m, and the resistivity surrounding the rocks was set to 1000 Ohm-m. A central loop style was adopted to collect data; the height of the receiving coil was 30 m and the radius of the transmitting coil was 10 m. In addition, the number of turns was 10 and the working current was 100 A. The transient electromagnetic responses were calculated and are shown in Figure 3. Because Western China is dry all the year round and the shallow surface of the earth is dry, the secondary electromagnetic field was very weak where it was far away from the thawing region, as shown in Figure 4. The curves in Figure 4 are the electromagnetic field of the z-component at the top part, and the rectangle filled with color is the geological model. As can be seen, the airborne transient electromagnetic response along the measured line from left to right is shown in Figure 4, where the electromagnetic anomalous response corresponded to the geo-electric model very well.

Numerical Simulation Two
As mentioned in introduction, there are severe terrains in the Qinghai-Tibet Plateau in Western China. To address this issue, the following numerical simulation employed a discretized geo-electrical model with a small sunken, as shown in Figure 5. The sunken is usually at a lower level than the surrounding area, and it is often submerged in water during the rainy seasons. Therefore, in winter, the sunken would become frozen soil to gradually form the permafrost. As shown in Figure 5, the ground surface is covered by a low resistivity overlay with a thickness of 5 m.

Numerical Simulation Two
As mentioned in introduction, there are severe terrains in the Qinghai-Tibet Plateau in Western China. To address this issue, the following numerical simulation employed a discretized geo-electrical Electronics 2020, 9, 1229 8 of 12 model with a small sunken, as shown in Figure 5. The sunken is usually at a lower level than the surrounding area, and it is often submerged in water during the rainy seasons. Therefore, in winter, the sunken would become frozen soil to gradually form the permafrost. As shown in Figure 5, the ground surface is covered by a low resistivity overlay with a thickness of 5 m. There is one low anomalous area of 30 m × 100 m. Similar to simulation model one, the low resistivity area was set to 10 Ohm-m, and the high resistivity surrounding earth was set to 1000 Ohm-m. The simulation calculation results are described in Figure 5. It can be seen in Figure 5 that three airborne transient electromagnetic responses were observed, and all values of the electromagnetic curves were positive, which was different from model one without the sunken. However, only the middle one of the airborne transient electromagnetic responses was true to simulation model two, whereas the other two were unreal responses. The reason is that the area of sunken was filled with air, whose high resistivity would influence the transient electromagnetic calculation. Compared with the medium of air, the earth medium, distributed in the left and right sides of the sunken, had a relatively low resistivity. As a result, there were two false anomalous responses. The middle anomalous response corresponds to a low resistivity anomalous geological body of permafrost distribution.
Electronics 2020, 9, x FOR PEER REVIEW 9 of 12 Figure 5. Geo-electrical model with terrain and the forward modelling results.

Airborne Electromagnetic Inversion
Airborne transient electromagnetic inversion continuously adjusted the iterative models to obtain the best responses to fit the forward data indicated by (k, P) where Δ is the observed data, such as receiving the data of ( ) by coil, and ( , ) represents the theoretical response calculated using Formula (22) to Formula (25) according to the geo-electrical model described by resistivity distribution P ( : , , , , ⋯ , ). K indicates the number of samples and M represents the total number of sampling points. Applying the Taylor expansion to ( , )at p 0 and neglecting the terms of the second and higher orders, it yields

Airborne Electromagnetic Inversion
Airborne transient electromagnetic inversion continuously adjusted the iterative models to obtain the best responses to fit the forward data indicated by f (k, P) and the observation data indicated by ∆z k . It firstly needed to construct a coupling function [39], as shown in Equation (28), Electronics 2020, 9, 1229 9 of 12 where ∆z k is the observed data, such as receiving the data of ∂H(t) ∂t by coil, and f (k, P) represents the theoretical response calculated using Formula (22) to Formula (25) according to the geo-electrical model described by resistivity distribution P (P : ρ 1 , ρ 2 , ρ 3 , ρ 4 , · · · , ρ N ). K indicates the number of samples and M represents the total number of sampling points. Applying the Taylor expansion to f (k, P) at p 0 and neglecting the terms of the second and higher orders, it yields where the superscript (k) means k-th iteration. Considering the misfit ϕ, Equation (28) should satisfy the following equation By introduction of damping operator [39], Equation (31) could be rewritten in a matrix as (∆z k − f (k, P (0) )) , ∆P = (∆p 1 , · · · , ∆p n ); I is an identity matrix; and λ is the damping factor.

Inversion Numerical Simulation One
Based on forward numerical simulation results of the geological model as described by Figure 4, we processed the forward numerical simulation data by inversion method as mentioned by Formula (32). In the procedure, the calculated result ∂H(t) ∂t can be regarded as ∇z k in Formula (28). Finally, we got the inversion result, as shown in Figure 6. Compared with the geo-electrical model as described by Figure 4, they match each other very well.
Electronics 2020, 9, x FOR PEER REVIEW 10 of 12 Formula (32). In the procedure, the calculated result ( ) can be regarded as ∇ in Formula (28).
Finally, we got the inversion result, as shown in Figure 6. Compared with the geo-electrical model as described by Figure 4, they match each other very well.

Inversion Numerical Simulation Two
A similar process as inversion numerical simulation one was done for the forward modeling results by geo-electrical model with terrain, and Figure 7 manifests the inversion result. Comparing the inversion result with the simulated geo-electrical model, one could note that, although terrain fluctuations exist, the airborne transient electromagnetic method was still effective for detecting the permafrost.

Inversion Numerical Simulation Two
A similar process as inversion numerical simulation one was done for the forward modeling results by geo-electrical model with terrain, and Figure 7 manifests the inversion result. Comparing the inversion result with the simulated geo-electrical model, one could note that, although terrain fluctuations exist, the airborne transient electromagnetic method was still effective for detecting the permafrost.

Inversion Numerical Simulation Two
A similar process as inversion numerical simulation one was done for the forward modeling results by geo-electrical model with terrain, and Figure 7 manifests the inversion result. Comparing the inversion result with the simulated geo-electrical model, one could note that, although terrain fluctuations exist, the airborne transient electromagnetic method was still effective for detecting the permafrost.

Conclusions
Railway construction in Western China is very important for economic development. However, the permafrost problem has always plagued the railway construction. It was imperative to explore the distribution of permafrost by geophysical methods to prevent potential geological disaster. According to the terrain characteristics of the Qinghai-Tibet Plateau, a new airborne transient electromagnetic method was developed to investigate two kinds of two geo-electric models: one was flat ground surface, and the other one was undulating ground surface. To achieve forward simulation, a new algorithm that calculated the airborne transient electromagnetic responses was proposed, and the detailed characteristics of the calculated electromagnetic responses of the geo-electrical models were presented. For the purpose of estimating the distribution of the permafrost using the forward modeling data, a coupling function based on the least squares was constructed to fit the observation data and the calculation data. The calculation results of the two simulation models demonstrated that it was feasible to solve the permafrost problem of Qinghai-Tibet Plateau in Western China by the proposed airborne transient electromagnetic method. The location and degree of the permafrost in the simulated models could be effectively recognized. In future plans, the proposed method will be further verified by field data.

Conclusions
Railway construction in Western China is very important for economic development. However, the permafrost problem has always plagued the railway construction. It was imperative to explore the distribution of permafrost by geophysical methods to prevent potential geological disaster. According to the terrain characteristics of the Qinghai-Tibet Plateau, a new airborne transient electromagnetic method was developed to investigate two kinds of two geo-electric models: one was flat ground surface, and the other one was undulating ground surface. To achieve forward simulation, a new algorithm that calculated the airborne transient electromagnetic responses was proposed, and the detailed characteristics of the calculated electromagnetic responses of the geo-electrical models were presented. For the purpose of estimating the distribution of the permafrost using the forward modeling data, a coupling function based on the least squares was constructed to fit the observation data and the calculation data. The calculation results of the two simulation models demonstrated that it was feasible to solve the permafrost problem of Qinghai-Tibet Plateau in Western China by the proposed airborne transient electromagnetic method. The location and degree of the permafrost in the simulated models could be effectively recognized. In future plans, the proposed method will be further verified by field data.