# Rigorous 2D Model for Study of Pulsed and Monochromatic Waves Propagation Near the Earth’s Surface

Authors: Seil S. Sautbekov; Yuriy K. Sirenko; Nataliya P. Yashina; Aleksey A. Vertiy
Journal: International Journal of Antennas and Propagation
 (2014)

---

## Abstract

A model problem considered in the paper allows solving rather complex 2D problems of the electromagnetic wave propagation with a required accuracy using conventional personal computers. The problems are of great importance for the theory and practical applications. The association of FDTD schemes with exact absorbing conditions makes up the basis for constructing models of the kind. This approach reduces the original open initial boundary value problems to the equivalent closed problems which can be solved numerically using the standard grid methods.


---

## Body

## 1. Introduction

The free-space radio wave propagation, with application of no special guiding structures, is widely used in communication and broadcast systems. In this condition the transmission line is represented by the medium (an assemblage of natural and artificial elements and objects) filling the space through which the radio waves travel from a transmitting antenna to a receiving one. In the case of not very long landlines such a medium includes all the objects located in the near-surface strata of the atmosphere and Earth.The Earth’s atmosphere and near-surface strata of the globe represent absorbing inhomogeneous media. Their relative permittivity and specific conductivity do not remain unchanged in space and time and at times are essentially dependent on the propagating electromagnetic wavelength. Usually, these media are regarded to be nonmagnetic since their relative permeability practically does not differ from unity.The literature sources (see, e.g., [1, 2]) present a comprehensive notion concerning the effects that should be taken into account in the calculations of the optimum propagation trajectories, loss, and possible distortions of the transmitted signals. Most importantly these are (i) the effects associated with electromagnetic wave scattering by local contrast irregularities of the space through which the signals propagate, (ii) interference effects which play an important role in the case of multipath propagation, and (iii) effects provoking gradual variations in the wave propagation direction and velocity in the case of smooth changes of the electric parameters of the medium (refraction). The suggested recommendations on the consideration of these effects are quite various [1, 2]. However, none of these can guarantee reliability of the propagation factor magnitudes obtained using these methods. Recall that this is a parameter to characterize attenuation of the radio wave field strength in the case of propagation in real conditions as compared with the respective magnitude corresponding to the free-space propagation.It is understandable why the capabilities of most approaches used in practice to calculate electromagnetic fields near the Earth’s surface are limited with respect to both the accuracy of the obtained results and the scope of the analyzable situations. The reason is that too many local centers of scattering with various geometrical and electrical parameters occur on the wave propagation path and too complex can be smooth variations of the electrical parameters of the air medium for each specific relief. Under these conditions the conventional deterministic and statistical methods (as a rule, approximate, or grounded on rough approximations of the real objects by objects of a regular geometrical shape) considering the relief and atmospheric inhomogeneities [2] are capable of providing, at the best, only a qualitatively true picture of what is going on. At the same time, application of the classical numerical methods, namely, the frequency domain approaches which are based on solving integral equations, for example, would require enormous computational resources and the real errors would occur far from always to be estimated with a sufficient accuracy.2D version of the problem which has been briefly described above, we suggest to solve rigorously using the familiar advantages of the time-domain approaches [3–6], specifically, their universality (the constitutive and geometrical parameters of the wave propagation medium can be practically whatever you like), computational efficiency (all the computations are carried out within explicit schemes not requiring inversions of any operators), and the possibility of the quick and accurate conversion of the obtained numerical data into usual amplitude-frequency characteristics. To take these advantages, it is necessary to correctly reduce (bound) the computational space of the primary initial boundary value problem describing space and time transformations of the electromagnetic waves which propagate near the Earth’s surface. In other words, it is necessary to reduce the open problem to an equivalent closed one. In the paper, this rather complex procedure is performed using the method of exact absorbing conditions [6, 7], which made it possible over the recent several years to comprehensively study a number of important problems of radio physics.The models of the method of exact absorbing conditions are quite universal and the spectrum of phenomena and situations analyzable within each of these are extremely wide. Thus, in this case we will analyze a minor complication of the standard problem, namely, the employment of instantaneous field sources collected in special way (or computed in advance), enables the study of the propagation features of directed electromagnetic waves and allows to “pass through” rather long channels, cutting them to the intervals that are much less exacting to computational recourses.


## 2. Formulation of the Model Problem

All the processes, the necessity of which has been discussed above, are described by the following 2D initial boundary value problem [6, 7]:
(1)
[
-
ε
(
g
)
∂
2
∂
t
2
-
σ
(
g
)
η
0
∂
∂
t
+
∂
2
∂
y
2
+
∂
2
∂
z
2
]
U
(
g
,
t
)
=
F
(
g
,
t
)
;
hhhhhhhhhhhhhhhhhhhhhhhhhhhhhhhh
t
>
0
,
g
∈
Ω
,
U
(
g
,
t
)
|
t
=
0
=
φ
(
g
)
,
∂
U
(
g
,
t
)
∂
t
|
t
=
0
=
ψ
(
g
)
;
g
=
{
y
,
z
}
∈
Ω
-
,
E
t
g
(
q
,
t
)
,
H
t
g
(
q
,
t
)
are
continuous
when
crossing
Σ
ε
,
σ
,
E
t
g
(
q
,
t
)
|
q
=
{
x
,
y
,
z
}
∈
Σ
=
0
;
t
≥
0
,
(
∂
/
∂
x
≡
0—each section by plane x
=
const has uniform geometry; fields and sources are x-independent). Here U
(
g
,
t
)
=
E
x
(
g
,
t
) in the case of waves of the E-polarization (E
y
=
E
z
=
H
x
≡
0, F
(
g
,
t
)
≡
η
0
·
∂
j
x
/
∂
t) and U
(
g
,
t
)
=
H
x
(
g
,
t
) in the case of H-polarized waves (H
y
=
H
z
=
E
x
≡
0, F
(
g
,
t
)
=
∂
j
y
/
∂
z
-
∂
j
z
/
∂
y); E
x, H
x, E
t
g, and so forth are components of the electric (E
→) and magnetic (H
→) field vectors; j
→
(
q
,
t
) is the vector of excitation current; ε
(
g
)
≥
1 and σ
(
g
)
≥
0 are the relative permittivity and specific conductivity of a nondispersive and nonmagnetic medium of wave propagation, respectively; η
0
=
(
μ
0
/
ε
0
)
1
/
2 is the free-space impedance; ε
0 and μ
0 are the electric constant and permeability of vacuum; g
=
{
y
,
z
} is a point in space R
2; and q
=
{
x
,
y
,
z
} is a point of space R
3. All the physical values figuring in the paper are presented in the SI system of units. An exception has been made for the “time” t which represents a product of the true time and the free-space velocity of light and hence is expressed in meters.The symbolΣ
=
Σ
x
×
[
-
∞
,
∞
] denotes surfaces of the perfect conductor, while the symbol Σ
ε
,
σ
=
Σ
x
ε
,
σ
×
[
-
∞
,
∞
] is used for the surfaces on which the constitutive parameters of the wave propagation medium, that is, piecewise smooth (in the E-polarization case) or piecewise constant (in the case of the H-polarization) functions ε
(
g
) and σ
(
g
), have a break. The domain of analysis Ω represents a part of the y
0
z-plane bounded by the contours Σ
x. The supports of the functions of the current (F
(
g
,
t
)) and instantaneous (φ
(
g
)
=
U
i
(
g
,
0
) and ψ
(
g
)
=
[
∂
U
i
(
g
,
t
)
/
∂
t
]
|
t
=
0; U
i
(
g
,
t
) are the exciting wave field sources bounded in Ω. The supports of the functions ε
(
g
)
-
1 and σ
(
g
) are bounded there as well; that is, all the scattering dielectric and metal irregularities of the wave propagation medium are concentrated within a bounded part Ω
int
⁡
=
{
g
=
{
y
,
z
}
:
y
∈
(
y
1
,
y
2
)
,
z
∈
(
z
1
,
z
2
)
} of the domain Ω.It is supposed that all the data of the problem (1), specifically the functions ε
(
g
), σ
(
g
), F
(
g
,
t
), φ
(
g
), and ψ
(
g
), are such that the problem is uniquely solvable in the Sobolev space W
2
1
(
Ω
T
), with Ω
T
=
Ω
×
(
0
,
T
) and T
<
∞ being the upper bound of the observation time interval (
0
,
T
) [6].The initial boundary value problem (1) is an open one. Its domain of analysis is unbounded and goes to infinity along two spatial directions, specifically, along the y- and z-axis. It can be rigorously solved by the net methods only for small values of T which case is of no interest for practice. We suggest replacing it by the following closed problem of analysis domain Ω
int
⁡ and exact boundary condition D
[
U
(
g
,
t
)
]
|
g
∈
Γ
,
t
≥
0
=
0 (see [8, 9]) in its virtual rectangular boundary Γ
=
{
g
:
y
∈
[
y
1
,
y
2
]
for
z
=
z
j
and
z
∈
[
z
1
,
z
2
]
for
y
=
y
j
,
j
=
1,2
}:
(2)
[
-
ε
(
g
)
∂
2
∂
t
2
-
σ
(
g
)
η
0
∂
∂
t
+
∂
2
∂
y
2
+
∂
2
∂
z
2
]
U
(
g
,
t
)
=
F
(
g
,
t
)
;
hhhhhhhhhhhhhhhhhhhhhhhhhhhhhh
t
>
0
,
g
∈
Ω
int
⁡
,
U
(
g
,
t
)
|
t
=
0
=
φ
(
g
)
,
∂
U
(
g
,
t
)
∂
t
|
t
=
0
=
ψ
(
g
)
;
g
=
{
y
,
z
}
∈
Ω
-
int
⁡
,
E
t
g
(
q
,
t
)
,
H
t
g
(
q
,
t
)
are
continuous
when
crossing
Σ
ε
,
σ
,
E
t
g
(
q
,
t
)
|
q
=
{
x
,
y
,
z
}
∈
Σ
=
0
,
D
[
U
(
g
,
t
)
]
|
g
∈
Γ
=
0
;
t
≥
0
.The problems (1) and (2) are equivalent in Ω
-
int
⁡ [7, 8] as exact absorbing conditions do not distort the simulated physical processes. The outgoing from domain Ω
int
⁡ wave U
(
g
,
t
) seems to be absorbed by domain Ω
ext
=
Ω
∖
Ω
-
int
⁡ or by its boundary Γ, and there is no reflection in domain Ω
int
⁡.


## 3. Computational Procedures

The standard discretization of the closed 2D initial boundary value problem (2) by the finite-difference (FD) method [4] using a uniform rectangular mesh referred to as the coordinates g
=
{
y
,
z
} leads to an explicit computational scheme with uniquely defined mesh function U
(
p
,
s
,
m
)
≈
U
(
y
p
,
z
s
,
t
m
). The approximation error is O
(
h
-
2
); however, it could be improved, for example, using higher-order schemes. Here h
- is the mesh width in the spatial coordinates; l
- is the mesh width with respect to the time variable t; y
p
=
p
h
-, z
s
=
s
h
-, and t
m
=
m
l
-. In order to achieve desired second-order accuracy, all integrals are computed using the composite trapezoid rule and all one-sided first-order derivatives are approximated using the FD operators [10]:
(3)
d
f
(
x
)
d
x
|
x
=
x
j
≈
B
±
[
f
(
x
j
)
]
=
[
∓
3
f
(
x
j
)
±
4
f
(
x
j
±
1
)
∓
f
(
x
j
±
2
)
]
2
h
-
.The range ofp
=
P
-, P
-
+
1, P
-
+
2
,
…
,
P
+, s
=
S
-, S
-
+
1, S
-
+
2
,
…
,
S
+, and m
=
0,1
,
…
,
M integers depends on the size of the Ω
int
⁡ area and the length of the observation time interval [
0
,
T
]; namely, g
p
s
=
{
y
p
,
z
s
}
∈
Ω
-
int
⁡ and t
m
∈
[
0
,
T
]. The condition
(4)
η
2
ξ
l
-
h
-
<
1
or
/
and
4
η
l
-
h
-
<
1
;
ξ
≤
ε
-
1
(
g
)
≤
η
,
g
=
{
y
,
z
}
∈
Ω
int
⁡
,
which ensures the uniform boundedness of the approximate solution U
(
p
,
s
,
m
) with decreasing h
- and l
- (see [11, formulas (10.35) and (10.49)]), is met. Hence the FD computational scheme is stable, and the mesh function U
(
p
,
s
,
m
) tends to the solution U
(
g
p
s
,
t
m
) of the original problem [11].The solutionU
(
g
,
t
) of the problem (2) which is constructed for the points g
∈
Ω
int
⁡ and values of t from the observation time interval [
0
,
T
] using the standard computational schemes of the finite-difference method [4] and which continued, if it is necessary, from the domain Ω
int
⁡ into the domain Ω
ext by the operator method (by the method of transport operators which determine space-and-time deformations of pulses along finite segments of their propagation in regular structures [6–9]) can be converted, using the integral transformation
(5)
f
~
(
k
)
=
∫
0
T
f
(
t
)
exp
⁡
(
i
k
t
)
d
t
,
into amplitude-frequency characteristics which are required for physical analysis. These are (i) distributions of the magnitudes U
~
(
g
,
k
) of the harmonically oscillating field in the domain Ω
int
⁡ and (ii) magnitudes V
(
g
,
k
)
=
|
U
~
(
g
,
k
)
|
/
|
U
~
free
(
g
,
k
)
| (or V
(
g
,
k
)
[
d
B
]
=
20
l
g
V
(
g
,
k
)) of attenuation factor characterizing variation of electromagnetic field of radio wave while propagating in real situation in comparison with the magnitude of the same variable if radio wave propagates in free-space. Here k
=
2
π
/
λ
>
0 is the wavenumber (frequency parameter or simply frequency); λ is wavelength in the free-space; T is the upper limit of the observation time interval 0
≤
t
≤
T; and for all t
>
T the function f
(
t
) in (5) is assumed to be equal to zero.


## 4. Some Numerical Results

Consider the case ofE-polarization and the channel (the wave propagation environment), comprising natural and artificial elements, shown in Figure 1. The points, where the attenuation factor has been computed, are fixed on the roofs of buildings (g
1 and g
5), inside dry soil at the depth one meter approximately (g
2), in woodland (g
3), and at the top of the hill that is out of line-of-sight (g
4).Figure 1
Geometry of the problem: 1: atmosphere (ε
=
1.0, σ
=
0), 2: dry soil (ε
=
4.5, σ
=
1
0
-
3), 3: fresh water (ε
=
90, σ
=
2.0
·
1
0
-
2), 4: forest (ε
=
1.2, σ
=
1.0
·
1
0
-
4), and 5: brick buildings (ε
=
3.0, σ
=
1
0
-
4). g
j, j
=
1,2
,
…
,
5-observation points.Current pulse(6)
F
(
g
,
t
)
=
G
(
g
)
P
(
t
)
;
G
(
g
)
=
χ
[
2.25
-
(
z
-
150
)
2
-
(
y
-
25
)
2
]
,
P
(
t
)
=
4
sin
[
Δ
k
(
t
-
T
~
)
]
cos
⁡
[
k
~
(
t
-
T
~
)
]
(
t
-
T
~
)
-
1
χ
(
T
-
-
t
)
,
k
~
=
1.0
,
Δ
k
=
0.5
,
T
~
=
100
,
T
-
=
200
,
(no instant sources) in computational domain Ω
int
⁡, generates the signal U
(
g
,
t
), covering the frequency range 0.5
≤
k
≤
1.5. Here χ
(
⋯
) is the Heaviside step function; k
~, T
~, and T
- are the central frequency, time of delay, and pulse P
(
t
) duration correspondingly.The principal peculiarities that are characteristic of the signal propagation (6) are illustrated by data shown in Figure 2 (patterns of E
x
(
g
,
t
), g
∈
Ω
int
⁡ at different observation moments) and Figure 3 (the values of attenuation factor in observation points g
j). The result presented in the last fragment of Figure 2 (t
=
1500) can be used for examination of following interval of signal propagation channel 1250
≤
y
≤
2750.Figure 2
Spatial patternsU
(
g
,
t
), g
∈
Ω
int
⁡ taken at various moments of time t for propagation of extra wideband signal.Figure 3
The magnitudes ofV
(
g
j
,
k
) obtained in different points g
j
∈
Ω
int
⁡.For that the spatial pattern ofU
(
g
,
t
)
|
t
=
1500 and [
∂
U
(
g
,
t
)
/
∂
t
]
|
t
=
1500 in the domain 0
<
z
<
250, 1250
<
y
<
1500 has to be used for assignment of instant sources φ
(
g
) and ψ
(
g
) in the problem (2) with new computational domain Ω
int
⁡
=
{
g
:
1250
<
y
<
2750,0
<
z
<
250
} and new current sources F
(
g
,
t
)
=
0.The investigation of spatial patterns ofU
(
g
,
t
) corresponding to quasi monochromatic signals propagating in the domain Ω
int
⁡ may be rather useful.For the results presented in Figure4, the signal is generated by super narrowband pulse of current F
(
g
,
t
)
=
G
(
g
)
cos
⁡
[
k
~
(
t
-
T
~
)
]
χ
(
T
-
-
t
), k
~
=
1.0, T
~
=
1.0, and T
-
=
2000. Here one can clearly see the local centers of signal scattering and regions of signals dying, which are caused by interference of multiple rays, generated of such scattering domains that are weakly or strongly shadowed due to the relief’s irregularity, and so forth.Figure 4
Spatial patternU
(
g
,
t
), g
∈
Ω
int
⁡ taken at various fixed moments t when the quasi monochromatic signal with central frequency k
~
=
1.0 is travelling along.


## 5. Conclusion and Future Research

In the paper, a novel algorithm has been applied for rigorously solving problems which allow studying the space-and-time and space-and-frequency transformations of electromagnetic waves propagating near the Earth’s surface. The algorithm is based on the closed model problems whose domain of analysis is confined by the exact absorbing conditions which do not distort physics of the modeled processes. Its capacity is supported by computational examples for rather complicated radio channels, comprising objects of natural and artificial origin.All computer codes, which are the implementation of algorithm developed, are authentic, developed by authors of the paper, and strongly oriented to the particular problem in the focus.The patterns of field strengthU
(
g
p
s
,
t
m
) in the points g
p
s
=
{
y
p
,
z
s
} of computational domain Ω
-
int
⁡ are calculated at each time stratum t
m and are displayed at monitor; some of them are presented in Figures 2 and 4. Such patterns at any required time moments t may be saved in.bmp-format. The sequences of such files are incorporated and are in.exe-files. Their activation allows visualization of the signal propagation within all possible time intervals of interest.The results of calculations can be treated in FD if necessary and arranged within computational domainΩ
-
int
⁡ as patterns |
U
(
g
,
k
)
|
=
const, arg
⁡
U
(
g
,
k
)
=
const, or V
(
g
,
k
)
=
const for each k from frequency band, covered by excitation pulse.The algorithms and computer codes developed open the opportunity for rather wide physical problems:(i)
the study of over-the-horizon wave propagation above both ground and sea surface;
(ii)
the study of tropospheric waveguides appearing due to vertical distortion in near surface atmosphere layer;
(iii)
the solution of the problems of subterranean and undersea surface communications.In simulation of waves and signal propagation the various types of souses may be used; besides the current type of souses generation omnidirectional radiation (see, e.g., (6)), the sources (current and instant) relevant to real life scenario and sources with given directivity of radiation field may be also considered.Natural limitations of this approach are basically connected with possible simulation errors (within 2D-models), with wave dimensions of computational domain and wave dimensions of particular scattering objects, filling the domain.The algorithms and computer codes presented in the paper, similar to other algorithms of the exact absorbing condition method [6–8], are thoroughly tested. The principal tasks for computer codes generalization and unification that we are about to resolve are the following: (i) the automation of the routine for trace “prolongation” by means of introduction of the new instant sources φ
(
g
) and ψ
(
g
) in the problem (2) with new computational domain Ω
int
⁡ and new current sources F
(
g
,
t
)
=
0; (ii) the construction of more efficient computational schemas for the solution to the problem (2) that can essentially reduce required computer resources (the usage of acceleration scheme based on the blocked fast Fourier transform [12] for calculation of the temporal convolutions in exact absorbing conditions of the problem (2), changing from computational schemes of the finite-difference method to schemas of the discontinuous Galerkin finite-element method [13]); (iii) the straightforward comparison (basing on criteria of universality, accuracy, and consumption of computer resources) of our approach with approaches of other authors, with those ones which are already considered to be classic [1, 2, 14, 15] and with relatively novel [16–24].As we can see the implementation of the task (iii) has to start with benchmark analysis of the results of (i) integral equation method or method of approximate boundary conditions (elementary scatterers, plane or spherical electrically inhomogeneous Earth’s surface, and homogeneous atmosphere), (ii) methods realizing approximation of small and sloping discontinuities (uneven or rough Earth’s surface and atmosphere), and (iv) perturbation method (turbulent atmosphere).The listed above method provides the possibility to simulate only idealized model traces or particular fragments of real traces of signal propagation [2, 14, 15]. But these methods are rather useful for description of various local effects (interference signal attenuation, obstacle gain, beyond-the-horizon propagation, rays diffusion, etc.), which contribution is still necessary to summarize correctly (that is rather a complicated problem) in order to obtain practically significant result.It is well known that field prediction models are fundamental tools for design, planning, and optimization of stationary and mobile radio systems. This is particularly true for arrangement of steady connection on the landscape with complicated relief and in an urban environment, where radio waves strongly interact with the multiple artificial and natural objects changing their amplitudes and direction of propagation and where the coverage prediction is no longer easy to perform. That is why the variety of new models and modifications of well-known models, enabling study and analysis of the process of wave propagation along traces of different complexity and filling, is permanently discussed in various scientific journals. Several of these models are planned to be implemented in research works of master and Ph.D. students at the department of radioengineering, electronics, and telecommunications of physic-technical faculty in L.N. Gumilyov Eurasian National University. It is planned to implement and to carry out the comparative analysis of various approaches (including approach briefly described in present paper) to the solution of fundamental and applied problems of the radio waves propagation theory. More likely we will consider approaches, basing on various models of the ray-tracing method [16–18], of the method utilizing the formalism of the integral geometry [19], of a novel dual-field time-domain finite-element domain-decomposition method [20], of the efficient integral and parabolic equation methods [21–23], and of Kirchhoff integral method [24].



---

*Source: 158297-2014-08-07.xml*