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

Brazilian Microwave and Optoelectronics Society-SBMO received 30 Sept 2018; for review 4 Oct 2018; accepted 22 Dec 2018 Brazilian Society of Electromagnetism-SBMag © 2019 SBMO/SBMag ISSN 2179-1074 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.


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 [1].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 -.
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

Ray Tracing and Applications to an Evaporation Duct Model Based on Data from
Oceanographic Buoy Sensors Leonardo L. Freitas and Emanoel Costa Centro de Estudos de Telecomunicações, Pontifícia Universidade Católica do Rio de Janeiro (CETUC PUC-Rio) Rua Marquês de São Vicente 225, 22451-900 Rio de Janeiro RJ, Brazil leonardofreitas@cetuc.puc-rio.br,epoc@cetuc.puc-rio.brvalid.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 [2], 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 [3].In this work, the Paulus-Jeske model [4], [5] 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 [4], [5] 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) [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 timeinvariant, 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 [7] - (1) where k 0 = 2 is the free space wave number,  is the wavelength, and the geometric wave fronts are defined by surfaces with constant .In the equation below, 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. ( The last equality in (2) results from the well-known eikonal equation and the orthogonality between and the wave front at the same point.Combining the eikonal equation with expression (2), one has 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 [7] (4) 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 [7] ( 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 In these expressions, is the angle of the beam axis with the horizontal (elevation) at the receiver, 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 [7] (8) 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) 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 , one arrives at [7] - (10) 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  and its initial values is obtained 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 [8], [9] 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.

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 [4], [5], briefly described below.The former author analyzed several semi empirical models of evaporation ducts and propagation measurements performed by the latter [5].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 [4].).An algorithm that estimates the evaporation duct height was developed from the conversion and adaptation of the original FORTRAN code described in [2], [4] 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 1 = 3.7 m was used, according to the height of Brazilian buoy sensors, instead of the fixed value z 1 = 6.0 m adopted by Paulus [4].Finally, the vertical profile of the modified refractivity M(h) due to an evaporation duct is calculated by where m (initial height) and M 0 = 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.

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: 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.

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.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 o C) are higher than those at Itajaí (22 o C).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.

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 (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 oRx ± 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 p p p p (12) 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 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). - Simulations of PL values were performed for the Canterbury [8] and English Channel [9] 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 [8], 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 [9], 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  , where 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 [8] and English Channel [9] 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 [8] link slightly exceed the corresponding value of FSL.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. [10], which compares ray tracing and Parabolic Equation Toolbox (PETOOLS) [11] propagation simulations for 1.5 GHz and 4.5 GHz and different ducts heights.Dinc and Akan [12] 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.

D. Irregular Terrain
The implemented ray-tracing algorithm also supports atmospheric environments above piecewiselinear irregular terrain.

V. CONCLUSION
The results from the ray-tracing model obtained for the regions of Canterbury [8] and English Channel [9] 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 [4], [5].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.

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

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

5 -
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 (that is, or ) and the average sea level.Considering the present measurements, the model estimates the evaporation duct height according to the following steps[4],[5]:1 -Ifthe wind speed u is less than 0.01 knots, = 0 m.Otherwise, the Richardson number Ri b is computed by: 2 -Determine Δ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: 3 -For her ally e ral a d able co d o (0 ≤ Ri b ≤ ), = 0 m for Δ p ≥ 0. If he la er inequality does not hold, the evaporation duct height is estimated by: Expressions for the Monin-Obukhov length will be provided below.If, from the above expression, ˂ 0 or 1, the result should be replaced by: 4 -On the other hand, for thermally unstable conditions (Ri b ˂ 0), the evaporation duct height should be estimated by: and the universal function is equal to: The Monin-Obukhov length is calculated by where has been determined in step 1 and: (1) Γ e =0.050, for Ri b ≤ −3. ; ( ) Γ e =0.065+0.004Ri b , for −3.˂ Ri b ≤ −0. ; (3) Γ e =0.109+0.367Ri b , for −0.˂ Ri b ≤ 0.14; and (4) Γ e =0.155+0.021Rib , for 0. 4 ˂ Ri b .

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

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

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

Fig. 10 .
Fig.10.Comparison between PL and FSL as functions of frequency: (a) for the Canterbury region[8], 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[9], 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.

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

Fig. 12 .
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.

Fig. 14 .
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.
Fig.16displays 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[5] in the above range interval.

Fig. 15 .
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 .
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.

TABLE I .
PERCENTAGE PROBABILITY OF MONTHLY OCCURRENCE OF EVAPORATION DUCTS HEIGHT HIGHER THAN 20 M ALONG THE BRAZILIAN COASTIt 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.