Measured parameters for different water types.

## Abstract

The radiative transfer equation (RTE) is a theoretical framework that can be used for predicting and interpreting underwater light fields in terms of the constituents of natural water bodies. However, the RTE is a complex integrodifferential equation and deriving exact solutions for it is a difficult task. In this chapter, we aim to present some details regarding Monte Carlo simulations and how this method may be applied to solve the RTE numerically. By solving the RTE, one may accurately predict the received power and estimate the channel bandwidth and several other measurable parameters with regard to multiple water conditions. Simulations will also be presented.

### Keywords

- RTE
- underwater
- optical
- wireless
- channel
- propagation
- photons
- UOWC

## 1. Introduction

When compared with free space optical (FSO) communication channel, the underwater optical wireless communication (UOWC) shows unique characteristics. Because of this, the channel models usually used in FSO are not suited for UOWC, and different channel models must be developed to describe the different photon interactions under play.

Mobley [1] classified the optical properties of water into two different groups, inherent optical properties (IOPs) and apparent optical properties (AOPs). The two major IOPs that will attenuate light propagation in UOWC are absorption and scattering. Hence, any attempt in modeling the UOWC channel will need to consider absorption and scattering effects within specific link configurations.

Several theoretical models were employed by researchers to model the UOWC channel, being the Beer-Lambert law the most employed due to its simplicity [2]. Despite its simplicity, the Beer-Lambert law contains two implicit assumptions: that the transmitter and receiver are perfectly aligned and that all scattered photons are lost. The assumption that all photons that undergo scattering are lost will severely underestimate the received optical power, especially in water environments where scattering is the dominant IOP.

Another theoretical model for modeling the underwater channel is the radiative transfer equation (RTE) [3]. The RTE simply says that as a beam of radiation travels through a medium, in this case water, it will lose energy to absorption, gain energy by emission, and redistribute energy by scattering. Yet, as the RTE is an integrodifferential equation involving different independent variables, an exact analytical solution is hard to find. In view of this, solving the RTE numerically is a preferred approach, being the most popular one the Monte Carlo simulation.

## 2. The radiative transfer equation

One important law of geometrical radiometry is the n-squared law for radiance. To derive it, we consider two mediums separated by a transparent surface; in that case, Snell’s law states that the angles of an incident ray from medium 1 to medium 2 are related by:

where *n* is the index of refraction of the medium. Considering that a ray is simply a narrow beam of photons traveling in almost the same direction, and if we consider two narrow rays, incident and refracted, that will have solid angles given by:

and as seen in [3], squaring each side of Eq. (1) and multiplying by the Azimuthal spread,

which is known as the Straubel’s invariant. Consider now the radiances of the two rays, incident and refracted, defined as:

where

from which we can obtain:

which is the n-squared law for radiance. For the case of

As seen in [3], this result is known as the fundamental theorem of radiometry and it states that the radiance divided by the refraction index squared is constant along any path; however, because all real substances will cause some absorption and scattering on the incident photons, the validity of theorem is restricted for paths in vacuum.

If we consider that a ray of radiance

where

Considering a beam of radiance traveling from point

being

where

The operator

and the gradient operator

From these developments, we can rewrite Eq. (11) as:

valid for photons traveling with speed

Writing the expression for the change in

In Eq. (17),

For UOWC systems, we are interested in the time-independent RTE in horizontally homogenous water bodies with a constant index of refraction. Because of this, the factor

It is also usual [3] to combine

and the path function for elastic scattering,

Combining all these terms in Eq. (18), we can come up with the standard form of the RTE. As a simple statement, the IOPs and boundary conditions will go into the RTE and the final radiance will come out as a result. The standard form of the RTE, as shown in [4], is given by:

where

### 2.1 Inherent optical properties

During the interaction of a photon with a water molecule, one of two things may happen: the photon may be absorbed, leaving the water molecule in a state with higher internal energy, or the photon may undergo scattering. Such scattering occurs when there is a change in the direction of the propagation of the photon or in the photon energy—or in both.

When scattering happens, if the molecule that interacted with the photon returns to its original internal energy state by emitting a photon of the same energy as the absorbed one, the scattering process is called elastic scattering; on the other hand, if the excited molecule emits a photon of lower energy, i.e., a photon with a longer wavelength, the incident photon goes through a process called inelastic scattering. As seen in [3], the absorption coefficient, that will be introduced later in this chapter, accounts for both the conversion of radiant energy into heat and the loss of power at wavelength

To derive a mathematical coefficient for scattering and absorption, we assume the model shown in Figure 1.

A volume of water

Using Eq. (22), we can define the ratio between the absorbed and incident power as the absorptance and the ratio between scattered power and incident power as scatterance. As already stated, when solving the RTE the IOPs usually employed are the absorption and the scattering coefficients which are the absorptance and scatterance per unit distance of the medium. Taking the limit of these terms as the water thickness

The beam attenuation coefficient present in the RTE (Eqs. (16) and (20)) is written as:

Another useful IOP showing in the RTE is the single scattering albedo,

In waters where the beam attenuation is due primarily to scattering,

In [3], it is possible to find the absorption, scattering, and the attenuation coefficients for several water types and wavelengths. The usual parameters considered for UOWC systems are the ones presented in Table 1.

Water type | ||
---|---|---|

Clear waters | 0.25 | |

Coastal waters | 0.55 | |

Harbor I waters | 0.83 | |

Harbor II waters | 0.83 |

As mentioned before, the ratio between the scattered power and incident power is referred as scatterance, that is the fraction of incident power scattered out of the beam through an angle

In [3], we can see that the spectral power scattered into the given solid angle

if the incident power falls on an area

As said in [5], the form of Eq. (28) suggests the name volume scattering function (VSF), that is, the scattered intensity per unit incident irradiance per unit volume of water. Integrating

Normalizing Eq. (26) with the scattering coefficient, we obtain the scattering phase function (SPF):

The scattering phase function can be interpreted as a probability density function (PDF) for scattering from an incident direction

## 3. Solving the radiative transfer equation: Monte Carlo simulation

The Monte Carlo name was coined by Nicholas Metropolis [7]. The modern version of the method was invented in the late 1940s and found use early use in the studies of neutron transport for the design of nuclear weapons [8, 9]. More recently, Monte Carlo techniques evolved and are now used to solve problems in physical sciences, economics, engineering, and pure mathematics.

Monte Carlo methods may be used to solve any problem that have a probabilistic interpretation, meaning that if the probability of occurrence of each separate event in a sequence of events is known, one can determine the probability that the entire sequence of events will occur. By tracing the fate of millions of photons according to statistical probabilities, a solution for the RTE is built. Each simulated photon path is randomly distinct from the others, as determined by the probabilities of absorption and scattering in the underwater channel.

The Monte Carlo simulations of photon propagation offer a flexible yet rigorous approach toward solving numerically the RTE. In the last decades, the exponential improvement in computer speeds attenuated the problem of the high computational cost of Monte Carlo simulations.

In the following subsections, a detailed description for building a solution for the RTE, step by step, will be shown.

As seen in [10], there are four main Monte Carlo methods for photon migration in turbid media, that is, four different ways to build photon’s trajectories. The four methods are referred as the albedo-weight method (AW), the albedo-rejection method (AR), the absorption-scattering path length rejection method (ASPR), and the microscopic Beer-Lambert law method (mBLL).

The AW method considers a source emitting packets of many photons and after each interaction a fraction of the photons in the packet are absorbed, i.e., lost. In this type of simulation, the packet is emitted with an initial weight of 1 and at each interaction the current weight is multiplied by the albedo to account for the loss of photons by absorption and the packet continues propagating in the scattered direction with a reduced weight. This method has in its favor that there are no wasted computations because all photon packets eventually reach the detector. This approach was favored by several research papers [11, 12, 13, 14, 15, 16].

In this chapter, we present a simulation using the AR method; the main advantage of this method is that it mimics what happens in nature by tracking individual photons emitted by some source. Naturally for this case, the computational time of the absorbed photons are wasted.

It is seen in the investigation conducted by Sassaroli that the total probability of detecting a photon is equivalent in the four methods and they have statistical equivalence, with some differences regarding the convergence time. When comparing the AW to the AR method, as one may expect, the AR will require, in most scenarios, more launched photons for convergence. However, it is also noted that in the AR method, the photon is detected or lost similarly to what happens in a time-resolved experiment in which photons are detected using the time-correlated single photon counting technique and the noise on the simulated response reproduces the Poisson statistics that typically characterize the noise in experiments. Therefore, as we already stated, the AR method reproduces a more physical simulation of propagation than the AW method, where the simulated photons are energy packets rather than physical photons. In Figure 2, we show a flow chart of the Monte Carlo algorithm used in this simulation, each step in this chart will be explained with further detail in the following sections.

### 3.1 Photon path length

The photon path length, or step size, determines how far the photon will travel before encountering a water molecule or particle, meaning the distance a photon travels between an absorption event and a scattering event. The path length is directly related to the water type and, consequently, to the beam attenuation coefficient introduced in Section 2.1, and shown for several water types in Table 1.

As stated in [4], from the derivation of the RTE, the radiance in a particular direction

where

putting it in terms of the photon path length

The exponential decay of radiance can be explained in terms of the fate of individual photons if the probability of any photon being absorbed or scattered out of the beam between photon path lengths

Mobley [3] notes that

with

selecting a random number

For homogenous water, where *c*(*r*) does not depend on *r*, where

### 3.2 Absorption events and scattering angles

Recalling that the photon path length, or step size, is the geometric distance a photon will travel between a scattering event and an absorption event, after determining the distance the photon will travel from origin, one must decide if the next interaction will be an absorption or a scattering event. For this, we recall the definition of the single scattering albedo,

From this description, one may clearly see where the name “Albedo-Rejection” comes from: if the photon is absorbed, it must be terminated, and its position is recorded for further analysis; if the photon was scattered, the scattering angles must be calculated, and its direction of propagation is updated.

#### 3.2.1 The Henyen-Greenstein phase function

Unlike the FSO channel, in an underwater environment, the photons will encounter a much larger number of particles, and thus multiple scattering events will play a more important role in the received optical power, when compared to other optical wireless communication channels.

One of the models proposed to describe the SPF is the Henyen-Greenstein phase function (HG) [12, 17, 18]. The HG was originally proposed by Henyen and Greenstein to describe the scattering of light by interstellar dust [19]:

where *g* is referred to as the anisotropy factor and is also the mean cosine of the scattering angle

To find the scattering angle, Eq. (41) must be solved for

where *q* is a uniformly distributed random number between 0 and 1. The complete method for obtaining Eq. (42) may be found elsewhere [4].

Because scattering is a three-dimensional process, besides the polar angle (

Comparing the HG phase function using Eq. (41) with the Petzold’s avg. Phase Function extracted from [3], we can see in Figure 3 that the HG is indeed an adequate approximation for the SPF.

The HG phase function found great use in ocean optics simulations and is heavily used in Monte Carlo simulations and other relevant channel modeling works by several researchers [20].

#### 3.2.2 The two-term Henyey-Greenstein phase function

Despite being an approximate fit, it can be seen in Figure 4 that the HG phase function shows noticeable differences from the Petzold’s measurements for small angles,

To overcome this limitation, Haltrin [21] proposed a two-term Heyen-Greenstein (TTHG) phase function that matched Petzold’s results better at small scattering angles and provided the elongation seen at larger angles. The TTHG is given by:

where

Haltrin also derived a relationship that relates the scattering coefficient, *b*, with the anisotropy factor,

After calculating the parameters using Eqs. (44)–(47) and computing the phase function for clear ocean waters with

#### 3.2.3 Haltrin empirical phase function

Despite the attempts from several researchers of developing analytical phase functions for natural waters, we can see that the ones described in the previous sections do not provide a good fit over the total range of

In [23], Haltrin developed an empirical fit using strong regression for all Petzold’s measurements. One of the strongest points of this model is that the Phase Function may be obtained only by using the scattering coefficient and the single scattering albedo as inputs:

and the scattering phase function:

with:

where *b* is the scattering coefficient and the remaining parameters are:

The strong regression given by these equations can be used for an empirical model for the phase functions. As seen from Figure 7, this model represents a better fit to the Petzold’s measurements than the phase functions presented before, especially in the regions where the phase function is the strongest.

In our simulations, we found that when comparing the CDF computed by directly interpolating the Petzold’s measurements, the CDF computed from the empirical phase function showed some divergence. Our CDF showed that for harbor waters with

The above presented divergence resulted in a small difference in the final received power; however, the empirical phase function largely overestimated the channel bandwidth, especially in harbor waters. This is due to the smaller scattering angles generated by this empirical phase function, causing the photons to travel in a straighter line when compared to the case of scattering with the CDF computed from Petzold’s data.

#### 3.2.4 The Fournier-Forand phase function

As seen the previous sections, several analytical and empirical formulae were used as empirical fits for ocean waters but none, specially the analytical ones provide a perfect fit for all scattering angles.

In [24], Fournier and Forand derived an analytic expression for the phase function of a set of particles that have a hyperbolic particle size distribution, with each particle scattering according to the anomalous diffraction theory. The Fournier-Forand (FF) phase function is given by:

where

Being *n* in Eq. (59) the real index of refraction of the particles,

It can be seen in [25] that with

When comparing the CDF computed from the FF phase function and the CDF from Petzold’s data for harbor waters, we get very close probabilities over all scattering angles.

#### 3.2.5 Sahu and Shanmugam (SS) phase function with flat back approximation

In [25], a semi-analytical model for a phase function was proposed by Sahu and Shanmugam (SS). The SS phase function predicts the light scattering properties for all forward angles, namely

where

with

Sahu found, in his most recent works [14, 26], that the set of values,

#### 3.2.6 Computing phase functions for Monte Carlo simulations

As seen in [27], phase functions satisfy the normalization condition:

and to determine the polar scattering angle:

where *q* is a uniformly distributed random number between 0 and 1.

However, given the shape of most phase functions, it is more efficient [27] to use the equation for the phase function to build up a table of

For the results presented in Table 2, we computed a CDF directly from Petzold average phase function extracted from [3] and used the same method described in [27] to obtain values for all angles between 0 and

P(θ _{PETZ}) | P(θ _{FF}) | P(θ _{SS}) | P(θ _{TTHG}) | |
---|---|---|---|---|

0.009094 | 0.012355 | 0.004059 | 0.002382 | |

0.031941 | 0.041314 | 0.025992 | 0.011982 | |

0.046044 | 0.057930 | 0.039860 | 0.020878 | |

0.062621 | 0.076605 | 0.056320 | 0.034634 | |

0.081934 | 0.097344 | 0.075422 | 0.055321 | |

0.104274 | 0.120254 | 0.097285 | 0.085451 | |

0.129946 | 0.145647 | 0.122155 | 0.127856 | |

0.158838 | 0.173561 | 0.149991 | 0.184165 | |

0.191217 | 0.204415 | 0.181067 | 0.254854 | |

0.352290 | 0.358677 | 0.335353 | 0.596932 | |

0.498215 | 0.505639 | 0.472517 | 0.788187 | |

0.659192 | 0.665537 | 0.630018 | 0.897428 | |

0.891501 | 0.882151 | 0.892093 | 0.974829 | |

0.958125 | 0.951993 | 0.960280 | 0.989702 | |

0.986110 | 0.983167 | 0.987188 | 0.996168 | |

0.994130 | 0.992771 | 0.994681 | 0.998674 | |

0.997280 | 0.996631 | 0.997562 | 0.999169 | |

1 | 1 | 1 | 1 |

Parameter | |
---|---|

Simulated beam | Gaussian beam |

Beam divergence at | 0, 10, 20 and 30 mrad |

Beam width | 1 mm |

Receiver diameter | 0.1 m |

Receiver FOV | |

Photons per simulation | |

Phase function | Fournier-Forand |

Parameter | |
---|---|

Simulated beam | Gaussian beam |

Beam divergence at | 0 |

Beam width | 1 mm |

Receiver diameter | 0.1 m |

Receiver FOV | |

Photons per simulation | |

Phase functions | Petzold average, FF, SS and HG |

In Table 2 and Figures 10 and 11, we compare the CDF obtained for three phase functions, the SS, FF, and the TTHG with

In this table, we can see that both the FF and SS phase functions are a great approximation for the Petzold’s average phase function, while the TTHG heavily overestimates the amount of scattering up to

### 3.3 Temporal dispersion and the underwater channel bandwidth

#### 3.3.1 Temporal dispersion

One of the most useful information regarding the underwater channel that may be provided by the Monte Carlo simulation, which is unavailable in most simple models like the Beer law, is the channel temporal dispersion. As seen before in this chapter, seawater is a dispersive medium in which light will suffer the effect of multiple scattering events. Scattering will cause the photons to arrive at the receiver at different time intervals causing temporal dispersion and a penalty in the channel bandwidth. In recent years, several researchers employed the double gamma functions to model the impulse response in underwater channels [15, 28]. Despite the accuracy of the double gamma function models, we find that the method described below consists in a simpler method by means of tracking the individual time of propagation of each photon.

By tracking the propagation distance of each photon (

where

with the index of refraction of sea water at

For a better visualization of the impulse response, the time delay may be normalized by the straight-line distance, or direct distance, from source to the receiver using:

Making a histogram of the received photons and distributing the received intensity in time slots defined as the bin widths, one can get a result for the channel impulse response (CIR) of the underwater channel. The results will be presented in Section 4.

#### 3.3.2 Underwater channel transfer function

When using laser diodes in UOWC systems, it is desirable to modulate the light source at very high frequencies. The CIR generated by the Monte Carlo simulation can be converted to frequency domain by the means of a Fourier Transform algorithm, from which we can obtain the channel transfer function.

Given that the impulse response is approximated by a discrete histogram of the time of arrival, it can be transformed to a frequency response, as proposed in [18, 19] by means of a discrete Fourier Transform (DFT), with the usual definition:

To numerically implement this, the fast Fourier transform (FFT) algorithm can be used to convert the

where

As expected, due to Nyquist sampling theorem, the maximum frequency that can be approximated in the frequency domain will be related to the time step by:

## 4. Simulation

In this section, we will provide some simulation results for several emitter parameters, underwater conditions, and different phase functions.

In the first part, we will simulate the received optical power for clear and coastal waters for different beam divergence values and compare it with the popular Beer-Lambert law. It can be seen in [14] that in these water conditions, the delay spread is so small that they can be considered as a nondispersive medium even for high data rate communications, and for this reason the channel impulse response (CIR) will be studied only for conditions where the medium is dispersive.

In the second part, we will study the effect of different phase functions in Harbor waters with a perfectly collimated beam. We believe that this is an interesting approach since in these waters’, due to the higher scattering coefficient

### 4.1 Clear and coastal waters

#### 4.1.1 Simulation results

In Figure 12, we compare the received power for clear ocean waters with varying beam divergence parameters of 0, 10, 20, and 30 mrad for an emitted beam with a FWHM width of 1 mm. For the phase function, we used the FF, since, as seen in Section 3.2.4, it provides a good fit for Petzold’s phase function. However, as mentioned in the “Introduction,” when discussing simulation, the impact of different phase functions, especially those which are a good approximation for the Petzold’s average phase function, is small for the clear and coastal waters due to the low contribution of scattered photons on the final received power.

Our simulation matches very well the Beer-Lambert law for the case of an emitted Gaussian beam with no divergence, since one of the main assumptions of the law is that the emitted beam is perfectly collimated, and as the beam divergence increases the power loss becomes more accentuated.

We note that for a beam with no divergence, after 10 m, the normalized received optical power would be

we may expect a beam spread of 20 cm after 10 m for a beam divergence of 10 mrad, potentially reducing by more than half the received optical power (due to the limited aperture of 10 cm). In fact, we can verify from Figure 12 that the received power for a 10 mrad divergence is

To further illustrate this, we show, in Figure 14, the geometrical losses corresponding to beam divergence. Given that for a perfectly collimated beam in clear ocean waters, power loss is mainly due to absorption and the additional power loss solely due to the spread of the Gaussian beam.

For the case of coastal waters with

As was the case for clear waters, the received power drops considerably with an increase in the beam divergence parameter.

Here, the contribution of scattered photons to the final received power increases: after 20 m, 25% of the received power came from scattered photons, and after 30 m, this contribution is elevated to almost 30%. This is further illustrated in Figure 16, where the scattering histograms for

Table 5 gives the received normalized power for clear ocean waters for various distances and beam divergence parameters. In the following table (Table 6), the same results are given for coastal waters.

Distance | ||||
---|---|---|---|---|

10 m | 30 m | 50 m | 70 m | |

0 | ||||

10 | ||||

20 | ||||

30 |

Distance | ||||
---|---|---|---|---|

10 m | 20 m | 30 m | 40 m | |

0 | ||||

10 | ||||

20 | ||||

30 |

### 4.2 Harbor waters

In this section, we compare the performance of different phase functions for Harbor I and Harbor II waters. First, we interpolated values for the Petzold average phase function and considered, as in Section 3.2.5, the benchmark for all other phase functions. The amount of unique scattering angles generated by this interpolation was enough for each photon to have a unique scattering angle.

The SS phase function, seen in Section 3.2.4, predicts scattering for all forward angles. For the angles between 90 and

where

with the values for

Using the values of

#### 4.2.1 Simulation results

For the received optical power for Harbor I and Harbor II waters, we can see in Figures 18 and 19 that the results obtained when using the FF and the SS phase functions are closely matched to the results obtained using the interpolated Petzold average PF, meanwhile the one term HG severely underestimates the received optical power. The main reason for this, as we have seen before, is that the HG is unable to reproduce the high amount of scattering that occur at small angles, as is the case in the Petzold average PF and matched by the other analyzed phase functions.

In Figure 20, the received beams are compared for Harbor I and Harbor II waters for

When looking at the CIR for Harbor I waters with *z* = 15 m, plotted in Figure 21, we can see that the results for the Petzold, FF, and the SS phase functions show similar behavior on the delay profile, while the CIR for the HG phase function shows a considerable number of delayed photons.

Figure 22 plots the transfer function of the channel, from which the 3-dB bandwidth can be extracted, for Harbor I waters corresponding to *z* = 15 and 20 m using the method described in Section 3.3.2. It is possible to verify from these results that the behavior of the SS phase function closely matches the one for the Petzold PF for both distances and the HG underestimates the bandwidth by a large margin in both scenarios. The FF phase function on the other hand, overestimates the bandwidth for the case *z* = 15 m. The 3-dB bandwidth is limited to ∼300 and 90 MHz, for 15 and 20 m, respectively.

The same analysis is now applied to Harbor II waters. The impulse response is shown in Figure 23. The interpretation is that the HG phase function exhibits a response where the power is much more distributed over time when comparing to what is seen in the other phase functions where most of the power is located at the instant

The effect of the delay spread is even more noticeable when we analyze the transfer function. In Figure 24, this analysis is done for *z* = 10 and 12 m, where we see that in these conditions, due to higher number of scattering events that each photon undergoes, the bandwidth is lower than what we found at Harbor I waters, even when considering smaller propagation distances. In this case, the phase functions FF and SS give similar results to the Petzold average PF, the bandwidth being limited to ∼150 and 100 MHz, for distances of 10 and 12 m, respectively.

In the following table, Table 7, we summarize the most relevant numerical results for Harbor I and Harbor II waters as we did for the case of clear and coastal waters in Section 4.1. For brevity, we only show the results obtained using the Petzold average phase function; however, as noted in the results presented previously, one may expect similar numerical results when using the FF or the SS phase function as they predict scattering angles similar to those predicted by Petzold.

Received power | 3-dB BW (MHz) | ||
---|---|---|---|

HI z = 15 m | 2.93 | 310 | |

HI z = 20 m | 6.18 | 80 | |

HII z = 8 m | 2.06 | 440 | |

HII z = 10 m | 3.51 | 185 | |

HII z = 12 m | 4.78 | 110 | |

HII z = 14 m | 7.87 | 53 |

### 4.3 Validation of simulation results

For further validation of the simulation model, we compared the values obtained in our simulation with the ones obtained by Sahu in [14]. For this, we used the same parameters used in his simulations, a Gaussian beam profile with FWHM beam width of 2 mm, a beam divergence parameter of 1.5 mrad, and a receiver aperture of 10 cm.

The comparison for clear and coastal waters between our simulation and [14] is shown in Figure 25, and we can see that our simulation closely matches the results obtained by Sahu.

For Harbor waters, Sahu only considered Harbor II waters with

Furthermore, Sahu defined a delay spread quantity which is the time period over which the CIR falls to

## 5. Conclusion

In this chapter, a simple yet powerful tool for modeling the propagation of photons in an underwater channel is presented, which provides more accurate results than the conventional Beer-Lambert law. The algorithm presented here highlights the fundamental processes of absorption and scattering, hopefully helping the reader to better understand the physics of photon propagation in an underwater environment.

In the first part of the simulation, we showed how important the collimation of a beam is and how the beam spread can cause losses up to −13 dB even after only 10 m. It was also possible to verify that in clear waters, a good signal to noise ratio is achievable for several tens of meters in both clear and coastal waters.

In the second part of the simulation, the impact of the water turbidity on the underwater link was assessed, where it was concluded that increasing turbidity limits not only the received optical power, but also the maximum bandwidth of the channel. An adequate choice of phase function is important as simpler phase functions, such as the Henyey-Greenstein, may wrongly predict the numerical results in these conditions.

For the first time results for the received power, channel impulse response and 3-dB bandwidth are obtained from Monte Carlo RTE simulations for different phase functions. The SS and FF are compared with the Petzold average PF and results are shown to agree well. Moreover, for validation purposes, our simulation results are further compared with those obtained by Sahu [14] and good agreement is shown to exist.

We believe that the approach presented here provides a general and flexible technique for numerically solving the radiative transfer equation, accessible to any reader with basic programming skills. Besides the conditions simulated in this chapter, the simulation program may be used to simulate innumerous other conditions, like misalignment between receiver and transmitter, smaller aperture sizes, limited FOVs, and other light sources other than Gaussian beams.

## Acknowledgments

This work is financed by the ERDF—European Regional Development Fund through the Operational Programme for Competitiveness and Internationalisation—COMPETE 2020 Programme, and by National Funds through the FCT—Portuguese Foundation for Science and Technology, I.P., within project POCI-01-0145-FEDER-031971.