An Underground Radio Wave Propagation Prediction Model for Digital Agriculture

Abstract

Underground sensing and propagation of Signals in the Soil (SitS) medium is an electromagnetic issue. The path loss prediction with higher accuracy is an open research subject in digital agriculture monitoring applications for sensing and communications. The statistical data are predominantly derived from site-specific empirical measurements, which is considered an impediment to universal application. Nevertheless, in the existing literature, statistical approaches have been applied to the SitS channel modeling, where impulse response analysis and the Friis open space transmission formula are employed as the channel modeling tool in different soil types under varying soil moisture conditions at diverse communication distances and burial depths. In this article, an electromagnetic field analysis is presented as an enhanced monitoring approach for subsurface radio wave propagation and underground sensing applications in the field of digital agriculture. The signal strength results are shown for different distances and depths in the subsurface medium. The analysis shows that the lateral wave is the dominant wave in subsurface communications. Moreover, the shallow depths are more suitable for soil moisture sensing and long-range underground communications. The developed paradigm leads to advanced system design for real-time soil monitoring applications.

1. Introduction

The subsurface radio wave propagation is vital for sending and receiving wireless signals in the next-generation precision agriculture Internet of Underground Things (IOUT) [1,2,3]. The Signals in the Soil (SitS) wireless communication approach has only focused mainly on the empirical methods [3,4,5,6,7,8,9,10,11,12,13,14,15,16], thereby ignoring the importance of the physical characteristics behind the fundamental electromagnetic subsurface physics [17,18,19,20,21,22].

The subsurface wireless underground communications involve two different material half-spaces (e.g., soil and air half-space). Moreover, both the transmitter and receiver can be placed in any of the half-spaces, which leads to three different types of subsurface wireless channels, which are different from Over-The-Air (OTA) wireless communications channels. These unique channels involve the soil-to-air half-space, air-to-soil half-space, and subsurface on-soil communications. Hence, the consideration of electromagnetic fields becomes very important in subsurface wireless communications. Therefore, the aim of this article is to present subsurface electromagnetic field analysis (an analytical solution) involving soil and air half-spaces, which is based on Maxwell’s physics. This article presents a subsurface radio wave propagation prediction model based on the fundamental EM analysis related to subsurface propagation modeling that considers the electromagnetic field of a unit vertical and horizontal electric dipole in the presence of a plane boundary.

This underground signals analysis is useful for subsurface radio wave propagation and underground antenna characterization in soil medium. It is also applicable to radio wave propagation in forests and other environments where there is a need for accurate wireless channel models that can be used to produce actual empirical results without the need for any linear and nonlinear regression. In addition, these channel models will also aid in the widespread deployment of agricultural wireless sensor networks, the development of improved efficiency protocols, and the improved subsurface communication system design. In this article, a direct approach based on electromagnetic field analysis has been applied to radio wave propagation in the stratified soil medium for smart farming applications in precision agriculture and the Internet of Underground Things.

The rest of the article is organized as follows: the background and related work is presented in Section 2. The electromagnetic field of a unit vertical electric dipole in the presence of a plane boundary is analyzed in Section 3. In Section 4, the electromagnetic field of a horizontal electric dipole in the presence of a plane boundary is discussed. Model evaluations and results are presented in Section 5. The applications of this work are discussed in Section 6. The paper concludes in Section 7.

2. Background and Related Work

The Per- and Poly-Fluoroalkyl Substances’ (PFAS) management techniques are important in order to establish a proper detection and analysis methods to achieve greater efficiency [23]. Timely PFAS management minimizes risks to humans and environment. PFAS treatment approaches that result from improper detection lead to detrimental environmental effects. Proper PFAS detection will help to reduce the potential of chemical leaching to the soil and water resources and also into the environment.

Among existing techniques, Granular Activated Carbon (GAC) is a growing technology in PFAS treatment in water [23,24,25]. However, there is a significant lack of data and procedure development in terms of fundamental understanding and quantification of the medium properties. The adsorptive and destructive technologies are considered for both soils and water [26,27,28]. Other remediation approaches are anion-exchange, ozofractionation, chemical oxidation, electrochemical oxidation, sonolysis, soil stabilization, and thermal technologies [29,30,31,32,33].

These treatment technologies are not best suited to provide PFAS management systems with almost real-time sensing data to facilitate fast decision making [34,35,36]. Thus, proper treatment cannot be applied at the right place of the waste system at the right time. Failure to consider the spatial variability of PFAS in management decisions results in inefficient management. Accordingly, the potential of chemical leaching from the landfill surface is increased. Timely information of temporal and spatial soil PFAS patterns can significantly aid PFAS decision makers in better managing their treatment operations to achieve higher treatment efficiency. The real-time knowledge of spatial PFAS demand can also further advance our understanding of variable physico-chemical and biological dynamics. Accordingly, more effective management and treatment strategies can be developed for PFAS removal.

The first generation of wireless underground communication systems was developed in the 2000s and was based on the application of the Over-The-Air (OTA) Friis transmission formula to the underground wireless channel [6,15]. In the early 2010s, the second generation of wireless underground communications was developed based on the impulse response analysis and advanced multi-carrier techniques [10,37]. The major development in the second-generation wireless underground communications is the statistical model of the underground channel, which has been developed based on the extensive empirical campaign both in field and testbed environments [9,10]. The model can predict the Root Mean Square (RMS) delay spread, coherence bandwidth, and path loss. Because the statistical model has been developed based on the statistical and information theoretic perspective, it lacks the insight into the physics of subsurface radio wave propagation, which can only be provided by Maxwell–Poynting theory-based signal processing. However, there is no underground channel model available based on the subsurface electromagnetic-based signal processing analysis, and only the numerical evaluations and results have been presented. Accordingly, to gain the physical insight into the propagation of wireless underground communication waves and to ensure that the wireless underground communication system performance is optimal in different subsurface media, it is vital to have an underground propagation prediction model based on electromagnetic using physics principals applicable to radiation efficiency, antenna design, and power transfer at different distances and depths in the soil medium.

Recently, many advanced analytical solutions have been reported in the literature for the Sommerfeld integral [38,39,40,41,42,43]. An alternative solution to the Sommerfeld integral based on the branch-cut integral has been discussed in [42]. An analysis of a two-layered structure using the longitudinal spectrum for space waves, surface waves, and Zenneck waves has been presented in [43]. A Zenneck wave-based UGchannel model has been developed in [44]. However, this model was developed using the approximation of the EM analysis and has higher error.

The aim of this paper is to present the wireless underground communications prediction model for different distances, and depths, under different soil moisture levels, capable of predicting in different soil textures. This work is the first work to develop a subsurface radio wave propagation prediction model completely based on the electromagnetic-based signal processing analysis.

3. The Electromagnetic Field of a Unit Vertical Electric Dipole in the Presence of a Plane Boundary

The change of the soil moisture impacts the signal propagation in the underground soil medium. The increase in Volumetric Water Content (VWC) leads to a decrease in the received signal strength at the underground receiver. In this section, we develop the field solutions by ascertaining the VWC percentage by using underground sensors and accordingly using Peplinski’s model to obtain the complex permittivity of the soil under consideration. The obtained complex permittivity value for a given VWC level is employed in the model along with other soil parameters (e.g., soil type, soil bulk density, burial depth in medium, and communication distance) to predict radio wave propagation for a particular agricultural field. Accordingly, the subsurface communication system can be tailored for site-specific and soil-specific applications. Overall, the developed model can be effectively used to predict propagation characteristics in the soil medium with the change in soil moisture for digital agriculture applications. These developments are vital to understanding the properties of subsurface waves and provide useful insights for system design.

A list of notations used in the paper is given in the “Abbreviations”. The electromagnetic field generated by a vertical electric dipole near the boundary [18] between two quite different material half-spaces is obtained from the explicit integrals for the component differentiation of the electromagnetic field. The vertical dipole with the unit electric moment (i.e., $I\mathsf{\Delta }l=1\mathrm{A}\mathrm{m}$) is located on the downward-directed z-axis at a distance d from the origin of the coordinates on the interface, i.e., the $xy$-plane. A drawing of the geometry of the underground antenna is shown in Figure 1. The electromagnetic field is to be determined at an arbitrary point $(x,y,z)$ in rectangular or $(\rho ,\varphi ,z)$ in cylindrical coordinates. The lower half-space ($z\ge 0$) is Region 1, and the upper half-space ($z\le 0$) is Region 2. The two regions are characterized by the complex wave numbers $k_j={\beta }_j+i{\alpha }_j=\omega {\left({\mu }_j{\widetilde{ϵ}}_j\right)}^{1/2}$, where ${\widetilde{ϵ}}_j={ϵ}_j+i{\sigma }_j/\omega$ and ${ϵ}_j={ϵ}_0{ϵ}_{jr}$ with $j=1,2$. It is assumed that both regions are nonmagnetic, so that ${\mu }_1={\mu }_2={\mu }_0=4\pi \times {10}^{-7}\mathrm{H}/\mathrm{m}$. The time dependence $e^{-i\omega t}$ is used.

From Maxwell’s equations, we obtain an ordinary differential equation for ${\overline{B}}_{jx}$:

(1)$\left(\frac{{\mathrm{d}}^2}{\mathrm{d}{z}^2}+{\gamma }_j^2\right){\overline{B}}_{jx}=-i\eta {\mu }_0\delta (z-d),$

with:

(2)${\gamma }_j={(k_j^2-{\xi }^2-{\eta }^2)}^{1/2},j=1,2.$

Other components can be expressed by ${\overline{B}}_{jx}$,

(3)$\begin{array}{ccc}\hfill {\overline{B}}_{jy}& =& -\frac{\xi }{\eta }{\overline{B}}_{jx},\hfill \end{array}$

(4)$\begin{array}{ccc}\hfill {\overline{E}}_{jx}& =& -\frac{i\omega }{k_j^2}\frac{\partial {\overline{B}}_{jy}}{\partial z}=-\frac{i\omega }{k_j^2}\frac{\xi }{\eta }\frac{\partial {\overline{B}}_{jx}}{\partial z},\hfill \end{array}$

(5)$\begin{array}{ccc}\hfill {\overline{E}}_{jy}& =& \frac{i\omega }{k_j^2}\frac{\partial {\overline{B}}_{jx}}{\partial z},\hfill \end{array}$

(6)$\begin{array}{ccc}\hfill {\overline{E}}_{jx}& =& \frac{\omega }{\eta k_j^2}\left(\frac{{\mathrm{d}}^2}{\mathrm{d}{z}^2}+k_j^2\right){\overline{B}}_{jx}.\hfill \end{array}$

The equations have been solved in [45], and the complete results for the two regions are summarized below. The following shorthand notation is used:

(7)$N\equiv k_1^2{\gamma }_2+k_2^2{\gamma }_1.$

Region 1, $z\ge 0$:

(8)$\begin{array}{ccc}\hfill {\overline{B}}_{1x}(\xi ,\eta ,z)& =& -{\mu }_0\eta (\frac{e^{i{\gamma }_1|z-d|}}{2{\gamma }_1}-\frac{e^{i{\gamma }_1(z+d)}}{2{\gamma }_1}\hfill \\ & & +\frac{k_2^2{e}^{i{\gamma }_1(z+d)}}{N}),\hfill \end{array}$

(9)$\begin{array}{ccc}\hfill {\overline{B}}_{1y}(\xi ,\eta ,z)& =& -\frac{\xi }{\eta }{\overline{B}}_{1x}(\xi ,\eta ,z),\hfill \end{array}$

(10)$\begin{array}{ccc}\hfill {\overline{B}}_{1z}(\xi ,\eta ,z)& =& 0,\hfill \end{array}$

(11)$\begin{array}{ccc}\hfill {\overline{E}}_{1x}(\xi ,\eta ,z)& =& -\frac{\omega {\mu }_0\xi }{k_1^2}(\mp \frac{e^{i{\gamma }_1|z-d|}}{2}+\frac{e^{i{\gamma }_1(z+d)}}{2}\hfill \\ & & -\frac{k_2^2{\gamma }_1{e}^{i{\gamma }_1(z+d)}}{N})\left\{\begin{array}{c}z>d\hfill \\ 0\le z\le d\hfill \end{array}\right.,\hfill \end{array}$

(12)$\begin{array}{ccc}\hfill {\overline{E}}_{1y}(\xi ,\eta ,z)& =& \frac{\eta }{\xi }{\overline{E}}_{1x}(\xi ,\eta ,z),\hfill \end{array}$

(13)$\begin{array}{ccc}\hfill {\overline{E}}_{1z}(\xi ,\eta ,z)& =& -\frac{\omega {\mu }_0({\xi }^2+{\eta }^2)}{k_1^2}(\frac{e^{i{\gamma }_1|z-d|}}{2{\gamma }_1}\hfill \\ & & -\frac{e^{i{\gamma }_1(z+d)}}{2{\gamma }_1}+\frac{k_2^2{e}^{i{\gamma }_1(z+d)}}{N}).\hfill \end{array}$

Region 2, $z\ge 0$:

(14)$\begin{array}{ccc}\hfill {\overline{B}}_{2x}(\xi ,\eta ,z)& =& -{\mu }_0\eta \frac{k_2^2{e}^{i{\gamma }_{1}d}{e}^{-i{\gamma }_{2}z}}{N},\hfill \end{array}$

(15)$\begin{array}{ccc}\hfill {\overline{B}}_{2y}(\xi ,\eta ,z)& =& -\frac{\xi }{\eta }{\overline{B}}_{2x}(\xi ,\eta ,z),\hfill \end{array}$

(16)$\begin{array}{ccc}\hfill {\overline{B}}_{2z}(\xi ,\eta ,z)& =& 0,\hfill \end{array}$

(17)$\begin{array}{ccc}\hfill {\overline{E}}_{2x}(\xi ,\eta ,z)& =& -\omega {\mu }_0\xi \frac{{\gamma }_2{e}^{i{\gamma }_{1}d}{e}^{-i{\gamma }_{2}z}}{N},\hfill \end{array}$

(18)$\begin{array}{ccc}\hfill {\overline{E}}_{2y}(\xi ,\eta ,z)& =& \frac{\eta }{\xi }{\overline{E}}_{2x}(\xi ,\eta ,z),\hfill \end{array}$

(19)$\begin{array}{ccc}\hfill {\overline{E}}_{2z}(\xi ,\eta ,z)& =& -\omega {\mu }_0({\xi }^2+{\eta }^2)\frac{e^{i{\gamma }_{1}d}{e}^{-i{\gamma }_{2}z}}{N}.\hfill \end{array}$

The electromagnetic field at an arbitrary position is the integrals of the components:

(20)$\begin{array}{ccc}\hfill \mathbf{E}(x,y,z)& =& \frac{1}{{\left(2\pi \right)}^2}{\int }_{-\infty }^{\infty }\mathrm{d}\xi {\int }_{-\infty }^{\infty }\mathrm{d}\eta e^{i(\xi x+\eta y)}\overline{\mathbf{E}}(\xi ,\eta ,z),\hfill \end{array}$

(21)$\begin{array}{ccc}\hfill \mathbf{B}(x,y,z)& =& \frac{1}{{\left(2\pi \right)}^2}{\int }_{-\infty }^{\infty }\mathrm{d}\xi {\int }_{-\infty }^{\infty }\mathrm{d}\eta e^{i(\xi x+\eta y)}\overline{\mathbf{B}}(\xi ,\eta ,z).\hfill \end{array}$

In order to evaluate the integrals for the components of the electromagnetic field at $(\rho ,z)$ in Region 1 when the vertical electric dipole is also in Region 1 at $(0,d)$, it is convenient to examine separately for different parts. These are the direct field, the reflected field due to an ideal image, the lateral-wave field, and a correction for the reflected field when this is not accurately that of an ideal image (Figure 2). Specifically, let:

(22)$\begin{array}{ccc}\hfill {\mathbf{E}}_1(\rho ,z)& =& {\mathbf{E}}_1^d(\rho ,z)+{\mathbf{E}}_1^i(\rho ,z)+{\mathbf{E}}_1^L(\rho ,z)+{\mathbf{E}}_1^c(\rho ,z),\hfill \end{array}$

(23)$\begin{array}{ccc}\hfill {\mathbf{B}}_1(\rho ,z)& =& {\mathbf{B}}_1^d(\rho ,z)+{\mathbf{B}}_1^i(\rho ,z)+{\mathbf{B}}_1^L(\rho ,z)+{\mathbf{B}}_1^c(\rho ,z).\hfill \end{array}$

The direct field is given by:

(24)$\begin{array}{ccc}\hfill B_{1\varphi }^d(\rho ,z)& =& -\frac{{\mu }_0}{4\pi }{e}^{i{k}_1{r}_1}\left(\frac{i{k}_1}{r_1}-\frac{1}{r_1^2}\right)\left(\frac{\rho }{r_1}\right),\hfill \end{array}$

(25)$\begin{array}{ccc}\hfill E_{1\rho }^d(\rho ,z)& =& -\frac{\omega {\mu }_0}{4\pi k_1^2}{e}^{i{k}_1{r}_1}\left(\frac{i{k}_1^2}{r_1}-\frac{3k_1}{r_1^2}-\frac{3i}{r_1^3}\right)\left(\frac{\rho }{r_1}\right)\hfill \\ & & \left(\frac{z-d}{r_1}\right),\hfill \end{array}$

(26)$\begin{array}{ccc}\hfill E_{1z}^d(\rho ,z)& =& \frac{\omega {\mu }_0}{4\pi k_1^2}{e}^{i{k}_1{r}_1}[\frac{i{k}_1^2}{r_1}-\frac{k_1}{r_1^2}-\frac{i}{r_1^3}\hfill \\ & & -{\left(\frac{z-d}{r_1}\right)}^2\left(\frac{i{k}_1^2}{r_1}-\frac{3k_1}{r_1^2}-\frac{3i}{r_1^3}\right)].\hfill \end{array}$

The ideal reflected field is the field at $(\rho ,z)$ due to an image antenna at $(0,-d)$. The electric moment of the image is the negative of that of the source.

(27)$\begin{array}{ccc}\hfill B_{1\varphi }^i(\rho ,z)& =& \frac{{\mu }_0}{4\pi }{e}^{i{k}_1{r}_2}\left(\frac{i{k}_1}{r_2}-\frac{1}{r_2^2}\right)\left(\frac{\rho }{r_2}\right),\hfill \end{array}$

(28)$\begin{array}{ccc}\hfill E_{1\rho }^i(\rho ,z)& =& \frac{\omega {\mu }_0}{4\pi k_1^2}{e}^{i{k}_1{r}_2}\left(\frac{i{k}_1^2}{r_2}-\frac{3k_1}{r_2^2}-\frac{3i}{r_2^3}\right)\left(\frac{\rho }{r_2}\right)\hfill \\ & & \left(\frac{z+d}{r_2}\right),\hfill \end{array}$

(29)$\begin{array}{ccc}\hfill E_{1z}^i(\rho ,z)& =& -\frac{\omega {\mu }_0}{4\pi k_1^2}{e}^{i{k}_1{r}_2}[\frac{i{k}_1^2}{r_2}-\frac{k_1}{r_2^2}-\frac{i}{r_2^3}\hfill \\ & & -{\left(\frac{z+d}{r_2}\right)}^2\left(\frac{i{k}_1^2}{r_2}-\frac{3k_1}{r_2^2}-\frac{3i}{r_2^3}\right)].\hfill \end{array}$

Given conditions:

(30)$|k_1|\ge 3|k_2|,\rho \ge 5|z|,\rho \ge 5d,|k_1\rho |\ge 3,$

the rest field can be written in the following compact form:

(31)$\begin{array}{ccc}\hfill B_{1\varphi }^L(\rho ,z)& =& -\frac{{\mu }_0{k}_2^2}{2\pi k_1^2}{e}^{i{k}_1(z+d)}{e}^{i{k}_2\rho }f(\rho ;k_1,k_2),\hfill \end{array}$

(32)$\begin{array}{ccc}\hfill B_{1\varphi }^c(\rho ,z)& =& \frac{i{\mu }_0{k}_2^2}{2\pi k_1^2}\left(\frac{z+d}{\rho }\right)\left(\frac{i{k}_1}{\rho }-\frac{3}{2{\rho }^2}\right)e^{i{k}_1{r}_1},\hfill \end{array}$

(33)$\begin{array}{ccc}\hfill E_{1\varphi }^L(\rho ,z)& =& -\frac{\omega {\mu }_0{k}_2^2}{2\pi k_1^3}{e}^{i{k}_1(z+d)}{e}^{i{k}_2\rho }f(\rho ;k_1,k_2),\hfill \end{array}$

(34)$\begin{array}{ccc}\hfill E_{1\varphi }^c(\rho ,z)& =& \frac{i\omega {\mu }_0{k}_2^2}{2\pi k_1^3}{e}^{i{k}_1{r}_2}\left(\frac{1}{{\rho }^2}+\frac{3i}{2k_1{\rho }^3}\right),\hfill \end{array}$

(35)$\begin{array}{ccc}\hfill E_{1z}^L(\rho ,z)& =& \frac{\omega {\mu }_0{k}_2^3}{2\pi k_1^4}{e}^{i{k}_1(z+d)}{e}^{i{k}_2\rho }g(\rho ;k_1,k_2),\hfill \end{array}$

(36)$\begin{array}{ccc}\hfill E_{1z}^c(\rho ,z)& =& -\frac{i\omega {\mu }_0{k}_2^2}{2\pi k_1^4}\left(\frac{z+d}{\rho }\right)\hfill \\ & & \left(\frac{i{k}_1^2}{\rho }-\frac{k_1}{2{\rho }^2}+\frac{7i}{8{\rho }^3}\right)e^{i{k}_1{r}_2},\hfill \end{array}$

where the radial functions $f(\rho ;k_1,k_2)$ and $g(\rho ;k_1,k_2)$ are defined as follows:

(37)$\begin{array}{ccc}\hfill f(\rho ;k_1,k_2)& =& \frac{i{k}_2}{\rho }-\frac{1}{{\rho }^2}-\frac{k_2^3}{k_1}{\left(\frac{\pi }{k_2\rho }\right)}^{1/2}{e}^{-ip\mathcal{F}\left(p\right)},\hfill \end{array}$

(38)$\begin{array}{ccc}\hfill g(\rho ;k_1,k_2)& =& \frac{i{k}_2}{\rho }-\frac{1}{{\rho }^2}-\frac{i}{k_2{\rho }^3}\hfill \\ & & -\frac{k_2^3}{k_1}{\left(\frac{\pi }{k_2\rho }\right)}^{1/2}{e}^{-ip\mathcal{F}\left(p\right)}.\hfill \end{array}$

Here:

(39)$p\equiv \frac{k_2^3\rho }{2k_1^2}$

is called the numerical distance. The function $\mathcal{F}\left(p\right)$ is:

(40)$\begin{array}{ccc}\hfill \mathcal{F}\left(p\right)& =& {\int }_p^{\infty }\frac{e^{it}}{{\left(2\pi t\right)}^{(}1/2)}\mathrm{d}t\hfill \\ & =& \frac{1}{2}(1+i)-C_2\left(p\right)-i{S}_2\left(p\right),\hfill \end{array}$

where $C_2+i{S}_2$ is the Fresnel integral. When the argument p is sufficiently great, specifically when:

(41)$\left|p\right|\ge 4\mathrm{or}|k_2\rho |\ge 8|\frac{k_1^2}{k_2^2}|,$

$f(\rho ;k_1,k_2)$ and $g(\rho ;k_1,k_2)$ can be approximated as:

(42)$f(\rho ;k_1,k_2)\sim g(\rho ;k_1,k_2)=-\frac{1}{{\rho }^2}\left(\frac{k_1^2}{k_2^2}+1\right)\sim -\frac{k_1^2}{k_2^2{\rho }^2}.$

If the conditions in (30) are not satisfied, the complete electromagnetic field in Region 1 is given by:

(43)$\begin{array}{ll}& E_{i\rho }(\rho ,z)\\ \hfill =& -\frac{\omega {\mu }_0}{2\pi k_1^2}\left\{\frac{k_2^2}{k_1}\right[f(\rho ;k_1,k_2)e^{i{k}_2\rho }{e}^{i{k}_1(z+d)}\\ & -i{e}^{i{k}_1{r}_2}\left(\frac{1}{{\rho }^2}+\frac{3i}{2k_1{\rho }^3}\right)]\\ & +\frac{e^{i{k}_1{r}_1}}{2}\left(\frac{i{k}_1^2}{r_1}-\frac{3k_1}{r_1^2}-\frac{3i}{r_1^3}\right)\left(\frac{\rho }{r_1}\right)\left(\frac{z-d}{r_1}\right)\\ & -\frac{e^{i{k}_1{r}_2}}{2}\left(\frac{i{k}_1^2}{r_2}-\frac{3k_1}{r_2^2}-\frac{3i}{r_2^3}\right)\left(\frac{\rho }{r_2}\right)\left(\frac{z+d}{r_2}\right)\},\end{array}$

(44)$\begin{array}{l}\hfill \\ & E_{iz}(\rho ,z)\\ \hfill =& \frac{\omega {\mu }_0}{2\pi k_1^2}\left\{\frac{k_2^2}{k_1^2}\right[k_{2}g(\rho ;k_1,k_2)e^{i{k}_2\rho }{e}^{i{k}_1(z+d)}\\ & -i{e}^{i{k}_1{r}_2}\left(\frac{z+d}{\rho }\right)\left(\frac{i{k}_1^2}{\rho }-\frac{k_1}{2{\rho }^2}+\frac{7i}{8{\rho }^3}\right)]\\ & +\frac{e^{i{k}_1{r}_1}}{2}[\frac{i{k}_1^2}{r_1}-\frac{k_1}{r_1^2}-\frac{i}{r_1^3}\\ & -{\left(\frac{z-d}{r_1}\right)}^2\left(\frac{i{k}_1^2}{r_1}-\frac{3k_1}{r_1^2}-\frac{3i}{r_1^3}\right)]\\ & -\frac{e^{i{k}_1{r}_2}}{2}[\frac{i{k}_1^2}{r_2}-\frac{k_1}{r_2^2}-\frac{i}{r_2^3}\\ & -{\left(\frac{z+d}{r_2}\right)}^2\left(\frac{i{k}_1^2}{r_2}-\frac{3k_1}{r_2^2}-\frac{3i}{r_2^3}\right)\left]\right\},\end{array}$

(45)$\begin{array}{l}\hfill \\ & B_{i\varphi }(\rho ,z)\\ \hfill =& -\frac{{\mu }_0}{2\pi }\left\{\frac{k_2^2}{k_1^2}\right[f(\rho ;k_1,k_2)e^{i{k}_2\rho }{e}^{i{k}_1(z+d)}\\ & -i{e}^{i{k}_1{r}_2}\left(\frac{z+d}{\rho }\right)\left(\frac{i{k}_1}{\rho }-\frac{3}{2{\rho }^2}\right)]\\ & +\frac{e^{i{k}_1{r}_1}}{2}\left(\frac{i{k}_1}{r_1}-\frac{1}{r_1^2}\right)\left(\frac{\rho }{r_1}\right)\\ & -\frac{e^{i{k}_1{r}_2}}{2}\left(\frac{i{k}_1}{r_2}-\frac{1}{r_2^2}\right)\left(\frac{\rho }{r_2}\right)\}.\end{array}$

For the complete field in Region 2, the following results are obtained.

(46)$\begin{array}{ll}\hfill & B_{2\varphi }(\rho ,z)\\ \hfill =& -\frac{{\mu }_0{k}_2^2}{2\pi k_1^2}\left\{e^{i{k}_{1}d}{e}^{i{k}_2{r}_0}\right[\left(\frac{i{k}_2}{r_0}-\frac{1}{r_0^2}\right)\left(\frac{\rho }{r_0}\right)\\ & -\frac{k_2^3}{k_1}{\left(\frac{\pi }{k_2\rho }\right)}^{1/2}{e}^{-iP\mathcal{F}\left(P\right)}]\\ & -i{e}^{k_{1}z}{e}^{i{k}_1{r}_d}\left(\frac{d}{\rho }\right)\left(\frac{i{k}_1}{\rho }-\frac{3}{2{\rho }^2}\right)\},\end{array}$

(47)$\begin{array}{ll}\hfill \sim & -\frac{{\mu }_0{k}_2^2}{2\pi k_1^2}[e^{i{k}_{1}d}{e}^{i{k}_2{r}_0}f(\rho ,k_1,k_2)\\ & -i{e}^{k_{1}z}{e}^{i{k}_1{r}_d}\left(\frac{d}{\rho }\right)\left(\frac{i{k}_1}{\rho }-\frac{3}{2{\rho }^2}\right)],\end{array}$

$\begin{array}{l}\hfill \\ & E_{2\rho }(\rho ,z)\\ \hfill =& -\frac{\omega {\mu }_0}{2\pi k_1^2}\left(\frac{k_2^2}{k_1}\right\{e^{i{k}_{1}d}{e}^{i{k}_2{r}_0}[\left(\frac{i{k}_2}{r_0}-\frac{1}{r_0^2}\right)\left(\frac{\rho }{r_0}\right)\\ & -\frac{k_2^3}{k_1}{\left(\frac{\pi }{k_2\rho }\right)}^{1/2}{e}^{-iP\mathcal{F}\left(P\right)}]\\ & -i{e}^{k_{1}z}{e}^{i{k}_1{r}_d}\left(\frac{1}{{\rho }^2}+\frac{3i}{2k_1{\rho }^3}\right)\}\\ & -e^{k_{1}z}{e}^{i{k}_1{r}_d}\left(\frac{i{k}_1^2}{r_d}-\frac{3k_1}{r_d^2}-\frac{3i}{r_d^3}\right)\left(\frac{\rho d}{r_d^2}\right)\end{array}$

(48)$\begin{array}{cc}& +e^{i{k}_{1}d}{e}^{i{k}_2{r}_0}\left(\frac{i{k}_2^2}{r_0}-\frac{3k_2}{r_0^2}-\frac{3i}{r_0^3}\right)\left(\frac{\rho z}{r_0^2}\right)),\end{array}$

(49)$\begin{array}{ll}\hfill \sim & -\frac{\omega {\mu }_0}{2\pi k_1^2}\left\{\frac{k_2^2}{k_1}\right[e^{i{k}_{1}d}{e}^{i{k}_2{r}_0}f(\rho ,k_1,k_2)\\ & -i{e}^{k_{1}z}{e}^{i{k}_1{r}_d}\left(\frac{1}{{\rho }^2}+\frac{3i}{2k_1{\rho }^3}\right)]\\ & -e^{k_{1}z}{e}^{i{k}_1{r}_d}\left(\frac{i{k}_1^2}{\rho }-\frac{3k_1}{{\rho }^2}-\frac{3i}{{\rho }^3}\right)\left(\frac{d}{\rho }\right)\\ & +e^{i{k}_{1}d}{e}^{i{k}_2{r}_0}\left(\frac{i{k}_2^2}{\rho }-\frac{3k_2}{{\rho }^2}-\frac{3i}{{\rho }^3}\right)\left(\frac{z}{\rho }\right)\},\end{array}$

(50)$\begin{array}{l}\hfill \\ & E_{2z}(\rho ,z)\\ \hfill =& \frac{\omega {\mu }_0}{2\pi k_1^2}\left\{e^{i{k}_{1}d}{e}^{i{k}_2{r}_0}\right[\frac{i{k}_2^2}{r_0}-\frac{k_2}{r_0^2}-\frac{i}{r_0^3}\\ & -\left(\frac{i{k}_2}{r_0}-\frac{3k_2}{r_0^2}-\frac{3i}{r_0^3}\right)\left(\frac{z^2}{r_0^2}\right)\\ & -\frac{k_2^4}{k_1}{\left(\frac{\pi }{k_2\rho }\right)}^{1/2}{e}^{-iP\mathcal{F}\left(P\right)}]\\ & -i{e}^{k_{1}z}{e}^{i{k}_1{r}_d}\left(\frac{d}{\rho }\right)\left(\frac{i{k}_1^2}{\rho }-\frac{k_1}{2{\rho }^2}+\frac{7i}{8{\rho }^3}\right)\},\end{array}$

(51)$\begin{array}{ll}\hfill \sim & \frac{\omega {\mu }_0}{2\pi k_1^2}[e^{i{k}_{1}d}{e}^{i{k}_2{r}_0}{k}_{2}g(\rho ,k_1,k_2)\\ & -i{e}^{k_{1}z}{e}^{i{k}_1{r}_d}\left(\frac{d}{\rho }\right)\left(\frac{i{k}_1^2}{\rho }-\frac{k_1}{2{\rho }^2}+\frac{7i}{8{\rho }^3}\right)],\end{array}$

where:

(52)$r_0={({\rho }^2+z^2)}^{1/2},\mathrm{and}{r}_d={({\rho }^2+d^2)}^{1/2}.$

The approximations are valid when the conditions:

(53)$\rho \ge 5\left|z\right|\mathrm{and}\rho \ge 5d$

are satisfied. In this case, $r_0\sim \rho$ and $r_d\sim \rho$ in amplitude, $P\sim R$, and the terms of the order $z^2/{\rho }^2$ are neglected.

4. The Electromagnetic Field of a Horizontal Electric Dipole in the Presence of a Plane Boundary

Consider the x-directed horizontal electric dipole [18] at the point $z=d$ on the downward-directed z-axis in the half-space (Region 1) defined by $z\ge 0$. Region 2 is defined by $z\le 0$. The wave numbers of the two regions are $k_j={\beta }_j+i{\alpha }_j=\omega {\left({\mu }_0{\widetilde{ϵ}}_j\right)}^{1/2}$, where ${\widetilde{ϵ}}_j=ϵ+i{\sigma }_j/\omega$ and $j=1,2$. It is assumed that both regions are nonmagnetic, so that ${\mu }_1={\mu }_2={\mu }_0$. Based on Maxwell’s equations, the following results are reached:

(54)$\begin{array}{ccc}\hfill {\overline{E}}_{jy}& =& \frac{1}{k_j^2-{\xi }^2}\left(-\xi \eta {\overline{E}}_{jy}+i\omega \frac{\partial }{\partial z}{\overline{B}}_{jx}\right),\hfill \end{array}$

(55)$\begin{array}{ccc}\hfill {\overline{E}}_{jz}& =& \frac{1}{k_j^2-{\xi }^2}\left(i\xi \frac{\partial }{\partial z}{\overline{E}}_{jx}+\eta \omega {\overline{B}}_{jx}\right),\hfill \end{array}$

(56)$\begin{array}{ccc}\hfill {\overline{B}}_{jy}& =& \frac{1}{k_j^2-{\xi }^2}\left(-\frac{i{k}_j^2}{\omega }\frac{\partial }{\partial z}{\overline{E}}_{jx}-\xi \eta {\overline{B}}_{jx}\right),\hfill \end{array}$

(57)$\begin{array}{ccc}\hfill {\overline{B}}_{jz}& =& \frac{1}{k_j^2-{\xi }^2}\left(-\frac{k_j^2\eta }{\omega }{\overline{E}}_{jx}+i\xi \frac{\partial }{\partial z}{\overline{B}}_{jx}\right).\hfill \end{array}$

Moreover, the following ordinary differential equations for ${\overline{E}}_{jx}$ and ${\overline{B}}_{jx}$ are also obtained:

(58)$\begin{array}{ccc}\hfill \left(\frac{{\mathrm{d}}^2}{\mathrm{d}{z}^2}+{\gamma }_j^2\right){\overline{E}}_{jx}& =& \frac{\omega {\mu }_0(k_j^2-{\xi }^2)}{i{k}_j^2}\delta (z-d),\hfill \end{array}$

(59)$\begin{array}{ccc}\hfill \left(\frac{{\mathrm{d}}^2}{\mathrm{d}{z}^2}+{\gamma }_j^2\right){\overline{B}}_{jx}& =& 0,\hfill \end{array}$

where:

(60)${\gamma }_j={(k_j^2-{\xi }^2-{\eta }^2)}^{1/2}.$

The six components of the electromagnetic field in the cylindrical coordinates are:

Region 1, $z\ge 0$:

(61)$\begin{array}{lll}\hfill E_{1\rho }& =& -\frac{\omega {\mu }_0}{4\pi k_1^2}cos\varphi \left[{\int }_0^{\infty }\right(k_1^2{J}_0\left(\lambda \rho \right)\hfill \\ & & -\frac{{\lambda }^2}{2}[J_0\left(\lambda \rho \right)-J_2\left(\lambda \rho \right)]){\gamma }_1^{-1}{e}^{i{\gamma }_1|z-d|}\lambda \mathrm{d}\lambda \hfill \\ & & +{\int }_0^{\infty }(\frac{{\gamma }_{1}Q}{2}[J_0\left(\lambda \rho \right)-J_2\left(\lambda \rho \right)]\hfill \\ & & -\frac{k_1^{2}P}{2{\gamma }_1}[J_0\left(\lambda \rho \right)+J_2\left(\lambda \rho \right)]\left)e^{i{\gamma }_1(z+d)}\lambda \mathrm{d}\lambda \right],\hfill \end{array}$

(62)$\begin{array}{lll}\hfill E_{1\varphi }& =& -\frac{\omega {\mu }_0}{4\pi k_1^2}sin\varphi \left[{\int }_0^{\infty }\right(k_1^2{J}_0\left(\lambda \rho \right)\hfill \\ & & -\frac{{\lambda }^2}{2}[J_0\left(\lambda \rho \right)+J_2\left(\lambda \rho \right)]){\gamma }_1^{-1}{e}^{i{\gamma }_1|z-d|}\lambda \mathrm{d}\lambda \hfill \\ & & +{\int }_0^{\infty }(\frac{{\gamma }_{1}Q}{2}[J_0\left(\lambda \rho \right)+J_2\left(\lambda \rho \right)]\hfill \\ & & -\frac{k_1^{2}P}{2{\gamma }_1}[J_0\left(\lambda \rho \right)-J_2\left(\lambda \rho \right)]\left)e^{i{\gamma }_1(z+d)}\lambda \mathrm{d}\lambda \right],\hfill \end{array}$

(63)$\begin{array}{lll}\hfill E_{1z}& =& \frac{i\omega {\mu }_0}{4\pi k_1^2}cos\varphi {\int }_0^{\infty }(\pm e^{i{\gamma }_1|z-d|}+Q{e}^{i{\gamma }_1(z+d)})\hfill \\ & & J_1\left(\lambda \rho \right){\lambda }^2\mathrm{d}\lambda ,\hfill \end{array}$

(64)$\begin{array}{lll}\hfill B_{1\rho }& =& -\frac{{\mu }_0}{4\pi }sin\varphi [\pm {\int }_0^{\infty }{J}_0\left(\lambda \rho \right)e^{i{\gamma }_1|z-d|}\lambda \mathrm{d}\lambda \hfill \\ & & {\int }_0^{\infty }(\frac{Q}{2}[J_0\left(\lambda \rho \right)+J_2\left(\lambda \rho \right)]\hfill \\ & & -\frac{P}{2}[J_0\left(\lambda \rho \right)-J_2\left(\lambda \rho \right)]\left)e^{i{\gamma }_1(z+d)}\lambda \mathrm{d}\lambda \right],\hfill \end{array}$

(65)$\begin{array}{lll}\hfill B_{1\varphi }& =& -\frac{{\mu }_0}{4\pi }cos\varphi [\pm {\int }_0^{\infty }{J}_0\left(\lambda \rho \right)e^{i{\gamma }_1|z-d|}\lambda \mathrm{d}\lambda \hfill \\ & & {\int }_0^{\infty }(\frac{Q}{2}[J_0\left(\lambda \rho \right)-J_2\left(\lambda \rho \right)]\hfill \\ & & -\frac{P}{2}[J_0\left(\lambda \rho \right)+J_2\left(\lambda \rho \right)]\left)e^{i{\gamma }_1(z+d)}\lambda \mathrm{d}\lambda \right],\hfill \end{array}$

(66)$\begin{array}{lll}\hfill B_{1z}& =& \frac{i{\mu }_0}{4\pi }sin\varphi {\int }_0^{\infty }(e^{i{\gamma }_1|z-d|}-P{e}^{i{\gamma }_1(z+d)}){\gamma }_1^{-1}\hfill \\ & & J_1\left(\lambda \rho \right){\lambda }^2\mathrm{d}\lambda .\hfill \end{array}$

Where two signs appear, the upper one is for $z>d$, and the lower is for $0\le z\le d$ (see Figure 3).

Region 2, $z\le 0$:

(67)$\begin{array}{lll}\hfill E_{2\rho }& =& -\frac{\omega {\mu }_0}{4\pi }cos\varphi {\int }_0^{\infty }(M^{-1}[J_0\left(\lambda \rho \right)+J_2\left(\lambda \rho \right)]\hfill \\ & & +\frac{{\gamma }_1{\gamma }_2}{N}[J_0\left(\lambda \rho \right)-J_2\left(\lambda \rho \right)])e^{i({\gamma }_{1}d-{\gamma }_{2}z)}\lambda \mathrm{d}\lambda ,\hfill \end{array}$

(68)$\begin{array}{lll}\hfill E_{2\varphi }& =& \frac{\omega {\mu }_0}{4\pi }sin\varphi {\int }_0^{\infty }(M^{-1}[J_0\left(\lambda \rho \right)-J_2\left(\lambda \rho \right)]\hfill \\ & & +\frac{{\gamma }_1{\gamma }_2}{N}[J_0\left(\lambda \rho \right)+J_2\left(\lambda \rho \right)])e^{i({\gamma }_{1}d-{\gamma }_{2}z)}\lambda \mathrm{d}\lambda ,\hfill \end{array}$

(69)$\begin{array}{ccc}\hfill E_{2z}& =& -\frac{i\omega {\mu }_0}{2\pi }cos\varphi {\int }_0^{\infty }\frac{{\gamma }_1}{N}{J}_1\left(\lambda \rho \right)e^{i({\gamma }_{1}d-{\gamma }_{2}z)}{\lambda }^2\mathrm{d}\lambda ,\hfill \end{array}$

(70)$\begin{array}{lll}\hfill B_{2\rho }& =& \frac{{\mu }_0}{4\pi }sin\varphi {\int }_0^{\infty }(\frac{{\gamma }_2}{M}[J_0\left(\lambda \rho \right)-J_2\left(\lambda \rho \right)]\hfill \\ & & +\frac{k_2^2{\gamma }_1}{N}[J_0\left(\lambda \rho \right)+J_2\left(\lambda \rho \right)])e^{i({\gamma }_{1}d-{\gamma }_{2}z)}\lambda \mathrm{d}\lambda ,\hfill \end{array}$

(71)$\begin{array}{lll}\hfill B_{2\varphi }& =& \frac{{\mu }_0}{4\pi }cos\varphi {\int }_0^{\infty }(\frac{{\gamma }_2}{M}[J_0\left(\lambda \rho \right)+J_2\left(\lambda \rho \right)]\hfill \\ & & +\frac{k_2^2{\gamma }_1}{N}[J_0\left(\lambda \rho \right)-J_2\left(\lambda \rho \right)])e^{i({\gamma }_{1}d-{\gamma }_{2}z)}\lambda \mathrm{d}\lambda ,\hfill \end{array}$

(72)$\begin{array}{ccc}\hfill B_{2z}& =& \frac{i{\mu }_0}{2\pi }sin\varphi {\int }_0^{\infty }{M}^{-1}{J}_1\left(\lambda \rho \right)e^{i({\gamma }_{1}d-{\gamma }_{2}z)}{\lambda }^2\mathrm{d}\lambda .\hfill \end{array}$

In these equations,

(73)$\begin{array}{ccc}\hfill M& \equiv & {\gamma }_1+{\gamma }_2,\hfill \end{array}$

(74)$\begin{array}{ccc}\hfill N& \equiv & k_1^2{\gamma }_2+k_2^2{\gamma }_1,\hfill \end{array}$

(75)$\begin{array}{ccc}\hfill P& \equiv & \frac{{\gamma }_2-{\gamma }_1}{M}=\frac{{\gamma }_2-{\gamma }_1}{{\gamma }_2+{\gamma }_1},\hfill \end{array}$

(76)$\begin{array}{ccc}\hfill Q& \equiv & \frac{k_1^2{\gamma }_2-k_2^2{\gamma }_1}{N}=\frac{k_1^2{\gamma }_2-k_2^2{\gamma }_1}{k_1^2{\gamma }_2+k_2^2{\gamma }_1},\hfill \end{array}$

and $J_n\left(\lambda \rho \right)$ is the integral form of the Bessel function. It is:

(77)$J_n\left(\lambda \rho \right)=\frac{i^{-n}}{2\pi }{\int }_0^{2\pi }{e}^{i\lambda \rho cos\theta }{e}^{in\theta }d\theta .$

With conditions:

(78)$\rho \ge 5d,\rho \ge 5\left|z\right|,\mathrm{and}|k_2\rho |\ge 8|\frac{k_1^2}{k_2^2}|,$

the expressions in Region 1 can be simplified.

(79)$\begin{array}{ll}& E_{1\rho }(\rho ,\varphi ,z)\\ \hfill =& -\frac{\omega {\mu }_0}{2\pi k_1^2}cos\varphi \{k_{2}g(\rho ;k_1,k_2)e^{i{k}_2\rho }{e}^{i{k}_1(z+d)}\\ & -\left(\frac{k_1}{r_1^2}+\frac{i}{r_1^3}\right)e^{i{k}_1{r}_1}+\left(\frac{z+d}{\rho }\right)\left(\frac{i{k}_1}{{\rho }^2}-\frac{3}{2{\rho }^3}\right)e^{i{k}_1{r}_2}\\ & -\frac{1}{2}[{\left(\frac{z-d}{r_1}\right)}^2\left(\frac{1k_1^2}{r_1}-\frac{3k_1}{r_1^2}-\frac{3i}{r_1^3}\right)\times e^{i{k}_1{r}_1}\\ & +{\left(\frac{z+d}{r_2}\right)}^2\left(\frac{i{k}_1^2}{r_2}-\frac{3k_1}{r_2^2}-\frac{3i}{r_2^3}\right)e^{i{k}_1{r}_2}\left]\right\},\end{array}$

(80)$\begin{array}{ll}& E_{1\varphi }(\rho ,\varphi ,z)\\ \hfill =& \frac{\omega {\mu }_0}{2\pi k_1^2}sin\varphi \{k_{2}h(\rho ;k_1,k_2)e^{i{k}_2\rho }{e}^{i{k}_1(z+d)}\\ & +\frac{1}{2}[\left(\frac{i{k}_1^2}{r_2}-\frac{3k_1}{r_2^2}-\frac{3i}{r_2^3}\right)e^{i{k}_1{r}_2}\\ & -\left(\frac{i{k}_1^2}{r_1}-\frac{k_1}{r_1^2}-\frac{i}{r_1^3}\right)e^{i{k}_1{r}_1}]\\ & +i\left(\frac{z+d}{r_2}\right)\left(\frac{i{k}_1^2}{\rho }-\frac{e{k}_1}{2{\rho }^2}-\frac{5i}{8{\rho }^3}\right)e^{i{k}_1{r}_2}\\ & -{\left(\frac{z+d}{r_2}\right)}^2\left(\frac{i{k}_1^2}{r_2}-\frac{3k_1}{r_2^2}-\frac{3i}{r_2^3}\right)e^{i{k}_1{r}_2}\},\end{array}$

(81)$\begin{array}{ll}& E_{1z}(\rho ,\varphi ,z)\\ \hfill =& \frac{\omega {\mu }_0}{2\pi k_1^2}cos\varphi \left\{\frac{k_2^2}{k_1}\right[f(\rho ;k_1,k_2)e^{i{k}_2\rho }{e}^{i{k}_1(z+d)}\\ & -\left(\frac{i}{{\rho }^2}-\frac{3}{2k_1{\rho }^3}\right)e^{i{k}_1{r}_2}]\\ & -\frac{1}{2}[\left(\frac{\rho }{r_1}\right)\left(\frac{z-d}{r_1}\right)\left(\frac{i{k}_1^2}{r_1}-\frac{3k_1}{r_1^2}-\frac{3i}{r_1^3}\right)e^{i{k}_1{r}_1}\\ & +\left(\frac{\rho }{r_2}\right)\left(\frac{z+d}{r_2}\right)\left(\frac{i{k}_1^2}{r_2}-\frac{3k_1}{r_2^2}-\frac{3i}{r_2^3}\right)e^{i{k}_1{r}_2}\left]\right\},\end{array}$

(82)$\begin{array}{ll}& B_{1\rho }(\rho ,\varphi ,z)\\ \hfill =& -\frac{{\mu }_0}{2\pi k_1}sin\varphi \{k_{2}h(\rho ;k_1,k_2)e^{i{k}_2\rho }{e}^{i{k}_1(z+d)}\\ & +\left[\frac{i{k}_1}{{\rho }^2}-\frac{7}{2{\rho }^3}+\left(\frac{z+d}{r_2}\right)\left(\frac{i{k}_1^2}{r_2}-\frac{5k_1}{r_2^2}-\frac{12i}{r_2^3}\right)\right]e^{i{k}_1{r}_2}\\ & \frac{1}{2}[\left(\frac{z-d}{r_1}\right)\left(\frac{i{k}_1^2}{r_1}-\frac{k_1}{r_1^2}\right)e^{i{k}_1{r}_1}\\ & +\left(\frac{z+d}{r_2}\right)\left(\frac{i{k}_1^2}{r_2}-\frac{k_1}{r_2^2}\right)e^{i{k}_1{r}_2}\left]\right\},\end{array}$

.

(83)$\begin{array}{ll}& B_{1\varphi }(\rho ,\varphi ,z)\\ \hfill =& -\frac{{\mu }_0}{2\pi k_1}cos\varphi \{k_{2}g(\rho ;k_1,k_2)e^{i{k}_2\rho }{e}^{i{k}_1(z+d)}\\ & +\frac{1}{2}[\left(\frac{2}{{\rho }^3}+\frac{3i}{k_1{\rho }^4}\right)e^{i{k}_1{r}_2}-\left(\frac{z-d}{r_1}\right)\left(\frac{i{k}_1^2}{r_1}-\frac{k_1}{r_1^2}\right)e^{i{k}_1{r}_1}\\ & -\left(\frac{z+d}{r_2}\right)\left(\frac{i{k}_1^2}{r_2}-\frac{3k_1}{r_2^2}-\frac{6i}{r_2^3}\right)e^{i{k}_1{r}_2}\left]\right\},\end{array}$

(84)$\begin{array}{ll}& B_{1z}(\rho ,\varphi ,z)\\ \hfill =& \frac{{\mu }_0}{2\pi }sin\varphi \{\frac{k_2^2}{k_1^2}\left(\frac{1}{{\rho }^2}+\frac{3i}{k_2{\rho }^3}-\frac{3}{k_2^2{\rho }^4}\right)e^{k_2\rho }{e}^{i{k}_1(z+d)}\\ & -\left(\frac{\rho }{r_2}\right)[\frac{1}{r_2^2}+\frac{3i}{k_1{r}_2^3}-\frac{3}{k_1^2{r}_2^4}\\ & +\left(\frac{z+d}{\rho }\right)\left(\frac{k_1}{\rho }+\frac{3i}{2{\rho }^2}-\frac{9}{8k_1{\rho }^3}\right)\\ & +{\left(\frac{z+d}{r_2}\right)}^2\left(\frac{i{k}_1}{r_2}-\frac{6}{r_2^2}-\frac{15i}{k_1{r}_2^3}\right)]e^{i{k}_1{r}_2}\\ & -\frac{1}{2}[\left(\frac{\rho }{r_1}\right)\left(\frac{i{k}_1}{r_1}-\frac{1}{r_1^2}\right)e^{i{k}_1{r}_1}\\ & -\left(\frac{\rho }{r_2}\right)\left(\frac{i{k}_1}{r_2}-\frac{1}{r_2^2}\right)e^{i{k}_1{r}_2}\left]\right\}.\end{array}$

In these formulas,

(85)$\begin{array}{ccc}\hfill f(\rho ;k_1,k_2)& =& \frac{i{k}_2}{\rho }-\frac{1}{{\rho }^2}-\frac{k_2^3}{k_1}\sqrt{\frac{\pi }{k_2\rho }}{e}^{-ip}\mathcal{F}\left(p\right),\hfill \end{array}$

(86)$\begin{array}{ccc}\hfill g(\rho ;k_1,k_2)& =& \frac{i{k}_2}{\rho }-\frac{1}{{\rho }^2}-\frac{i}{k_2{\rho }^3}-\frac{k_2^3}{k_1}\sqrt{\frac{\pi }{k_2\rho }}{e}^{-ip}\mathcal{F}\left(p\right),\hfill \end{array}$

(87)$\begin{array}{ccc}\hfill h(\rho ;k_1,k_2)& =& \frac{2}{{\rho }^2}+\frac{2i}{k_2{\rho }^3}+\frac{i{k}_2^2}{k_1\rho }\sqrt{\frac{\pi }{k_2\rho }}{e}^{-ip}\mathcal{F}\left(p\right),\hfill \end{array}$

where:

(88)$\mathcal{F}\left(p\right)=\frac{1}{2}(1+i)-C_2\left(p\right)-i{S}_2\left(p\right),p\equiv \frac{k_2^3\rho }{2k_1^2},$

and:

(89)$C_2\left(p\right)+i{S}_2\left(p\right)={\int }_0^p\frac{e^{it}}{\sqrt{2\pi t}}\mathrm{d}t$

is the Fresnel integral.

In the simplified expressions, the direct field is given by all terms multiplied by $e^{i{k}_1{r}_1}$; the complete reflected field includes all terms with the factor $e^{i{k}_1{r}_2}$; and the lateral-wave field consists of the terms multiplied by $e^{i{k}_2\rho }{e}^{i{k}_1(z+d)}$. The quantities $e^{i{k}_1{r}_1}$ and $e^{i{k}_1{r}_2}$ represent spherical waves traveling radially outward from the dipole at $(0,0,d)$ or the image dipole at $(0,0,-d)$. The quantity $e^{i{k}_2\rho }{e}^{i{k}_1(z+d)}$ suggests a plane wave $e^{i{k}_{1}d}$ that travels upward in Region 1 from the source at $(0,0,d)$ to the boundary at $(0,0,0)$, then radially outward in Region 2 as a cylindrical wave $e^{i{k}_2\rho }$ to $(\rho ,0,0)$, and finally, vertically downward from the boundary in Region 1 as a plane wave $e^{i{k}_{1}z}$ to the point of observation at $(\rho ,0,z)$.

5. Model Evaluations and Results

In order to show the path loss in subsurface wireless underground communications, model evaluations were done, and results are shown for a horizontal electric dipole in the presence of a plane boundary. The Volumetric Water Content (VWC) was 20% (obtained through the complex permittivity $10.8+2.4i$ [46,47]). Due to the lossy characteristics of the soil medium, it is represented through complex permittivity. Moreover, because the soil propagation medium is non-magnetic, the permeability of the free space (${\mu }_0$) was used (Section 3). A total of five transmitter burial depths depths of 10 cm, 20 cm, 30 cm, 40 cm, and 50 cm were considered for the evaluation of an x-directed horizontal electric dipole. The lateral wave power density magnitude is shown as a function of distance. The communication distances up to 5 m were evaluated at the frequency of 433 MHz. From all six components of the lateral signal (e.g., three electric and three magnetic components), the results are only shown for the components $E_z$ and $E_x$ ($E_{\rho }$ with $\varphi =0$) in the subsurface region.

The comparison of the power density of $E_x$ and $E_z$ components at different burial depths is shown in Figure 4. It can be observed that at a 10-cm depth (Figure 4a), the power density of the horizontal component $E_x$, which is the lateral wave component $E_{\rho }$($\rho$,0,z), was on par with the $E_z$ component. However, for distances greater than 1 m, the horizontal component became 11 dB weaker as compared to the $E_z$. This happened because of the higher attenuation of the $E_x$ component in the medium. This became clearer when the depth increased to 30 cm (Figure 4b). At a 30-cm depth, the lateral component had a strong power density at shorter distances (<1 m). With the increase in distance, the horizontal component decreased gradually. Similar results have also been obtained in [10], and those empirical results showing the dominance of the lateral wave confirm the accuracy of this underground propagation prediction model. Moreover, the model accurately captured the physical propagation differences due to the propagation of different components.

The results of the numerical evaluations of the signal strength for different depths for both the $E_z$ and $E_z$ components are shown in Figure 5. The abrupt drops in the signal strength can be observed for a distance value of less than 2 m for both components. This happened because when the wave traveled in the soil medium, then a cross-sectional plane was formed at a $\lambda =2\pi /k$ wavelength. At this cross-sectional plane, the intensity of the signal was at the positive maximum value. However, half-way, the reversed direction of the electric field vector was at the negative maximum value, which led to these drops in the power density. At the higher distances, this effect became less significant because of the low signal strength caused by the higher attenuation of the wave traveling in the subsurface medium. Moreover, it can also be observed that the signal strength of the dominant lateral wave decreased with the increase in the depth and distance. The developed model can also capture the water sensitivity of the soil medium.

6. Applications

For an underground communications’ system design, the propagation characteristics of the subsurface soil environment should be determined accurately. This knowledge is also important for underground network design in the Internet of Underground Things (IOUT) in precision agriculture. For accurate soil moisture sensing, the location of the placement of underground nodes plays an important role in node deployment. The analysis presented in the paper can be used to determine the path loss between the underground Transmitter-Receiver pair (T-R), which helps to achieve the maximum coverage area in the crop field with optimal deployment. This model also provides information about the ideal depths in the underground medium. The model output also aids in selecting optimal communication parameters (e.g., modulation scheme and transmit power) for the wireless channel in a particular environment. Without the use of the model, this path loss information can only be obtained through empirical evaluations in the field. These field measurements are site specific, harder to implement in complex field environments, and lead to inefficient and delayed system deployment.

7. Conclusions

In this article, a physical electromagnetic analysis related to subsurface propagation modeling, which considers the electromagnetic field of a unit vertical and horizontal electric dipole in the presence of a plane boundary, has been carried out. It has been shown that the lateral wave is the dominant wave in subsurface communications, and shallow depths are more suitable for long-range communications and sensing. This analysis forms the basis of subsurface radio wave propagation modeling in the stratified soil medium for smart farming applications in precision agriculture [13,20,48,49,50,51] and urban underground infrastructure monitoring [52,53].

Furthermore, we highlight that the empirical evaluations of UG solutions under realistic application requirements is essential due to the interplay between the physical soil medium and radios. The main contribution of this paper is to present the physics-based subsurface propagation prediction model. To this end, we have evaluated the results and have shown the agreement with the empirical channel measurement results obtained through extensive impulse response measurements in an indoor and field testbed. The change of the soil moisture impacts the signal propagation in the underground soil medium [37]. The increase in VWC leads to a decrease in the received signal strength at the underground receiver [10]. The developed approach works by ascertaining the VWC percentage by using underground sensors and accordingly using Peplinski’s model [46] to obtain the complex permittivity of the soil under consideration. The obtained complex permittivity value for a given VWC level can accordingly be inserted into the model along with other soil parameters (e.g., soil type, soil bulk density, burial depth in medium, and communication distance) to predict radio wave propagation for a particular agricultural field. Accordingly, the subsurface communication system can be tailored for site-specific and soil-specific applications. Overall, the developed model can be effectively used to predict propagation characteristics in the soil medium with a change in soil moisture for digital agriculture applications. For future work, we will perform the evaluation of the developed solutions through both unit-level experimentation and simulations and application-level evaluations. The wide variety of applications of UG communications including agriculture, soil contaminants detection, mining, transportation, and border patrol results in a diverse set of requirements that can benefit from this approach.

Funding

This research received no external funding.

Acknowledgments

Publication of this article was funded in part by the Purdue University Libraries Open Access Publishing Fund.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

Notation	Description
$k_s$	wavenumber in the soil
$k_0$	wavenumber in the air
${\lambda }_s$	wavelength in the soil
${\beta }_s$	phase constant
${\alpha }_s$	attenuation constant
${ϵ}_s$	permittivity of soil
${\mu }_0$	permeability of free space
l	half length of antenna
$\omega$	angular frequency
$Z_a$	antenna impedance
$Z_s$	transmission line impedance
$E_r$	reflected electric field
$I_r$	reflected current
$I_m$	current amplitude
$Z_r$	reflected impedance
$Z_0$	characteristics impedance
$R_a$	radiation resistance
$I_0\left(\zeta \right)$	current distribution along the antenna
$Z_a^u$	soil-air interface adjusted impedance
$RL$	return loss
$f_r$	resonant frequency
$f_0$	over-the-air frequency
$BW$	bandwidth
h	burial depth
${\eta }_a$	refractive index of air
${\eta }_s$	refractive index of soil
$m_v$	volumetric water content
${\rho }_b$	bulk density
${\rho }_s$	particle density
$\delta$, ${\nu }^{\prime }$, and ${\nu }^{″}$	empirical constant
${ϵ}_s^{{}^{\prime }}$,${ϵ}_s^{{}^{″}}$	real and imaginary part of the permittivity of soil
S and C	sand and clay particle percentage
${ϵ}_{fw}^{\prime }$ and ${ϵ}_{fw}^{″}$	relative permittivity of free water, real and imaginary
${\sigma }_{eff}$	effective conductivity of soil
${ϵ}_{w\infty }$,	limit of permittivity of water
${ϵ}_{w0}$	static permittivity of water
${\tau }_w$	relaxation time of water

Figures

Figure 1. The drawing of the geometry of the underground antenna in the soil propagation medium: (a) schematic; (b) coordinate system.

Figure 2. The coordinates for a lateral wave at the surface, ϕ=0.

Figure 3. The burial depth and the z-axis.

Figure 4. The comparison of the power density of Ex and Ez components at different burial depths: (a) 10-cm depth; (b) 30-cm depth.

Figure 5. The power density for different burial depths: (a) Ex component; (b) Ez component.