Model of time-distance curve of electromagnetic waves diffracted on a local feature in the georadar study of permafrost zone rock layers

Introduction

One of the geophysical problems to be solved by using GPR (georadar) method is the study of physical and mechanical properties of rocks. However, accumulation of experimental data and development of methodological support for GPR in this area is much slower [1, 2] than in other areas [3] that leads to underestimation of the capabilities of the georadar method.

The reasons that led to this state of affairs in georadar application can be different. One of them is the incorrect use of the procedure for determining the velocity v of electromagnetic wave (EMW) propagation using hyperbolic time-distance curves (lineups of georadar impulses). This procedure is the most common way to estimate EMW velocity [1], based on which the material part of the relative complex dielectric permittivity ε′, depending on moisture, density, and cryogenic state of rocks, is calculated [4]. In training manuals (both domestic [1, 2] and foreign [5, 6]), as well as in data processing manuals of georadar manufacturers (GSSI, GEOTECH) and American Society for Testing and Materials (ASTM) standard¹, justification of application of EMW propagation velocity determination by hyperbolic time-distance curves is given for the case when the host medium is homogeneous.

Currently, in the practice of georadar measurements [3, 7, 8], as well as in scientific works devoted to the automation of the search for hyperbolic time-distance curves in GPR data [9–11], including on a real-time basis [12], the study subject is, as a rule, an inhomogeneous medium. As a consequence, the EMW propagation velocity determined by a hyperbola located in some layer is an averaged (integral) characteristic of all overlying layers, as mentioned in the work of one of the GPR classics [13]. When conducting georadar studies in a permafrost zone, it is possible to incorrectly assess the cryogenic state of rocks and, correspondingly, their physical and mechanical properties in the presence of a layer of thawed rocks, whose effect on the shape of the hyperbolic time-distance curve may be significant but insufficient for the result of v determination by the method of approximation of the time-distance curve by a hyperbola [2] proved to be within the range of values characteristic of thawed rocks. Thus, in a layered medium it is possible to determine the true velocity v of EMW propagation directly from the hyperbolic time-distance curve only in the first layer of rocks, and for the correct use of v values in the practice of georadar works it is necessary to establish the regularities of formation of EMW time-distance curves diffracted at a local feature in a layered rock mass. In order to achieve the above goal, the following tasks need to be accomplished:

• develop a model of a hyperbolic time-distance curve of impulses obtained in the study of a layered rock mass;
• establish the dependence of v and ε′ determined from the hyperbolic time-distance curve on the values of v and ε′ of the overlying layers;
• determine the influence of a thawed layer in a frozen rock mass on the v value calculated from the hyperbolic time-distance curve;
• verify the validity of the obtained theoretical results using the data of computer simulation.

1 ASTM D6432-11, Standard guide for using the surface ground penetrating radar method for subsurface investigation, ASTM International, West Conshohocken, PA; 2011. https://doi.org/10.1520/D6432-11

Model of a hyperbolic time-distance curve of georadar impulses obtained from probing of a layered rock mass

EMW emitted by a georadar located at point x propagates in a layered rock rock mass according to Fermat’s principle (gray solid line in Fig. 1), but in this paper we consider a model of EMW propagation along a raypath (black dashed line in Fig. 1) in a rock mass consisting of n layers of hi thickness with given values of ε′_i and v_i with i ranging 1 to n. A local feature, indicated by the black circle in Fig. 1, is located in the last layer at depth h0. The distance h_r, traveled by the ray from georadar to the local feature, will be equal to:

When moving the georadar along the profile, the x coordinate will increase, while hr will correspondingly decrease, forming the left branch of the hyperbola and reaching a minimum at the point х = х₀, where the vertex of the hyperbola will be located. At х > х₀, the values of h_r will increase and correspond to the right branch of the hyperbola. In the intermediate layer numbered i, the ray travels a distance h_ri, which is greater than the thickness of the layer h_i at all points except х₀:

where

Since in the last layer the ray travels a distance smaller than h_n, then h_rn will be equal to:

Fig. 1. Schematic diagram of the model of electromagnetic wave propagation in a layered rock mass

The ray propagation time in layer i:

The total time t_r of the ray propagation from the georadar to a local feature will be:

Then the averaged ray velocity will be equal to:

Equation of a hyperbolic time-distance curve in a homogeneous medium [2]:

For a horizontally layered medium, equation (1) will take the following form:

Based on the basic trigonometrical identity and the positivity of the cosine function in the area of arcsine values, we transform the denominator of the fraction:

Then the equation of the hyperbolic time-distance curve of GPR impulses reflected from a local feature located in a layered rock mass can be represented as (2). When substituting parameters for a single-layer medium into equation (2), it coincides with expression (1):

When processing georadar data, the hyperbolic time-distance curve having the form corresponding to expression (2) is approximated by a hyperbola having the form according to equation (1). The values of EMW propagation velocity calculated as a result of such approximation are not true, but apparent (v_app), and represent some integral value of EMW velocities in all overlying layers. To determine the dependence of v_app on the values of v in the layers overlying the local feature, we equate equations (1) and (2) to each other:

And let’s express v_app:

According to the known dependence v = c/√ε′ (с = 300 000 km/s = 0.3 m/ns) [2], the apparent dielectric permittivity will be equal to:

It is impossible to estimate the EMW propagation velocity in the rocks of a particular layer from the value of v_app, and it can be assumed that v of each layer is in the range of v_app ± Δv. For example, if in GPR measurements of a permafrost rock mass with v in the narrow range of 100–150 m/μsec (the average value v_avr ≈ 125 m/μsec) [14] v_app will be higher than 100 m/μsec, then the rock mass as a whole can be characterized as frozen. However, in such rock masses there may be a layer of rocks in thawed state, which can be detected by georadar data, but there is a problem with its identification [15]. In this connection there is a problem of determining the influence of a layer of rocks with a low value of vavr on v_app.

To solve this problem, we used formula (3) and the fact that v of frozen rocks varies within narrow limits, and the velocity of EMW propagation in unfrozen (thawed) rocks v_th is much lower [2]. Let’s represent a part of a frozen rock mass as layers of the same thickness h_avr, in each of which the EMW propagation velocity will be equal to the average value v_avr. The frozen part of the rock mass can be divided into layers arbitrarily, since at the same v their number and thickness do not affect the traveltime of GPR impulses that make up the hyperbolic lineup. The thickness h_th and velocity v_th in and thawed rock layer are represented as proportional to h_avr and v_avr:

h_th = k_hv_avr,  v_th= k_vv_avr.

After substitution into (3), we obtain an expression for the apparent velocity v_app.th of EMW propagation in a frozen rock mass containing a layer of thawed rocks:

At k_h = 0, according to formula (5), v_app.th = v_avr, i.e. it will correspond to frozen rock mass. In order to determine how v_app will change compared to v_avr in the presence of a low-velocity layer of thawed rocks, let us divide formula (3) with the frozen rock parameters havr and v_avr for all layers, denoted by v_app.fr, by expression (5):

Verification of the obtained theoretical expressions using computer simulation data

A simulation of GPR data was performed in the gprMax system [16], which has positively proved itself in studies devoted to the determination and analysis of hyperbolic lineups of GPR impulses [17–19]. The following parameters were used in the simulation: probing impulse – Ricker pulse with Fourier spectrum center frequency of 400 MHz, time-base sweep of 150, baseline of 0 mm. The data for the simulation are presented in the table below.

Input file text for gprMax for model #1:

#domain: 4 9.1 0.002
#dx_dy_dz: 0.002 0.002 0.002
#time_window: 150e-9
#material: 6 0 1 0 sloi1
#material: 4 0 1 0 sloi2
#waveform: ricker 10 0.4e9 my_ricker
#hertzian_dipole: z 0.1 9 0 my_ricker
#rx: 0.1 9 0
#src_steps: 0.01 0 0
#rx_steps: 0.01 0 0
#box: 0 0 0 4 9 0.002 sloi1
#box: 0 0 0 4 7 0.002 sloi2
#cylinder: 2 5 0 2 5 0.002 0.01 pec

Table

Parameters of rock mass models

The simulation result (out-file) was exported to the format of the GeoScan32 software program (of SPC GEOTECH manufacturer), in which the depth scale origin and the Baseline parameter equal to 1 were set. Fig. 2, a shows the results of the simulation (Model #1), the values of ε′_app, calculated using the “Hyperbola” procedure, and v_app. The traveltimes of impulses reflected from the lower boundary of layer 1 t₁ and from the local feature, the top of the hyperbola, t_h, according to formula (1) are equal to:

t₁ = 32.7 ns; t_h = t₁ + 26.7 = 59.4 ns.

According to formulas (3) and (4), v_app and ε′_app are equal to:

v_app = 0.1349 m/ns = 134.9 m/µs; ε′_app = 4.9484.

Verification of the values obtained:

Thus, the relative error of the values of v_app and ε′_app calculated in the GeoScan32 program was 0.07 and 0.37%, respectively. For model No. 2, from the calculations by formulas (3) and (4), the values of v_app = 123.7 m/μsec, ε′_app = 5.8853 were obtained.

To verify formula (6), we first calculate v_app.fr for the model of frozen rock mass with thickness h₀ = 5 m, consisting of n = 5 layers with h = 1 m, v = 122.5, 150, 122.5, 134.2, 150 m/μs and ε′ = 6; 4; 6; 5; 4. According to expression (3), v_app.fr = 134.7 m/µs.

Further, instead of the fourth layer, we introduce a low-velocity layer (model No. 3 in the Table) with such parameters h_th, v_th so that the value of v_app.th is in the range characteristic of thawed rocks. Let us set the thickness of the thawed layer h_th = 0.5 m (thus h₀ decreases to 4.5 m) as half (k_h = 0.5) of the thickness of the averaged layer h_avr = 1 m, and the EMW propagation velocity in it is two times less (k_v = 0.5) than the average value (at ε′_avr = 5) of v_avr = 134.2 m/µs, i.e. v_th = 67.1 m/µs, ε′_th = (avr/v_th)² = 20. Other parameters of the simulated rock mass are presented in the Table. According to formula (3), v_app.th = 121.2 m/μsec, which is confirmed by the result of the calculation of v_app based on the computer simulation data (Fig. 2, c), the relative error of which, compared to the exact value, was less than 0.5%.

Fig. 2. Synthetic radargrams of models Nos. 1(a), 2(b) and 3(c)

Now we can calculate how the apparent EMW propagation velocity in the rock mass model has changed with the introduction of a low-velocity layer (into the model):

That is, in the presence of a low-velocity layer with the above parameters, v_app decreases by ≈10%. In order to obtain this result, we had to carry out the whole set of calculations to calculate v_app both for a fully frozen rock mass and for the case with a layer of unfrozen (thawed) rocks. Such calculations can be substantially simplified if we use formula (6), which allows to obtain the same result with accuracy to thousandths:

To confirm that formula (6) is correct when the frozen part of a rock mass is subdivided into an arbitrary number of layers, calculations were performed for the model with 3 and 9 layers. Only one parameter, k_h, changes, for instance, at n = 3 (h = 2, 0.5, 2 m), k_h = 0.25, at n = 9 (thickness of each layer is 0.5 m), k_h = 1:

In general, to determine whether v_app.th will be within the range of the values characteristic of frozen rocks, for example, in Central Yakutia, at arbitrary k_h, k_v, we substitute into formula (5) the averaged values characteristic of this region, v_avr = 125 m/µs and v_app.th = 100 m/µs [14], then with a very significant difference between v_app.fr and v_app.th:

we obtain from (5):

from where

When v_app.th increases or v_app.fr decreases, expression (7) will be greater than 0.8, and the ratio v_app.fr/v_app.th will decrease, which will lead to difficulties in interpreting GPR data to identify the presence of a layer of thawed rocks on the basis of analyzing the hyperbolic time-distance curve of GPR impulses. The estimation of the feasibility of detecting the layer of thawed rocks based on the values of parameters k_h, k_v according to formula (7) was performed for the region of Central Yakutia. For other regions the calculation should be accomplished with the corresponding values of EMW propagation velocities in frozen and thawed rocks.

Conclusions

The performed study allowed to develop a model of hyperbolic time-distance curve of GPR impulses reflected from a local feature located in a rock mass with an arbitrary number of layers. On the basis of the developed model, the expressions for the apparent values of electromagnetic wave propagation velocity and the material part of the relative complex dielectric permittivity calculated from the hyperbolic time-distance curve of georadar impulses were obtained. The obtained expressions allowed us to determine how the velocity of electromagnetic wave propagation in a rock mass containing a layer of unfrozen (thawed) rocks decreases compared to a fully frozen rock mass. The findings of the theoretical studies were confirmed by comparison with the results of the analysis of the data of the georadar measurements computer simulation in the gprMax system (the relative error was less than 0.5%).

The results obtained in the course of the study are of great importance for the development of methodological support of GPR for determining electrophysical properties of rocks, which will increase the reliability of the assessment of their physical and mechanical properties, especially in the area of permafrost occurrence. Practical application of the obtained results in studies aimed at automated determination of electrophysical properties of rocks and soils by hyperbolic time-distance curves will allow to create a database with up-to-date information on dielectric permittivity of rocks.