# Long-Range Three-Dimensional Ground Wave Propagation Modeling Over Flat, Irregular Terrain

J. Vincent, P. Borderies, Jean-René Poirier, V. Gobin

► To cite this version:

J. Vincent, P. Borderies, Jean-René Poirier, V. Gobin. Long-Range Three-Dimensional Ground Wave Propagation Modeling Over Flat, Irregular Terrain. IEEE Transactions on Antennas and Propagation, 2016, 64 (5), pp.1900-1906. 10.1109/TAP.2016.2540647. hal-01721742

HAL Id: hal-01721742

https://hal.science/hal-01721742

Submitted on 15 Oct 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Long-Range Three-Dimensional Ground Wave Propagation Modeling Over Flat, Irregular Terrain

J. Vincent, P. Borderies, J. R. Poirier, and V. Gobin

*Abstract*—Ground wave propagation of low-frequency (LF) electromagnetic waves is well known for the canonical case of a flat, azimuthally homogeneous, soil. Moreover, it is well investigated for realistic propagation problems, some of which include Earth's curvature, irregular and lossy terrain profiles, and mixed paths. In this paper, kinds of environment are considered in the far-field of the radiating antenna with an original methodology. It rests on the use of hybridization of finite-difference in time-domain (FDTD) with the direct numerical integration of Sommerfeld-type integrals which is also presented and validated. This hybridization uses domain decomposition and Huygens' surfaces in a stratified dielectric medium in three-dimensional (3-D) space. The approach is validated with respect to a reference solution both in the canonical flat ground case and in a hilly one.

*Index Terms*—Ground wave propagation, hybridization, low-frequency wave, problem of Sommerfeld, three-dimensional (3-D) FDTD.

## I. INTRODUCTION

USUALLY the electromagnetic problem of low-frequency (LF) radiating antennas is decomposed into two parts: 1) the near-field radiation emitted by the antenna and 2) the propagation of this near-field in the far range for the ground wave [1]. LF electromagnetic field propagation over the ground is well known for the canonical case of a flat homogeneous soil [1], [2], and in the case of a stratified one [3] through approximate analytical approaches which are widely used. An appropriate numerical method for solving the propagation problem over long distances is the solution of parabolic equations (PE) above stratified media [4], [5]. Three-dimensional (3-D) vector PE approach exists [6], but PE method is mostly used in two dimensions [7] and it is not a full-wave method [8]. Perturbation methods using the compensation theorem treat typical forms of radial heterogeneity such as mixed paths and irregular ground [9], [10].

However, these techniques do not take into account the 3-D topography near the propagation path. Then, there is still a need of rigorous numerical computation of the radiated field when the distant irregularity does not correspond to the domain of validity of the usual methods or when there is a need of a reference solution for the interaction of the ground wave with a very general discontinuity.

Finite-difference in time-domain (FDTD) method is a well-known method used for computational electromagnetics [11]. Among its advantages is the ability to treat very general geometries involving heterogeneous complex media in volume and surface [12]. In [13], a time-domain wave propagator based on two-dimensional (2-D) FDTD computes the electromagnetic field propagation over Earth's surface. In [14] and [15], bi-dimensional FDTD is applied to obtain the electric field both in the near and the far-field over an irregular terrain, with the constraint of a radial cylindrical geometry. FDTD in three dimensions has been applied to study the near-field radiation of wire antennas in a complex environment [16]. However, FDTD is not applicable in three dimensions at large distance because of the huge amount of computational resources which are required and also because of the numerical problems arising at very long distances (numerical dispersion, precision, boundary conditions, etc.).

In this paper, an original hybrid approach using the FDTD method is proposed. It rests on the use of the direct numerical integration of Sommerfeld's integrals, which is also presented, and its hybridization with 3-D FDTD applied to the surroundings of the irregularity under consideration, yielding an approach free of the former numerical problems. This hybridization technique aims at treating complex kinds of environment located in the far-field.

Section II is dedicated to the implementation of the hybridization in the 3-D FDTD method. First the Sommerfeld's half-space problem is solved with a numerical integration along the positive real axis using an adaptive algorithm called NISP (numerical integration of Sommerfeld's problem). Then, the new hybrid NISP—3-D FDTD method is presented. The validation and illustration of the method are presented in Section III.

## II. PRINCIPLES OF HYBRIDIZATION

### A. Sommerfeld Half-Space Problem

In [17], Sommerfeld computed the radiation of an infinitesimal vertical electric current element $I = I_0 dlz$ located at a height $h$ above a homogeneous dielectric half-space. The XOY plane separates the two media, the medium 1 is free space in

J. Vincent was with ONERA, 31000 Toulouse, France. He is now with the Université de Toulouse; INPT, UPS, CNRS; LAPLACE; ENSEEIHT, 31000 Toulouse, France (e-mail: vincent@laplace.univ-tlse.fr).

P. Borderies and V. Gobin are with ONERA, 31000 Toulouse, France (e-mail: pierre borderies@onera.fr; vincent.gobin@onera.fr).

J. R. Poirier is with the Université de Toulouse; INPT, UPS, CNRS; LAPLACE; ENSEEIHT, 31000 Toulouse, France (e-mail: poirier@laplace.univ-tlse.fr).
<figure><img src="image_2.png" /><figcaption>Fig. 1. Vertical z-directed current element $I_0 dl$ at height $h$ above a planar lossy ground with electrical parameters ($\mu_2$, $\epsilon_2$, $\sigma_2$).</figcaption></figure>

the upper half-space with constants $\mu_0$, $\epsilon_0$ and soil is the medium 2 in the lower half-space with constants $\mu_2 = \mu_0$, $\epsilon_2$, and $\sigma_2$. The geometry is presented in Fig. 1.

Using Hertz vector potentials and the boundary conditions at height $z = 0$, one obtains in (1) the exact formulas of all components of the electromagnetic field **E** and **H** in both media 1 and 2 [18]. The wavenumbers in free space $k_0$ and soil $k_2$ are given by $k_i = [-j\omega\mu_i(\sigma_i + j\omega\epsilon_i)]^{1/2}$. [See (1), shown at the bottom of the page.]

The free space components in (1) are expressed as the sum of three terms: the first represents a wave traveling directly from the current element $I_0 dl$; the second can be seen as radiated from a negative image current in the ground at $-h$; and the third is related to a surface wave term. The distance between the image current at distance $h$ below the plane and the observation point is $R' = [r^2 + (z+h)^2]^{1/2}$. $J_0$ and $J_1$ are Bessel functions of the first kind for orders 0 and 1, respectively [19].

## B. Numerical Integration

The integrands in all components in both media are the product of three factors. The first element is the quotient whose denominator is $(k_0^2(\lambda^2 - k_2^2)^{1/2} + k_2^2(\lambda^2 - k_0^2)^{1/2})$. An oscillating term depending on the distance $r$ of the observation point is created by the spherical Bessel functions $J_0$ and $J_1$. Finally, there is an exponential attenuation factor depending on the height of the dipole $h$ and the observation point $z$. The system has a pair of complex poles given by $s^2 = k_0^2 k_2^2 / (k_2^2 + k_0^2)$. There is no singularity on the integration path and the integrand is convergent but very oscillating along the positive real axis. It can be computed by a numerical integration with a small error. To determine the best strategy, the behavior of the integrand functions is detailed. In the following, the z-directed electric field integral in (2) is only studied but the approach is identical for the other five components of the electromagnetic field

$$ 
\begin{aligned}
N(r, z) &= \int_{0}^{+\infty} n(\lambda, r, z) d\lambda \\
&= \int_{0}^{+\infty} J_{0}(\lambda r) \frac{2k_{2}^{2}\lambda^{3}d\lambda e^{-(z+h)\sqrt{\lambda^{2}-k_{0}^{2}}}}{k_{2}^{2}\sqrt{\lambda^{2}-k_{0}^{2}}+k_{0}^{2}\sqrt{\lambda^{2}-k_{2}^{2}}}. \quad (2)
\end{aligned}
$$

For a current element located at $h = 25$ m above a medium wet soil ($\epsilon_r = 15$, $\sigma = 10^{-2} \text{ S} \cdot \text{m}^{-1}$) and an observation point located at $z = 1$ m and $r = 1$ km, the real and imaginary values of the integrand are plotted in Fig. 2 as a function of the integration parameter $\lambda$. The four main parts of the integrand are bounded by the vertical lines.

$$ 
\begin{aligned}
H_{\phi}^{2}(r, z) &= \frac{k_{2}^{2}I_{0}dl}{2\pi} \int_{0}^{+\infty} J_{1}(\lambda r) \frac{\lambda^{2}d\lambda e^{z\sqrt{\lambda^{2}-k_{2}^{2}}-h}\sqrt{\lambda^{2}-k_{0}^{2}}}{k_{2}^{2}\sqrt{\lambda^{2}-k_{0}^{2}}+k_{0}^{2}\sqrt{\lambda^{2}-k_{2}^{2}}} \\
E_{\rho}^{2}(r, z) &= \frac{j\omega\mu_{0}I_{0}dl}{2\pi} \int_{0}^{+\infty} J_{1}(\lambda r) \frac{\lambda^{2}d\lambda e^{z\sqrt{\lambda^{2}-k_{2}^{2}}-h}\sqrt{\lambda^{2}-k_{0}^{2}}}{k_{2}^{2}\sqrt{\lambda^{2}-k_{0}^{2}}+k_{0}^{2}\sqrt{\lambda^{2}-k_{2}^{2}}} \\
E_{z}^{2}(r, z) &= -\frac{j\omega\mu_{0}I_{0}dl}{2\pi} \int_{0}^{+\infty} J_{0}(\lambda r) \frac{\lambda^{3}d\lambda e^{z\sqrt{\lambda^{2}-k_{2}^{2}}-h}\sqrt{\lambda^{2}-k_{0}^{2}}}{k_{2}^{2}\sqrt{\lambda^{2}-k_{0}^{2}}+k_{0}^{2}\sqrt{\lambda^{2}-k_{2}^{2}}}
\end{aligned}
\quad (1) $$

$$ 
\begin{aligned}
H_{\phi}^{1}(r, z) &= \frac{I_{0}dl}{4\pi} \left[ -\frac{\partial}{\partial r} \frac{e^{-jk_{0}R}}{R} + \frac{\partial}{\partial r} \frac{e^{-jk_{0}R'}}{R'} + \int_{0}^{+\infty} J_{1}(\lambda r) \frac{2k_{2}^{2}\lambda^{2}d\lambda e^{-(z+h)\sqrt{\lambda^{2}-k_{0}^{2}}}}{k_{2}^{2}\sqrt{\lambda^{2}-k_{0}^{2}}+k_{0}^{2}\sqrt{\lambda^{2}-k_{2}^{2}}} \right] \\
E_{\rho}^{1}(r, z) &= -\frac{j\omega\mu_{0}I_{0}dl}{4\pi K_{0}^{2}} \left[ \frac{\partial^{2}}{\partial z \partial r} \frac{e^{-jk_{0}R}}{R} - \frac{\partial^{2}}{\partial z \partial r} \frac{e^{-jk_{0}R'}}{R'} + \int_{0}^{+\infty} J_{1}(\lambda r) \frac{2k_{2}^{2}\lambda^{2}d\lambda e^{-(z+h)\sqrt{\lambda^{2}-k_{0}^{2}}}}{k_{2}^{2}\sqrt{\lambda^{2}-k_{0}^{2}}+k_{0}^{2}\sqrt{\lambda^{2}-k_{2}^{2}}} \right] \\
E_{\rho}^{1}(r, z) &= -\frac{j\omega\mu_{0}I_{0}dl}{4\pi K_{0}^{2}} \left[ \left( \frac{\partial^{2}}{\partial Z^{2}} + k_{0}^{2} \right) \frac{e^{-jk_{0}R}}{R} - \left( \frac{\partial^{2}}{\partial Z^{2}} + k_{0}^{2} \right) \frac{e^{-jk_{0}R'}}{R'} + \int_{0}^{+\infty} J_{0}(\lambda r) \frac{2k_{2}^{2}\lambda^{3}d\lambda e^{-(z+h)\sqrt{\lambda^{2}-k_{0}^{2}}}}{k_{2}^{2}\sqrt{\lambda^{2}-k_{0}^{2}}+k_{0}^{2}\sqrt{\lambda^{2}-k_{2}^{2}}} \right].
\end{aligned}
$$
<figure><img src="image_3.png" /><figcaption>Fig. 2. Real (top) and imaginary (bottom) parts of the integrand in (2) as a function of $λ$, with a medium wet soil, $h = 25$ m, $z = 1$ m, and $r = 1$ km. The four main parts are bounded by the three vertical lines.</figcaption></figure>

The first interval value is small; the function is regular and has low oscillations. The second is the most important contributor because the function is less regular. In fact, it contains a singularity which corresponds to the positive real part of the complex pole s. The third part is longer to evaluate than the previous two because both real and imaginary values of the integrand oscillate due to the Bessel function. The last part is dominated by the exponential attenuation term and may not be calculated because it is negligible despite the oscillating behavior. These intervals have different properties and numerical integration does not have the same convergence behavior. Integration functions are available in many software products; however, changes in the integration scheme are complicated; it is therefore difficult to implement an optimal strategy.

Indeed the previous analysis shows that different low- and high-order integration methods have to be used for each part of the integrand. Third degree Newton–Cotes method is chosen for parts one and four integrations and Gauss–Legendre quadrature with more than ten points is chosen for the integrations of parts two and three.

An adaptive algorithm for numerical integration is realized to ensure a good compromise between computational time and numerical accuracy [20]. The algorithm performs the integration independently on the four intervals with minimal discretization $Δ$ min. At each iteration, the step is halved for a new integration until the fixed convergence or the maximum number of iterations is reached. The four integration intervals are heuristically defined as follows:

$$ 
\begin{aligned}
I_1 &= [0; 0.9 \cdot \Re e(s)]; I_2 = [0.9 \cdot \Re e(s); 1.1 \cdot \Re e(s)] \\
I_3 &= [1.1 \cdot \Re e(s); \lambda_{\eta}]; I_4 = [\lambda_{\eta}; 3 \cdot \lambda_{\eta}]
\end{aligned} 
\quad (3) 
$$

where $\Re e(s)$ is the real part of complex pole s and $\lambda_{\eta}$ is determined from 1% of the maximal value of the integrand and the

<table><caption>TABLE I
DIELECTRIC CONSTANTS AND CONDUCTIVITY FOR DIFFERENT SOILS</caption><thead><tr><th>Soils</th><th>Dry</th><th>Medium Wet</th><th>Wet</th></tr></thead><tbody><tr><td>ε<sub>r</sub></td><td>5</td><td>15</td><td>30</td></tr><tr><td>σ (S/m)</td><td>10<sup>-4</sup></td><td>10<sup>-3</sup></td><td>10<sup>-2</sup></td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Comparison of the vertical electric field component in Sommerfeld’s problem solved with two different methods: NISP and numerical approximation of Norton. A $dl = 25$ m long current element at frequency $f = 100$ kHz is at height $h = 25$ m above different soils, the magnitude of the vertical electric field is computed a long distance $r$ at height $z = 1$ m.</figcaption></figure>

exponential attenuation factor given in the following:

$$ k_0 \le \lambda_{\eta} = \frac{1}{(z+h)} \cdot \ln \left( \frac{100}{\max\{n(\lambda, r, z)\}} \right). \quad (4) $$

The algorithm NISP has been tested for vertical component $E_z$ radiated at 100 kHz by a $dl = 25$ m long current element located at $h = 25$ m over different soils, along the direction of propagation $r$ at height $z = 1$ m. The soil is characterized with the dielectric constant $ε_r$ and the conductivity $σ$, which depend on the humidity, see Table I, and can be found in [21].

The minimal discretization of NISP is set to $Δ$ min = $2π/(10r)$ because of the $2π$-periodic Bessel function. The convergence is computed with an absolute relative error, its value is fixed to $cv = 10^{-3}$, and the maximal number of iteration is max = 4. Fig. 3 shows the comparison between two different methods: the numerical integration and compared with the analytical formulas given by Norton [2]. We can see that there is a good agreement between these two approaches on the upper half-space.

Consequently, NISP may be used to compute the radiated field whatever the distance beyond a wavelength. It has to be noted that it yields the radiated field also in the ground: it is then relevant to initialize the hybrid method presented in the next section.

### C. Hybrid NISP—3-D FDTD Method

The main objective of this section is to present a method able to compute numerically the effects of a perturbation in/on the
<figure><img src="image_5.png" /><figcaption>Fig. 4. SF/TF representation for 1-D FDTD. The Huygens' surface separates the computational space into regions: a region where the total field is considered and the scattered field is considered in the other one.</figcaption></figure>

soil at large distance from the radiating antenna. It rests on the hybridization of the approach presented in Section II-B and the 3-D FDTD algorithm. The scattered-field/total-field (SF/TF) formulation has been first proposed in [22]: along a Huygens' surface surrounding the scatterer, the incident plane wave is introduced on the six faces with the goal of dividing the computational domain into a total-field zone and a scattered-field one. In [23] and [24], still in free space, the incident plane wave is substituted by the field radiated by a remote source, leading to multidomain approaches. In [25], again with an incident plane wave, free space is substituted by a stratified medium. As a simple example, the SF/TF one-dimensional (1-D) formalism is represented in Fig. 4. Considering the total electromagnetic field as a sum of the incident and the scattered electromagnetic field (5), one can write (6) and (7): the nth time step of the FDTD scheme applied on the whole Huygens' surface located at grid point $i_S$ in a dielectric medium ($\mu_0, \epsilon = \epsilon_r \epsilon_0, \sigma$) with space step $\Delta x$ and time step $\Delta t$. All components of the electromagnetic field on the whole 3-D Huygens' surface are easily deduced from (5)–(7), shown at the bottom of the page.

In this paper, the incident field is the ground wave created by the remote source which is computed by NISP and the computation is done in a medium which is stratified.

The computational domain includes soil and free space above, both of them bounded by adapted absorbing boundaries. NISP algorithm allows to identify the values of the electromagnetic field at any point in space in the frequency domain: it is used to compute at all times the incident field both in free space and in the soil along the six faces of the Huygens' surface with a sine excitation at the frequency of interest. It is

<figure><img src="image_6.png" /><figcaption>Fig. 5. NISP—3-D FDTD Hybridization: a vertical z-directed current element source $I_0 dl$ creates the electric and magnetic fields ($E_{inc}, H_{inc}$) which are incorporated as incident wave sources in the SF/TF FDTD formulation.</figcaption></figure>

expected that when the soil is flat, the total field is the incident one: the unperturbed ground wave is retrieved. When there is some irregularity either in the soil composition or on the relief (Fig. 5), then a scattered field is generated.

### III. RESULTS

#### A. Validation of the Method

In order to validate the method, various scenarios were implemented. First the hybridization technique is used with a homogeneous flat ground (medium wet) and compared with the analytical formulas given by Norton [2]. The antenna radiating at $f = 100$ kHz is $dl = 25$ m long and at height $h = 25$ m. The electromagnetic field is computed at height $z = 76$ m along the distance $r$ from 1 to 100 km with several reduced domains. The computation space for each hybridization implementation is a 100 cells side cube with the mesh size $\Delta x = \Delta y = \Delta z = 75$ m and the ground surface ($z = 0$ m) is located at the middle of the z-directed cells. Thus the total computational volume is $7500 \times 7500 \times 2500$ m. The volume inside the Huygens' surfaces is a 60 cell side cube centered in the domain which means that $60 \times 60 \times 30$ cells mesh the ground inside the volume. Uniaxial-PML is used to terminate the 3-D computational domain [26] and the NISP parameters are set to $\Delta \min = 2\pi/(10r)$ and $\max = 4$. Fig. 6 shows the results of the electromagnetic field computations. It validates the method when applied to an infinite flat ground since it agrees with

$$E_{TOT} = E_{INC} + E_{SCAT}, \quad H_{TOT} = H_{INC} + H_{SCAT} \qquad (5)$$

$$E_{TOT} \begin{vmatrix} n+1 \\ i_S \end{vmatrix} = \frac{2\epsilon - \sigma\Delta t}{2\epsilon + \sigma\Delta t} E_{TOT} \begin{vmatrix} n \\ i_S \end{vmatrix} + \frac{2\Delta t}{(2\epsilon + \sigma\Delta t) \Delta x} \left( H_{TOT} \begin{vmatrix} n + \frac{1}{2} \\ i_S + \frac{1}{2} \end{vmatrix} - H_{SCAT} \begin{vmatrix} n + \frac{1}{2} \\ i_S - \frac{1}{2} \end{vmatrix} \right) - \frac{2\Delta t}{(2\epsilon + \sigma\Delta t) \Delta x} H_{INC} \begin{vmatrix} n + \frac{1}{2} \\ i_S - \frac{1}{2} \end{vmatrix} \qquad (6)$$

$$H_{SCAT} \begin{vmatrix} n + \frac{1}{2} \\ i_S - \frac{1}{2} \end{vmatrix} = H_{SCAT} \begin{vmatrix} n - \frac{1}{2} \\ i_S - \frac{1}{2} \end{vmatrix} + \frac{\Delta t}{\mu_0 \Delta x} \left( E_{TOT} \begin{vmatrix} n \\ i_S \end{vmatrix} - E_{SCAT} \begin{vmatrix} n \\ i_S - 1 \end{vmatrix} \right) - \frac{\Delta t}{\mu_0 \Delta x} E_{INC} \begin{vmatrix} n \\ i_S \end{vmatrix} \qquad (7)$$
<figure><img src="image_7.png" /><figcaption>Fig. 6. Solution comparison between the NISP—3-D FDTD hybridization and Norton’s analytic formulas. The electromagnetic field along the distance r from 1 to 100 km at height z = 75 m is emitted by a 25-m long vertical current element located at height h = 25 m along distance r at height z = 1 m.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. 3-D representation of the computed hill with the hybridization. A stack of 10 dielectric blocks with different widths are located on the flat ground.</figcaption></figure>

the reference method for a wide range of distances to the source.

Then the hybridization technique in the vicinity of the source is compared to 3-D FDTD method including the source when both methods are relevant (FDTD is here the reference solution since the distance of observation is short enough). A 250-m high hill is located at 3.4 km of the infinitesimal antenna previously described, over a medium wet soil. It is modeled with a stack of 10 dielectric blocks which are all 25 m high (Fig. 7). Horizontally, the largest block is 1.5 km × 3 km; the second is 1.35 km × 3 km, up to the top which is 0.15 km × 3 km. For the FDTD, the current element and the hill are computed in a 240 cells side cube with Δx = Δy = 75 m and Δz = 25 m. The previous monopole antenna is centered inside the computed volume. The same geometry of the hill and the space steps are used with a 100 cells side cube for the hybridization. The volume inside the Huygens’ surfaces is still a 60 cell side cube centered in the domain. The normalized magnitude of the electric field at constant elevation 175 m over the soil and the hill is displayed in Fig. 8. Both methods are in agreement better than 0.3 dB and a relative error of 3.39%, which demonstrates the accuracy of the hybridization technique when an obstacle is present.

<figure><img src="image_9.png" /><figcaption>Fig. 8. Solution comparison between the NISP—3-D FDTD hybridization and the 3-D FDTD in this vicinity of the source and at a constant elevation 175 m over the soil. A hill is located at 3.4 km, electric field values are plotted in dBV·m⁻¹ unit.</figcaption></figure>

## B. Study of Ground Wave Propagation With Obstacles

To complete the study of obstacles along the ground wave propagation path, several natural elements are remotely located from the current element. The FDTD parameters, the source, the ground, the distance, and the size of the Huygens’ surface are the same in all these computations. The current element at f = 100 kHz is dl = 25 m long and at height h = 25 m. The FDTD mesh for each hybridization implementation is a 100 cells side cube with the mesh size Δx = Δy = 75 m and Δz = 25 m; the ground surface (z = 0 m) is located at the middle of the z-directed cells. The volume inside the Huygens’ surface is a 60 cells side cube and is located at 23 km from the source. Since the volume bounded by the Huygens’ surfaces is the same for all the following studies, the incident electromagnetic field computed by NISP is stored in a file. Thus, only the nature of the obstacle changes between each computation which will save time. The chosen obstacles in this section are a hill, a valley, a lake, and a forest.

First, the geometry of the hill previously used for the validation is translated along x-axis (Fig. 7). Thus the total electric field around the same hill, remotely located at 24.5 km from the current element is computed in the XOZ vertical plane at y = 0 (Fig. 9) and XOY horizontal plane at z = 25 m (Fig. 10). Then these values of the total electric field in the air are compared to the case of a flat homogeneous soil with (8) applied in the XOZ vertical plane at y = 0 and z > 0

$$ \eta = |E_{\text{natural element}}|_{\text{dBV} \cdot \text{m}^{-1}} - |E_{\text{flat soil}}|_{\text{dBV} \cdot \text{m}^{-1}} \quad (8) $$

Results are shown in Fig. 11, in which a field enhancement brought about by the natural element is determined by the positive values of η and conversely a decrease by a negative value. On the top of the hill, a 3-dB field enhancement is observed and some interferences appear in front of it (1.5 dB). Further than the obstacle, after a shaded area there is no visible effect beyond the distance of 27 km as η ≈ 0 dB. Now other natural elements such as valley, a lake, and a forest are computed at the same distance of 23 km from the current element. The
<figure><img src="image_10.png" /><figcaption>Fig. 9. Magnitude of the electric field, computed in the XOZ vertical plane, with the hybridization at 23 km of the antenna. A hill is located at 24.5 km and the values are plotted in dBV·m⁻¹.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10. Magnitude of the electric field, computed in the XOY horizontal plane at z = 25 m, with the hybridization at 23 km of the antenna. A hill is located at 24.5 km and the values are plotted in dBV·m⁻¹.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. Difference between the value of the electric field in the air with the hill and the homogeneous flat soil. It is computed with (5) in the XOZ vertical plane at z > 0 at 23 km of the antenna. The natural element is located at 24.5 km and the values of η are plotted in dB.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12. Difference between the value of the electric field in the air with the valley and the homogeneous flat soil. It is computed with (5) in the XOZ vertical plane at z > 0 at 23 km of the antenna. The natural element is located at 24.5 km and the values of η are plotted in dB.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 13. Difference between the values of the electric field in the air with a lake and the homogeneous flat soil. They are computed with (5) in the XOZ vertical plane at z > 0 at 23 km of the antenna. The natural element is located at 24.5 km and the values of η are plotted in dB.</figcaption></figure>

dimensions of the valley are the same as those of previous hill, the dielectric constants and the z-axis are reversed in order to build the hollow in the ground. In Fig. 12, contrasts are opposite to the previous case. The lake is a square-shaped dielectric block with dielectric constants ($ε_r$ = 80; $σ$ = $5 \cdot 10^{-2} S·m^{-1}$). It is located at the surface of the ground and its dimensions are 1.5-km large and 100-m deep (Fig. 13).

The forest is a 25-m height square-shaped dielectric block with dielectric constants [27] ($ε_r$ = 1.065; $σ$ = $10^{-3} S·m^{-1}$). It is located at the surface of the ground and its dimensions are 1.5-km large. Edge enhancements are presented in Fig. 14.

The results shown do not exceed 3 dB of contrast in all cases. One can conclude from these observations that all natural elements have an effect on the electromagnetic field which is in their vicinity: from 0.3 dB with the lake, 1 dB with the forest, to 3 dB with the hill and the valley. A little further on the propagation path of the ground wave, η approaches the zero value for all natural elements, the electric field is the same as the case of a flat homogeneous half-plane. Therefore, the presence
<figure><img src="image_15.png" /><figcaption>Fig. 14. Difference between the value of the electric field in the air with the forest and the homogeneous flat soil. It is computed with (5) in the XOZ vertical plane at $z > 0$ at 23 km of the antenna. The natural element is located at 24.5 km and the values of $\eta$ are plotted in dB.</figcaption></figure>

of these obstacles does not affect the long-range ground wave propagation.

## IV. CONCLUSION

In this paper, a method to study the propagation of ground wave at low frequencies with the presence of an irregular half-space has been studied with a 3-D approach.

First, the NISP algorithm was proposed to compute the electromagnetic field emitted by a vertical current element over an infinite flat ground and validated for the relevant values of soils parameters.

Then, an hybrid NISP—3-D FDTD method was introduced to compute the field in the vicinity of remote irregularities and validated both in the case of an infinite half-space and in the presence of a hill. This study includes the exploitation of this technique to treat various kinds of ground irregularities with emphasis on the three dimensions disturbance effects. These irregularities are computed over simplified models, located at 24.5 km from the current element and computed with the hybridization technique. The obtained values of the total electric field are compared to the case of a flat homogeneous soil. The observed natural elements create a field enhancement up to 3 dB in its vicinity but they do not affect the ground wave propagation beyond them.

## REFERENCES

[1] N. DeMinco, "Propagation prediction techniques and antenna modeling (150 to 1705 kHz) for Intelligent Transportation Systems (ITS) broadcast applications," *IEEE Trans. Antennas Propag.*, vol. 42, no. 4, pp. 9-34, Aug. 2000.

[2] K. A. Norton, "The propagation of radio waves over the surface of the earth and in the upper atmosphere," *Proc. IRE*, vol. 25, no. 9, pp. 1203-1236, Sep. 1937.

[3] J. Wait, *Electromagnetic Waves in Stratified Media*. New York, NY, USA: Pergamon, 1962.

[4] M. Levy, *Parabolic Equation Methods for Electromagnetic Wave Propagation*. London, U.K.: The Institution of Electrical Engineers, 2000.

[5] O. Ozgun, G. Apaydin, M. Kuzuoglu, and L. Sevgi, "PETOOL: MATLAB-based one-way and two-way split-step parabolic equation tool for radiowave propagation over variable terrain," *Comput. Phys. Commun.*, vol. 182, no. 12, pp. 2638-2654, 2011.

[6] R. Janaswamy, "Path loss prediction in the presence of buildings on flat terrain: A 3D vector parabolic equation approach," *IEEE Trans. Antennas Propag.*, vol. 51, no. 8, pp. 1716-1728, Aug. 2003.

[7] F. Akleman and L. Sevgi, "A novel MoM and SSEP-based groundwave propagation field-strength prediction simulator," *IEEE Trans. Antennas Propag. Mag.*, vol. 49, no. 5, pp. 69-82, Oct. 2007.

[8] V. Fabbro, P. F. Combes, and N. Guillet, "Apparent radar cross section of a large target illuminated by a surface wave above the sea," *Progr. Electromagn. Res.*, vol. 50, pp. 41-60, 2005.

[9] T. S. M. Maclean and Z. Wu, *Radiowave Propagation Over Ground*. London, U.K.: Chapman & Hall, 1993.

[10] G. Apayadin and L. Sevgi, "MATLAB-based FEM—Parabolic equation tool for path-loss calculations along multi-mixed-terrain paths," *IEEE Trans. Antennas Propag. Mag.*, vol. 56, no. 3, pp. 221-236, Jun. 2014.

[11] A. Taflove and S. Hagness, *The Finite-Difference Time Domain Method*, 3rd ed. Norwood, MA, USA: Artech House, 2005.

[12] H. Chanal, J. P. Segaud, P. Borderies, and M. Saillard, "Homogenization and scattering from heterogeneous media based on finite-difference-time-domain Monte Carlo computations," *J. Opt. Soc. Amer. A*, vol. 23, pp. 370-381, 2006.

[13] F. Akleman and L. Sevgi, "Realistic surface modeling for a finite-difference time-domain wave propagator," *IEEE Trans. Antennas Propag.*, vol. 51, no. 7, pp. 1675-1679, Jul. 2003.

[14] L. L. Zhou, X. L. Xi, J. F. Liu, and N. M. Yu, "LF ground-wave propagation over irregular terrain," *IEEE Trans. Antennas Propag.*, vol. 59, no. 4, pp. 1254-1260, Apr. 2011.

[15] L. L. Zhou, X. L. Xi, and Y. X. Du, "Application in low-frequency ground-wave propagation of parallel FDTD based on GPU," in *Proc. 10th Int. Symp. Antennas Propag. EM Theory (ISAPE)*, Oct. 2012, pp. 917-920.

[16] J. Vincent, P. Borderies, J. R. Poirier, and V. Gobin, "Simulations and effects of natural environments on low frequency antennas with 3D FDTD method," *Progr. Electromagn. Res. M*, vol. 38, pp. 45-52, Aug. 2014.

[17] A. Sommerfeld, "Propagation of waves in wireless telegraphy," *Annalen der Physik*, vol. 28, pp. 665-736, Mar. 1909.

[18] L. Robin and P. Poincelot, "Rayonnement d'un Dipôle Électrique en Présence de Deux Milieux Séparés par un Plan," *Annales des Télécommun.*, vol. 22, pp. 243-248, Sep. 1967.

[19] G. N. Waston, *A Treatise on the Theory of Bessel Functions*. Cambridge, U.K.: Cambridge Univ. Press, 1995.

[20] A. Quarteroni, R. Sacco, and F. Saleri, "Numerical mathematics," in *Texts in Applied Mathematics Series*, vol. 37, 2nd ed. New York, NY, USA: Springer, 2007.

[21] A. R. von Hippel, *Dielectric Materials and Applications*. Cambridge, MA, USA: MIT Press, 1954.

[22] G. Mur, "Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic-field equations," *IEEE Trans. Electromagn. Compat.*, vol. EMC-23, no. 4, pp. 377-382, Nov. 1981.

[23] A. Reineix, F. Jecko, and B. Jecko, "Time domain modelling techniques improvement for new EMC problems," in *Proc. IEEE Int. Symp. Electromagn. Compat.*, Aug. 1997, pp. 231-235.

[24] V. Mouysset, P. A. Mazet, and P. Borderies, "A multidomain decomposition method to solve electromagnetic scattering problems in time domain," *Radio Sci.*, vol. 42, no. 4, Aug. 2007. Online. Available: http://onlinelibrary.wiley.com/doi/10.1029/2005RS003417/full

[25] K. Demarest, R. Plumb, and Z. Huang, "FDTD modeling of scatterers in stratified media," *IEEE Trans. Antennas Propag.*, vol. 43, no. 10, pp. 1164-1168, Oct. 1995.

[26] S. D. Gedney, "An anisotropic perfectly matched layer absorbing medium for the truncation of FDTD lattices," *IEEE Trans. Antennas Propag.*, vol. 44, no. 12, pp. 1630-1639, Dec. 1996.

[27] R. K. Tewari, S. Swarup, and M. N. Roy, "Evaluation of relative permittivity and conductivity of forest slab from experimentally measured data on lateral wave attenuation constant," *Int. J. Electron.*, vol. 61, no. 5, pp. 597-605, May 1986.