Scattering and Attenuation in 5G Electromagnetic Propagation (5 GHz and 25 GHz) in the Presence of Rainfall: A Numerical Study

Abstract

Rainfall has always been a concern for wireless communications systems. As 5G technology relies on high-frequency bands, it is fundamental to model and simulate the interaction of such radio waves with rainfall, as the deployment of large-scale infrastructure for 5G is highly expensive. This research presents a reformulation of the Maxwell equations for a bi-dimensional space in a transverse electric propagation mode, for a linear, inhomogeneous, and isotropic propagation medium with its magnetic and electric properties dependent on time. This reformulation was solved using the Finite Differences in Time Domain (FDTD) method with the Convolutional Perfectly Matched Layer (CPML) boundary condition. Two main frequency propagation scenarios were studied: 5 GHz (corresponding to Wi-Fi in the 802.11n standard as well as to the lowest bands of 5G) and 25 GHz (corresponding to 5G), within a $10\mathrm{m}\times 3\mathrm{m}$ rectangular domain in air and with rain. The rainfall was simulated using a parallel Ziggurat algorithm. According to the findings, while 5 GHz waves experience scattering processes, 25 GHz waves experience substantial dispersion and attenuation throughout the domain in low- to moderate-intensity rain.

1. Introduction

In the current context of telecommunications, the evolution towards a fully networked environment is becoming more and more noticeable; an example of this is the concept known as IoT (Internet of Things) [1,2] and its interoperability with traditional communications systems of any scale, which day after day is leading us to the so-called smart world [3].

With the emergence of autonomous systems reliant on real-time transmission, communication between systems, and the digitalisation of services, among others, massive data to be transmitted at high transfer data rates arise. This set of demanding requirements is being met by the leap to the fifth generation of mobile communications or 5G [4,5,6], which consists of applying new data transmission technologies at high power waves between 2.4 GHz and 300 GHz in wireless telecommunication networks [5,7].

As the deployment of large-scale 5G infrastructure for the purpose of validation of the technology is extremely expensive, modelling and simulation play a fundamental role in the analysis of multiple implementation issues and potential project solutions, mitigating costs and developing time in prototypes and physical tests, as well as improving the development of the necessary infrastructure.

In this sense, much of the current research is devoted to the modelling and simulation of different aspects mostly related to the technology’s performance. Initially, as a 5G propagated wavelength is very short compared with the propagated distances, the ray tracing method became popular to model 5G paths. For instance, Hou et al. [8] applied it to analyse the effects between a 5G wireless communication system and the industrial environment, including the present equipment, while Hsiao et al. [9] found that since the attenuation and blockage of millimetre-wave propagation in urban environments are severe, Line-of-Sight (LoS) and Non-Line-of-Sight (NLoS) reflected rays should be the dominant propagation mechanisms.

Further research has specialised in the study of propagation path loss models for 5G frequencies indoors and outdoors. For indoor environments, Samad et al. [10] studied the propagation properties of LoS links at 3.7 and 28 GHz within long corridors in Korea to model and simulate the path loss, while Alabdullah et al. [11] studied the performance of 28, 39, 60, and 73 GHz waves in both LoS and NLoS scenarios in indoor environments over Tx-Rx separations of 1.5 m to 62 m with both measurements and software. Muttair et al. [12] simulated wave propagation including path loss, delay spread, and received power for outdoor 28, 39, 60, and 73 GHz LoS and NLoS links through Wireless InSite^® program, finding that LoS paths have a high receiving capacity and fewer path losses than NLoS ones, while high frequencies such as 73 GHz have greater effects on propagation than low frequencies. Muttair et al. [13] studied 10, 17, 30, and 60 GHz links for outdoor to indoor antennas, using measurement analysis and ray tracing simulation, and found that an increase in frequency yields an increase in path loss and a decrease in the received signal strength (RSS), the delay spread, and the received power. Oladimeji et al. [14] provided a comprehensive review of propagation path loss prediction in enclosed environments for 5G networks. Casillas-Aviña et al. [15] implemented RF path loss models for cellular 5G links in a practical free link budget calculator. A further review on advanced simulation methods for 5G antennas and propagation can be found in [14,16]. Other methods, such as those based on machine learning, have become popular for 5G [17,18,19,20,21,22,23].

Of particular concern to telecommunications developers is that several atmospheric conditions can affect free space wave propagation, which can complicate and degrade wireless communication links. Some of the effects with the most significant impact are the presence of aerosols, the concentration of gases, cloudiness, and precipitation [24].

Although the impacts of these factors on electromagnetic propagation have been studied for many years [25], it has been determined that aerosols primarily affect the optical atmospheric links in cases where the trajectory of the particles is noticeable [26]. In addition, the repercussions of aerosols and airborne particles are highly complex because of their great diversity, size, geometry, mass, chemical interactions, etc., and require complex and costly mathematical algorithms for their modelling [27,28]. In particular, fine dust often affects low-frequency propagation due to the change in the refractive index of the atmosphere, causing power loss and multi-paths. However, these effects are representative only in extreme cases such as dust storms or sandstorms [29].

Moreover, gas concentrations in the atmosphere can adversely affect propagation. Since nitrogen and oxygen make up a high percentage of the total gases in the troposphere, the information sources of radio-based communication systems are designed away from the absorption bands of mono-atomic nitrogen and oxygen, so the gases in concentration in the atmosphere that can affect wireless links are dry air (molecular oxygen) and water vapour [30]. These components cause wave fluctuations that tend to increase with frequency. The peak absorption for water vapour (cloudiness) is found at 22.235 GHz and for molecular oxygen or dry air at 61.100 GHz [31,32]. However, both frequencies were declared extremely important for radio astronomy at the 21st General Assembly of the International Astronomical Union (IAU), stating the need to protect these bands from anthropogenic emissions, mainly from space transceivers. In this way, radio communication systems are automatically shielded from these effects [31].

On the other hand, precipitation has one of the most notable effects on terrestrial wireless communications networks. Precipitation, in its broadest sense, is any concentration of water or hydrometeor produced in the atmosphere and falling to the planet’s surface [33]. Depending on the water concentration and temperature, such hydrometeors may take the form of rain, snow, hail, or fog. In the case of free space electromagnetic propagation, precipitation in the form of rain has the most significant impact [34], as it causes an attenuation directly proportional to the frequency; so in a wave with a frequency above 10 GHz, it can have a very representative effect [35,36]. In addition to the rainfall intensity rate, the impact depends on the droplets’ shape, size, and distribution. In this manner, high-frequency electromagnetic waves sustain attenuation and dispersion problems of considerable magnitude in the presence of hydrometeor phenomena [35,36] due to diffraction in the medium and shorter penetration lengths in the materials, increasing both the importance of direct transmission to the user (LoS propagation), as well as the development of a more extensive infrastructure to add redundancy and integrity to the network, concerning the case of 4G [7]. Thus, the analysis and study of the atmospheric effects on the medium for the propagation of electromagnetic waves for 5G technology have become fundamental for its successful local and global implementation [37].

The response of an electromagnetic field to a time-dependent atmospheric environment, as that represented by the presence of hydrometeors, is fundamental for understanding the magnitude of the reflection, refraction, scattering, and attenuation processes present in propagation at these operating frequencies. Studies suggest that the most significant effect is attenuation due to hydrometeors being the most disruptive to high-frequency telecommunication systems [38].

In the case of low-frequency wireless systems such as 5 GHz (Wi-Fi and low band of 5G), modelling of the behaviour of electromagnetic waves in different environments such as rain, snow, sand, and forests shows a substantial impact on the transmissibility of electromagnetic waves, which is observed to cause scattering [36,39], while more recent studies show that 5G is particularly vulnerable to rain scattering [40,41] for ultra-high frequencies (20.2, 39.4, 73.0, and 83.0 GHz).

Although different models have been developed to study the impact of precipitation on wireless telecommunications [42,43], the results are insufficient in a field with a rapidly evolving technology as seen currently. Thus, the core part of this work is to identify the effects of precipitation in the form of rain in the main transmission rate band for the emerging 5G technology. In this work, we seek to provide a quantitative solution to this problem by employing a computational simulation of these phenomena based on Maxwell’s electromagnetic theory [44], considering an appropriate treatment of the constitutive relations that allow us to simulate the hydrometeoric phenomena, which are time-dependent. Among the different numerical algorithms used for the solution of hyperbolic partial differential equations, such as Maxwell’s equations, the Finite Difference Time Domain (FDTD) method, based on the application of Yee’s step cell [45], stands out for its ease of implementation as well as for the reliability of its results.

Regarding 5G modelling and simulation, the FDTD method has recently been applied. For instance, Gorniak [46] developed an effective FDTD method in order to be able to simulate stochastic electromagnetic fields in the frequency band of 5G, while Asif et al. [47] combined the FDTD and the Finite Element Method (FEM) to design a MIMO antenna for 5G in the context of cellular phones. Moreover, a large amount of concern has arisen regarding the impact of 5G on human health, which has been studied with (mostly parallel) FDTD methods. Jariyanorawiss and Chongburee [48] reported the effects of the 2.6 GHz Mid-Band of 5G on different human head exposures, while Yoshida et al. [49] studied the shadowing generated by the human body for 5G indoor propagation. Moreover, Yoshida et al. [50], estimated the propagation loss generated by the human body in 5G frequencies.

This research tackles the problem of electromagnetic wave propagation in two different frequencies in the presence of hydrometeors, in particular rainfall, by developing a numerical algorithm based on the second-order of precision in the space and time FDTD, with Convolutional Perfectly Matched Layer (CPML)-type absorbing boundary conditions [51], to solve Maxwell’s equations considering constitutive relations for linear, isotropic, and inhomogeneous time-dependent media. To simulate rain, a parallel version of the Ziggurat algorithm was used to generate pseudo-random numbers. In this way, the effects of rain on electromagnetic wave propagation were studied for both sources at 5 GHz and 25 GHz. Furthermore, the generated code was accelerated at a high level by introducing directives in the serial code, which allow its parallelisation. In this case, the FDTD algorithm was accelerated with OpenACC directives (in GPU), while the Ziggurat problem was accelerated through OpenMP. The main motivation of this work is to provide a numerical study of the impact of sweet water rainfall on the propagation of the 5 and 25 GHz bands of 5G technology for 5G developers to be aware of the associated drawbacks, and to be able to provide a more robust and extensive infrastructure.

This paper is structured as follows: Section 2 presents the development of Maxwell’s equations for the time-dependent constitutive relations and their development in finite differences, the implementation of the CPML boundary condition, the development of the case studies, and the designed computational performance tests. Section 3 presents the multiplicity of results obtained from the simulations for both propagation frequencies in free space and the presence of rain, while Section 4 discusses the results, both physically and computationally. Finally, some conclusive aspects of this work are presented in Section 5, while Appendix A briefly reviews the development of the Ziggurat pseudo-random number sampling algorithm.

2. Materials and Methods

2.1. Maxwell’s Equations for Time-Dependent Media

The classical electromagnetic phenomena in material media are governed by Maxwell’s equations, which in their most general differential or strong form in the International System of Units take the form [44]

$$\overrightarrow{\nabla }·\overrightarrow{D}={\rho }^{(\mathrm{free})},$$

(1)

$$\overrightarrow{\nabla }·\overrightarrow{B}=0,$$

(2)

$$\overrightarrow{\nabla }\times \overrightarrow{E}+\frac{\partial \overrightarrow{B}}{\partial t}=0$$

(3)

and

$$\overrightarrow{\nabla }\times \overrightarrow{H}-\frac{\partial \overrightarrow{D}}{\partial t}=\overrightarrow{J},$$

(4)

where $\overrightarrow{D}$ represents the vector of electric displacement that provides a measure of the polarisation of the medium; $\overrightarrow{H}$ represents the magnetic response vector of the material; $\overrightarrow{B}$ and $\overrightarrow{E}$ are the magnetic and electric field, respectively. Moreover, ${\rho }^{\left(free\right)}$ and $\overrightarrow{J}$ are the electric charge and current densities, respectively. However, the current density in the presence of conducting media is divided into the contribution of free and induced current in the conductor, ${\overrightarrow{J}}^{\left(\mathrm{free}\right)}$ and ${\overrightarrow{J}}^{\left(\mathrm{ind}\right)}$, respectively, so the electric current density takes the form [44]

$$\overrightarrow{J}={\overrightarrow{J}}^{\left(\mathrm{free}\right)}+{\overrightarrow{J}}^{\left(\mathrm{ind}\right)}.$$

(5)

The set of non-closed Equations (1)–(4) need to be provided with constitutive relations to exhibit a unique solution. Constitutive relations determine the response of the medium according to the applied fields, and in a certain approximation, they take the form $\overrightarrow{D}=\overrightarrow{D}\left(\overrightarrow{E}\right)$, $\overrightarrow{H}=\overrightarrow{H}\left(\overrightarrow{B}\right)$, and ${\overrightarrow{J}}^{\left(\mathrm{ind}\right)}={\overrightarrow{J}}^{\left(\mathrm{ind}\right)}\left(\overrightarrow{E}\right)$. In this paper, we proposed constitutive relations for linear, homogeneous and isotropic time-dependent media, which are given by

$$\overrightarrow{D}(\overrightarrow{r},t)=ϵ(\overrightarrow{r},t)\overrightarrow{E}(\overrightarrow{r},t),$$

(6)

$$\overrightarrow{H}(\overrightarrow{r},t)=\frac{1}{\mu (\overrightarrow{r},t)}\overrightarrow{B}(\overrightarrow{r},t)$$

(7)

and

$${\overrightarrow{J}}^{\left(ind\right)}(\overrightarrow{r},t)=\sigma (\overrightarrow{r},t)\overrightarrow{E}(\overrightarrow{r},t),$$

(8)

where the proportionality coefficients $ϵ(\overrightarrow{r},t)$, $1/\mu (\overrightarrow{r},t)$, and $\sigma (\overrightarrow{r},t)$ are the permittivity, the inverse of the permeability, and the electrical conductivity of the propagation media, respectively, and are functions of position and time. Then, by substituting such constitutive relations to the set of Maxwell’s equations in their vector components, we obtain the set of Equations (9)–(16).

From Equation (1), we obtain

$$\frac{\partial \left[ϵ(\overrightarrow{r},t)E_x\right]}{\partial x}+\frac{\partial \left[ϵ(\overrightarrow{r},t)E_y\right]}{\partial y}+\frac{\partial \left[ϵ(\overrightarrow{r},t)E_z\right]}{\partial z}={\rho }^{\left(free\right)};$$

(9)

from Equation (2),

$$\frac{\partial \left[\mu (\overrightarrow{r},t)H_x\right]}{\partial x}+\frac{\partial \left[\mu (\overrightarrow{r},t)H_y\right]}{\partial y}+\frac{\partial \left[\mu (\overrightarrow{r},t)H_z\right]}{\partial z}=0;$$

(10)

from Equation (3),

$$\frac{\partial H_x}{\partial t}=\frac{1}{\mu (\overrightarrow{r},t)}\left[\frac{\partial E_y}{\partial z}-\frac{\partial E_z}{\partial y}-H_x\frac{\partial \left(\mu (\overrightarrow{r},t)\right)}{\partial t}\right],$$

(11)

$$\frac{\partial H_y}{\partial t}=\frac{1}{\mu (\overrightarrow{r},t)}\left[\frac{\partial E_z}{\partial x}-\frac{\partial E_x}{\partial z}-H_y\frac{\partial \left(\mu (\overrightarrow{r},t)\right)}{\partial t}\right],$$

(12)

and

$$\frac{\partial H_z}{\partial t}=\frac{1}{\mu (\overrightarrow{r},t)}\left[\frac{\partial E_x}{\partial y}-\frac{\partial E_y}{\partial x}-H_z\frac{\partial \left(\mu (\overrightarrow{r},t)\right)}{\partial t}\right];$$

(13)

and finally, from Equation (4),

$$\frac{\partial E_x}{\partial t}=\frac{1}{ϵ(\overrightarrow{r},t)}\left[\frac{\partial H_z}{\partial y}-\frac{\partial H_y}{\partial z}-E_x\frac{(\partial ϵ(\overrightarrow{r},t))}{\partial t}-{\overrightarrow{J}}_x^{\left(free\right)}-\sigma (\overrightarrow{r},t)E_x\right],$$

(14)

$$\frac{\partial E_y}{\partial t}=\frac{1}{ϵ(\overrightarrow{r},t)}\left[\frac{\partial H_x}{\partial z}-\frac{\partial H_z}{\partial x}-E_y\frac{\partial \left(ϵ(\overrightarrow{r},t)\right)}{\partial t}-{\overrightarrow{J}}_y^{\left(free\right)}-\sigma (\overrightarrow{r},t)E_y\right],$$

(15)

and

$$\frac{\partial E_z}{\partial t}=\frac{1}{ϵ(\overrightarrow{r},t)}\left[\frac{\partial H_y}{\partial x}-\frac{\partial H_x}{\partial y}-E_z\frac{\partial \left(ϵ(\overrightarrow{r},t)\right)}{\partial t}-{\overrightarrow{J}}_z^{\left(free\right)}-\sigma (\overrightarrow{r},t)E_z\right].$$

(16)

For reasons of computational capacity, the system represented by Equations (9)–(16) is reduced to two dimensions. For this purpose, if we let the propagation plane be the XY plane, we restrict the electromagnetic energy flux to such a plane, which requires that the $\widehat{z}$ component of the Poynting vector be zero,

$${\overrightarrow{S}=\overrightarrow{E}\times \overrightarrow{H}=\left[E_y{H}_z-E_z{H}_y\right]\widehat{x}-\left[E_x{H}_z-E_z{H}_x\right]\widehat{y}+\overline{)\left[E_x{H}_y-E_y{H}_x\right]}\widehat{z}.}^0$$

(17)

Thus, the system is obliged to comply with the condition

$$E_x{H}_y=E_y{H}_x.$$

(18)

If, additionally, the electric transverse mode (TE) is imposed, then the electric field components in the direction of field propagation are zero, so Equations (9)–(16) take the form

$$\frac{\partial H_x}{\partial t}=-\frac{1}{\mu (\overrightarrow{r},t)}\left[\frac{\partial E_z}{\partial y}+H_x\frac{\partial \mu (\overrightarrow{r},t)}{\partial t}\right],$$

(19)

$$\frac{\partial H_y}{\partial t}=\frac{1}{\mu (\overrightarrow{r},t)}\left[\frac{\partial E_z}{\partial x}-H_y\frac{\partial \mu (\overrightarrow{r},t)}{\partial t}\right]$$

(20)

and

$$\frac{\partial E_z}{\partial t}=\frac{1}{ϵ(\overrightarrow{r},t)}\left[\frac{\partial H_y}{\partial x}-\frac{\partial H_x}{\partial y}-E_z\frac{\partial ϵ(\overrightarrow{r},t)}{\partial t}-J_z^{\left(free\right)}-\sigma (\overrightarrow{r},t)E_z\right],$$

(21)

which is the final set of equations to be implemented in the simulations of this work.

2.2. Numerical Modelling

Equations (19)–(21) are discretised using the classical FDTD method of second order in space and time, obtaining the algebraic Equations (22)–(24). In Figure 1, the calculation molecule of $E_z$, $H_x$, and $H_y$ in space and time is depicted.

$$H{x}_{(i,j-1/2)}^{n+1/2}=H{x}_{(i,j-1/2)}^{n-1/2}\left(\frac{{\mu }_{(i,j-1/2)}^{n-1/2}}{{\mu }_{(i,j-1/2)}^{n+1/2}}\right)-\frac{\Delta t}{{\mu }_{(i,j-1/2)}^{n+1/2}}\left(\frac{E{z}_{(i,j)}^n-E{z}_{(i,j-1)}^n}{\Delta y}\right),$$

(22)

$$H{y}_{(i-1/2,j)}^{n+1/2}=H{y}_{(i-1/2,j)}^{n-1/2}\left(\frac{{\mu }_{(i-1/2,j)}^{n-1/2}}{{\mu }_{(i-1/2,j)}^{n+1/2}}\right)+\left(\frac{\Delta t}{{\mu }_{(i-1/2,j)}^{n+1/2}}\right)\left(\frac{E{z}_{(i,j)}^n-E{z}_{(i-1,j)}^n}{\Delta x}\right)$$

(23)

and

$$E{z}_{(i,j)}^{n+1}=E{z}_{(i,j)}^n\left(\frac{4{ϵ}_{(i,j)}^n-\Delta t\left({\sigma }_{(i,j)}^{n+1}+{\sigma }_{(i,j)}^n\right)}{4{ϵ}_{(i,j)}^{n+1}+\Delta t\left({\sigma }_{(i,j)}^{n+1}+{\sigma }_{(i,j)}^n\right)}\right)+\frac{\Delta t}{{ϵ}_{(i,j)}^{n+1}+\frac{\Delta t}{4}\left({\sigma }_{(i,j)}^{n+1}+{\sigma }_{(i,j)}^n\right)}$$

(24)

$$\left[\frac{H{y}_{(i+1/2,j)}^{n-1/2}-H{y}_{(i-1/2,j)}^{n-1/2}}{\Delta x}-\frac{H{x}_{(i,j+1/2)}^{n-1/2}-H{x}_{(i,j-1/2)}^{n-1/2}}{\Delta y}-\frac{J{z}_{(i,j)}^{n+1}+J{z}_{(i,i)}^n}{2}\right],$$

where $i,j,n$ denote the discrete position corresponding to the coordinates $x,y$, and t, respectively, while $\Delta x$, $\Delta y$, and $\Delta t$ are the FDTD grid size in the x, y, and t variables, respectively.

The time step must satisfy the Courant–Friedrichs–Lewy (CFL) criterion $\Delta t\le {\displaystyle \frac{\Delta x}{c_0}}\sqrt{{\displaystyle \frac{{\epsilon }_r{\mu }_r}{2}}}$ to ensure numerical stability [52,53], where $c_0$ is the speed of light in a vacuum, and ${\epsilon }_r$ and ${\mu }_r$ are the propagating medium’s relative permittivity and permeability, respectively. Several numerical experiments in the context of electromagnetic propagation have shown that a relation of $\Delta x=\Delta y=\lambda /20$, where $\lambda$ is the propagated wavelength, allows the quality of the numerical solution to be maintained [53]. That is, when setting the computational grid, there should be at least 20 spatial grid cells of the FDTD method for an accurate resolution of one wavelength of the signal.

The absorbing boundary condition implemented is the Convolutional PML (CPML) [51] in its simplified version by Martin and Komatisch [54]; this version of the Convolutional PML saves memory and is relatively straightforward to code. To implement the CPML boundary conditions, it is necessary to introduce the Luebbers convolutional memory variables [55] for each spatial derivative. The convolutional term $\mathrm{\Psi }$ is introduced in Equations (22)–(24) as

$${\mathrm{\Psi }}_i^n=\sum _{m=0}^{n-1}\frac{\partial }{\partial x}{|}_{i\Delta x}^{\left[n-(m+\frac{1}{2})\right]\Delta t}Z\left(m\right),$$

(25)

where

$$Z\left(m\right)=\left(\frac{d}{\kappa d+\alpha {\kappa }^2}\right)\left\{e^{-\left(\frac{d}{\kappa }+\alpha \right)(m+1)\Delta t}-e^{-\left(\frac{d}{\kappa }+\alpha \right)m\Delta t}\right\},$$

(26)

and where $\alpha \left(x\right),\kappa \left(x\right)$, and $d\left(x\right)$ are the CPML parameters. If, additionally, we define the functions

$$b\left(q\right)=e^{-\{\frac{d\left(q\right)}{\kappa }\left(q\right)+\alpha \left(q\right)\}\Delta \left(t\right)}$$

(27)

and

$$a\left(q\right)=\left(\frac{d\left(q\right)}{d\left(q\right)+\alpha \left(q\right)\kappa {\left(q\right)}^2}\right)\{b\left(q\right)-1\},$$

(28)

the convolutional variable in Equation (25) for the fields $E_z$ and $H_{\nu }$, with $\nu =x,y$ in Equations (22)–(24), takes the form

$${\mathrm{\Psi }}_i^{n+1}\left[H_{\nu }\right]=b{\mathrm{\Psi }}_i^n\left[H_{\nu }\right]+a\left(\frac{{\left(H_{\nu }\right)}_{i+\frac{1}{2}}^{n+1/2}-{\left(H_{\nu }\right)}_{i-\frac{1}{2}}^{n+1/2}}{\Delta x}\right)$$

(29)

and

$${\mathrm{\Psi }}_{i-\frac{1}{2}}^{n+\frac{1}{2}}\left[E_z\right]=b{\mathrm{\Psi }}_{i-\frac{1}{2}}^{n-\frac{1}{2}}\left[E_z\right]+a\left(\frac{{\left(E_z\right)}_i^n-{\left(E_z\right)}_{i-1}^n}{\Delta x}\right).$$

(30)

Following the numerical experiments’ development in [53], where $\alpha \left(x\right),\kappa \left(x\right)$, and $d\left(x\right)$ are studied as polynomial functions in the domain of the CPML boundary, and $\mathcal{N}$ is the number of layers, thus $T_h=\mathcal{N}\Delta x$ is the thickness of the CPML. Then, the following equations are defined as

$$x{\left(q\right)}_{\mathrm{norm}}=T_h-q\frac{\Delta x}{\mathcal{N}}=\left(\mathcal{N}-\frac{q}{\mathcal{N}}\right)\Delta x,$$

(31)

$$\mathrm{d}\left(q\right)\equiv {\mathrm{d}}_{\max }{x}_{\mathrm{norm}}^{\eta },$$

(32)

$$\kappa \left(q\right)\equiv 1+({\kappa }_{\max }-1)x_{\mathrm{norm}}^{\eta }$$

(33)

and

$$\alpha \left(q\right)\equiv {\alpha }_{\max }{(1-x_{\mathrm{norm}})}^{\xi },$$

(34)

where $q=1,\cdots ,\mathcal{N}$, $\eta$ and $\xi$ are the polynomial degrees of $d\left(x\right),\kappa \left(x\right)$, and $\alpha \left(x\right)$, respectively, and ${\mathrm{d}}_{\max }=0.75([0.8\eta +1]/\left[{\Delta }_x\sqrt{{\mu }_0/{ϵ}_0}\right])$ [53]. Thus, Equations (22)–(24) with the implemented CPML are shown below. Equation (22) is rewritten as

$$\begin{array}{c}\hfill H{x}_{(i,j-1/2)}^{n+1/2}=H{x}_{(i,j-1/2)}^{n-1/2}\left(\frac{{\mu }_{(i,j-1/2)}^{n-1/2}}{{\mu }_{(i,j-1/2)}^{n+1/2}}\right)-\left\{\left(\frac{\Delta t}{{\mu }_{(i,j-1/2)}^{n+1/2}}\right)\right[\frac{1}{{\kappa }_{(i,j-\frac{1}{2})}}\left(\frac{E{z}_{(i,j)}^n-E{z}_{(i,j-1)}^n}{\Delta y}\right)\\ \hfill +\mathrm{\Psi }{y}_{(i,j-\frac{1}{2})}^{n+\frac{1}{2}}\left(E_z\right)\left]\right\},\end{array}$$

(35)

Equation (23) takes the form

$$\begin{array}{c}\hfill H{y}_{(i-1/2,j)}^{n+1/2}=H{y}_{(i-1/2,j)}^{n-1/2}\left(\frac{{\mu }_{(i-1/2,j)}^{n-1/2}}{{\mu }_{(i-1/2,j)}^{n+1/2}}\right)+\left\{\left(\frac{\Delta t}{{\mu }_{(i-1/2,j)}^{n+1/2}}\right)\right[\frac{1}{{\kappa }_{(i-\frac{1}{2},j)}}\left(\frac{E{z}_{(i,j)}^n-E{z}_{(i-1,j)}^n}{\Delta x}\right)\\ \hfill +\mathrm{\Psi }{x}_{(i-\frac{1}{2},j)}^{n+\frac{1}{2}}\left(E_z\right)\left]\right\},\end{array}$$

(36)

whereas, as Equation (24) has two derivatives of the fields $Hx$ and $Hz$, it is necessary to introduce two convolutional terms that will take the final form

$$\begin{array}{c}\hfill E{z}_{(i,j)}^{n+1}=E{z}_{(i,j)}^n\left[\frac{4{ϵ}_{(i,j)}^n-\Delta t\left({\sigma }_{(i,j)}^{n+1}+{\sigma }_{(i,j)}^n\right)}{4{ϵ}_{(i,j)}^{n+1}+\Delta t\left({\sigma }_{(i,j)}^{n+1}+{\sigma }_{(i,j)}^n\right)}\right]+\{\frac{\Delta t}{{ϵ}_{(i,j)}^{n+1}+\frac{\Delta t}{4}\left({\sigma }_{(i,j)}^{n+1}+{\sigma }_{(i,j)}^n\right)}\\ \hfill [\frac{1}{\kappa (i,j)}\frac{H{y}_{(i+1/2,j)}^{n-1/2}-H{y}_{(i-1/2,j)}^{n-1/2}}{\Delta x}+\mathrm{\Psi }{x}_{(i,j)}^{n+1}\left(H_y\right)-\frac{1}{\kappa (i,j)}\frac{H{x}_{(i,j+1/2)}^{n-1/2}-H{x}_{(i,j-1/2)}^{n-1/2}}{\Delta y}+\\ \hfill \mathrm{\Psi }{y}_{(i,j)}^{n+1}\left(H_x\right)-\frac{J{z}_{(i,j)}^{n+1}+J{z}_{(i,i)}^n}{2}\left]\right\}.\end{array}$$

(37)

To guide the search of the optimum CPML parameters for this study in order to avoid spurious waves within the domain, we start with the values of the parameters given in [53], which proved to work for electromagnetic propagation in different frequency domains. Such parameter values are shown in

Table 1. Initial CPML parameters used in [53].

Parameter	Value
${\mu }_0$ (H/m)	$4\pi \times {10}^{-7}$
${ϵ}_0$ (F/m)	$8.85\times {10}^{-12}$
$\eta$	3.0
$\xi$	1.0
${\kappa }_{max}$	3.0
${\alpha }_{max}$	0.08
$\mathcal{N}$	20

.

To validate the numerical implementation, the electromagnetic energy flux is evaluated at each time step by calculating the magnitude of the Poynting vector at each mesh centre point as

$$S_{(i,j)}^{n+1}=\sqrt{{\left(E{z}_{(i,j)}^{n}H{y}_{(i,j)}^{n+1/2}\right)}^2+{\left(E{z}_{(i,j)}^{n}H{x}_{(i,j)}^{n+1/2}\right)}^2};$$

(38)

consequently, the electromagnetic energy flux at time $n+1$ in the entire domain is given by

$$S^{n+1}=\sum _{i,j}{S}_{(i,j)}^{n+1}.$$

(39)

2.3. Numerical Experiments

For the numerical experiments, a 2D dimensional domain of rectangular geometry of size 10 × 3 m bordered by the CPML boundary is set up, as shown in Figure 2.

To analyse the behaviour and quality of the solution, twelve uniformly distributed viewers are included in the domain, whose coordinates can be checked in

Table 2. Position of numerical viewers within the computational domain for recording electric field values $E_z$.

Viewer	01	02	03	04	05	06	07	08	09	10	11	12
$\mathit{x}$ (m)	2.00	4.00	6.00	8.00	2.00	4.00	6.00	8.00	2.00	4.00	6.00	8.00
$\mathit{y}$ (m)	0.75	0.75	0.75	0.75	1.50	1.50	1.50	1.50	2.25	2.25	2.25	2.25

, and whose locations can also be observed in Figure 2.

The sinusoidal propagation source is located in the centre of the domain, which is given by

$$J(x,y,t)=sin\left[2\pi f_0(t-t_0)\right]\delta (x-x_c)\delta (y-y_c)\mathrm{\Theta }(t-t_0),$$

(40)

where $f_0$ is the frequency, $t_0$ the initial time, and $(x_c,y_c)$ is the location at the centre of the computational domain. The source discretised in the computational domain is expressed as

$$J_{i,j}^n=\left\{\begin{array}{ccc}\sin \left[2\pi f_0(n\Delta t-t_0)\right],\hfill & \mathrm{if}(i,j)\hfill & =(i_c,j_c)\hfill \\ 0,\hfill & \mathrm{if}(i,j)\hfill & \ne (i_c,j_c)\hfill \end{array}\right.,$$

(41)

where $(i_c,j_c)$ are the corresponding nodes to the centre of the domain.

Two frequencies of the 5G band are investigated [56]: $f_0=5$ GHz (also corresponding to Wi-Fi in its 802.11n standard [37]) as well as $f_0=25$ GHz. The experiments were carried out in the presence of precipitation and free space for comparison purposes. The parameters of each simulation can be found in

Table 3. Parameters of the four numerical experiments.

Characteristic	5 GHz		25 GHz
	Free Space	Rainfall	Free Space	Rainfall
$\Delta x=\Delta y$ (m)	0.0030		$5.9958491600000004\times {10}^{-4}$
$\Delta t$ (s) ${}^{†}$	$7.07106781186547\times {10}^{-12}$		$1.4000714267493641\times {10}^{-12}$
Nodes (1)	3,339,336		83,461,716
Iterations (1)	33,480	134,100	133,920	334,800

${}^{†}$ This time step is calculated through the CFL criterion, for air (${ϵ}_r=1$ and ${\mu }_r=1$), which yield the smallest $\Delta t$ possible, ensuring quality of solution for other media.

.

Precipitation Simulation

The rainfall simulation was carried out using the parameters of conventional water. Nevertheless, different water parameters could be used, such as those shown in

Table 4. Electromagnetic properties of water in the various forms in which it is found on Earth, as well as of conventional air.

Material	${\mathit{ϵ}}_{\mathit{r}}$ (1)	$\mathit{\sigma }$ (mS/m)	${\mathit{\mu }}_{\mathit{r}}$ (1)
Air	1	0	1.0003
Distilled water		0.01
Fresh water	80–81	0.1–10	1
Salt water (and sea water)	81–88	4000
Polar snow	1.4–3	-	1
Polar ice	3.3–1.5	0.02–0.003
Warm ice	3.2	5.10${}^{-4}$–8.10${}^{-6}$
Pure ice	3.2
Freshwater frozen lake	4
Sea ice	2.5–8
Permafrost	1–8	1–0.1	1

. For this research, we used the specific values ${ϵ}_r=80$, ${\mu }_r=1$, and $\sigma =10$ mS/m corresponding to the representative rainfall of lake areas in the continental interior [57].

The geometry of the droplets is considered to be that of oblate ellipsoids [58], and a fully vertical drop trajectory is considered. For the droplet terminal drop velocities and sizes, we use limit values reported in the literature [59,60], which can be seen in

Table 5. Physical parameters of the water droplets for the experiments shown in Table 3 [59].

Drop Properties	5 GHz	25 GHz
Semi-minor axis (mm)	$2\Delta x_{\mathbf{5}\mathbf{GHz}}$–$3\Delta x_{\mathbf{5}\mathbf{GHz}}$	$5\Delta x_{\mathbf{25}\mathbf{GHz}}$–$10\Delta x_{\mathbf{25}\mathbf{GHz}}$
Semi-major axis (mm)	$4\Delta x_{\mathbf{5}\mathbf{GHz}}$–$8\Delta x_{\mathbf{5}\mathbf{GHz}}$	$8\Delta x_{\mathbf{25}\mathbf{GHz}}$–$20\Delta x_{\mathbf{25}\mathbf{GHz}}$
Eccentricity (1)	0.7071–0.7906	0.6124–0.8660
Terminal velocity (m/s)	7.61–9.90	8.72–9.40

. The terminal velocity of droplets is used as the simulation is intended to be at ground level. Drop sizes and velocities at values shown in Table 5 were generated pseudo-randomly using the Ziggurat method [61], which can be consulted in Appendix A.

To obtain realistic precipitation, the pseudo-random droplet generation process is coupled to the propagation velocities of the electromagnetic waves so that the droplet layers will only be generated every time $t=M\Delta t$, where M is a specific number of iterations for each propagation case. The numerical FDTD algorithm coupled to the Ziggurat-based rain simulation can be observed in the flowchart at Figure 3.

The code was implemented in FORTRAN 2018 standard [62], and the results were visualised manually, generating ppm graphics files. OpenACC pragmas were added to the source code to port the program to the GPU because OpenACC has demonstrated a better performance than OpenMP for a bi-dimensional FDTD numerical scheme [53,63,64]. Nevertheless, OpenMP was used to generate the pseudo-random numbers in a parallel way, through the mentioned Ziggurat algorithm (Appendix A).

The computer’s capabilities on which the experiments were carried out are:

• Intel^® processor Xeon^® CPU E5-2630 V4.
• Number of cores: 10.
• Processor frequency: 2.2 GHz.
• Total RAM memory: 128 GiB.
• Operating System: Linux Ubuntu 22.04.
• NVIDIA RTX 3060 12 GB.
• NVIDIA TITAN RTX 24 GB.

3. Results

Before proceeding with the simulations, it is necessary to mention that the boundary condition efficiency was tested, and the $\xi$ parameter had to be modified to obtain optimum absorption results with respect to the initial values proposed by [53] (see Table 1). While the other CPML parameters reported in Table 1 remained fixed, $\xi$ took the value 2.0.

3.1. Five GHz Experiments

For the 5 GHz experiments, Figure 4 shows the snapshots of the $E_z$ component propagation at different times in air (free space) without the presence of precipitation. The distribution of the energy in the computational domain is shown in Figure 5, and the behaviour of the $E_z$ electric field in the twelve numerical viewers is depicted in Figure 6.

Analogously, the snapshots of the $E_z$ component, the energy distribution, and the viewer’s behaviour for the rain scenario are shown in Figure 7, Figure 8 and Figure 9, respectively.

3.2. Twenty-Five GHz Experiments

For the experiments at 25 GHz, the snapshots of the $E_z$ component propagation at different times in the air (free space) case are shown in Figure 10; energy and viewers in the computational domain are depicted in Figure 11 and Figure 12, respectively.

Finally, and in an analogous manner, the snapshots of the $E_z$ component, the energy distribution, and the viewer’s behaviour for the rain scenario for the 5G propagation are shown in Figure 13, Figure 14, and Figure 15, respectively.

4. Discussion

4.1. Propagation Results

4.1.1. Free Space Propagation

In the snapshot $n=7000$ of the 5 GHz experiment, shown in Figure 4, it can be observed how the source was stopped and the propagation continued until the time iteration $n=8500$, where it can be noticed how the CPML efficiently absorbed the waves, and there was no electromagnetic field remaining visible in the domain. No spurious propagations can be seen within the domain coming from the region of the absorbing boundaries. The efficient absorption can be corroborated by observing the electric field viewers’ behaviour $E_z$ depicted in Figure 6. It is observed that once the source emits the last wavefront and passes through the viewers, the electric field abruptly drops to zero, and no spurious propagation is recorded within the domain, even though the simulation continued for a long time. Additionally, to validate the absorption performance, the electromagnetic energy flow in the domain increases until the 5 GHz source is turned off, as can be observed in Figure 5. After such a moment, the electromagnetic energy flux remains constant and then begins to decrease as it is absorbed by the CPML boundary, reaching the zero value, with no spurious energy appearing in the domain in the additional time.

Therefore, the effectiveness of the absorbing boundary technique implemented in this work, with the optimised parameters reported at the beginning of Section 3, is demonstrated.

Analogous results can be found for the case of a 25 GHz scenario in free space. The great similarity of the boundary condition performance with the 5 GHz experiment is remarkable since there are no remnant or spurious waves coming from the absorbing boundaries. This can be observed in the behaviour of the propagation snapshots (see Figure 10).

The $E_z$ electric field viewers for the 25 GHz case in free space show a completely analogous behaviour to the 5 GHz case (see Figure 12). In the same way, the energy flux behaviour is similar, as can be observed in Figure 11.

Thus, the successful application of the CPML absorbing boundary conditions is efficient in the frequency range of 5G herein studied. It is necessary to note that although the initial values of the parameters in Table 1 [53] provided good results with 5 GHz propagation, their performance was not good enough with 25 GHz source propagation, so the CPML parameter optimisation performed in this work was fundamental in order to have a unique and coherent set of CPML parameters for the frequency range herein explored.

4.1.2. Propagation in the Presence of Rain

For the 5 GHz experiment, the snapshots in Figure 7 show the wave scattering effect produced by the rain, which can be better appreciated in the zoom presented in Figure 16. In the snapshot at $n=6000$, it can be observed that the source is switched off, and the wavefronts are absorbed by the CPML in step $n=8000$, while the rain continues on the right side of the computational domain. It is very interesting to note the presence of residual energy at $n=$ 50,000 within the domain, which is not completely absorbed by the boundary conditions generated by the scattering caused by the rain interacting with the wavefronts. The latter can be corroborated by analysing the electromagnetic energy flux in the entire domain (see Figure 8).

The scattering effect is associated with the wavelength $\lambda =0.0599584916$ m, and the simulated droplets, in this case, are of the order of $0.012$–$0.024$ m, so they essentially act as electromagnetic field scattering particles in the short wavelength limit [44].

For the 25 GHz propagation experiment under the presence of rainfall, which is the core experiment of this work, Figure 13 shows how the wavefront reaches the region where the rain occurs, specifically, at snapshot $n=$ 50,000, where its interaction with the raindrops can be observed, proving a strong absorption process. This effect can be better appreciated in the zoom made of such a snapshot, which can be observed in Figure 17. At the same time, perfectly formed wavefronts arrive at the left boundary (see Figure 13d); a visual inspection of the right side of the domain where rain is present reveals that the vast majority of the wave has been absorbed (Figure 17). Only small, extremely thin wave beams (which have not interacted with random drops) remain, propagating towards the right boundary of the domain.

Later, in the snapshot at $n=$ 70,000, the wave source has been switched off, and the last wavefronts can be observed to be spreading within the domain (Figure 13c). Finally, for the snapshot at $n=$ 80,000, all the wavefronts have been absorbed by the CPML, leaving no remaining energy in the domain (see Figure 13d). In fact, the energy flow in the whole domain can be observed in Figure 14, where it can be appreciated that the energy cannot remain constant when propagation is turned on due to the interactions between the electric field and the water droplets, generating a continuous process of energy absorption. Subsequently, when the source is switched off, the energy is rapidly absorbed by the CPML boundaries, confirming that there is no residual electromagnetic field in the domain, as can be witnessed in the domain viewers shown in Figure 15. $E_z$ remains zero throughout the simulation, proving the CPML absorption.

Analysing the results concerning the 25 GHz propagation in the presence of rain, it can be concluded that the simulated rain generates strong absorption processes of the electromagnetic waves in question. In this particular case, the characteristic length of the droplets (see Table 5) coincides with the propagated wavelength. In this way, the droplets essentially act as obstacles. These results are consistent with those reported in the scientific literature, where it has been reported that electromagnetic waves with frequencies higher than 10 GHz suffer strong attenuation in the presence of precipitation such as rain and snow [35,36].

Finally, it should be stated that the CPML boundary condition becomes fully effective when simulating the absorption of the energy flow for both propagation scenarios, free space and with rain, considering the found set of optimised CPML parameters. To extend this set of optimum CPML parameters to propagation frequencies farther than 25 GHz would require further numerical experiments.

4.2. Computational Performance

For this type of application, code execution acceleration is crucial for achieving results in moderate execution times, particularly when the size of the domain increases exponentially. Therefore, OpenACC and OpenMP directives were used to enable the code for GPU and multi-core CPUs, using the methodology described in [63].

The implementation consists of incorporating OpenACC to accelerate the wave propagation, while using OpenMP to parallelise the rainfall (based on the Ziggurat algorithm). The NVIDIA compilers with OpenACC v23.1 were used, and the code was tested in two cards: RTX Titan and RTX 3060. The flow diagram of the parallel kernels’ execution sequence is shown in Figure 18.

In general, bi-dimensional FDTD code is of low computational intensity, which means that the number of float operations is not much more significant than the number of memory accesses. Thus, OpenMP is unsuitable for accelerating wave propagation in this case. Nevertheless, it is very useful to accelerate the Ziggurat algorithm.

The main limitation of using the GPU is the amount of available NVRAM. For this reason, an efficient implementation of the application is performed, yielding a negligible difference in the use of NVRAM between the free space and the presence of rain scenarios.

Table 6. NVRAM use for the propagation scenarios herein studied.

Scenario	5 GHz	25 GHz
Free Space	137.2 Mb	11.2 Gb
Rainfall	141.5 Mb	11.5 Gb

shows the NVRAM consumption for the four scenarios analysed in this work.

An optimised implementation of the convolutional variables used in the boundary condition is the key to reducing memory usage, following the algorithm presented in [53]. It should be stated that the memory usage in the rain simulations occupied only 2.28% more on average compared to the free space simulations, which is a consequence of the efficiency in coupling the FDTD method with the Ziggurat algorithm.

On the other hand, execution times were measured by taking as the reference one CPU core (Xeon(R) CPU E5-2630 v4 @ 2.20 GHz); the results for the 5 GHz scenario are depicted in Figure 19. As can be observed, OpenACC reduces the computing serial time by 5.06X and 5.39X for the rain scenario (for RTX Titan and RTX 3060, respectively), as well as between 5.48X and 6.93X for the free space scenario (for RTX Titan and RTX 3060, respectively), considering standard double precision variables.

For the 25 GHz experiments, a similar reduction in execution times was found, as can be observed in Figure 20. In this case, OpenACC directives reduce the execution time by 5.00X and 5.39X for the free space propagation scenario (for RTX Titan and RTX 3060, respectively), and 5.48X and 6.93X for the precipitation case (for RTX Titan and RTX 3060, respectively).

Although the rain scenario introduces computational overload, it only represents between 12.4805% and 12.5490%. Due to the large number of pseudo-random numbers used in each stage of the numerical method, this percentage is reasonable. In other words, implementing the Ziggurat method using OpenMP is highly parallel.

5. Conclusions

In this work, a formulation for Maxwell’s equations in two dimensions for the TE mode is developed, considering constitutive relations for linear, inhomogeneous, and isotropic time-dependent propagation media, in their electrical, magnetic, and conducting properties. This formulation is numerically solved by means of a conventional second-order in space and time FDTD method, coupled to Convolutional PML-type absorbing boundary conditions. Two sinusoidal sources were studied: $f_0=5$ GHz (corresponding to Wi-Fi in its 802.11n standard as well as the lowest 5G band) as well as $f_0=25$ GHz (5G), for a rectangular domain of $10\mathrm{m}\times 3\mathrm{m}$, both in free space (air) and in the presence of rain. The rain was simulated by means of a highly efficient pseudo-random numbering generator based on the parallel Ziggurat algorithm.

Despite the precipitations being idealised by considering completely vertical rainfall and by testing only one rain intensity, the main result of proving the intense absorption of 25 GHz waves by precipitation holds. Moreover, more realistic rainfall simulations can be straightforwardly performed by adding an angle for precipitation to simulate wind, as well as further values of the M parameter to test different rain intensities. Furthermore, although we present the results of electromagnetic propagation considering raindrops with sweet water properties, the developed model is suitable for studying the interaction of electromagnetic waves with other hydrometeors.

Simulations performed on $f_0=5$ GHz propagation reveal that in the presence of rain, the droplets act as scatterers of the electromagnetic field, generating quasi-stationary states (electromagnetic noise), hindering the propagation and absorption of the waves by the CPML boundary for extended periods after the rain stops. This effect occurs because the droplet size is smaller than the propagated wavelength.

In the case of $f_0=25$ GHz propagation, on the other hand, the propagated wavelength is of the same order of magnitude as the characteristic length of droplets, causing a strong attenuation phenomenon and yielding complete absorption. This result is essential for the implementation of 5G infrastructure in its more promissory operating band in the following years, according to 5G standards, as such vulnerabilities in this emerging technology must be addressed by developers and manufacturers through protection and backup measures in the presence of hydrometeors, allowing for higher service integrity and operability.