Bidimensional Ray Tracing Model for the Underwater Noise Propagation Prediction

: An increasing attention has recently been paid to the effect of the underwater noise ﬁeld generated by ship activities on the marine environment. Although this problem is widely discussed in international treaties and conventions, it has not yet found a consolidated technical-scientiﬁc treatment capable of quantifying the level of underwater noise emissions produced by naval systems. As part of a national research collaboration, a novel code has been developed to predict noise propagation according to the Ray Tracing approach. Such optical geometry-based technique allows for calculating the Transmission Loss (TL) trend in its respective contributions: geometrical loss (due to the distance between the source and receiver), dissipation loss (due to the characteristics of the propagation environment), and reﬂection loss (due to the surfaces that delimit the ﬁeld). The simulation requires as input parameters the source info as spatial position, frequency, and sound pressure level (SPL) as well as the sea properties like seabed depth, the speed of sound proﬁle, the layers thickness the water column is divided into, the sea salinity, temperature, and pH. The simulation code provides the SPL spatial distribution useful as a fast industrial tool in the future studies addressed to identify the emission limits for the protection of marine wildlife.


Introduction
A topic of scientific interest for many years now is the study of noise propagated in the sea. As pioneering contributions in this sense, it is worth recalling the works of Knudsen [1] and Wenz [2], published at the beginning of the twentieth century. The sound disturbances emitted by naval systems are certainly among the main components of noise pollution-both sonic and ultrasonic-negatively affecting the marine ecosystem [3][4][5][6][7]. Awareness perspectives are supported by binding regulations of certain states or of EU Directives on the care of the natural environment, such as to impose limits on underwater noise levels [8]. The United Nations Convention on the Law of the Sea (UNCLOS) established underwater noise produced by both coastal and offshore anthropogenic activities as a source of pollution of the marine environment [9]. Article 1 defines the "pollution of the marine environment". It is classified on the basis of continuous emissions (such as drilling, ship traffic) and intermittent (such as pile drivers, geoseismic surveys), which are added to the noise normally present in the marine environment. The characteristics of ship noise are a function of various parameters such as the type of vessel, its size, type of propulsion, and speed performance. The noise produced by the propeller in particular could have an emission range of up to hundreds of kilometers. In 2008, underwater noise was defined as a qualitative parameter of the state of health of the marine environment, thus requiring Member States to consider the problem and take precautionary measures aimed at containing cross-border pollution emissions (European Community, through the "Marine Strategy Framework Directive" [10]). In order to assess the extent of the environmental impact, national and international organizations were the first to issue procedures

Theoretical Model of Ray Tracing
The development of predictive tools can also prove extremely useful in order to optimize the noise levels emitted at sea: the mathematical models are based on the phenomenological laws that govern the physics of wave propagation in fluids. The "Ray Tracing theory" represents a possible modeling approach for this purpose. This type of modeling is used in the optical field to analyze the light diffusion, but its applicability can also be extended to the acoustic field. This technique is a high-frequency approximation of the solution of the wave equation. It allows to define the propagation of sound through the emission of rays emitted by a point source in a medium at a constant speed of sound. The rays emitted by the source are perpendicular to the wave fronts, which are instead spherical surfaces radiated by the source itself. The reference theory is due to the Snell-Descartes law, which is used in acoustic geometry to describe the interference phenomena that occur Fluids 2021, 6,19 3 of 14 between two homogeneous fluids having different speeds. Let us consider two fluids of which the first has density ρ 1 and speed c 1 and the second which has density ρ 2 and speed c 2 (with c 2 > c 1 ), as represented in Figure 1.
Fluids 2021, 6, x FOR PEER REVIEW 3 of 14 emission of rays emitted by a point source in a medium at a constant speed of sound. The rays emitted by the source are perpendicular to the wave fronts, which are instead spherical surfaces radiated by the source itself. The reference theory is due to the Snell-Descartes law, which is used in acoustic geometry to describe the interference phenomena that occur between two homogeneous fluids having different speeds. Let us consider two fluids of which the first has density ρ1 and speed c1 and the second which has density ρ2 and speed c2 (with c2 > c1), as represented in Figure 1. Snell-Descartes Theorem. The ratio of the cosines of the angles of incidence and refraction is equivalent to the ratio of phase velocities in the two media, or equivalent to the reciprocal of the ratio of the indices of refraction.
The difference in the velocity of the fluids means that the incident ray (inclined by the angle β1 with respect to the horizontal line) partially reflects upwards symmetrically and partially refracts in the second fluid with an angle β2 that is obtained from the following relationship (1): This equation is valid only if β2 ≤ 1 and therefore cosβ1 ≤ c1/c2. From this relation, it is possible to introduce another angle βcritic called "critical angle", equal to Formula (2): = arccos For incidence angles β1 smaller than βcritic, there is total reflection. The same law can also be extended to study sound propagation phenomena when velocity varies linearly with depth.
Considering the law of variation of speed with depth as per Equation (3) (see Figure 2): In which z is the vertical dimension and g is the speed gradient. Substituting the Eq. (3) in the Snell-Descartes law Equation (1), the Equation (4) is obtained: Given the radius RC = c0/gcosβ0, the relationship (4) is reduced as follows: It is observed, therefore, that in the case in which the speed varies with depth, the rays do not maintain straight trajectories but bend defining curved sections, unlike the case in which the speed remains constant and the rays propagate according to rectilinear sections. Snell-Descartes Theorem. The ratio of the cosines of the angles of incidence and refraction is equivalent to the ratio of phase velocities in the two media, or equivalent to the reciprocal of the ratio of the indices of refraction.
The difference in the velocity of the fluids means that the incident ray (inclined by the angle β 1 with respect to the horizontal line) partially reflects upwards symmetrically and partially refracts in the second fluid with an angle β 2 that is obtained from the following relationship (1): This equation is valid only if β 2 ≤ 1 and therefore cosβ 1 ≤ c 1 /c 2 . From this relation, it is possible to introduce another angle β critic called "critical angle", equal to Formula (2): For incidence angles β 1 smaller than β critic , there is total reflection. The same law can also be extended to study sound propagation phenomena when velocity varies linearly with depth.
Considering the law of variation of speed with depth as per Equation (3) (see Figure 2): In which z is the vertical dimension and g is the speed gradient. Substituting the Equation (3) in the Snell-Descartes law Equation (1), the Equation (4) is obtained: Given the radius R C = c 0 /gcosβ 0 , the relationship (4) is reduced as follows: It is observed, therefore, that in the case in which the speed varies with depth, the rays do not maintain straight trajectories but bend defining curved sections, unlike the case in which the speed remains constant and the rays propagate according to rectilinear sections. The method therefore consists in determining of the propagation of sound through linear segments and arcs that are connected along the water column, Figure 2. Relying upon this physical principle, the numerical prediction has been developed in this paper. The method therefore consists in determining of the propagation of sound through linear segments and arcs that are connected along the water column, Figure 2. Relying upon this physical principle, the numerical prediction has been developed in this paper.

Numerical Prediction Code: Objectives and Functionality of the Program
The simulation program developed and described within this paper implements the basic foundations of this theory. The main purpose was to evaluate the Transmission Loss (TL) in water, and therefore the calculation of the sound pressure level (SPL) of a given point-like source that radiates a sound wave. The phases necessary to achieve this goal will be explained below, giving the reader the opportunity to get the program steps and then describing the applicative example about the noise analysis of the Adriatic Sea (Italy). From a general standpoint, the simulation code is structured in two distinct parts in order to guarantee greater speed and fluidity. The first part is dedicated to the input of the characteristic parameters of the source (the minimum angle, the maximum angle, and the step of emission of the rays) and of the surrounding environment (depth, temperature, salinity, pH). These data are fundamental to evaluate the water absorption coefficient α (Appendix A). Subsequently, the program, through a series of iterative cycles, will provide the different propagation paths of the rays emitted by the source. In the second part of the program, the user has the possibility to set up a matrix of virtual receivers through which to map the trend of the Transmission Loss (TL) and the SPL sound pressure level. The global layout is represented in Figure 3.

Numerical Prediction Code: Objectives and Functionality of the Program
The simulation program developed and described within this paper implements the basic foundations of this theory. The main purpose was to evaluate the Transmission Loss (TL) in water, and therefore the calculation of the sound pressure level (SPL) of a given point-like source that radiates a sound wave. The phases necessary to achieve this goal will be explained below, giving the reader the opportunity to get the program steps and then describing the applicative example about the noise analysis of the Adriatic Sea (Italy). From a general standpoint, the simulation code is structured in two distinct parts in order to guarantee greater speed and fluidity. The first part is dedicated to the input of the characteristic parameters of the source (the minimum angle, the maximum angle, and the step of emission of the rays) and of the surrounding environment (depth, temperature, salinity, pH). These data are fundamental to evaluate the water absorption coefficient α (Appendix A). Subsequently, the program, through a series of iterative cycles, will provide the different propagation paths of the rays emitted by the source. In the second part of the program, the user has the possibility to set up a matrix of virtual receivers through which to map the trend of the Transmission Loss (TL) and the SPL sound pressure level. The global layout is represented in Figure 3.
segments and arcs that are connected along the water column, Figure 2. Relying upon this physical principle, the numerical prediction has been developed in this paper.

Numerical Prediction Code: Objectives and Functionality of the Program
The simulation program developed and described within this paper implements the basic foundations of this theory. The main purpose was to evaluate the Transmission Loss (TL) in water, and therefore the calculation of the sound pressure level (SPL) of a given point-like source that radiates a sound wave. The phases necessary to achieve this goal will be explained below, giving the reader the opportunity to get the program steps and then describing the applicative example about the noise analysis of the Adriatic Sea (Italy). From a general standpoint, the simulation code is structured in two distinct parts in order to guarantee greater speed and fluidity. The first part is dedicated to the input of the characteristic parameters of the source (the minimum angle, the maximum angle, and the step of emission of the rays) and of the surrounding environment (depth, temperature, salinity, pH). These data are fundamental to evaluate the water absorption coefficient α (Appendix A). Subsequently, the program, through a series of iterative cycles, will provide the different propagation paths of the rays emitted by the source. In the second part of the program, the user has the possibility to set up a matrix of virtual receivers through which to map the trend of the Transmission Loss (TL) and the SPL sound pressure level. The global layout is represented in Figure 3.

Propagation Paths
In evaluating the propagation paths, it was assumed that the velocity varies linearly with the depth through a law of this type (6): where c 1 is the surface speed, g 1 the speed gradient, and z n the various levels in which the water column is subdivided. The variation of the incident ray along the water column is evaluated as a function of the velocity along the different levels.
According to the Snell-Descartes law at each level: -If the angle of the i-th ray is greater than the critical angle, the latter will continue to propagate according to its direction until it reaches the bottom or the surface of the sea and then will reverse its direction. -If the angle of the i-th ray is greater than the critical angle, the latter will obviously reverse its direction of propagation before reaching the bottom.
Finally, one set of rays will reverse their direction when they interfere with the sea floor and surface while other rays will reverse their direction before reaching the bottom, tracing arch-shaped trajectories, Figure 4a. Once the different propagation paths have been obtained, the second part of the program requires the definition of the map of receivers such as to discretize the fluid domain, Figure 4b.

Propagation Paths
In evaluating the propagation paths, it was assumed that the velocity varies linearly with the depth through a law of this type (6): where c1 is the surface speed, g1 the speed gradient, and zn the various levels in which the water column is subdivided. The variation of the incident ray along the water column is evaluated as a function of the velocity along the different levels.
According to the Snell-Descartes law at each level: -If the angle of the i-th ray is greater than the critical angle, the latter will continue to propagate according to its direction until it reaches the bottom or the surface of the sea and then will reverse its direction. -If the angle of the i-th ray is greater than the critical angle, the latter will obviously reverse its direction of propagation before reaching the bottom.
Finally, one set of rays will reverse their direction when they interfere with the sea floor and surface while other rays will reverse their direction before reaching the bottom, tracing arch-shaped trajectories, Figure 4a. Once the different propagation paths have been obtained, the second part of the program requires the definition of the map of receivers such as to discretize the fluid domain, Figure 4b.

Calculation of the Loss of Propagation (TL)
When a wave propagates from a source, a reduction in intensity due both to the absorption of acoustic energy and to a loss of geometric diffusion (Appendix B) could be considered. Absorption was taken into account for the initial stage of the program, aiming to calculate the transmission loss ratio as follows in the Equation (7): In the case of rays incident on the sea floor, an absorption coefficient equal to 0.5 has been assumed, determining, therefore, an increase of 3 dB in the propagation loss according to Equation (8): where n is the number of reflections on the sea floor.

Calculation of the Sound Pressure Level (SPL)
The SPL magnitude with respect to the acting pressure p is defined as the following quantity (9):

Calculation of the Loss of Propagation (TL)
When a wave propagates from a source, a reduction in intensity due both to the absorption of acoustic energy and to a loss of geometric diffusion (Appendix B) could be considered. Absorption was taken into account for the initial stage of the program, aiming to calculate the transmission loss ratio as follows in the Equation (7): In the case of rays incident on the sea floor, an absorption coefficient equal to 0.5 has been assumed, determining, therefore, an increase of 3 dB in the propagation loss according to Equation (8): where n is the number of reflections on the sea floor. The SPL magnitude with respect to the acting pressure p is defined as the following quantity (9): L P = 10 log 10 p p re f 2 [dB] (9) in which p ref represents the reference pressure evaluated at 1 m from the source and equal to 1 µPa. In order to evaluate the SPL due to the various sound fields, the linear sum of the pressures is automatically performed in Matab ® (The MathWorks, Inc. Natick, MA, USA), assuming that: where p i is the i-th sound field pressure calculated at the receiver. The total sound field is obtained from the sum of the contributions of all pressures as per Equation (11):

Case Study Overview
As part of a collaboration between the University of Naples "Federico II" and Fincantieri technical department, the following program was developed to evaluate the propagation of sound in the Adriatic Sea, from a source located at a depth of 10 m from the sea surface. Five case studies were carried out in the same hydrogeological site, considering different depths with the purpose of making a final comparison of the different modes of transmission. The parameters used during the application and the expected results are outlined below.
• SOURCE: A point source located at -10m from the sea surface and with an SPL equal to 100 dB was assumed as a source of disturbance. With a range of rays varying between −90 • and +90 • with a step of 0.5 • .
• SEA AND CHARACTERISTICS: The Adriatic Sea, as mentioned, was the reference marine environment for this research. This sea develops for a length of 800 km and a width of 200 km. This sea can be divided into three characteristic areas, Figure   For the second part of the program, a matrix of 6 × 6 receivers distributed at intervals of equal width was set up as reported in Figure 4b.

Results
Following the logic of the previous paragraphs and exploiting the input data shown, the program will be able to graphically provide the different sound fields. At the end of processing, the program will graphically provide a series of sound fields. The first to be provided are: the direct sound field, the sound field after reflection on the free surface, the sound field after reflection on the seabed. Since sound can be composed, it is observed that the compound sound fields are obtained from the sum of direct and reflected sound fields. In fact, the program provides: the sound field composed between direct and reflected on the free surface, the sound field composed between direct and reflected on the bottom. Finally, the total one is obtained from the composition of the three sound fields, that is, total sound field composed of direct, reflected on the free surface, and reflected on the bottom. This reasoning was carried out for each case study and therefore at each depth. The simulations have been performed using an Intel ® Core™ i7-6500U processor @2.50 GHz with 2 cores. These sound fields (direct, reflected, and combined) obtained for the particular case of maximum height of 250 m are represented in Figure 6. The composed acoustic contours representative of the other heights from the seabed, i.e., 500 m, 750 m, 1000 m, and 1250 m, are reported in Figure 7. The simulations show the importance of considering the reflected fields for the global propagation loss quantification. The global sound field in Figure 6f actually follows almost the same wave pattern as the direct field in Figure 6a but higher in amplitude: the resulting field reaches up to 55 dB against the 50 dB of the direct field only. As the coverage range increases, the SPL decays with a logarithmic law and as a function of the attenuation parameters illustrated above in paragraph 2: the acoustic level reaches a maximum of 40 dB at the distance of 1250 m, Figure 7d.

Results
Following the logic of the previous paragraphs and exploiting the input data shown, the program will be able to graphically provide the different sound fields. At the end of processing, the program will graphically provide a series of sound fields. The first to be provided are: the direct sound field, the sound field after reflection on the free surface, the sound field after reflection on the seabed. Since sound can be composed, it is observed that the compound sound fields are obtained from the sum of direct and reflected sound fields. In fact, the program provides: the sound field composed between direct and reflected on the free surface, the sound field composed between direct and reflected on the bottom. Finally, the total one is obtained from the composition of the three sound fields, that is, total sound field composed of direct, reflected on the free surface, and reflected on the bottom. This reasoning was carried out for each case study and therefore at each depth. The simulations have been performed using an Intel ® Core™ i7-6500U processor @2.50 GHz with 2 cores. These sound fields (direct, reflected, and combined) obtained for the particular case of maximum height of 250 m are represented in Figure 6. The composed acoustic contours representative of the other heights from the seabed, i.e., 500 m, 750 m, 1000 m, and 1250 m, are reported in Figure 7. The simulations show the importance of considering the reflected fields for the global propagation loss quantification. The global sound field in Figure 6f actually follows almost the same wave pattern as the direct field in Figure 6a but higher in amplitude: the resulting field reaches up to 55 dB against the 50 dB of the direct field only. As the coverage range increases, the SPL decays with a logarithmic law and as a function

Discussion
In order to assess the Ray Tracing trend accuracy, the transmission loss has been compared with the Pekeris waveguide-based model; the purpose of such investigation was to deepen the basic understanding of the propagating mechanisms of underwater noise with respect to well-known theoretical schemes, especially for far-field cases. A maximum range of 25 km has been assumed. Below are the graphs representing a final comparison among the case studies, Figure 8.

Discussion
In order to assess the Ray Tracing trend accuracy, the transmission loss has been compared with the Pekeris waveguide-based model; the purpose of such investigation was to deepen the basic understanding of the propagating mechanisms of underwater noise with respect to well-known theoretical schemes, especially for far-field cases. A maximum range of 25 km has been assumed. Below are the graphs representing a final comparison among the case studies, Figure 8.
The Pekeris model represents the simplest one to characterize shallow water acoustics, considering calm sea, constant depth, homogeneous medium non-absorbent seabed. The analyses showed an acceptable level of correlation between the Ray Tracing model, whose decay is a logarithmic function, and the transcended Pekeris model based instead on the propagation of normal modes; an average error of about 3 dB was estimated between the two approaches, which is acceptable in a preliminary design phase. For Pekeris, the rate of dispersion and waves attenuation in the waveguide depends on the iterative search of mode numbers. The wave equation is solved to determine the normal mode solutions which prevail at large distances from the source, according to Lamb procedure. For the present calculation, four normal modes have been supposed: the pressure field is calculated according to a two-dimensional modal extension (Kirchhoff equation). The discrepancy between two approaches could be found, therefore, in the different schematization of reflection and attenuation mechanisms, with the TL increasing with mode number [41]. Due to the broadband nature of the problem, the model is very computationally demanding. The calculation times essentially depend on the pitch of the rays and the distance to be covered. Performance can certainly be optimized by consolidating the simulation code structure. However, for a preliminary phase, these times are still advantageous with respect to finite element models, which required much more computational cost and mesh modeling steps as well. At the end of the processing, the time required for each simulation was evaluated: it is the function of two main parameters: • HEIGHT OF THE SEABED The higher the height of the bottom, the more the processing time required for the simulation increases.

Discussion
In order to assess the Ray Tracing trend accuracy, the transmission loss has been compared with the Pekeris waveguide-based model; the purpose of such investigation was to deepen the basic understanding of the propagating mechanisms of underwater noise with respect to well-known theoretical schemes, especially for far-field cases. A maximum range of 25 km has been assumed. Below are the graphs representing a final comparison among the case studies, Figure 8. The Pekeris model represents the simplest one to characterize shallow water acoustics, considering calm sea, constant depth, homogeneous medium non-absorbent seabed. The analyses showed an acceptable level of correlation between the Ray Tracing model, whose decay is a logarithmic function, and the transcended Pekeris model based instead on the propagation of normal modes; an average error of about 3 dB was estimated between the two approaches, which is acceptable in a preliminary design phase. For Pekeris,  The code allows for inserting the interval between the emission rays of the source; the smaller the step, the more detailed the noise map will be. The computational cost varies quadratically with the depth of the sea, Figure 9. In Table 1, the times required by the program in the five cases.

Conclusions
This paper is part of the macro-topic inherent in the current need to study forecasting systems for optimizing the sources of underwater noise pollution generated by ship traffic. In particular, the activity arises from a collaboration between the University of Naples "Federico II" (industrial engineering department) and the technical office of Fincantieri Group with the purpose of creating a fast predictive code for a preliminary evaluation of the acoustic propagation fields in submarine environments, with special attention to the basin of the Adriatic Sea. The simulation program in Matlab ® is essentially based on the optical-acoustic technique of Ray Tracing, which allows to characterize the sound propagation path also as a function of the characteristic parameters of the source and of the diffusion medium. Methodology is very suitable for preliminary calculations for two reasons: • For the identification of the sound field radiated even at a great distance and, therefore, of the impact this may have on the marine environment. • For the interpretation of the possible experimental values detected by the hydrophones in a field relatively close to the ship for an effective characterization of the disturbance source.
Future developments conducted by the research group aim at expanding the functionality of the program, taking into account also the third spatial dimension, while always limiting the computational cost, which, for applications of this kind, can be excessively expensive. Experimental campaigns of field measurements are also underway to validate this first part of the program.  Quadratic: y = 7e-05*x 2 + 2.10*e-02*x + 8.02e-01 R 2 = 1 Figure 9. CPU time estimation as function of sea depth.

Conclusions
This paper is part of the macro-topic inherent in the current need to study forecasting systems for optimizing the sources of underwater noise pollution generated by ship traffic. In particular, the activity arises from a collaboration between the University of Naples "Federico II" (industrial engineering department) and the technical office of Fincantieri Group with the purpose of creating a fast predictive code for a preliminary evaluation of the acoustic propagation fields in submarine environments, with special attention to the basin of the Adriatic Sea. The simulation program in Matlab ® is essentially based on the optical-acoustic technique of Ray Tracing, which allows to characterize the sound propagation path also as a function of the characteristic parameters of the source and of the diffusion medium. Methodology is very suitable for preliminary calculations for two reasons: • For the identification of the sound field radiated even at a great distance and, therefore, of the impact this may have on the marine environment. • For the interpretation of the possible experimental values detected by the hydrophones in a field relatively close to the ship for an effective characterization of the disturbance source.
Future developments conducted by the research group aim at expanding the functionality of the program, taking into account also the third spatial dimension, while always limiting the computational cost, which, for applications of this kind, can be excessively expensive. Experimental campaigns of field measurements are also underway to validate this first part of the program. Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The theory used to evaluate the water absorption coefficient α is the model proposed by Francois-Garrison, in which the absorption coefficient is expressible into three terms that correspond, respectively, to the contribution of boric acid, of the sulphate of magnesium, and water viscosity: Factors taking into account the contribution of boric acid B(OH) 3 :

•
Factors that take into account the contribution of water viscosity:

Appendix B
When a source emits a sound, it propagates in all directions in the form of surfaces whose sound levels are constant on each surface but decrease as the distance from the source increases. In fact, the sound wave propagating outwards causes a decrease in intensity proportional to the inverse of the surface. Near the source, the sound surfaces are spherical; as the radius increases, the diffusion will not always be spherical since there will be on one side the sea surface and on the other the seabed and, therefore, from a spherical diffusion, we pass to a cylindrical diffusion. For the two types of diffusion, the following (TL) transmission loss can be considered: The transmission loss can be evaluated also taking into account the attenuation coefficient α. In this case, the relationship becomes: