Ray Tracing and Applications to an Evaporation Duct Model Based on Data from Oceanographic Buoy Sensors

Abstract

The two-dimensional ray tracing method allows an easy and fast modeling of tropospheric propagation in the microwave frequency range. A version of this method that determines ray trajectories, amplitudes and delays of the electromagnetic field, as well as the propagation loss in a two-dimensional inhomogeneous environment will be described. The implemented algorithm may consider generalized maps of modified refractivity (or refractive modulus), including not only the vertical gradients, but also their horizontal variations along the path between transmitter and receiver. Next, the present paper will discuss the results from the application of a model of evaporation duct height to data from instruments installed in sea buoys located along the Brazilian coast. Finally, the results from the application of the ray-tracing model to evaporation ducts will be presented, to analyze the propagation of microwave signals in the maritime environment. These results will also be compared with corresponding ones from the software Advanced Refractive Effects Prediction System 3.6, based on the numerical solution of a parabolic equation.

Index Terms
amplitude; evaporation duct; microwave; parabolic equations; ray tracing; troposphere

I. INTRODUCTION

The troposphere is the lowest layer of Earth's atmosphere. Its upper limit varies between 8 km (at the poles) and 18 km (at the equator), and its temperatures varies regularly, usually decreasing as altitude increases [¹]. The troposphere refractive index n is related to the total air pressure, absolute temperature and relative humidity, varying with position. Since the refractive index n is very close to 1, its variations are better described by the refractivity $N={10}^6(n-1)$.

Refractivity gradients are very important for the analysis of propagation of radio waves in the troposphere. The vertical components of these gradients are more significant than the horizontal ones. Stratified media with single vertical profiles of the refractivity are widely used in the analysis of propagation phenomena such as sub-refraction, super-refraction and trapping by ducts, generally considering a horizontally stratified medium. That is, it is generally assumed that the vertical profile of the refractivity N(h), where h represents the height, does not change along the path between transmitter and receiver. However, the present work will show that this simplification is not always valid. In cases to be analyzed, it is necessary to consider the horizontal variations of the refractivity vertical profiles along the path between transmitter and receiver.

As an application of the ray-tracing model, this article will consider the propagation of radio waves in the presence of evaporation ducts, in particular those that occur near the sea surface. It will be seen that, over the oceans, evaporation ducts occur during high percentages of time [²], due to a decrease with height of the partial pressure of water vapor from its saturation value at the sea surface. This behavior is responsible for sufficiently negative vertical gradients of the refractivity. There are mathematical models that use meteorological measurements performed in these environments to estimate the evaporation duct heights [³]. In this work, the Paulus-Jeske model [⁴], [⁵] based on measurements of the sea surface and air temperatures, wind speed, relative humidity, and atmospheric pressure, was used with the immediately above purpose.

In section II, the main results of ray tracing theory will be presented and initially applied to a medium with a constant refractivity gradient (dN/dh) above the curved surface of the Earth. It will then be demonstrated that the results of this theory can be: (1) expressed in terms of ray tracing through a medium with a modified refractivity gradient (dM/dh) above a flattened Earth; and (2) generalized to consider horizontal variations of the modified refractivity vertical profile along the path between transmitter and receiver. In Section III, the two-dimensional modified refractivity maps available in the literature will be described. This section will also review the Paulus-Jeske model [⁴], [⁵] that estimates the heights of evaporative ducts, as well as its application to measurements made by instruments installed in nine oceanographic buoys of the Programa Nacional de Boias (PNBOIA), located along the Brazilian coast. Section IV will show the results from the application of the raytracing model to the two-dimensional modified refractivity maps and to a scenario with an irregular ground characterized by a triangular obstacle. This section will also present and discuss results from the application of the developed model to simulate the propagation of radio waves in the presence of evaporation ducts, in comparisons with corresponding results from the software Advanced Refractive Effects Prediction System (AREPS 3.6) [⁶], based on the numerical solution of a parabolic equation. The conclusions of the paper will be presented in section V.

II. Ray Tracing

A. Curved Earth

A two-dimensional (2D) environment with the troposphere (linear, isotropic, lossless and time-invariant, without variations across the great circle containing the link) immediately above a curved Earth will be considered. Assuming a time-harmonic dependence exp(-iωt), the electric field along a ray can be represented by [⁷]

(1) $$\overrightarrow{E}(\overrightarrow{r},t)=\overrightarrow{e}(\overrightarrow{r})e^{i{k}_0\zeta (\overrightarrow{r})-i\omega t}$$

where $k_0=2\pi /\lambda$ is the free space wave number, λ is the wavelength, and the geometric wave fronts are defined by surfaces with constant $\zeta (\overrightarrow{r})$. In the equation below, $\overrightarrow{r}(s)$ is the position vector of a ray in function of its length s measured from the source. Therefore, the direction of the ray is given by (2) at any of its points.

(2) $$\widehat{s}=\frac{d\overrightarrow{r}}{ds}=\frac{\nabla \zeta }{n}$$

The last equality in (2) results from the well-known eikonal equation and the orthogonality between $\widehat{s}$ and the wave front at the same point. Combining the eikonal equation with expression (2), one has

(3) $$\frac{d}{ds}(n\widehat{s})=\frac{d}{ds}\left(n\frac{d\overrightarrow{r}}{ds}\right)=\nabla n$$

Discrete versions of equations (2) and (3) can be used to determine the ray trajectories in generic 2D media. The phase along the trajectory is determined from [⁷]

(4) $$\frac{d\zeta }{ds}=\widehat{s}\cdot \nabla \zeta =\frac{\nabla \zeta }{n}\cdot \nabla \zeta =n\to \zeta (s)={\displaystyle {\int }_o^{s}n(s\text{'})ds\text{'}}$$

The right-hand side of the above equation (optical path) is integrated along the trajectory. Finally, power conservation along any narrow beam of rays, as seen in in the left plot of Fig. 1, provides the following relation [⁷]

(5) $$I_{1}d{S}_1=I_{2}d{S}_2$$

Fig. 1
(Left) narrow beam of rays, representing power conservation; (Center and Right) upper and lateral beam views.

This expression indicates that the product of the power flux I(W/m²) by the transverse elementary section dS(m²) remains constant along the beam axis. Applying this result to a beam whose top and side views are shown in the center and right plots of Fig. 2, one obtains

(6) $$n_T|E_T{|}^2[(S_0\delta {\psi }_0)\Delta z]=n_R|E_R{|}^2\left[(\delta h\mathit{cos}\psi )\Delta z\frac{s}{s_0}\right]$$

(7) $${\left|\frac{E_R}{E_T}\right|}^2=\frac{n_T}{n_R}\frac{s_0}{s}{\left(\frac{\delta h}{\delta {\psi }_o}\frac{\mathit{cos}\psi }{s_0}\right)}^{-1}$$

Fig. 2
Trajectory and angles of an infinitesimal section of a ray.

In these expressions, ψ is the angle of the beam axis with the horizontal (elevation) at the receiver, δψ_o is the elevation difference between two rays at the transmitter and is their height difference, measured along the radial that contains the origin and the receiver.

B. Flattened Earth Approximation

The application of equation (3) to a 2D horizontally stratified medium provides the well-known result [⁷]

(8) $$[n(h)\left(1+\frac{h}{a}\right)]\cos \Psi =\left[n_o\left(1+\frac{h_o}{a}\right)\right]\cos {\Psi }_o=C$$

Where h is the height of the trajectory above the Earth's surface, a is its radius and the subscript o indicates values at the transmitter. It is observed that the format of equation (8) holds for two configurations: (1) in the Flat Earth (a → ∞) limit, it can be expressed as n(h) cos ψ = C; and (2) in the Flattened Earth approximation, it is interpreted as n_m(h)cosψ = C, where the modified refractive index

n_m(h) is defined by $n_m(h)=n(h)(1+h/a)\approx n(h)+h/a$. Therefore, a flattened Earth approximation is introduced, by simultaneously replacing the refractive index n(h) with the modified refractive index n_m(h). This simplifies the ray tracing analysis by allowing the use of the Cartesian coordinate system (distance and height). By similarity with the relation $n(h)=1+{10}^{-6}N(h)$, the modified refractivity (or refractivity modulus) M(h) is defined in such a way that $n_m(h)=1+{10}^{-6}M(h)$. From the last expressions, it is immediate that

(9) $$M(h)=N(h)+156,8h_{km}=M_o+\alpha (h-h_o)$$

The right side of equation (9) was introduced to indicate that the following development will be performed for a medium with a modified refractivity that varies linearly with height (with gradient α).

Fig.2 shows an infinitesimal ray section ds. Since ds, dh, and dθ have infinitesimal values, the triangle OAC is isosceles, the side AC is horizontal and the angles OAC and OCA are approximately equal to π/2. Exploring the relations of triangle ABC, in combination with the expression n_m(h)cosψ = C, one arrives at [⁷]

(10) $$x=a\theta =\pm C{\displaystyle {\int }_{h_0}^h\frac{dh\text{'}}{\left(1+\frac{h\text{'}}{a}\right)\sqrt{n_m^2(h\text{'})-C^2}}}$$

Remembering that a » h in the region of interest, replacing the right side of expression (9) into equation (10), integrating and inverting the result, the following expression for the height h of the trajectory as a function of the distance x from the source (measured along the surface of the Earth), the gradient a and its initial values is obtained

(11) $$h=h_0+\mathit{tan}{\Psi }_{0}x+\frac{{10}^{-6}\alpha }{2{\mathit{cos}}^2{\Psi }_0}{x}^2$$

Equations (4) and (7) describe the phase and amplitude of the electric field component that propagates in the medium already characterized along the path defined by equation (11).

III. Modified Refractivity Maps and Evaporation Duct Height Model

A. Modified Refractivity Map

Two maps of modified refractivity available in the literature [⁸], [⁹] are reproduced in Fig.3. These maps clearly show horizontal variations of the vertical profile of the modified refractivity. A regular grid (x_i, h_j) was defined for each one and values of M (x_i, h_j) between 2 level curves were estimated by the quadratic mean of vertical and horizontal linear interpolations, generating a complete map of M. vertical and horizontal M gradient maps were found by applying central differences to the values of the complete maps. However, the horizontal gradient in each element was neglected, since its value was on average 525 times smaller than the vertical gradient. Therefore, only the vertical gradient of M in each element was considered.

Fig. 3
Modified refractivity maps available in literature, resulting from measurement performed in the: (a) Canterbury Region, reproduced from [⁸]; and (b) English Channel Region, reproduced from [⁹].

B. Evaporation Ducts in the Brazilian Sea

Results from measurements carried out by meteorological instruments installed in nine oceanographic buoys of the Programa Nacional de Boias (PNBOIA) located along the Brazilian coast are available at https://www.marinha.mil.br/chm/data-do-goos-brasil/pnboia-map. The positions of the buoys are shown in Fig.4. These data were accessed on February 4, 2018 and processed using the Paulus-Jeske model [⁴], [⁵], briefly described below. The former author analyzed several semi empirical models of evaporation ducts and propagation measurements performed by the latter [⁵]. The result from this analysis is a practical and well-accepted model that uses data from sensors located in oceanographic buoys or ships to estimate the evaporation duct height [⁴].

Fig. 4
Locations of the PNBOIA buoys.

The model uses the following input data, restricted to the corresponding limits: (1) wind speed (between 0 knots and 50 knots); (2) air temperature (between −20°C and 50°C, mapped onto the associated absolute temperature in the interval 253.2 K <T_ak< 323.2 K); (3) relative air humidity RH (between 0% and 100%); (4) air total pressure P_o(mbar); and (5) sea water temperature (between 0°C and 40°C, mapped onto the associated absolute temperature in the interval 273.2 K <T_sk< 313.2 K). As output, the model calculates the evaporation duct height, with limits between 0 m and 40 m, associated with the critical value of the potential refractive gradient ∇N_pot=−125 N units/km (∇N=−157 N units/km) or ∇M= 0 M units/km). An algorithm that estimates the evaporation duct height was developed from the conversion and adaptation of the original FORTRAN code described in [²], [⁴] for MATLAB language. Two modifications were made to the original code: (1) the atmospheric pressure P_o measured by the sensor installed on the buoy was used instead the standard value P_oo = 1000 mbar; and (2) the reference height z₁ = 3.7 m was used, according to the height of Brazilian buoy sensors, instead of the fixed value z₁ = 6.0 m adopted by Paulus [⁴].

The model output is the evaporation duct height δ(m) from 0 m to 40 m, measured by the difference between that at the critical value of the potential refractive gradient ∇N_p = – 125 N units/km (that is, ∇N = – 157 N units/kmor ∇M = – 0 M units/km) and the average sea level. Considering the present measurements, the model estimates the evaporation duct height δ according to the following steps [⁴], [⁵]:

1. If the wind speed u is less than 0.01 knots, δ = 0 m. Otherwise, the Richardson number Ri_b is computed by:

   $$R{i}_b=1365.3\cdot (T_{ak}-T_{sk})/(T_{ak}\cdot u^2)$$

2. Determine $\Delta N_p=N_A-N_S$, where the potential refractivity of the air N_A at 3.7 m (the height of the air temperature sensor) and the potential refractivity of the air immediately in contact with the sea surface N_S are estimated by:

   $$\begin{array}{c}{N}_A=\frac{77.6}{T_{ak}}\left\{P_o+RH\frac{293.65}{T_{ak}}\cdot e^{\left[25.22\frac{T_{ak}-273.2}{T_{ak}}-5.31\text{Ln}\left(\frac{T_{ak}}{273.2}\right)\right]}\right\}\\ N_S=\frac{77.6}{T_{sk}}\left\{P_o+\frac{29365}{T_{sk}}\cdot e^{\left[25.22\frac{T_{sk}-273.2}{T_{sk}}-5.31\text{Ln}\left(\frac{T_{sk}}{273.2}\right)\right]}\right\}\end{array}$$

3. For thermally neutral and stable conditions (0 ≤ Ri_b ≤ 1),δ = 0 m for ΔNp ≥ 0. If the latter inequality does not hold, the evaporation duct height is estimated by:

   $$\delta =-\frac{\Delta N_p}{1263.75+\frac{2405+5.2\Delta N_p}{L\text{'}}}$$

   Expressions for the Monin-Obukhov length L’ will be provided below. If, from the above expression, δ < 0 or δ/L’ > 1, the result should be replaced by:

   $$\delta =-\frac{2405+6.2\Delta N_p}{1264.15}$$

4. On the other hand, for thermally unstable conditions (Ri_b < 0), the evaporation duct height should be estimated by:

   $$\begin{array}{ccc}\delta ={(D^4-\frac{18}{L\text{'}}{D}^3)}^{-1/4}& \text{where}& D=-125\frac{10.11-\Psi (3.7/L\text{'})}{\Delta N_p}\end{array}$$

   and the universal function Ψ(p) is equal to: (1) Ψ(p) = — 4.500 p, for p ≥ — 0.010; (2) Ψ(p) = 4.898 (—p)^1.020, for — 0.026 ≤ p ≤ —0.010; (3) Ψ(p) = 2.023 (—p)^0.776, for —0.100 ≤ p < —0.026; (4) Ψ(p) = 1.445 (—p)^0.630, for — 1.000 ≤ p < — 0.100; (5) Ψ(p) = 1.445 (—p)^0.414, for —2.200 ≤ p < — 1.000; and (6) Ψ(p) = 2.0 0 0, for p < — 2.200.

5. The Monin-Obukhov length is calculated by L’ = 3 7 Γ_e/Ri_b, where Ri_b has been determined in step 1 and: (1) Γ_e =0.050, for Ri_b ≤ −3.75; (2) ${\Gamma }_e=0.065+0.004{\mathit{Ri}}_b$, for −3.75 < Ri_b ≤ −0.1 2; (3) ${\Gamma }_e=0.109+0.367{\mathit{Ri}}_b$, for −0.12 < Ri_b ≤ 0.14; and (4) ${\Gamma }_e=0.155+0.021{\mathit{Ri}}_b$, for 0.14 < Ri_b.

Finally, the vertical profile of the modified refractivity M(h) due to an evaporation duct is calculated by

$$M(h)=M_0+0.125h-0.125\delta \text{Ln}(h/h_0)$$

where h₀ = 1.5 · 10^-4 m (initial height) and M₀ = 340 M units (initial refractivity). This expression immediately shows that dM/dh = 0 at h = δ, as defined above.

The algorithm also divides the measurements into 4 periods of the day: dawn (0 h to 6 h), morning (6 h to 12 h), afternoon (12 h to 18 h) and night (18 h to 24 h). Simple arithmetic means of each physical quantity (temperatures, humidity, total pressure and wind speed) in each period are calculated and then used to simulate the ducts heights. Thus, four values of each parameter (averages of the dawn, morning, afternoon, and night values) and, consequently, four evaporation ducts heights are associated with each measurement day. Table I shows the probability of occurrence of ducts with a height higher than 20 m in the nine geographic coverage areas of the PNBOIA, namely: Rio Grande - RS, Itajaí - SC, Santos - SP, Niterói - RJ, Cabo Frio - RJ, Vitória - ES, Porto Seguro - BA, Recife - PE and Fortaleza - CE. These results were obtained by applying the algorithm described above to 101,004 measurements of each physical parameter (wind speed, air temperature, sea temperature, relative humidity and atmospheric pressure), available between 2009 and 2017.

Table I
Percentage probability of monthly occurrence of evaporation ducts height higher than 20 m along the Brazilian coast

Month\Buoy	Rio Grande	Itajaí	Santos	Niterói	Cabo Frio	Vitória	Porto Seguro	Recife	Fortaleza
JANUARY	42	30	1	11	14	15	60	60	85
FEBRUARY	29	45	13	9	14	36	68	80	63
MARCH	60	63	42	0	27	28	61	74	51
APRIL	69	44	34	*	42	50	52	78	61
MAY	56	38	48	3	25	34	44	91	85
JUNE	44	27	23	1	0	18	31	76	50
JULY	46	17	21	2	17	39	28	65	93
AUGUST	37	19	28	2	29	*	60	90	98
SEPTEMBER	22	22	35	3	20	*	38	84	92
OCTOBER	22	20	21	2	17	7	42	97	80
NOVEMBER	27	34	21	4	18	12	9	96	94
DECEMBER	31	33	9	18	16	30	50	71	85
MONTHLY MEAN	40	33	25	5	20	27	45	80	78

• ^*
  Data not avalaible

It is observed that the northeast region has a very high probability of occurrence of evaporation ducts thicker than 20 m, with low variations between the months of the year. The average probabilities of occurrence of evaporation ducts thicker than 20 m in the two southernmost buoys (Rio Grande and Itajaí) are also relatively high. It can also be observed that Niterói and Cabo Frio results are substantially different, despite their geographic proximity, and that the Niterói buoy has the lowest probability values among the nine buoys. This occurs because the Niterói buoy is the only one not located in the open sea. It is located inside the Guanabara Bay, being less influenced by winds.

C. Analysis of Evaporation Duct Using the Itajaí Buoy Data

The Itajaí buoy is located off the coast of the state of Santa Catarina, at the following coordinates: latitude 27.4° S and longitude 47.26° W. The total of 30,947 measurements of each physical quantity available by PNBOIA from year 2009 to 2017 were distributed as follows: 4,922 measurements in 2009; 7,618 in 2011; 6,647 in 2012; 1,846 in 2013; 671 in 2014; 4,346 in 2015; 2,233 in 2016; and 2,664 in 2017. The distributions of evaporation ducts heights for most of the above years are plotted in Fig. 5.

Fig. 5
Occurrences of evaporation ducts heights estimated from the measurements performed by the Itajaí buoy during the indicated years.

It is clearly observed that the distributions obey a relatively regular pattern in all years, with higher occurrences between 10 m and 25 m. This behavior of the distribution of the evaporation ducts height has been observed for all years and geographic positions.

Fig. 6 shows the average values of the evaporation duct heights for each month, considering all the 30,947 Itajaí buoy measurements. The average evaporation duct heights are small during the winter and highest in March.

Fig. 6
Variation of the monthly averages of the evaporation duct heights estimated from the Itajaí buoy measurements.

D. Evaporation Duct Behaviors at Two Different Buoys

This section compares the behaviors of the meteorological parameters and the corresponding evaporation duct heights estimated in the same winter month (July 2017) by the Itajaí and Fortaleza buoys, located in the south and the northeast of Brazil, respectively. Figs. 7 and 8 show the average dawn, morning, afternoon, and night values of the air temperature, humidity, total pressure, sea temperature, and wind speed, as well as the associated evaporation duct heights, corresponding to the two buoys, respectively.

Fig. 7
Meteorological data and evaporation duct height from the Itajaí buoy in July 2017.

Fig. 8
Duct and meteorological data of Fortaleza buoy in July 2017.

Figs. 7 and 8 show that the air and sea temperatures remain relatively constant and equal for most of the days at both locations. However, the Fortaleza temperatures (28 ^oC) are higher than those at Itajaí (22 ^oC). On the other hand, it is observed that, on average, the Itajaí humidities and wind speeds are smaller and more variable than those at Fortaleza. In some cases, the variability of the parameters occur simultaneously.

These variabilities are mapped into the ones in the evaporation duct height. At both locations, it is observed that the evaporation duct height tend to increase when the humidity decreases.

Finally, comparing Fig. 7 and 8, one concludes that the difference between the average evaporation duct heights is caused by the significant difference between the air and sea temperatures of the two locations during the winter (July 2017) of the southern hemisphere.

IV. Numerical Simulations

A. Regions with Horizontal Variations of ∇M

Initially, the element of each map described in section III.A that contains the transmitter is identified. The background of Fig. 9a displays the rectangular mesh that defines one such map. Then, the ray tracing procedure is performed on this element, based on equation (11), considering the appropriate values of h_o, Ψ_o and α. The edge of the element intercepted by the ray is identified, as well as the new element in which the procedure should be performed next (again, with appropriate parameters). This partial procedure is repeated until the ray reaches the desired distance. The ray tracing procedure is also repeated for all desired initial angles Ψ_oq (q = 1, …,Q), and the corresponding heights h_q (q = 1, …, Q) at the desired distance are stored.

Next, all elevation intervals (Ψ_oq, Ψ_oq+1) associated with a height interval (h_q, h_q+1) that contains the receiver height h_Rx at the desired distance are identified. Using the bisection method in combination with the procedure of the previous paragraph, the rays in these elevation intervals that reach the receiver within a prescribed precision are determined. Their delays ${\tau }_q=\zeta (s)/c$ (where c is the velocity of light in vacuum) are then estimated, with the aid of equation (4). Let Ψ_oRx represent the initial elevations of the rays that reach the receiver. Rays with initial elevations ${\Psi }_{\mathit{oRx}}\pm \delta {\Psi }_o/2$ are then traced and their height differences at the desired distance used to estimate the central ray amplitude, through equation (7).

In addition to the ray tracing trajectory, the model determines the value of the Transfer Function H(f) of the channel defined by the transmitter (fixed) and each observation point at the desired frequency, considering the effects of multipath. That is, it calculates

(12) $$H(f)={\displaystyle {\sum }_{p=0}^{p=P}{\left|\frac{E_R}{E_T}\right|}_p{e}^{-i2\pi f{\tau }_p}}$$

considering the amplitude and phase of P rays that arrive at the observation point, determined as described above.

The propagation factor PF is defined by PF=H (f) and the path loss PL is represented by equation (13), where r(km) is the distance between the transmitter and receiver and f(MHz) is the operation frequency. In this equation, the term within parentheses represents the Free Space Loss (FSL).

(13) $$PL=(\mathrm{32,44}+20{\mathit{log}}_{10}r+20{\mathit{log}}_{10}f)–20{\mathit{log}}_{10}PF$$

Simulations of PL values were performed for the Canterbury [⁸] and English Channel [⁹] modified refractivity maps. The following configurations were used: frequencies from 1 GHz to 20 GHz, antenna heights of 15 m, 1000 rays traced with equally-spaced initial elevations in the interval (-1°, +1°), Gaussian transmitter antenna with a 13° half-power beam width and horizontal polarization. The ground was considered flat, with relative permittivity ε_r = 75 and conductivity σ = 5 S/m (sea water) for the range between 0 km and 235 km, and ε_r = 30 and σ = 0.01 S/m (wet soil) for the range between −10 km and 0 km. For the Canterbury simulation [⁸], the transmitting antenna was positioned at – 10 km (10 km on land) and the receiving antenna at 50 km into the sea. Fourteen rays arrived at the receiving antenna. For the English Channel simulation [⁹], the transmitting antenna was positioned at 0 km, the receiving antenna at 90 km and two rays were received.

Figs. 9a and 9b show the effects from ducts on ray tracing in the two regions. These effects allow Beyond Line-of-Sight (B-LoS) communications (that is, beyond 32 km for the corresponding antenna heights, estimated with basis on the well-known approximation $d_{hor}=\sqrt{2a_{eff}}\left(\sqrt{h_t}+\sqrt{h_r}\right)$, where a_eff is the effective Earth's radius for a reference atmosphere). As shown in Fig.9b, the ray tracing model also identified: (1) shadow regions, such as the one limited by the range interval (120 km, 180 km) and altitudes less than 100 m; and (2) multipath regions, such as the one limited by the range and altitude intervals (140 km, 160 km) and (250 m, 300 m), respectively. Fig.10 shows that the path loss PL in the Canterbury [⁸] and English Channel [⁹] links are generally less than the free-space loss FSL. Only for frequencies in the band (1.0 GHz, 2.0 GHz) does PL in the Canterbury [⁸] link slightly exceed the corresponding value of FSL.

Fig. 9
Ray tracing with colored background ∇M for the regions of: (a) Canterbury [⁸]; and (b) English Channel [⁹],

Fig. 10
Comparison between PL and FSL as functions of frequency: (a) for the Canterbury region [⁸], the transmitting antenna was positioned at −10 km (on land) and the receiving antenna at 50 km (into the sea); (b) for the English Channel region [⁹], the transmitting antenna was positioned at 0 km and the receiving antenna at 90 km. All antennas are located 15 m above the ground or sea.

B. Transmitting and Receiving Antennas Inside an Evaporation Duct

In this section, the ray tracing model is applied to a simple 20-m high evaporation duct environment, with antennas 15 m above the sea surface (ε_r = 7 5 and σ = 5 S/m), using horizontal polarization and maximum range of 60 km. The left plot of Fig. 11 displays the assumed vertical profile of the modified refractivity. A grid of 1000 rays was traced from the transmitting antenna with equally-spaced initial elevations in the interval (-5°, +5°). These initial rays, displayed in the right plot of Fig.11, were used to determine the two rays (direct and reflected at the sea surface) that reach each observation point, as explained in Section IV.A. Fig. 12 displays the ray-tracing predictions of the path loss PL dependence on range, for the frequencies 3 GHz, 5 GHz, and 10 GHz. The results are compared with the corresponding ones provided by the free-space model FSL and the AREPS 3.6 software [⁶], which is based on a numerical solution of a parabolic equation.

Fig. 11
Ray tracing in the presence of a 20-m high evaporation duct, for15-m high transmitting and receiving antennas.

Fig. 12
Comparison between the path loss PL values estimated by the ray tracing and AREPS 3.6 models for 15-m high transmitting and receiving antennas within the 20-m high evaporation duct and the frequencies: (a) 3GHz; (b) 5 GHz; and (c) 10 GHz.

Note that the corresponding ray tracing and AREPS 3.6 results become increasingly closer together as the frequency increases. Differences similar to those observed in Fig. 12 can also be seen in Dahman et al. [¹⁰], which compares ray tracing and Parabolic Equation Toolbox (PETOOLS) [¹¹] propagation simulations for 1.5 GHz and 4.5 GHz and different ducts heights. Dinc and Akan [¹²] also used the AREPS software to observe that the confinement of transmitted power by evaporation ducts is more effective at higher frequencies (10.5 GHz, 15.0 GHz, and 20.0 GHz). Indeed, Fig. 12(c) shows that the path loss PL is smaller than the free-space loss FSL for all ranges greater than 20 km.

This interval includes transmitter-receiver paths that would ordinarily be obstructed by the Earths curvature under reference atmospheric conditions. Thus, diffraction effects would increase the corresponding path losses beyond the free-space values [¹³].

C. Transmitting and Receiving Antennas Above an Evaporation Duct

Next, the ray tracing model was applied to an evaporation duct environment characterized by an M profile with null gradient at the height of 20 m, as displayed in the left plot of Fig.13. Both transmitting and receiving antennas (with 13° half-power beam width and horizontal polarization) were located 35 m above the sea surface (ε_r = 75 and σ = 5S/m).

Fig. 13
Ray tracing in the presence of a 20-m high evaporation duct, for 35-m high transmitting and receiving antennas.

Again, 1000 rays were traced from the transmitting antenna with equally-spaced initial elevations in the interval (-5°, +5°). These initial rays, displayed in the right plot of Fig.13, were used to determine the two rays (direct and reflected at the sea surface) that reach each observation point. This plot shows rays arriving at approximately 82 km, well beyond the line-of-sight limit (49 km, for a link with 35-m high antennas). Fig. 14 displays the ray-tracing predictions of the path loss PL dependence on range, for the frequencies 1 GHz and 3 GHz. The results are compared with the corresponding ones provided by the free-space model FSL, the one combining FSL with diffraction effects due to the Earth's curvature [¹³], and the AREPS 3.6 software [⁶]. At 1 GHz, the results from the ray-tracing model and the AREPS 3.6 software are similar for ranges in the interval (1 km, 35 km). At 3 GHz, the results are similar for ranges in the interval (1 km, 73 km). Additionally, Fig.14(a) shows a difference between the predictions of the two models for the 1 GHz frequency: while the ray-tracing model predicts nearfree space propagation losses for distances less than 60 km, the AREPS 3.6 software foresees a transition from this regime to that of intermediate propagation losses between the FSL and FSL+Diffr curves. On the other hand, for the 3 GHz frequency, the two models agree with a regime of propagation losses always close to those of free space for distances less than 70 km, as shown by Fig. 14(b).

Fig. 14
Comparison between the path loss PL values estimated by the ray tracing and AREPS 3.6 models for 35-m high transmitting and receiving antennas above the 20-m high evaporation duct and the frequencies: (a) 1 GHz; and (b) 3 GHz.

D. Irregular Terrain

The implemented ray-tracing algorithm also supports atmospheric environments above piecewise-linear irregular terrain. Figs. 15(a) and 15(b) show the results from the application of the ray-tracing model to an environment characterized by the standard atmosphere (∇M = 118 M units/km) above a 20 m high triangular terrain. Both transmitting and receiving antennas (with half-power beam width 13° and horizontal polarization) were placed 20 m above the horizontal baseline. The following electrical parameters were assumed for the terrain: ε_r = 75 and σ = 1mS/m. Initially, 500 rays were traced from the transmitting antenna with equally-spaced initial elevations in the interval (-4°, +4°), as shown in Fig. 15(a). Then, to minimize difficulties in the calculation of the path loss PL for observation points near the top of triangular terrain (for ranges near 10 km), 30 additional rays were similarly traced in the interval (-0.05°, 0°), as shown in Fig. 15(b). Combining these grids of initial rays, the algorithm was able to determine two received rays (direct and ground-reflected) at all observation points with ranges less than 10 km, which is the limit of geometric optics for the selected antenna heights.

Fig. 16 displays the ray-tracing predictions of the path loss PL dependence on range, for 1 GHz. A very good agreement is observed between the ray-tracing results and those provided by the AREPS 3.6 software [⁵] in the above range interval.

Fig. 15
Ray tracing in an environment characterized by the standard atmosphere (∇M = 118 M units/km) above a 20-m high triangular terrain, for transmitting and receiving antennas located 20 m above the horizontal baseline: (a) 500 rays with equally-spaced initial elevations in the interval (-4°, +4°); and (b) 30 additional rays were similarly traced in the interval (-0.05°, 0°).

Fig. 16
Path loss dependence on the distance between transmitting and receiving antennas located 20 m above the horizontal baseline and operating at 1 GHz in an environment characterized by the standard atmosphere (∇M = 118 M units/km) above a 20-m high triangular terrain.

V. Conclusion

The results from the ray-tracing model obtained for the regions of Canterbury [⁸] and English Channel [⁹] have shown the importance of considering horizontal variations of the vertical profile of the modified refractivity for propagation predictions.

Evaporation duct heights have been estimated using data recorded by oceanographic buoys of the PNBOIA located along the Brazilian coast, applied to the Paulus-Jeske model [⁴], [⁵]. The results have shown a high percentage of occurrence of evaporation ducts with heights greater than 20 m, particular in the southern and northeastern regions. These results confirm the importance of considering evaporation ducts in communications in marine environments of the Brazilian coast.

For the above reason, the propagation electromagnetic fields have been simulated by means of the ray-tracing model in the presence of a hypothetical evaporation duct of 20 m height. Two scenarios were considered, with antennas located above (35 m) or inside (15 m) of the evaporation duct. Differences have been observed between the ray tracing and AREPS 3.6 results at 1 GHz. However, both models have always agreed that, in the presence of evaporation ducts, the field at great distances is more intense than the one predicted by the association of the free space loss with the additional attenuation due to the diffraction by the curvature of the Earth. Therefore, it may be interesting to further investigate the use of propagation in the presence of evaporation ducts for maritime communication applications.

Finally, it is important to note that the ray-tracing model can provide information on broadband channels in a natural and relatively easy way. Among the parameters that could be predicted by the model, one would list arrival and departure angles, channel transfer functions and impulse responses, as well as average delays and rms delay spreads.

REFERENCES

• ^[1]
  International Telecommunication Union - Radiocommunication Sector, Handbook on Radiometeorology, Geneva, CH, 2013.
• ^[2]
  R. A. Paulus, “Practical aplication of an evaporation duct model,” Radio Science, vol. 20, pp. 887-896, 1985.
• ^[3]
  S. M. Babin, G. S. Young, and J. A. Carton, “A new model of the oceanic evaporation duct,” Journal of Applied Meteorology, vol. 36, pp. 193-204, 1997.
• ^[4]
  R. A. Paulus, Specification for evaporation duct height calculations, Technical Document 1596, Naval Ocean Systems Center, San Diego, CA, USA, 1989.
• ^[5]
  H. Jeske, State and limits of prediction methods of radar wave propagation conditions over sea. In Zancla A. (eds) Modern Topics in Microwave Propagation and Air-Sea Interaction. NATO Advanced Study Institutes Series (Series C - Mathematical and Physical Sciences), vol.5, pp. 130-148, Springer, Dordrecht, NL, 1973.
• ^[6]
  W. L. Patterson, Advanced refractive effects prediction system (AREPS), 2007 IEEE Radar Conference, pp. 891-895, Boston MA, april 2007.
• ^[7]
  M. Born and E. Wolf, Principles of Optics: Electromagnetic Theory of Propagation Interference and Diffraction of Light, 7^th ed., Cambridge University Press, Cambridge, UK, 2003.
• ^[8]
  B. R. Bean and E. J. Dutton, Radio Meteorology, National Bureau of Standards Monograph 92, Boulder, CO,USA, 1966.
• ^[9]
  L. W. Barclay, Propagation of Radiowaves, 3^rd ed., The Institution of Engineering and Technology, London, UK, 2013.
• ^[10]
  G. Dahman, F. Gagnon, and G. Poitau, Ship-to-ship beyond line-of-sight communications: a comparison between ray tracing simulations and PETOOL, XXXII URSI GASS, Montreal, Canada, August, 2017.
• ^[11]
  O. Ozgun, G. Apaydin, M. Kuzuoglu, L.Sevgi,“PETOOL: MATLAB-based one-way and two-way split-step parabolic equation tool for radiowave propagation over variable terrain,” Computer Physics Communications, vol. 182, pp. 2638-2654, 2011.
• ^[12]
  E. Dinc and O. B. Akan, “Beyond-line-of-sight communications with ducting layer,” IEEE Communications Magazine, vol. 52, pp. 37-43, 2014.
• ^[13]
  International Telecommunication Union - Radiocommunication Sector, Propagation by diffraction, Recommendation ITU-R P.526-14, Geneva, CH, 2018.

Publication Dates

• Publication in this collection
  Mar 2019

History

• Received
  30 Sept 2018
• Reviewed
  04 Oct 2018
• Accepted
  22 Dec 2018