Channel Modeling for Underwater Acoustic Network Simulation

Simulation forms an important part of the development and empirical evaluation of underwater acoustic network (UAN) protocols. The key feature of a credible network simulation model is a realistic channel model. A common approach to simulating realistic underwater acoustic (UWA) channels is by using specialised beam tracing software such as BELLHOP. However, BELLHOP and similar modeling software typically require knowledge of ocean acoustics and a substantial programming effort from UAN protocol designers to integrate it into their research. This paper is a distilled tutorial on UWA channel modeling with a focus on network simulation, providing a trade-off between the flexibility of low level channel modeling via beam tracing and the convenience of automated channel modeling, e.g. via the World Ocean Simulation System (WOSS). The tutorial is accompanied by our MATLAB simulation code that interfaces with BELLHOP to produce channel data for UAN simulations. As part of the tutorial, we describe two methods of incorporating such channel data into network simulations, including a case study for each of them: 1) directly importing the data as a look-up table, 2) using the data to create a statistical channel model. The primary aim of this paper is to provide a useful learning resource and modeling tool for UAN protocol researchers. Initial insights into the UAN protocol design and performance provided by the statistical channel modeling approach presented in this paper demonstrate its potential as a powerful modeling tool for future UAN research.


I. INTRODUCTION
R ECENT developments in underwater acoustic modem capabilities [1]- [4] will make large scale underwater acoustic networks (UANs) feasible in the near future.Such large scale UAN deployments will have a wide range of applications, e.g.water quality monitoring [5], seismic monitoring [6], marine animal tracking [7], off-shore asset monitoring [8], and ocean exploration using autonomous underwater vehicles (AUVs) [9].However, compared with terrestrial radio systems, the performance of UANs is severely limited by the adverse characteristics of the underwater acoustic (UWA) communication medium [10]: extremely slow propagation (sound speed through water is approximately 1500 m/s), low available bandwidth (typically on the order of several kHz), large multipath delay spread and Doppler effect.These challenging channel characteristics necessitate the design of networking protocols dedicated specifically to UANs [11] [12].
The development, testing and validation of UAN protocols involve two principal steps: simulations and sea experiments.In addition to circumventing the high cost and logistical challenges involved in performing sea experiments, the major advantage of simulation-based studies is that they enable researchers to test their network protocols under controlled, reproducible conditions, and obtain more comprehensive, statistically valid results, e.g.via parameter sweeps, Monte Carlo simulations etc.In contrast, implementing and testing the network protocols at sea is more suitable as a validation step to prove that they work in a real deployment.It is usually not logistically feasible at sea to perform parameter sweeps, benchmark comparisons, and obtain large statistical samples of the network protocol performance.Instead, a UAN sea experiment is usually a demonstration of the network operating in a specific environment.Therefore, simulation is of particular importance in performing a thorough empirical evaluation.
One of the key challenges in developing a credible network simulation model is a realistic representation of the UWA channel characteristics.Generally, the channel models found in the UAN protocol literature can be split into three categories: • Binary range-based model.The simplest way to model a UAN communication environment is to derive a binary connectivity pattern among the nodes based on a fixed connection range (e.g. if the distance between any two nodes is less than the maximum connection range, there is a link between them) and to assume a fixed propagation speed of 1500 m/s, e.g.[13] [14].Although this is a simple and intuitive approach that is useful for theoretical UAN protocol development, it oversimplifies the behaviour of a realistic UWA channel.• Analytical transmission loss model (often referred to as the Urick model [15]).This model takes the above approach a step further and calculates the transmission loss on every link using mathematical expressions for distance-related spreading loss and frequency-related absorption loss [16].In contrast with the range-based model, it gives a measure of the received signal strength, allowing the researchers to estimate the Signal-to-Noise Ratio (SNR) and the Signal-to-Interference-plus-Noise Ratio (SINR).However, this model still omits many typical features of UWA channels, e.g.shadow zones due to acoustic wave refraction, delay spread and frequency selective fading due to multipath.• Specialized channel modeling software.In order to model more advanced characteristics of the UWA channel listed above, specialized simulation models are required, e.g. based on ray/beam tracing or normal mode calculations [15].A popular open source platform for this is BELLHOP [17] [18], which employs beam tracing to predict acoustic pressure fields in specified underwater environments.There are multiple extensions to BELL-HOP that enable the researchers to adopt it in their studies, e.g.VirTEX [19] for simulating time-varying UWA channels, or the World Ocean Simulation System (WOSS) [20] for simulating a UAN in an environment representing a specified geographical location (based on real measurements).
The use of channel modeling software is a common approach to obtaining realistic representations of UWA channels.However, the UAN protocol researchers, especially those coming from the terrestrial wireless communications background, face a steep learning curve in ocean acoustics when learning how to model the UWA channel correctly, e.g.setting up the environment using BELLHOP and interpreting the beam tracing results.To alleviate this problem, the WOSS simulation platform [20] abstracts the user from the low level BELLHOP channel modeling process and enables them to simply specify the desired geographical location of the nodes and allow WOSS to set up BELLHOP automatically with the right environmental parameters measured in sea experiments.WOSS can be integrated with any C++ based network simulator, e.g.ns2-MIRACLE [21] or ns-3 [22], and is widely used as part of the well-established underwater network simulation/emulation suites, e.g.DESERT [23], SUNSET [24].However, we argue that learning about the key characteristics of UWA propagation via a more hands-on channel modeling process provides the UAN protocol researchers with valuable insights into the communication environment that they are investigating.
In this paper, we aim to bridge the gap between low level channel modeling via BELLHOP beam tracing (or similar software) and automated channel modeling via WOSS, by providing a detailed tutorial with MATLAB simulation code [25], that focuses on several key characteristics of the UWA channel most relevant for networking protocol design -signal attenuation, propagation delay, multipath fading and delay spread.As such, our proposed simulation framework does not aim to replace the established fully integrated platforms, such as WOSS, nor to replace the standard BELLHOP beam tracing interface designed more widely for ocean acoustics research.Rather, the main purpose of the simulation framework proposed in this paper is to make beam tracing accessible for the underwater networking research community.The main contributions of this paper can be summarized as follows: • Survey of existing channel simulators -we provide an overview on the features, capabilities and relative merits of the state-of-the-art in UWA channel simulation with the focus on networking research; • Tutorial on UWA propagation -we give a detailed tutorial on the UWA communication environment, focusing on the features most relevant for network simulations; • BELLHOP-based channel simulation platform -the tuto-rial is accompanied by our user-friendly MATLAB code that creates channel models from basic (for a simple introduction) to more advanced UWA environments using BELLHOP; • Integration of the channel data into network simulatorswe also propose a framework for processing our channel simulator data and integrating it into network simulations, including the demonstration of this approach in two case studies.The rest of the paper is organized as follows: Section II surveys the state-of-the-art in UWA channel and network simulation; Section III gives an introduction on UWA propagation; Section IV provides a more detailed description of UWA communication links and how to model them using beam tracing; Section V describes how this UWA link model can be efficiently incorporated into network simulations: 1) as a direct look-up table, 2) via statistical channel modeling; Section VI presents two network simulation case studies; finally, Section VII concludes the paper.

II. STATE-OF-THE-ART IN UNDERWATER ACOUSTIC CHANNEL SIMULATION
A widely used method of channel modeling in the UWA communications research community is by using BELLHOP [17], e.g.see [26]- [32].BELLHOP is a beam/ray tracing model for predicting acoustic pressure fields in the underwater environment [17] [18], which is publicly available as part of the Acoustics Toolbox [33], originally developed by M. Porter and currently maintained by the Woods Hole Oceanographic Institute.Beam/ray tracing is based on ray theory which approximates the propagation of acoustic waves as rays travelling along particular spatial paths from the source to the receiver [34].The difference between a beam and a ray is that the former adds an intensity profile (e.g.Gaussian) normal to the ray trajectory, thus allowing more accurate calculations of the total acoustic intensity at a given point in space [35]- [37].The beam tracing approach is considered an accurate approximation of acoustic wave propagation in cases where the curvature of the ray trajectory and the change in the acoustic pressure amplitude within a single wavelength are negligible [15].A more appropriate way of calculating the acoustic intensity at low frequencies is by solving the wave equation using normal mode theory [15].In these cases the KRAKEN simulation program [38] can be used instead of BELLHOP.However, in most cases considered in UAN research the carrier frequencies are significantly higher than 1.5 kHz, i.e. the wavelengths are shorter than 1 m (given 1500 m/s propagation speed), which comfortably satisfies the high frequency criterion of the beam tracing approach.
A common approach to channel modeling in simulationbased UAN research is to use the outputs of BELLHOP beam tracing to synthesize realistic impulse responses of UWA multipath channels, and calculate characteristics of received signals, e.g.signal amplitude and delay, using these simulated channel realizations.For example, Yildiz et al. [26] propose a framework for jointly optimizing the packet size and transmit power in UANs and use BELLHOP to simulate a PHY layer UAN simulation/emulation suite based on ns2-MIRACLE • Designed to facilitate easy transition between simulations and at-sea testing (more reliably than DESERT [41]) • Includes an interface with WOSS for channel modeling (same as DESERT) • More complex than DESERT (for the transition from simulation to at-sea testing) A more advanced and accurate method of modeling the UWA channel is to simulate a full virtual PHY layer transmission using a time-varying channel impulse response, e.g.via the Virtual Timeseries EXperiment (VirTEX) program [19] [45].It takes into account the received signal distortion due to the Doppler effect, that is not captured by simulating a single BELLHOP channel realization.Instead, VirTEX performs a series of BELLHOP beam tracing evaluations taking into account the motion of the source, the receiver and the sea surface during a signal transmission.Furthermore, the original VirTEX was modified to include new platform and sea surface motion algorithms that significantly reduce the computation time [34].Similarly to VirTEX, the Waymark model [39] [46] simulates a virtual underwater acoustic transmission assuming a specified trajectory of the relative source-receiver motion.However, while VirTEX is based on BELLHOP, the Waymark model can incorporate any propagation modeling tool (including normal mode models) that produces a channel frequency/impulse response given a set of environmental parameters.Although such simulators provide a much more detailed insight into the behaviour of the channel, they are much more suitable for single point-to-point link PHY layer research.In most cases it is not computationally feasible to simulate a full virtual signal transmission for an entire network consisting of many pointto-point links.
There are multiple open-source simulation suites that have been developed specifically for underwater network simulation.For example, the World Ocean Simulation System (WOSS) [20] [47] is one of the earlier and most well-known UAN simulation platforms.It binds BELLHOP beam tracing outputs to the physical layer of C++ based network simulation platforms, e.g.ns2-MIRACLE [21] or ns-3 [22].WOSS provides a highly integrated solution for UAN modeling, where the user can specify the time of the year and geographic locations of the nodes, and the simulator automatically queries the relevant databases, fetches the corresponding sea bottom characteristics and the SSPs, uses them as environment parameters for BELLHOP beam tracing, and integrates the BELLHOP outputs into the network simulation.Similarly to WOSS, the Aqua-Sim simulator [40] [48] combines the ns-2 network simulation suite with a UWA channel model to produce an integrated UAN simulation tool.However, the channel model used in Aqua-Sim is based on a simple analytical signal attenuation model [16] and a constant 1500 m/s propagation speed.Therefore, Aqua-Sim captures the characteristics of the UWA channel in less detail compared with WOSS.
There are also multiple UAN simulation suites that focus on providing a seamless transition between testing the network performance in simulation and testing the developed protocols at sea using real hardware.Two notable examples of such simulation suites are DESERT [23] [49] and SUNSET [24].Both of these simulators are based on the ns2-MIRACLE network simulation platform and both have been verified to successfully facilitate the transition from simulation to at-sea testing using real acoustic modems [47] [24]; however, an investigation by Petroccia and Spaccini [41] showed that SUNSET provides a more mature and efficient solution for transitioning from simulation to real-time at-sea implementation.To incorporate a realistic UWA propagation model in the simulation mode, both DESERT and SUNSET include an interface to WOSS, thus allowing them to simulate BELLHOP-based multipath channels.UnetStack [42] is another increasingly popular simulation platform that was developed to streamline the process of UAN protocol development and testing, similarly to DESERT and SUNSET, by enabling the users to port their simulation code onto UnetStack-compatible acoustic modems [1], e.g. the Subnero modems [50].It includes the in-built options to simulate a simple range-based channel model, or a basic acoustic channel model consisting of the commonly used analytical transmission loss model [16] and a Rayleigh or Rician fading model in [42].However, it is also possible to integrate a custom channel model into the UnetStack simulations, e.g. by specifying the per-link detection and decoding probabilities.The two in-built examples of this UnetStack functionality include the channel models based on the real measurements from the MISSION 2012 [51] and MISSION 2013 [52] experiments.
Table I reviews the capabilities, advantages and disadvantages of the UWA channel and network simulation tools discussed in this section.

III. THE UNDERWATER ACOUSTIC CHANNEL
This section gives a brief introduction of key characteristics of the UWA channel related to UAN protocol design.To summarize, in this paper we look at the following UWA channel features: • slow propagation of acoustic waves (typically in the range of 1450-1550 m/s), • multipath scattering due to reflections off the sea surface and bottom, • long channel delay spread due to slow propagation and the refraction and reflection of acoustic waves, • signal attenuation due to spreading and absorption.There are other significant challenges stemming from slow propagation of acoustic signals investigated by the PHY layer and signal processing researchers, such as rapid channel variability and Doppler distortion [53]- [55].However, in this paper we focus on the basics of UWA propagation necessary for network simulations, assuming appropriate acoustic modem design that is able to deal with the PHY layer.

A. Sound Speed
The dominant physical property affecting the performance of UAN protocols is the low sound propagation speed.In contrast with terrestrial radio networks with a propagation speed of 3×10 8 m/s, acoustic waves propagate through water at approximately 1500 m/s, i.e. slower by a factor of 2×10 5 .For example, if an acoustic link length is 1.5 km, it will take roughly 1 second for the signal to propagate from transmitter to the receiver.Furthermore, the sound speed depends on the temperature, pressure and salinity of the water and is, therefore, variable in space and time [56].Fig. 1 shows an example of a depth-dependent sound speed profile (SSP) derived by Dushaw [57] from the 2009 World Ocean Atlas temperature, pressure and salinity data in summer at (56.5 o N, 11.5 o W), i.e. in the North Atlantic Ocean off the coast of the UK and Ireland.
The depth-dependent SSP causes refraction of the acoustic waves, which in turn results in curved wave propagation trajectories as shown in Fig. 2.These plots were obtained using the BELLHOP ray tracing program [33] based on the SSP data shown in Fig. 1b.
The ray trajectories illustrated in Fig. 2 demonstrate that calculating propagation delays based on a Euclidean distance between two communication nodes, a method often used in UAN research [59] [60], is not necessarily valid, since the signal arriving at the receiver may not travel in a straight line.There also may not be a direct path between two nodes, but only a path reflected off the sea surface or bottom.Using a single value of the propagation speed could also be inaccurate, e.g. a typical 1500 m/s approximation [13] [59] [60], since typical sound speed values can vary between 1450 and 1550 m/s depending on location and time of the year.Furthermore, curved trajectories of the acoustic waves can result in acoustic shadow zones with no coverage, and challenging multipath channel conditions, where several refracted copies of the same signal arrive at the receiver at different times and with different amplitudes, in addition to the echoes reflected off the surface and bottom of the sea.

B. Multipath Propagation
Fig. 3a shows a ray trace from the same scenario as for Fig. 2, but with a specific receiver location at 5 km range and 250 m depth.There are three distinct paths (marked 1-3) between the source and receiver due to refraction caused by the SSP from Fig. 1b.The acoustic waves tend to refract towards the lower sound speed region, sometimes forming waveguides at particular depths.Fig. 3a gives an example of a waveguide, where paths 2 and 3 depart upwards until they reach a steep positive sound speed gradient causing them to refract downwards, whereas path 1 starts propagating towards the sea bottom and gradually refracts upwards.However, in addition to the three refracted signal paths in Fig. 3a, there is another possibility (path 4) for the signal to reach the receiver -by reflecting off the sea surface and/or bottom.This would result in several different arrivals of the transmitted signal as shown in Fig. 3b, with the most direct path taking 8 ms less to propagate to the receiver than the other three paths.Such large differences in the arrival times of different multipath components present a challenge for the receiver design, and are in stark contrast with typical terrestrial RF networks, where, for example, only a 5 µs cyclic prefix is sufficient in OFDMbased 4th generation cellular networks to avoid inter-symbol interference (ISI) due to multipath [61].

C. Spreading and Absorption Loss
The attenuation of the acoustic signal power underwater is caused by two phenomena -geometric spreading and absorption, and can be computed as follows [16]: where L(d, f ) = 10 log(P rx /P src ) is the power loss in dB, defined as the ratio between the received power P rx and the The spreading loss is determined as: where k ∈ [1,2] is the exponent that describes the propagation geometry (equivalent to the pathloss exponent in terrestrial RF propagation models [62]).Setting k = 1 describes cylindrical spreading, where water depth is significantly smaller than the horizontal communication range, whereas k = 2 describes spherical spreading and is equivalent to the free space path loss in terrestrial radio systems.Equation (2) divides the distance by a reference distance d ref , thus expressing the spreading loss relative to the signal strength at distance d ref away from the source.The unit commonly used to describe acoustic "signal power" is dB relative to 1 µPa r.m.s.pressure 1 m away from the source (dB re 1 µPa @ ( For frequencies above a few hundred Hz, the absorption loss is often computed using Thorp's empirical formula derived from ocean measurement data [16] [63] [64]: (4) Although the Thorp formula is the most widely used model for calculating the absorption loss, other empirical models have been proposed in the literature [65], including the Francois-Garrison model [66] that was validated by field measurements in many locations across the globe, e.g.North Pacific Ocean, Atlantic Ocean, Mediterranean Sea.
Fig. 4 shows a contour plot of the total transmission loss (comprising geometric spreading and Thorp absorption) obtained via BELLHOP simulations of a source at 200 m depth transmitting at 24 kHz.It highlights the waveguide formed around the 100 m depth caused by the SSP in Fig. 1b.

IV. MODELING AN UNDERWATER ACOUSTIC LINK
This section describes the UWA multipath propagation environment in more detail and how an acoustic communication link can be modeled using beam tracing.Appendix A describes how the simulation code linked with this tutorial [25] can be used to follow and replicate the BELLHOP beam tracing results discussed in this section.

A. Sea Surface
The sea surface has a significant impact on the multipath structure of the UWA channel [67] [68].The acoustic waves are reflected off the sea surface with negligible loss and a 180 o phase shift, thus potentially resulting in destructive multipath interference.Therefore, it is important to include a realistic shape of the sea surface in the beam tracing simulations to obtain more realistic scattering patterns of reflected signal paths, compared with the perfectly flat sea surface shown in the plots in the previous section.Fig. 5 gives an example of a rough sea surface synthesized using the Pierson-Moskowitz spectral model for fully developed wind seas [69] [70] (depicted in Fig. 6), with the power spectral density (PSD) given by: where α = 0.0081 and β = 0.74 are empirically derived, g = 9.82 m/s 2 is the acceleration of gravity, U is the wind speed in m/s at 19.5 m height above the sea surface (the only parameter of the Pierson-Moskowitz spectrum), and k = 2π/λ is the angular spatial frequency in rad/m (spatial wave frequency multiplied by 2π, λ -wavelength in m).
To obtain a random surface wave realization such as those shown in Fig. 5, we follow the spectral method described by Mobley et al. [71] and shown in Fig. 6.Taking an Inverse Fast Fourier Transform (IFFT) of the resulting spectrum yields a surface wave in the spatial domain seen in Fig. 5.The spectra in Fig. 6 show that at higher wind speeds, the spectrum extends further into the lower frequencies with higher PSD, which results in longer wave components with greater peakto-peak wave elevation as depicted in Fig. 5. Appendix A-B describes how such sea surface realizations can be generated and incorporated into BELLHOP ray tracing simulations.
The key point of generating such surface wave realizations is the simulated "roughness" of the sea surface that would appropriately scatter the reflected UWA waves.For example, Fig. 7 shows the results of a ray tracing simulation equivalent to that in Fig. 3 but with the random surface waves at 10 m/s wind speed shown in Fig. 5. Fig. 7 shows the change in the multipath structure of the channel caused by the rough sea surface.The three refracted signal paths remain identical, whereas there are now two sea surface reflections at the receiver which are different from that shown in Fig. 3 because these rays are reflected off the sea surface at different, random angles.

B. Bathymetry
Modeling the bathymetry, i.e. characteristics of the sea bottom, plays a similar role in providing a more realistic multipath scattering pattern, as the surface waves discussed in the previous subsection.The interaction of acoustic waves with the sediment at the bottom of the ocean is highly complex and is a standalone topic of many research projects and publications, e.g.[72] [73].In general terms, acoustic waves partly penetrate the sediment layer, which introduces an attenuation and phase shift varying with the angle of incidence.The degree of absorption and reflection of the acoustic waves by the sediment layer depends on such factors as grain size, porosity, grain density and gas content [72], that are specific to many types of sediment around the world.
The shape of the sea bottom has a direct impact on the angles of incidence and departure of the reflected signal paths and may reduce the coverage in some areas by obstructing the line-of-sight.Global bathymetry data of the ocean is freely available to the public, e.g.via the British Oceanographic Data Centre [74].However, the spatial resolution of the data in [74] is 30 arc-seconds, i.e. approximately 1 km.This is a good source for large scale bathymetry shapes over long distances, but does not give enough granularity to simulate multipath scattering due to a rough sea bottom (small-scale bathymetry variations).It is possible to obtain datasets with significantly more detailed bathymetry, including the physical parameters such as density, reflectivity etc., e.g. the SWellEx-96 experiment data [75].However, such detailed bathymetry data including sub-bottom properties is difficult to find and interpret, and is beyond the scope of more generally applicable underwater acoustic simulations presented in this paper.
A more generic approach to simulate the small-scale roughness of the sea bottom is to assume a sinusoid shaped bathymetry [76] [77].Here, the exact shape of the sea bottom is not as important as the fact that different rays get reflected at different angles due to the variable slope of the sea bottom, which would yield a generally more realistic multipath scattering pattern.
For example, in this paper we synthesize a generic sinusoidal topology of the sea bottom with random elevation of the hills z(x) as follows: x is the horizontal range, z max is the maximum hill elevation, and L hill is the length of a single hill, equal to the distance between two adjacent peaks.We shift the sinusoid by −π/2 to align the base of the first hill with the zero range; this does not have to be the case, but makes it easier to derive R(x).R(x) ∈ (0, 1] is a scaling function that returns a uniform .Rough sea surface and uneven seabed result in a significant increase in the number of ray-traced multipath components at the receiver, compared with the flat seabed used to obtain the results in Fig. 3 and 7. random number at different ranges, but is constant across a single hill length between two adjacent minima, thus scaling the hill elevation randomly between 0 and z max .Fig. 8 shows an example of a bathymetry with z max = 20 m and L hill = 200 m, generated using the code described in Appendix A-B.Fig. 9 shows the results of ray tracing with the addition of such a randomly uneven seabed (generated using the method described in Appendix A-B).Here, we use the default BELLHOP characterisation of a generic sea bottom layer as an acousto-elastic half-space with 1600 m/s sound speed (representative of sand-silt [64]) and 1 g/cm 3 density [17].
Fig. 9 shows that there are significantly more multipath components arriving at the receiver due to the increase in the number of possible paths reflected from the rough sea surface and uneven sea bottom.The number of multipath components in Fig. 9b is a more realistic representation of challenging underwater acoustic channels encountered in practice.However, this is a relatively extreme example in terms of the multipath propagation, chosen by us for illustrative purposes.If the communication range was shorter, the sea was shallower, and/or the nodes were placed near the sea surface or seabed instead of the middle of the water column, the number of ray-traced multipath components and their delay spread would likely be significantly smaller.

C. Gaussian vs Geometric Beams
There are two types of beams that are typically used for modeling the UWA multipath propagation: • Geometric (hat-shaped) -the beams are separated at the point of departure by linear boundaries half-way between neighbouring rays, and only those rays whose "hat-shaped" boundary encloses the receiver location are recorded [17].• Gaussian -the energy of every beam spreads more broadly using a Gaussian intensity profile normal to the ray [35].Geometric beams are a better option for graphical ray tracing (e.g. in BELLHOP) because it restricts the resulting plots to only include the signal paths that arrive in very close vicinity of the receiver, e.g.Fig. 3a, 7a, 9a.However, for more advanced simulations, Gaussian beam spreading is considered a more accurate approach for estimating the total acoustic intensity at the receiver by calculating a superposition of multiple Gaussian beams in the vicinity of the receiver [35]- [37].
For example, an equivalent beam tracing simulation to that in Fig. 9 but with Gaussian instead of geometric beams produced the set of arrivals shown in Fig. 10.A lot more echoes are traced to the receiver due to broader Gaussian beams.The relative amplitude of the multipath components is also different from the previously discussed geometric beam simulations due to the Gaussian spreading of the beam energy.This set of amplitudes, phases and delays will yield a more accurate calculation of the total acoustic intensity, described in Subsection IV-D.

D. Calculating Wideband Received Signal Power
This subsection explains how the total received signal power can be calculated using the channel impulse response data (i.e. the attenuation, phase and delay of each multipath component) generated via beam tracing.The key feature of the modeling approach described in this paper is to enable the calculation of the wideband received power, i.e. across a frequency bandwidth that is not negligible compared with the central frequency, which is often the case in UWA communications.
First, consider the two separate factors of transmission loss discussed in Subsection III-C -geometric spreading and absorption.While absorption loss depends on both the distance and the frequency, the spreading loss is only distancedependent.Therefore, if we use beam tracing (e.g. via BELL-HOP) to compute the spreading loss of every signal path, but calculate absorption loss separately for any specified frequency, we can calculate the overall channel gain G = P Rx /P Tx , i.e. the received power P Rx relative to the transmitted power at the source P Tx , by integrating across a frequency bandwidth [f min , f max ] as follows: where: • G -channel gain (linear scale), • f min and f max -minimum and maximum frequencies in the simulated channel, • N -total number of multipath components, • A spr [n] -spreading loss of the n th path, • θ[n] -a phase shift of the n th path due to reflections, • τ [n] -propagation delay of the n th path, • τ 0 -reference time, e.g.propagation delay of the first received signal path, • A abs (n, f ) -absorption loss of the n th path at frequency f .Fig. 11 compares the received power of narrowband and wideband signals, simulated on a grid of receiver locations spanning 500 m depth and 5 km range, for a source located at 200 m depth, with rough sea surface and seabed introduced in this section.Fig. 11a shows the result with 1 Hz bandwidth that demonstrates the sensitivity of a narrowband signal to multipath interference due to the phase of the multipath components at a given geographical location and frequency.In contrast, Fig. 11b shows the result of the same beam tracing simulation, but post-processed using 7.2 kHz bandwidth (acoustic modem frequency specifications taken from [78]), and as a result significantly smoother due to a decreased sensitivity to the phase of the multipath components.
Appendix A-D gives details of our implementation of the wideband UWA channel model described in this subsection, including the code to replicate the plots in Fig. 11.

E. Calculating Ambient Noise Power
The effect of the noise on UWA communications in realistic environments is an ongoing research topic due to the complex spatially and temporally variable noise environment underwater, e.g.generated by propellers, hydraulic pumps, snapping shrimp etc. [79] [80].In order to provide a generic noise  model, not specific to a particular location in the ocean, we can approximate the common sources of noise using Gaussian statistics and a continuous PSD as described by Stojanovic and Preisig [16] [81].The PSDs of turbulence, shipping, surface wave and thermal noise can be calculated, respectively, using the following empirical formulae [16]: where the PSDs are in dB re 1µPa @ 1m per Hz, s ∈ [0, 1] is the shipping activity factor (0 -low, 1 -high), and w is the wind speed in m/s that causes noise due to the surface waves.
Fig. 12 shows the PSD of the individual noise sources and the total noise PSD between 1 Hz and 1 MHz.The plot shows that particular noise sources are dominant at particular frequencies, e.g. the turbulence noise at very low frequencies, the thermal noise at very high frequencies, the noise due to shipping activity at tens of Hz, and the surface wave noise at 100 Hz -100 kHz.Frequency, Hz Power spectral density of the ambient acoustic noise due to turbulence, shipping, wave and thermal noise sources, shipping activity -0.5, wind speed -10 m/s (reproduced using the model from [16]).

F. Signal-to-Noise Ratio
The ambient noise model described in the previous subsection can be combined with the received power calculations from Subsection IV-D to compute the Signal-to-Noise Ratio (SNR) of wideband signals at specified receiver locations.The SNR is defined as the ratio of the received signal power and the integral of the noise PSD from Fig. 12 across the bandwidth [f min , f max ], as follows: where P Tx G = P Rx is the received signal power on linear scale, and S noise (f ) is the combined noise PSD from Fig.
converted from dB to the linear scale.Fig. 13 shows a plot of the SNR that combines the wideband received power data from Fig. 11b with the intergal of the noise PSD from Fig. 12 in the 20.4-27.6 kHz frequency band.This SNR pattern is useful for estimating the overall coverage area and potential coverage holes for a particular source depth, in this case 200 m.For example, if we assume that the receiver requires a minimum SNR of 0 dB to decode the signal, the approximate communication range of the signal transmitted at 170 dB re 1µPa @ 1m source level is 3.5 km (not accounting for internal receiver noise characteristics), with propagationdependent shadowing patterns extending the range at some depths and reducing it at others.

V. CHANNEL MODELLING FOR NETWORK SIMULATION
In this section we propose a computationally efficient method of incorporating the UWA link model described in the previous section into network simulations with potentially hundreds or thousands of links that must be modeled (the maximum number of links is N (N − 1)/2, i.e. proportional to N 2 , where N is the number of nodes).The key idea of our approach is to separate the channel simulation from the network simulation as depicted in the block diagram in Fig. 14.In this way, the channel data is generated separately via an extensive series of beam tracing simulations, but is then used in the network simulations via the pre-generated look-up table (e.g.saved as a CSV file) at a negligible computational cost.In particular, we propose the following channel metrics, particularly relevant to the network protocol design, to be saved in a look-up table for every link in the network: • Channel gain -overall channel gain for a wideband signal between the source and receiver i.e. the difference in delays between the first and last significant multipath arrivals.In our model, we consider the strongest multipath components constituting 95% of the total received energy, i.e. ignoring the negligible signal paths with longer propagation delays such as those seen in Fig. 10.

Channel simulator (BELLHOP)
Network simulator (e.g.UnetStack, ns-  This approach can also be extended to include other relevant channel metrics as additional columns in the link look-up tables.Alternatively, the raw channel impulse response data can be stored for later processing as described in Appendix B.
The key additional parameter required for SNR calculations within the network simulator is the ambient noise power for the given frequency bandwidth, that can be calculated using the noise model described in Subsection IV-E.

A. Generating Link Look-Up Tables
Fig. 15 shows the flowchart that describes the key steps of our proposed method of generating a UWA channel look-up table, taking an arbitrary 3D network topology as the input and simulating the channel between every pair of nodes using BELLHOP beam tracing, as described in Section IV.It iterates through all node positions, selects one node as the source, and maps all other node positions onto the 2D range-depth plane by calculating their horizontal ranges relative to the source node, as depicted in Fig. 16.
Conceptually, a limitation of our 3D-2D mapping approach is the inconsistency in the surface wave and bathymetry shapes, that are randomly generated for every node acting as the source.An internally consistent alternative to this approach would be to generate a 3D sea surface and bottom and perform 3D BELLHOP ray tracing directly.Another alternative is to simulate the link between every pair of nodes separately via 2D BELLHOP using a vertical cross-section of the 3D environment connecting the two nodes; however, this would increase the number of required ray tracing runs, and therefore the computation time, by a factor of N , where N is the number of nodes (receivers) in the network.These approaches would dramatically extend the simulation time compared with  our proposed approach without necessarily providing benefits for the evaluation of communication protocols.In reality the UWA channel characteristics between every pair of nodes will vary in time with every transmission, mainly due to the small scale motion of the nodes and the sea surface.Therefore, we argue that simulating a specific "frozen" 3D surface wave and bathymetry pattern would not provide a more valid channel model than generating these patterns randomly, unless the given simulation study is specifically focused on dealing with obstructed paths due to underwater objects in 3D space.The proposed approach, in addition to being significantly more computationally efficient, preserves the consistency in direct signal paths but introduces stochasticity into the reflected/scattered multipath components, thus resembling the behaviour of a real UWA channel.
Our MATLAB implementation of the proposed channel look-up table generation approach is described in Appendix B.

B. Statistical channel modeling
An alternative approach to modeling the stochastic nature of the UWA communication channel is to simulate the link between a given pair of nodes many times and build a statistical model of its behaviour.Despite previous efforts in statistical characterization of UWA channels validated using real measurement data [82]- [84], the UWA research community still lacks widely accepted statistical models used for simulations.This is partly due to the large variability of the UWA environments that have been observed to follow different probability distributions, e.g.Rician, Rayleigh, lognormal, K-distribution, ranging from highly time-variant to almost static channels [81].In this subsection we offer a tool for statistical modeling of UWA channels based on BELLHOP simulations that takes into account the key factors affecting the time variability of the UWA channel: random small scale motion of the source, the receiver and the sea surface.
Instead of simply simulating a link between every transmitter and receiver in the network, we select one pair of the source and receiver location, generate matrices with random small scale perturbations in their coordinates (e.g.within several metres of their average location) and simulate the channel for every combination of the randomly perturbed source-receiver locations.For example, in this subsection the locations of both the source and the receiver are randomly varied (50 times each) within a 10 m radius sphere (uniform random azimuth, elevation and radius).This enables the statistical representation of the UWA channel between two quasi-static nodes based on 2500 realizations (50×50 combinations of the source and receiver displacements).As an alternative to the random node displacement model described above, a database of node displacement measurements from real UAN experiments in [85] can also be used as the basis for creating a statistical UWA channel model in the same way.Fig. 17 shows the results generated by this script for two nodes spaced 4 km apart (horizontally), 500 m sea depth, rough sea surface and bottom, SSP from Fig. 1b.We performed three separate sets of simulations: 1) The source is near the sea surface (20 m depth) and the receiver is near the seabed (480 m depth); 2) Both the source and the receiver are in the middle of the water column (250 m depth); 3) Both the source and the receiver are near the seabed.Firstly, Fig. 17a shows a considerable statistical spread of channel gain values caused by variable multipath scattering.Secondly, it shows that the channel gains are visibly different for the three scenarios considered despite roughly the same distance between the source and the receiver.For example, the crucial factor negatively affecting the performance of the seabed to seabed acoustic links is the upward refraction of acoustic waves caused by the sound speed gradient (Fig. 1b), thus steering the direct signal paths away from the receiver, often resulting in the sea surface reflections being the only received signal paths.
The stepped shape of the cumulative distribution function (CDF) of the propagation delay in Fig. 17b for the seabed scenario reveals that in approximately 75% of the cases the first received signal path is reflected off the sea surface.In contrast, Fig. 17a shows that the presence of at least one direct signal path between two mid-column nodes increases the average channel gain and reduces its variability, compared with the seabed scenario.Another interesting observation from Fig. 17a is that the propagation between a node near the sea surface and a node near the sea bottom is better than that between two mid-column nodes, despite the slightly increased propagation distance.Due to the proximity of a node to the sea surface, a lot of the reflected acoustic energy travels a very similar distance as the direct signal path, thus forming additional "quasi-direct" signal paths and statistically increasing the received signal strength (resembling cylindrical spreading).Whereas the case where both nodes are located in the middle of the water column is closer to spherical spreading with no reflective surface in close proximity of the nodes.
Fig. 17c gives a valuable insight into a typical multipath delay spread in a UWA channel.It shows that a UAN protocol designer for this scenario should accommodate at least a 200 ms delay spread (to cover most links), for example, by separating the scheduled packet reception times by a 200 ms guard interval.These important features of the UWA channel behaviour are typically not captured by the simplified analyt- Channel gain, linear Example of statistical channel modeling -fitting lognormal distributions to the empirical linear channel gain data.Horizontal range -4km; water depth -500 m; SSP from Fig. 1b.ical propagation models.For example, the widely used Urick model [16] based on the Euclidean distance between the nodes would not incorporate any of the channel gain variability, direct path refraction effects or multipath delay spread observed in Fig. 17.
One way of using this empirical data for statistical channel modeling is by randomly selecting one of these channel realizations for every transmission between the corresponding two nodes.In this example, drawing a uniform random integer between 1 and 2500 to pick a channel gain, delay and multipath spread from the look-up table would, in the long run, result in the same statistical channel properties as those shown in Fig. 17.However, simulating many channel realizations for every pair of nodes in the network may not be computationally feasible, especially if the network size is in the order of hundreds or possibly thousands of links.
A more flexible and widely applicable way of using such empirical data is to characterize the observed stochastic channel behaviour using analytical probability distributions.Fig. 18 gives an example of such statistical channel modeling for the three scenarios from Fig. 17.Here, the histograms of the linear channel gain data can be approximated by lognormal probability density functions, where µ is the mean and σ is the standard deviation.The challenge in this approach is to identify mathematical relationships between the environment parameters, e.g.source and receiver depths, range, frequency band, SSP etc., and the probability distribution type and its parameters, in order to generalize these models and eliminate the need to simulate every link many times using BELLHOP.We do not propose a statistical channel model, as this would require an extensive study and as such is beyond the scope of this paper, but rather offer a tool for researchers to generate their own statistical models tailored to the UWA environment parameters most relevant to them, e.g.deep/shallow water, short/long range etc.

VI. NETWORK SIMULATOR CASE STUDIES
In this section we present two case studies of integrating the proposed beam tracing based channel modeling approach into network simulators.The first case study investigates the effects of custom beam tracing channel data on the Riverbed Modeler [86] simulations of the ALOHA protocol [87] in a singlehop UAN.The second case study applies statistical channel modeling described in Subsection V-B to investigate its effects on custom MATLAB simulations of Spatial TDMA (STDMA) [88] in a linear UAN scenario.

A. Riverbed Modeler Case Study
Riverbed Modeler (formerly known as OPNET) is a discrete-event packet-level network simulation platform.It provides a customizable broadcast medium to model wireless communications via the Radio Transceiver Pipeline (RTP), which allows the user to model every transmitter-receiver link in the network.This pipeline consists of fourteen stages executed on a per-receiver basis whenever a packet is transmitted.These stages, shown in Fig. 19, use a number of Transmission Data Attributes (TDAs) offering a standard set of values to support the implementation of communication links.Each stage is defined by a module written in C and saved in an editable file with the extension ".ps.c".

UWA Channels in Riverbed Modeler
By modifying a number of these pipeline stages, Riverbed Modeler can be used to simulate other types of wireless channels including UWA channels.To this end, at least three stages highlighted as shaded blocks in Fig. 19 must be customized: the propagation delay, the received power, and the background noise.In this section we compare three different methods of modeling UWA links using the Riverbed Modeler's RTP: • Basic binary collision model -the link connectivity is defined by a fixed connection range and any temporal overlap in received packets results in a collision and loss of both packets.The propagation delay is calculated using the Euclidean distance between the nodes and a fixed 1500 m/s propagation speed.This model does not consider the noise or the received signal power (abstracted by the distance based connection range).• Urick propagation model -the received power is calculated using the geometric spreading loss and the Thorp absorption formula described in Subsection III-C.This enables the calculation of the Signal-to-Interference-plus-Noise Ratio (SINR), the resulting bit error rate (BER), and the probability of packet error computed by the Riverbed Modeler.The propagation delay is calculated using the Euclidean distance between the nodes and a fixed 1500 m/s propagation speed.• BELLHOP-based channel model -The channel gain and propagation delay values are precomputed using BELLHOP beam tracing, and are directly imported as a look-up table.The received power is then calculated using the imported channel gains and used by the Riverbed Modeler to compute the probability of packet error based on the SINR in the same way as when using the Urick propagation model.The key customization steps of the Riverbed Modeler's RTP are described in more detail below.
a) Propagation Delay: In this stage, a default propagation model, called dra_propdel, is used to compute the propagation delay of each transmitted packet (i.e. each link) based on a pre-defined propagation speed and transmission distance.For an acoustic-link scenario, this pipeline model can be used to set the desired speed of sound in each link.In the default propagation model dra_propdel.ps.c, the speed of sound can be set as a fixed value to provide a single value to all links (the usual assumed speed of UWA propagation is 1500 m/s).Alternatively, the propagation speed on every link can be set using the delay look-up table produced by our BELLHOP channel model via the dra_propdel.ps.c source file.
b) Rx Power: The default model for this stage is called dra.power, which takes into account the transmitted power, path loss and Rx/Tx antenna gains to compute the received power for every link.For the UWA link scenario, this pipeline model can be used to calculate the received power by inserting an empirical model (e.g. the Urick model) into the dra.power.ps.c file to calculate the propagation loss and estimate the received power.Another approach is to import the channel gain values produced by our BELLHOP-based model into the dra.power model in order to compute the received power.c) Ambient Noise: This task is defined by a default model called dra.bkgnoise, which is a simple procedure taking into account fixed ambient noise and thermal noise power levels.For UWA simulations, the empirical formulae for the ambient noise due to turbulence, shipping, surface waves and thermal noise described in Subsection IV-E can be inserted into the dra.bkgnoise.ps.c file.This noise model is used in both the Urick propagation model and our proposed BELLHOP-based channel model.d) Packet Receiving Stages: in the SNR and BER stages, Riverbed Modeler works out the SNR and BER values respectively.Following this, the probability of bit errors for each packet segment is obtained.This is based on the SINR value and a built-in look-up table for a given modulation scheme (e.g.BPSK).Next, in the Error Allocation stage, the number of bit errors in a packet is calculated.Then, it is determined whether the arriving packet can be accepted at the destination node.The acceptability test of a packet at the receiver can be customised to reflect one of the three different methods of modeling UWA links listed above in this subsection.For the Urick and BELLHOP-based models, the acceptability test is based on the comparison between the instantaneous SINR value of the arriving packet and a predefined SINR threshold.The instantaneous SINR value is calculated based on the outcome of the Rx Power, Noise and Interference stages.For the basic binary collision model, the Error Correction stage is adjusted to reject all packets involved in an overlap, if a non-zero-length overlap between successive arriving packets is detected in the Interference stage.
e) Inactive Stages: some stages, with dashed line boundaries in Fig. 19, are specific to the internal Riverbed Modeler simulation setup.They are concerned with creating an ini-tial potential receiver group for each transmitter, computing Rx/Tx antenna gains and determining the closure between the transmitter and the receiver (i.e. the ability of physically establishing a link with regard to the intersections of this link with the earths surface).These stages are outside the scope of this case study and have no effect on UWA link modeling, but they must be executed on a per-receiver basis.

Simulation Setup
The effects of using different channel models in Riverbed Modeler simulations, as described above, are investigated in this case study using a single-hop UAN network depicted in Fig. 20, where a number of underwater sensor nodes communicate with a single surface node.In our simulations 50 sensor nodes were placed randomly in a 6×6 km coverage area at uniformly distributed random depths between 20 and 480 m, with the surface node located at the centre of the coverage area at 10 m depth.The simulation parameters are summarized in Table II.

Simulation Results
In Fig. 21 the network performance is evaluated in terms of the overall network throughput and the packet loss recorded at every individual node.Firstly, the results show that the channel model has a visible effect on the simulated network performance.For example, the well-known ALOHA throughput curve in Fig. 21a peaking at 50% offered traffic and 18% throughput, obtained via the simplistic binary collision model, is in fact a pessimistic estimate compared with more detailed channel models which consider the received signal and interference power.This is because some nodes, typically located closer to the receiver, are able to transmit their packets successfully despite interference from more distant nodes, due to a high SINR, whereas the binary collision model discards all packets involved in a collision.The difference in the packet loss between the binary collision model and the two SINRbased models is shown in more detail in Fig. 21b.
A comparison between the analytical Urick propagation model and the BELLHOP-based channel model reveals that the former is more optimistic in terms of the received power values, thus resulting in less packet loss and slightly higher network throughput.Furthermore, the packet loss distribution in Fig. 21b shows that one of the 50 nodes experienced complete outage due to low signal strength under the BELLHOP-based channel model, thus highlighting that an important network topology feature is not captured by the more basic channel models.For example, if the link quality is too poor for one or more nodes to communicate with the surface node in the simulated environment, this should inform the network topology and protocol design, e.g.multi-hop connectivity should be considered.

B. Statistical Channel Modelling Case Study
In the second case study, we evaluate the effects of statistical channel modeling on the performance of Spatial TDMA (STDMA) applied to a UAN with a line topology, representative of a subsea asset monitoring scenario shown in Fig. 22.Here, the network consists of multiple underwater sensor nodes arranged in a line such that every node has two connections -a node one hop closer to the sink node (up the chain) and a node one hop further down the chain.The job of a sensor node is to transmit its own packets up the chain and forward data packets from the nodes down the chain.
The inherent sparsity of linear network topologies is wellsuited for STDMA, since it can be exploited by assigning TDMA slots to several spatially separated transmissions simultaneously without collision [88]- [90], thus reducing the number of slots in the TDMA frame.In fact, Chitre et al. [91] show that it is theoretically possible to design packet schedules for networks with long propagation delays (UANs are a typical example of this) that exceed the throughput of networks

Simulation Setup
In this case study we simulate the scenario from [92], where 10 sensor nodes and one sink node are arranged in a linear topology near the seabed with 1 km spacing between the adjacent nodes.The simulation parameters are summarized in Table III.
The statistical channel modeling approach from Subsection V-B with random 10 m radius node displacement is used to generate 2500 UWA channel realizations for every possible hop distance, i.e. from 1-hop node separation (adjacent nodes) up to 10-hop separation, resulting in 10 separate channel lookup tables.A full network model is then generated by selecting a random channel realization (channel gain, delay and delay spread) for every link in the network from a corresponding look-up table.
After a channel realization is assigned to every link in the network, a binary N × N interference matrix I with the elements defined as: where I[i, j] indicates if there is a link between nodes i and j based on the 0 dB SNR threshold, i.e. if a signal from node i is received at node j with ≥0 dB SNR, they are considered interfering nodes.G[i, j] is the channel gain between nodes i and j in dB; and P tx and P n are the source power and the ambient noise power in dB re 1 µPa @ 1m, respectively.
The STDMA schedule is derived by computing an N sn × N slots matrix, where N sn = 10 is the number of transmitting sensor nodes and N slots is the number of time slots, indicating which node transmits in which time slot, such that N slots is minimized subject to no collisions according to I. In this way, the interference matrix I dictates the efficiency of the spatial reuse pattern and the STDMA frame length achievable in a given network realization.Furthermore, in the classical contention-free TDMA the slot duration must incorporate the propagation delay and delay spread into the guard interval in order to avoid inter-slot interference.The TDMA slot duration τ slot for a given network realization is determined as follows: where τ dp is the packet duration, T p [i, j] is the propagation delay between nodes i and j, and T spr [i, j] is the multipath delay spread on the link between nodes i and j.
The frame length N slots and the slot duration τ slot can be used to compute the network throughput in packets per second under full buffer traffic conditions as follows: where N packets is the total number of packets transmitted within a single frame.In the scenario considered in this case study, where 10 sensor nodes transmit packets to the sink node in a line topology, the number of packets per frame is: since every node transmits its own packet and forwards the packets from all other nodes down the chain.

Simulation Results
Fig. 23a shows the statistical distribution of the slot duration calculated using (14) in 10000 network realizations for every simulated source level.The slot duration is the shortest at 160 dB re 1 µPa @ 1m source level because in the vast majority of cases the maximum interference range is limited to 2 hops, thus eliminating the need to extend the guard interval to accommodate propagation delays to the nodes further away.However, this source level results in approximately 10-11 dB SNR of the intended transmissions which may not leave a sufficient margin for reliable communication if the ambient noise increases or the channel experiences increased fading.However, Fig. 23a shows that increasing the source level by 5 or 10 dB in most cases extends the maximum interference range (and with it the slot duration) to 3 hops or even 4 hops, thus providing a trade-off between the idle guard time added to the slots and the reliability of transmissions.Fig. 23b shows that statistical variations in the channel gain have a direct impact on the frame length of the STDMA protocol due to the differences in the spatial reuse patterns governed by the interference matrix in a given network realization.For example, at 165 dB re 1 µPa @ 1m source level the STDMA frame length varies between 27 and 34 slots only due to the channel gain variability among different network realizations, demonstrating the effect of statistical channel modeling in this scenario.Furthermore, the significant variability in the STDMA frame length (up to a factor of two) observed across all simulations at the three source levels in Fig. 23b gives researchers a valuable insight for MAC protocol design, that typically would not be captured by simplified interference models often used in the literature.
Finally, Fig. 23c quantifies the variability in the expected STDMA network throughput.For example, it reveals that the network throughput is superior at a lower source level (assuming no packet loss) due to the combined effect of shorter slot duration and more efficient spatial reuse patterns.However, the main conclusion from Fig. 23c is to reiterate the considerable effect of statistical channel variations on the network performance, which should be taken into account when designing UAN protocols to be deployed in real-world environments.

VII. CONCLUSION
This paper has presented a detailed tutorial on modeling multipath UWA channels, primarily aimed at UAN protocol researchers.The tutorial was particularly focused on modeling the channel gain, propagation delay and multipath delay spread, as the key parameters affecting the performance of network protocols.We described two methods of incorporating the beam tracing channel data into network simulations, including a case study for each of them: 1) directly importing the data as a look-up table, 2) using the data to create a statistical channel model.The Riverbed Modeler case study revealed that a simple binary collision model provided a pessimistic estimate of the packet loss and throughput performance of ALOHA, compared with the BELLHOP-based channel model.In contrast, a widely used analytical UWA propagation model provided an optimistic estimate of the network performance by omitting the multipath structure of the UWA channel captured by our proposed model.The second case study showed that the slot duration, frame length and network throughput of STDMA can be greatly affected by the variability captured by a statistical channel model, demonstrating the importance of considering such statistical channel variability when designing UAN protocols to be deployed in real-world environments.tic propagation physics via beam/ray tracing.However, for most researchers with network protocol background it requires learning the basics of ocean acoustics and a substantial programming effort before they can start simulating underwater acoustic networks.This appendix describes how a UWA communication link can be modeled using BELLHOP with our provided codebase [25], replicating the plots discussed in Section IV.
Fig. 24 shows how BELLHOP operates in terms of reading the user input and producing a beam/ray tracing output.It reads a plain text file with a .envextension which follows a pre-defined format that specifies the environment to be simulated.There are also several optional input files a user can create to customize surface waves (.ati), bathymetry (.bty), range-dependent SSP (.ssp), top/bottom reflection coefficients (.trc, .brc)and source directivity (.sbp).Examples of these environment files generated during the BELLHOP experiments in this section can be found in the data directory of the provided codebase; the reader is also encouraged to study the BELLHOP manual [17] for other simple examples.
Depending on the type of simulation specified by the user in the .envfile, BELLHOP produces one of the following types of output files: • .ray-coordinates of the rays for a graphical ray tracing output (e.g.Fig. 2, 3a), • .shd-transmission loss data for a specified 2D rangedepth area (e.g.Fig. 4), • .arr-attenuation, phase and delay of every signal path traced to all specified receiver locations (e.g.Fig. 3b).

A. Ray Tracing Using BELLHOP
As the first simple example, the ray trace plot in Fig. 2 can be produced by running the code in Listing 1.The key variables and functions there are the following: • pars -structure of the environment and simulation parameters, i.e. a structure containing the information to be written into the .envfile.The comprehensive list of the fields of this structure is described in default_sim_pars.the custom altimetry and/or bathymetry files.3) Run bellhop(pars.filename).4) Process the output file generated by BELLHOP.The single_sim script in the provided codebase allows the user to try different types of BELLHOP simulation, other than simple graphical ray tracing.For example, setting the pars.simtypefield to 'eray' simulates a large number of rays and only plots those that arrive near the specified receiver location; it will reproduce the plot in Fig. 3a.It also shows the user how to change other simulation parameters, for example, the source depth and receiver depths and ranges.The amplitude-delay plot of the multipath arrivals in Fig. 3b can be reproduced by setting pars.simtype to 'arr'.Similarly, setting pars.simtype to 'loss' will configure the single_sim script to produce the transmission loss plot in Fig. 4.

B. Simulating the Altimetry and Bathymetry
Random wind-induced sea surface waves, such as those shown in Fig. 5, can be incorporated into the BELLHOP simulations using our interface by including the code in Listing 2 before executing the bellhop function.First, the pars.use_altimetryflag must be set to true to instruct BELLHOP to use a custom altimetry.Next, the spatial wave sampling frequency and the wind speed for % Set up altimetry parameters in the pars structure pars.use_altimetry= true; % use custom altimetry pars.wave_resolution= 5; % sampling at 5m pars.wind_speed= 10; % 10 m/s wind % Create .ATI file with random sea surface create_sea_surface_file(pars); Listing 2. Matlab code which generates random surface waves at the specified wind speed as shown in Fig. 5.The filename by default is same as the env file, and the length of the area is automatically set to maximum range.the Pierson-Moskowitz PSD need to be specified.While choosing the spatial wave sampling frequency (referred to as wave resolution), a trade-off between the level of detail and simulation speed needs to be determined; better wave resolution will include higher frequency components in the wave realization but will cause BELLHOP to run more slowly due to an increased number of altimetry sampling points.The create_sea_surface_file function then creates a .atifile that follows a format defined in BELL-HOP, which specifies the depth of the sea surface at fixed pars.wave_resolution range increments.
Furthermore, the reader can reproduce the wave spectra and the sea surface realizations from Fig 5 and 6, or experiment with other wind speed and wave resolution values, using the wave_modelling script.
The code for generating the bathymetry depicted in Fig. 8 is given in Listing 3, which follows the same pattern as the altimetry code in Listing 2. By default, the common BELLHOP characterisation of a generic sea bottom layer is used -an acousto-elastic half-space with 1600 m/s sound speed (representative of sand-silt [64]) and 1 g/cm 3 density [17].

C. Compressed Set of Multipath Arrivals
A core feature of our approach to building channel models for network simulators is to obtain the data for a set of multipath arrivals via BELLHOP for every pair of nodes, and save them all into a large look-up table (e.g. in the Comma-Separated Values (CSV) file format).We can then import this look-up table into any network simulator to characterize the channel for every link in the network.
An important step in creating such a look-up table for large networks with hundreds or thousands of links is the option to compress the amount of information we store about any individual link without compromising the accuracy of the stored model.For example, there are 78 multipath components in Fig. 10, most of which have near-zero amplitude and a negligible effect on the overall channel properties, but which would make the file size of our look-up table unmanageable.strongest signal paths that constitute 99% or 95% of the total received energy.This reduces the number of multipath components from 78 to 20 and 10, respectively, whilst preserving the vast majority of the important information about the channel.This compression is done using the dedicated compress_arr_set function, in this example as part of the single_sim script.

D. Wideband Channel Gain and SNR
The geometric spreading loss, phase shift and propagation delay values are generated by BELLHOP for every signal path as the first three columns of the output .arrfile, which can be processed by the process_arr_file function in the provided codebase.To calculate the absorption loss across a given frequency bandwidth instead of a single frequency, BELLHOP needs to be instructed not to incorporate Thorp absorption into its calculations by setting the pars.thorpabsorbparameter to false.Then, the absoption loss of the n th signal path at frequency f is computed as: where d km [n] is the length of the n th signal path in km, and L abs (f ) is the Thorp absorption loss in dB given by ( 4) with f specified in kHz.If the lengths of individual signal paths are not known from beam tracing (e.g. by default BELLHOP does not output them), they can be approximated using the average sound speed c[n] and the propagation delay τ [n] of every path as follows: Since the sound speed is variable with depth, e.g. as shown in in Fig. 1b, we approximate c[n] as the mean value of the of using this function.There, the user only needs to specify a set of node positions in 3D Cartesian coordinates (which are used both as the source positions and as the receiver positions) and the name of the output CSV file.This node position format provides the data compatible with 2D BELLHOP simulations using an irregular grid of receivers (configured by setting pars.regulargrid=false),i.e. simulating the UWA propagation for receivers at an arbitrary set of (depth, range) pairs.
The example_3d_channel_gen script provides a more specific example of how to use this function.It creates a channel look-up table for 30 nodes randomly placed within a 6km × 6km × 500m area and a surface node at 10 m depth in the centre of the area.
The create_3d_channel_lut creates a CSV file containing a table with the following columns: 1) Source node index -index of the source node in the array of node coordinates, 2) Receiver node index -index of the receiver node in the array of node coordinates, 3) Channel gain -overall channel gain for a wideband signal between this source and receiver [dB], 4) Channel delay -propagation delay of the first received path [sec], 5) Delay spread -channel delay spread [sec], considering the strongest multipath components constituting 95% of the total received energy.An optional input of create_3d_channel_lut is a binary flag asking the user whether the raw data for every multipath arrival should be saved instead of the processed wideband channel model described above.If this flag is set to true, the column format of the output CSV file is changed to include the amplitude in dB, propagation delay and phase shift of each multipath component, i.e. 20 log(A spr [n]), τ [n] and θ[n] from Equation (7).While the source and receiver index columns are identical in both formats, the number of subsequent columns in a given row is variable depending on the number of multipath components traced for the given sourcereceiver pair.In this way, the channel data is independent of the centre frequency and bandwidth of the signals, e.g. it is more generally applicable and not limited to a particular frequency band specification.However, if necessary, a CSV file with the raw data can be post-processed using our wideband channel gain model described in Subsection IV-D via the process_raw_ch_imp function.
As an example of a more advanced approach to UWA channel modeling, the example_stat_channel_model script shows how the create_3d_channel_lut

Fig. 1 .
Fig. 1.Example of an SSP in the North Atlantic Ocean based on average summer temperature, pressure and salinity data at (56.5 o N, 11.5 o W) [58].

Fig. 2 .
Fig. 2. Underwater acoustic signal propagation with refraction due to variable sound speed from Fig. 1b, and with reflections off the sea surface and bottom; generated using BELLHOP at 200 m source depth.

Fig. 3 .
Fig. 3. Underwater acoustic multipath channel with the sound speed profile from Fig. 1b; source depth -200 m, receiver at 5 km range and 250 m depth.

Fig. 4 .
Fig. 4. Incoherent transmission loss of an acoustic signal due to spreading and absorption at 24 kHz frequency; source depth -200 m

Fig. 7 .
Fig.7.Underwater acoustic multipath channel with randomly generated surface waves at 10 m/s wind speed; two distinct surface-reflected paths are traced to the receiver compared with the flat surface case in Fig.3.

Fig. 8 .Fig. 9
Fig. 8. Example of a sinusoidal bathymetry with 200m long hills and random hill height between 0 and 20 m

Fig. 10 .
Fig.10.Set of multipath arrivals in an underwater acoustic channel equivalent to Fig.9but with broader Gaussian beams, as opposed to geometric beams, resulting in more multipath components traced to the receiver, with more accurately estimated amplitude.
Wideband signal, 24 kHz centre frequency with 7.2 kHz bandwidth.

Fig. 13 .
Fig. 13.Signal-to-Noise Ratio analysis for the source at 200 m depth; source level -170 dB re 1 µPa @ 1m, 24 kHz centre frequency, 7.2 kHz bandwidth.Blank parts of the plot indicate the areas with SNR < -10 dB.

Fig. 14 .
Fig.14.Our proposed simulation framework, where the channel for every combination of source and receiver location is simulated using BELLHOP beam tracing.A network simulator then uses the channel data (e.g. a CSV file) to characterize the link between every pair of nodes.

Fig. 15 .Fig. 16 .
Fig.15.Generating BELLHOP channel data given a set of node positions, by iterating over every node as the source, and mapping the other nodes onto a range-depth set of receivers for 2D BELLHOP simulations Fig.18.Example of statistical channel modeling -fitting lognormal distributions to the empirical linear channel gain data.Horizontal range -4km; water depth -500 m; SSP from Fig.1b.

Fig. 20 .
Fig. 20.Single-hop UAN scenario where the sensor nodes send their data directly to a gateway node located on the sea surface.

Fig. 21 .
Fig. 21.The UWA link model has a visible effect on the throughput and packet loss performance of a single-hop UAN simulated in Riverbed Modeler.

Fig. 25 Fig. 25 .
Fig. 25.Compressing the set of multipath arrivals in an underwater acoustic channel to include only 99% or 95% of total received energy by eliminating weak, negligible echoes; source depth -200 m, receiver at 10 km range and 250 m depth.
function is used to create the statistical channel model presented in Subsection V-B.Likewise, the create_stat_ch_model_linnet script utilizes the create_3d_channel_lut function to generate a statictical channel model for the linear UAN case study in Subsection VI-B.Paul Mitchell (M'00-SM'09) received the MEng and PhD degrees from the University of York in 1999 and 2003, respectively.His PhD research was on medium access control for satellite systems, which was supported by British Telecom.He has been a member of the Department of Electronic Engineering at York since 2002, and is currently Senior Lecturer.He has gained industrial experience at BT and QinetiQ.Current research interests include medium access control and networking for underwater and mobile communication systems, and the application of artificial intelligence techniques to such problems.He is currently a lead investigator on EPSRC USMART (EP/P017975/1) and a co-investigator of H2020 MCSA 5G-AURA.Dr Mitchell is an author of over 100 refereed journal and conference papers and he has served on numerous international conference programme committees.He was General Chair of the International Symposium on Wireless Communications Systems which was held in York in 2010, and a Track Chair for IEEE VTC in 2015.He is an Associate Editor of the IET Wireless Sensor Systems journal and the Sage International Journal of Distributed Sensor Networks.Dr Mitchell is a Senior Member of the IEEE, a member of the IET and a Fellow of the Higher Education Academy.Yuriy V. Zakharov (M'01) received the M.Sc.and Ph.D. degrees in electrical engineering from the Power Engineering Institute, Moscow, Russia, in 1977 and 1983, respectively.From 1977 to 1983, he was with the Special Design Agency in the Moscow Power Engineering Institute.From 1983 to 1999, he was with the N. N. Andreev Acoustics Institute, Moscow.From 1994 to 1999, he was with Nortel as a DSP Group Leader.Since 1999, he has been with the Communications Research Group, University of York, U.K., where he is currently a Reader in the Department of Electronic Engineering.His research interests include signal processing, communications, and acoustics.

TABLE I .
UNDERWATER ACOUSTIC CHANNEL AND NETWORK SIMULATION PLATFORMS Fig. 19.Using a custom UWA channel model in the Riverbed Modeler's Radio Transceiver Pipeline.

TABLE III
BELLHOPFig. 24.BELLHOP simulation setup: environment parameters are read from .ENV, .ATI, .BTY files, results are stored in .RAY/.SHD/.ARR and .PRT files Structure with default simulation parameters 2 pars = default_sim_pars; 3 % Specify the name for files generated by BELLHOP 4 pars.filename= 'example_ray_trace'; 5 % Create the ENV file using the given parameters 6 create_bellhop_env_file(pars); 7 % Run BELLHOP using this ENV file 8 bellhop(pars.filename);9 % Plot the ray trace produced by BELLHOP 10 plotray(pars.filename);Listing 1. Minimal example of Matlab code that uses our proposed interface to run BELLHOP and produce a ray trace plot in Fig. 2. • bellhop -the main BELLHOP function that invokes the FORTRAN executable; the input is a string specifying the name of the .envfile.• plotray -Acoustics Toolbox function that plots the rays saved by BELLHOP in a .rayfile; the input is a string specifying the name of the .rayfile, which by default is the same as that of the .envfile.The sequence of steps in Listing 1 describes in general the setup of any BELLHOP simulation using our interface: 1) Populate the pars structure with environment and simulation parameters.2) Create the BELLHOP input .envfile, and (if needed) 1 %