# Under-Sea Ice Diffusing Optical Communications

Abdallah S. Ghazy, Haitham S. Khallaf, Steve Hranilovic, Mohammad Ali Khalighi

► To cite this version:

Abdallah S. Ghazy, Haitham S. Khallaf, Steve Hranilovic, Mohammad Ali Khalighi. Under-Sea Ice Diffusing Optical Communications. IEEE Access, In press, pp.1-1. 10.1109/ACCESS.2021.3131276. [hal-03464499](https://hal.science/hal-03464499)

HAL Id: hal-03464499

https://hal.science/hal-03464499

Submitted on 3 Dec 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Date of publication xxxx 00, 0000, date of current version xxxx 00, 0000.

Digital Object Identifier 10.1109/ACCESS.2017.DOI

# Under-Sea Ice Diffusing Optical Communications

ABDALLAH S. GHAZY¹, HAITHAM S. KHALLAF², (Senior Member, IEEE) STEVE HRANILOVIC³, (Senior Member, IEEE) and MOHAMMAD-ALI KHALIGHI⁴, (Senior Member, IEEE)

¹,²,³Electrical and Computer Engineering School, McMaster University, Ontario, Canada (e-mail: ghazya@mcmaster.ca, khallafh@mcmaster.ca, hranilovic@mcmaster.ca).

⁴Aix-Marseille University, CNRS, Centrale Marseille, Institut Fresnel, Marseille, France (e-mail: Ali.Khalighi@fresnel.fr).

Corresponding author: Steve Hranilovic (e-mail: hranilovic@mcmaster.ca).

This work has been supported by the Natural Sciences and Engineering Research Council of Canada and the France-Canada Research Fund.”

**ABSTRACT** In this paper, we propose a novel approach to establish a reliable high-speed broadcast communication link between a group of autonomous underwater vehicles (AUVs) swarm under-sea ice. We utilize the fact that sea ice exists above the AUVs to diffuse the optical beam sent from AUV transmitter. We model this channel using a new seawater-sea ice cascaded layers (SSCL) model in which the vertical channel is divided into multiple layers based on their optical characteristics. The diffusing pattern of the SSCL model is computed using Monte Carlo numerical ray-tracing technique. We derive a quasi-analytic equation for the channel impulse response (CIR) which is valid for AUV receivers with different configurations, locations and orientations. The communication performance of underwater sea ice diffusing system is quantified via bit error rate performance, power penalty and maximum achievable bit rate. Our results reveal that, for a snow-covered sea ice sheet with thickness of 36 cm and bare sea ice sheet with thickness 12 cm, the proposed system can achieve a broadcast communication rate of 100 Mbps with ranges up to 3.5 meters and 3 meters, respectively, with BER less than $10^{-3}$ and average transmitted power of 100 mW.

**INDEX TERMS** underwater wireless optical communication, diffusing communication, AUVs, channel impulse response.

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

## I. INTRODUCTION

Sea ice regions are key zones as they play an important role in climate change and ecosystems of the Earth [1]. They cover roughly 7% and 15% of the earth and the sea-waters, respectively. To understand this rapidly changing environment, researchers have been working on measuring campaigns such as mapping thickness of the sea ice sheets and measuring ice characteristics (e.g., the temperature and salinity) [2], [3]. In addition, it is important to observe, monitor, and protect this ecosystem (e.g., detecting and removing oil spills) [4]. Due to their reliability, cost-effectiveness, and ability to improve human safety, autonomous underwater vehicles (AUVs) are commonly used in under sea ice measurements [2]–[4].

A reliable communication link between the AUVs is essential in order for them to work collaboratively to tackle complex tasks, such as the case of cooperative agents in AUV swarms [5]–[7]. For any technologies deployed on AUVs, limitations on the size, weight and power consumption are critical [8]. Acoustic, radio frequency, and optical communi-

cations are the three main wireless communication systems used in underwater purposes. Compared to both acoustic and radio frequency systems, optical wireless communication (OWC) systems achieve higher transmission data-rate, better power efficiency, and smaller size on the order of cubic centimeters [3], [9].

The mobility of the AUVs, the nature of the sea ice terrain, and presence of the marine groups (e.g., bears, seals, pENGUINS) can degrade the performance of line-of-sight (LOS) OWC systems because of high misalignment and blockage probabilities [10], [11]. While, non-line-of-sight (NLOS) links based on omni-directional sources such as light emitted diodes (LEDs) offer relatively higher reliability, they provide relatively lower speed communication due to their limited modulation bandwidth. In indoor environments, broadcast OWC systems have been proposed where the ceiling, walls and floor have been employed as diffuse reflectors of the optical signal [12]. There has been much work on the optimizing of the diffusing pattern, and the system complexity

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38
of such indoor systems to reduce the effects of interference and background noise [13]–[16]. Recently, the application of diffusing communication links to underwater scenarios has started to be considered [9]. Arnon et al. [17] and Liu et al. [18] proposed using seawater-air interface as reflective surface and turbidity seawaters as scattering mediums for NLOS communications, respectively. Anous et al. [19] modeled a vertical underwater link taking into account the in-homogeneous nature of the seawater environment with the depth for both LOS and NLOS scenarios. Anous et al. used the concept of the layering to discretize the vertical variation in the temperature and salinity profiles of the seawaters. This discretized modelling approach of using multi-layers representing the vertical variation in the temperature, salinity and pressures profiles is commonplace in such systems and widely used in the geoscience literature (e.g., [20]–[24]).

In this paper, we propose the concept of *sea ice diffusing optical communications* (SDOCs) where the sea ice is utilized as a diffusing surface with a LD source to establish high-speed short-distance broadcast communication links between the AUVs. Link reliability is improved due to the multiple reflections/scattering from the sea ice and thanks to high impurities contaminating ice mediums and snow caps covering the sea ice sheets. To the best knowledge of the authors, this is the first introduction of this approach in the literature. The main contributions in this paper are summarized as follows:

* For the first time, we introduce a new approach in which the ice sheet is utilized as a diffusing surface to establish reliable diffusing-based broadcast link between underwater AUVs.

* The channel is presented using a seawater-sea ice cascaded layers (SSCL) model where the ice and snow are divided into layers according to the variations in their temperature and salinity profiles.

* In order to obtain transmitter to receiver channel impulse response (CIR), we propose a new simulation methodology consisting of two-steps. In the first step, a Monte Carlo numerical ray tracing (MCNRT) method is used to numerically obtain the ice sheet diffusing pattern. In the second step, the CIR is derived analytically considering the configuration, position and orientation of the AUVs. This methodology reduces the computation time of the CIR, where the first step is computed once, regardless of the number of the receivers, while the second step is only repeated for each receiver.

* An appropriate transceiver is proposed by which the SDOC system achieves a high speed and longer communication range with low bit error rate (BER).

* We numerically investigate the CIR for different sea ices, seawater, and receiver configurations. As well, the performance of the SDOC system is evaluated considering the BER, normalized optical power penalty (NOPP), and maximum achievable bit rate.

The balance of the paper is organized as follows. In Section

II, we introduce the SDOC approach and the SSCL channel model. In Section III, we use the MCNRT method to model the upward transmission, then derive a quasi-analytic equation for the CIR. We introduce and model the proposed transceiver architecture in Section IV. In Section V, we numerically investigate the channel characteristics and system performance. Finally, conclusions are given in Section VI.

## II. PROPOSED SEA ICE DIFFUSING OPTICAL COMMUNICATION (SDOC) APPROACH

In this section we introduce the SDOC link as a new approach to establish communication between AUVs operating under sea ice. We discuss the temperature and salinity profile of the sea ice. Then, we introduce a new approach to model optical characteristics of the sea ice.

### A. SDOC ARCHITECTURE

As shown in Fig. 1, we consider a group of AUVs, for example an AUVs-swarm<fnref n="1" />, navigating several meters beneath a sheet of sea ice. The AUVs move together in the coordinated fashion with a separation of a few meters. In the proposed approach a broadcast communication link between the AUV transmitter (AUV-Tx) and the AUV receivers (AUV-Rxs) is accomplished in two steps: upward and downward transmissions. In the upward transmission, the AUV-Tx sends a narrow collimated laser beam toward the sea ice. Due to impurities (particles)<fnref n="2" />, the transmitted beam is subject to intense scattering at the surface and during propagation in the interior of the ice sheet. Inside the sea ice, a portion of the power will be transmitted through the sheet and lost to the atmosphere. Alternatively, the transmitted light may be trapped in the interior of the sheet where it is absorbed. Finally, a portion of the incident light will be diffused back from the ice sheet into the water. This diffused light which escapes the ice sheet is the useful signal which is used to establish the broadcast communication link. Given that the light is diffused inside the sheet, as shown in the green ellipse in Fig. 1, a wide coverage area is possible. The AUV-Tx can control the position of the diffusing spot by adjusting the direction of the laser beam, i.e., polar and azimuthal launching angles. For instance, if the AUV-Rxs are distributed symmetrically around the AUV-Tx, the beam should be vertically oriented toward the ice sheet to offer a fair coverage for all AUV-Rxs, as shown for the case in Fig. 1. However, if the AUV-Rxs are biased to one side, the AUV-Tx can orient its beam toward the direction of the AUV-Rxs to improve link quality. In the downward transmission from the ice sheet, the diffused beam propagates in the seawater and covers the AUV-Rxs with a large spot. Regardless of the position and orientation of the AUV-Rxs, each AUV-Rx receives a portion of this diffused beam, and the AUV-Tx establishes a broadcast communication with the AUV-Rxs.

<footnote n="1">¹Such a swarm typically employs a number of AUVs, however for simplicity, just five AUVs are shown in Fig. 1.</footnote>

<footnote n="2">²In this paper, the term of impurity and particle refer to any of: solid matter, dissolved matter, brine pockets, solid salt, air bubbles or air gaps.</footnote>
<figure><img src="image_2.png" /><figcaption>FIGURE 1: A topology for the SDOC approach: AUVs navigate underneath a sea ice and communicate with up and downward transmissions.</figcaption></figure>

The intensity of the diffused optical signal that emanates from the sea ice to the seawater depends on the density of impurities which contaminate the ice sheet as well as the sea ice surface roughness. The optical characteristics (e.g., absorption and scattering coefficients) of the ice sheet are highly affected by changes in impurity density which depend on ice sheet temperature and salinity [25], [26]. Temperature and salinity affect the freezing process of the sea ice which can introduce contaminants such as brine pockets, solid salts, and air bubbles. Given the high values for the temperature and salinity, the ice is most likely contaminated by particles and air bubbles [25]. For sea ice covered by snow, the optical properties will be impacted by temperature changes as well as the gaps between snowflake particles [26].

An example of the measured temperature and salinity profiles shown in Fig. 2. This figure represents the temperature and salinity of a 36 cm snow-covered sea ice sheet with 3 cm of snow cap and 33 cm of ice. The shown profiles are measured between November 2007 and June 2008 in the southern Beaufort Sea–Amundsen Gulf, Canadian Arctic [27]³. As shown in Fig. 2, the temperature $T(z)$ and salinity $S(z)$ change with the depth $z$ inside the ice sheet. The two curves in Fig. 2 can be well fitted by the following equations

$$T(z) = 0.2668 z - 10.74, \qquad (1)$$

$$S(z) = -3.24 \times 10^{-7} z^6 + 3.58 \times 10^{-5} z^5 - 1.47 \times 10^{-3} \\ \times z^4 + 2.74 \times 10^{-2} z^3 - 0.205 z^2 - 0.095 z + 13.63, \qquad (2)$$

where $T$ is the temperature in Celsius (°C), $S$ is the salinity in part per thousand (ppt), and $0 \le z \le 36$ cm. The equations are shown in the figure, and there is good agreement between the measured and the fitted profiles⁴.

Another example is a 12 cm bare-sea ice sheet whose temperature and salinity profiles are shown in [29, Fig. 3].

³Although the given profiles are for specific ice sheet, they hold the common linear relationship and C-shape for temperature $T$ and salinity $S$, respectively [28].

⁴The corresponding goodness of the fit criteria are; R-square= {0.9916, 0.9931} for the temperature and salinity curves, respectively.

<figure><img src="image_3.png" /><figcaption>FIGURE 2: The temperature and salinity profiles versus the sea ice depth for a snow-covered sea ice sheet as measured by [27].</figcaption></figure>

The sheet is young laboratory-grown saline sea ice. The two profiles of the sheet can be well fitted in $T(z)$ and $S(z)$ functions as⁵

$$T(z) = 1.176 z - 15.61, \quad 0 \le z \le 12 \text{cm} \qquad (3)$$

$$S(z) = 0.05003 z^2 - 0.7432 z + 8.203. \qquad (4)$$

These two ice sheet examples will be used later in the numerical results as case studies.

As shown in Fig. 2, the top surface of the sea ice is lower in the temperature than the bottom due to a cooling of the atmosphere and a warming of the seawater. As well, the salinity at the top and bottom is much higher than at the middle of the sea ice sheet. The vertical variations in the temperature and salinity with the thickness of the ice sheet result in changes in particle densities, which impact the channel optical characteristics. Given that the scattering inside the ice sheet is extensive and varies through the thickness of the ice sheet, channel modeling is challenging. In the following we introduce a simplified channel model.

## B. SEAWATER-SEA ICE CASCADED LAYERS (SSCL) CHANNEL MODEL

In this subsection, inspired by the geoscience literature [25], [26], we propose a SSCL channel model for upward trans-

⁵The corresponding goodness of the fit criteria are; R-square= {0.9997, 0.9788} for the temperature and salinity curves, respectively.
<figure><img src="image_4.png" /><figcaption>FIGURE 3: The SSCL channel model.</figcaption></figure>

mission in the SDOC approach. By SSCL, the vertical upward transmission link is modelled using cascaded layers of the seawater, ice, snow, and the air as shown in Fig. 3. Each of the layers of seawater and air are presented using a single layer since in the scale of few meters range, the particle densities do not change greatly with the depth [21], [22], [30], Ch. 3]. However, as mentioned in the previous subsection, the optical characteristics inside the ice sheet change continuously with the depth. Thus, we divide the ice sheet and snow layers into $m_i$ and $m_s$ of cascaded layers, respectively, over which the temperature and salinity are approximated as being uniform and presented using the average temperature $T(m)$ and average salinity $S(m)$. The thickness of each layer (and consequentially the number of layers) depends on the thickness of the sea ice sheet<fnref n="6" /> and the rates of change of the temperature and salinity profiles with the depth. Each layer in the SSCL model is characterized by thickness $d(m)$, lengths of $L_x(m)$ and $L_y(m)$ in x and y axes, respectively, and two rough interfaces between layer and the adjacent ones. By considering a constant temperature and salinity inside each layer, the particle density and the optical characteristics i.e., absorption coefficient $a(m)$, scattering coefficient $b(m)$, and effective refractive index $n_e(m)$ are also constant for each layer in the SSCL model.

As shown in Table 1, each layer in the SSCL model is composed of a mixture of particles, i.e., a hosting medium with additional impurities. For instance, ice layers are composed of the pure ice as a hosting medium with a mixture of particles (e.g., brine pockets, air bubble, solid salt, algae

and soot). However, snow layers are composed of air as a hosting medium with a fewer numbers of mixture particles (e.g., snow grains, algal and non-algal particles and soot). Due to these particles, the optical ray propagating inside the $m^{th}$ layer of the SSCL model suffers from absorption and scattering effects. The absorption coefficient, $a(m)$, is the weighted summation of the contribution from the mixture components as [56]

$$a(m) = f_{v_o} a_o(m) + \sum_{j=1}^{J_m} f_{v_j} a_j(m), \quad (5)$$

where $a_o$ and $f_{v_o}$ are the absorption coefficient and the volume fraction of the hosting medium, respectively. As well, $a_j$ and $f_{v_j}$ are the absorption coefficient and the volume fraction associated with the $j^{th}$ particle, respectively, where $f_{v_o} + \sum_{j=1}^{J_m} f_{v_j} = 1$. Symbol $J_m$ is the number of mixture particles in layer $m$, and the value of $J_m$ depends on the hosting medium of the layers and its surrounding environment. The hosting medium does not contribute to the scattering effect, thus, the scattering coefficients for each layer, $b(m)$, are weighted summations of the contribution from the impurity components only as [56]

$$b(m) = \sum_{j=1}^{J_m} f_{v_j} b_j(m) \quad (6)$$

where $b_j$ is the scattering coefficient associated with the $j^{th}$ particle.

Based on the assumptions given in [26] and [54], the one term Henyey-Greensteen (OTHG) function is a good approximation to the phase scattering function [57]

$$p_{\theta_s}(\theta_s, m) = \frac{1}{4\pi} \frac{1 - g(m)^2}{(1 + g(m)^2 - 2g(m) \cos(\theta_s))^{3/2}}, \quad (7)$$

where $g(m)$ is the asymmetry factor and $\theta_s$ is a scattering angle. The asymmetry factor is obtained using the weighted sum as [56]

$$g(m) = \frac{1}{b_m(m)} \sum_{j=1}^{J_m} b_j(m) g_j(m), \quad (8)$$

where $g_j$ is the asymmetry factor of the $j^{th}$ particle. The effective refractive index of the layer is computed using the volume fraction $f_{v_j}$ as [58]

$$n_e(m) = f_{v_o} n_o(m) + \sum_{j=1}^{J_m} f_{v_j} n_j(m), \quad (9)$$

where $n_o$ is the refractive index of the hosting medium, and $n_j$ is the refractive index of the $j^{th}$ particle.

The interfaces between the adjacent layers are assumed to be rough surfaces which leads to optical surface scattering at

<footnote n="6">⁶The thickness of the sheet depends on the climate and the location of the sea ice. For instance, Worby et al. [31] reported the mean and standard deviation of the ice and snow thickness in Arctic, e.g., 0.87 ± 0.91 and 0.16 ± 0.2 metres, respectively, with a correlation length in kilometre range.</footnote>
<table><caption>TABLE 1: Summary of references used to quantify the surface roughness and optical parameters of the SSCL model.</caption><thead><tr><th colspan="5">Surface Scatter Model</th></tr><tr><th>Medium</th><th>Interface Roughness</th><th>Fixed Particle Density with z?</th><th>Thickness</th><th>Layers Numbers</th></tr></thead><tbody><tr><td>Seawater and atmosphere</td><td rowspan="2">[32]–[35]</td><td>Yes ([21], [22], [30, Ch. 3])</td><td>Meters range [21], [30]</td><td>Single layer</td></tr><tr><td>Ice and snow</td><td>No ([21], [22])</td><td>Centimeter range [31], [36]</td><td>Multi layers [20]–[22]</td></tr><tr><th colspan="5">Particle Scatter Model</th></tr><tr><th>Mixture Particles<br>j = {1, 2, ..., J<sub>m</sub>}</th><th>Absorption Coefficient<br>a(m)</th><th>Scattering Coefficient<br>b(m)</th><th>Asymmetry Parameter<br>g(m)</th><th>Refractive Index<br>n(m)</th></tr><tr><td colspan="5" style="text-align: center;">Seawater Layer (m = 1)</td></tr><tr><td>Pure water</td><td>[37, Tab. 3]</td><td>[38, First term in Eq. (19)]</td><td>[39]</td><td>[40, Tab. 1]</td></tr><tr><td>Chlorophyll-a</td><td>[41, Second term in Eq. (16)]</td><td>[38, Second term in Eq. (19)]</td><td>[39]</td><td>[42]</td></tr><tr><td>Yellow substance</td><td>[43, Eq. (1)]</td><td>≈ Zero</td><td>Null</td><td>[42]</td></tr><tr><td>Equivalent Seawater</td><td>Eq. (5)</td><td>Eq. (6)</td><td>Eq. (8)</td><td>Eq. (9)</td></tr><tr><td colspan="5" style="text-align: center;">m<sub>i</sub> Ice Layers (m = [2, 3, ..., m<sub>i</sub> + 1])</td></tr><tr><td>Pure ice</td><td>[44, Fig. 3]</td><td>≈ Zero</td><td>Null</td><td>[40, Tab. 1]</td></tr><tr><td>Brine pocket</td><td>≈ Zero</td><td>[26], [25], [26]</td><td>[45, Fig. 7 (b)]</td><td>[45, Fig. 7]</td></tr><tr><td>Air bubble</td><td>≈ Zero</td><td>[26], [25], [20]</td><td>[45, Fig. 7 (b)]</td><td>[45, Fig. 7 a]</td></tr><tr><td>Solid salt</td><td>≈ Zero</td><td>[26], [25], [46]</td><td>[45, Fig. 7 (b)]</td><td>[45, Fig. 7 a]</td></tr><tr><td>Algae</td><td>[26, Fig. 7]</td><td>≈ Zero</td><td>Null</td><td>Ref. [42]</td></tr><tr><td>Soot</td><td>[42, Fig. 13]</td><td>≈ Zero</td><td>Null</td><td>[47]</td></tr><tr><td>Equivalent ice</td><td>Eq. (5)</td><td>Eq. (6)</td><td>Eq. (8)</td><td>Eq. (9)</td></tr><tr><td colspan="5" style="text-align: center;">m<sub>s</sub> Snow Layers (m = [m<sub>i</sub> + 2, ..., m<sub>i</sub> + m<sub>s</sub> + 1])</td></tr><tr><td>Air</td><td>≈ Zero</td><td>≈ Zero</td><td>Null</td><td>≈ 1</td></tr><tr><td>Snow grains</td><td>[44, Fig. 3]</td><td>Eqs. [26], [26]</td><td>[48, Fig. 4]</td><td>[40, Tab. 1]</td></tr><tr><td>Algal impurities</td><td>[26, Fig. 9]</td><td>≈ Zero</td><td>Null</td><td>[42]</td></tr><tr><td>Non-algal impurities</td><td>[26, Fig. 9]</td><td>≈ Zero</td><td>Null</td><td>[42]</td></tr><tr><td>Soot</td><td>[42, Fig. 13 (a)]</td><td>≈ Zero</td><td>Null</td><td>Ref. [47]</td></tr><tr><td>Equivalent snow</td><td>Eq. (5)</td><td>Eq. (6)</td><td>Eq. (8)</td><td>Eq. (9)</td></tr><tr><td colspan="5" style="text-align: center;">Atmosphere Layer (m = m<sub>i</sub> + m<sub>s</sub> + 2)</td></tr><tr><td>Gases (i.e. free space)</td><td>≈ Zero</td><td>≈ Zero</td><td>Null</td><td>1</td></tr><tr><td>Snowflakes</td><td>[49, Eq. (7)]</td><td>[50, Eq. (13)]</td><td>[51, Eq. (7)]</td><td>[40, Tab. 1]</td></tr><tr><td>Rain drops</td><td>[49, Eq. (6)]</td><td>[52, Eq. (7)]</td><td>[53], [54]</td><td>[40, Tab. 1]</td></tr><tr><td>Fog droplets</td><td>[49, Eq. (4)]</td><td>[55, Eq. (9)]</td><td>[53], [54]</td><td>[30, Eq. (27.6)]</td></tr><tr><td>Equivalent atmosphere</td><td>Eq. (5)</td><td>Eq. (6)</td><td>Eq. (8)</td><td>Eq. (9)</td></tr></tbody></table>

the entrance of each layer. The surface roughness of the interface is presented with the random height in the z direction for each point (x, y), which can be well described in the x and y directions using the two-dimensional Gaussian distribution as measured in [32], [33]. To generate a realization of the ice surface, a two dimensional Gaussian random variable is generated with independent components in x and y according to [59]

where z is the height at (x, y) point, and σx(m) and σy(m) are the RMS values in x and y directions<fnref n="7" />, respectively. As measured in [32], [33], the correlation between heights over the surface is well approximated using the two-dimensional generalized power-law function. Thus, to represent the correlation in space of the surface, the Gaussian realization can be filtered by a generalized power-law function. This function is given with one dimension in [34] and can be generalized to

two dimensions $p_{\rho_m}(\rho_x, \rho_y)$ as

$$p_{\rho_m}(\rho_x, \rho_y) = \exp\left(-\left[\left(\frac{\rho_x}{l_x(m)}\right)^\xi + \left(\frac{\rho_y}{l_y(m)}\right)^\xi\right]\right), \quad (11)$$

where $\rho_x$ and $\rho_y$ are the distances between correlated points in x and y directions, respectively, $l_x(m)$ and $l_y(m)$ are the correlation lengths in x and y directions, respectively. The value of $\xi$ depends on the geographical location of the sea ice sheet, and is equal to 1 and 2 in cases of exponential-correlated and Gaussian-correlated surfaces, respectively. Note that, the surface roughness includes parts of the ice suspended in seawater. Due to the low density of these parts, they typically float up toward the ice sheet and settle on its bottom surface [1].

For the reader convenience, a summary of the equations and parameter values needed to quantify surface and optical parameters of the SSCL layers are given in Table 1. The compositions of each layer in the SSCL model are given in the table with references and equations needed to calculate the optical characteristics of each material.

### III. THE SDOC LINK MODEL

In this section, we obtain an expression for the CIR of links between the AUV-Tx and the AUV-Rxs considering the effects of scattering, attenuation, as well as AUV-Rxs

<footnote n="7">⁷The experimental measurements in the Arctic and Antarctic regions revealed that the roughness parameters, RMS and correlation length, are in the millimetre and centimetre ranges, respectively [34], [35].</footnote>
<figure><img src="image_5.png" /><figcaption>FIGURE 4: A model for the AUV-Tx to an AUV-Rx link (note that, the scattering in the seawater is only shown in the downward transmission for the sake of illustration simplicity.).</figcaption></figure>

configuration, position and orientation. Here, we introduce a new methodology that consists of two steps to obtain the CIR. In the first step, due to dense scattering occurring in the interior of the sea ice sheet, the upward transmission is evaluated numerically using an MCNRT approach. The MCNRT method obtains the diffusing pattern (e.g., the green ellipse in Figs. 1 and 4) that feeds the AUV-Rxs. In the second step, the downward transmission from the bottom of the sea ice sheet to an AUV-Rx is modeled analytically under a single scattering assumption in the seawater layer. This two-step methodology reduces computational complexity where the upward transmission is evaluated once regardless the number of the AUV-Rxs. As well, the CIR equation is a function of the configuration, position and orientation of the AUV-Rxs.

Figure 4 shows a link model between AUV-Tx and an AUV-Rx. The distances and angles are measured relative to the diffusing axes, ($X_d, Y_d, Z_d$), which is centred at the bottom of the sea ice sheet. Relative to these axes, we assume that the AUV-Tx is located at ($x_o, y_o, z_o$) position and with aperture orientation polar and azimuth angles ($\theta_o, \phi_o$). While the AUV-Rx is located at ($\Delta_x, \Delta_y, \Delta_z$) position with

aperture polar and azimuthal inclination angles ($\theta_{in}, \phi_{in}$). Thus, the AUV-Rx position can be described using the position and orientation (PO) vector (5 × 1) as $\Delta_r := [\Delta_x; \Delta_y; \Delta_z; \theta_{in}; \phi_{in}]$. The AUV-Rx is equipped with a lens with diameter $D_r$ and field of view (FOV) of $\theta_{FOV}$.

A. UPWARD TRANSMISSION MODEL

As shown in Fig. 4, the AUV-Tx launches an optical beam with profile $I_o$, power $P_o$, wavelength $\lambda_o$, and beam width $W_o$ toward the sea ice. The center of the beam is presented by a ray $\vec{e}_o$ with directions ($\theta_o, \phi_o$) and a photon packet weight $w_o$ (equivalent to optical intensity). Angles ($\theta_o, \phi_o$) correspond to intended and non-intended orientation for the optical beam. An intended orientation when is the AUV-Tx directs the optical beam with a specific direction toward the ice sheet. A non-intended orientation occurs disturbances in the environment such as sea currents and waves. Without loss of generality, we assume the spot of the beam on the bottom of the sea ice is centered at the origin, i.e., (0, 0, 0). Thus, the position ($x, y$) of the AUV-Tx is obtained as $x_o = z_o \sin(\theta_o) \cos(\phi_o)$ and $y_o = z_o \sin(\theta_o) \sin(\phi_o)$. The depth
335 and orientation of the AUV-Tx are noted in a PO vector (3×1) as $\Delta_t = [z_o; \theta_o; \phi_o]$.

336 Given the challenge of using analytic approaches to obtain the diffusing pattern produced from the ice sheet in the upward transmission, an MCRT method is used instead. In MCNRT, many optical rays $\vec{e}_o$ are launched from the AUV-Tx to ensure the reliability of the result. The launched rays are diffused due to the surface and particle scattering taken place between and in the layers of the SSCL channel, respectively. The seawater, sea ice, snow and atmosphere layers contribute in producing the diffusing pattern, however, the sea ice and snow layers are the dominant contributors. The surface and particle scattering are simulated using geometric equations and numerical random process with associated PDFs, respectively, as given in the following subsections.

340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
1) Surface Scattering
Surface scattering occurs when the optical ray strikes the rough interface between the $m^{th}$ layer with refractive index $n_e(m)$ and the $(m+1)^{th}$ layer with refractive index $n_e(m+1)$ in upward propagation. Since $n_e(m) \neq n_e(m+1)$, the optical ray $\vec{e}_o$ incident on the interface with an angle $\theta_i(m)$ is split into a reflected ray $\vec{e}_1$ to the $m^{th}$ layer with an angle $\theta_r(m)$, where $\theta_i(m) = \theta_r(m)$, and a transmitted (i.e., refracted) ray $\vec{e}_2$ to the $(m+1)^{th}$ layer with an angle $\theta_t(m+1)$. The angle of transmitted ray between the $m^{th}$ and $(m+1)^{th}$ layers is given by [60]

$$ \theta_t(m+1) = \arcsin \left( \frac{n_e(m)}{n_e(m+1)} \sin(\theta_i(m)) \right), \quad (12) $$

362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
where, $\theta_i(m)$, $\theta_t(m+1)$, and $\theta_r(m)$ are measured relative to the local normal of the incident point which has a random direction due the randomness of the surface roughness. The reflection coefficient is computed for non polarized-light<fnref n="8" /> using angles $\theta_i(m)$ and $\theta_t(m+1)$ as [60]

$$ R_s(m) = 0.5 \left[ \left( \frac{\sin(\theta_i(m) - \theta_t(m+1))}{\sin(\theta_i(m) + \theta_t(m+1))} \right)^2 + \left( \frac{\tan(\theta_t(m+1) - \theta_i(m))}{\tan(\theta_i(m) + \theta_t(m+1))} \right)^2 \right], \quad (13) $$

367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
and the corresponding transmission coefficient is obtained as $T_s(m+1) = 1 - R_s(m)$. Accordingly, the reflected and transmitted rays, $\vec{e}_1$ and $\vec{e}_2$, propagate in $m^{th}$ and $(m+1)^{th}$ layers with packet weights $w_1(m) = w_o \times R_s$ and $w_1(m+1) = w_o \times T_s$, respectively.

372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
401
402
403
404
405
406
407
408
409
410
2) Particle Scattering
After the optical ray $\vec{e}_1$ enters the $m^{th}$ layer, it will propagate a random distance $\mu_{u_o}(m)$ with a likelihood of particle scattering $p_\mu(\mu_{u_o}, m)$ given as [61]

$$ p_{\mu}(\mu_{u_0}, m) = c(m) \exp[-c(m) \mu_{u_0}(m)]. \quad (14) $$

377
378
379
380
381
382
383
384
385
and the random distance is generated as [61]

$$ \mu_{u_o}(m) = -\frac{\log(1 - u_{\mu})}{c(m)}, \quad (15) $$

where $u_\mu$ is a uniform random variable, $u_\mu \sim U[0, 1]$, and $c(m)$ is the extinction coefficient of the $m^{th}$ layer representing the loss in the power of the ray. The value of the extinction coefficient $c(m)$ is computed as

$$ c(m) = a(m) + b(m). \quad (16) $$

When a scattering event occurs, the weight of the photon packet is dropped to [61]

$$ w_3(m) = w_1(m) \frac{b(m)}{c(m)}. \quad (17) $$

Upon scattering, the optical ray arriving from the direction $\vec{e}_1$ will have a new direction $\vec{e}_3$ determined randomly according to polar and azimuthal scattering angles $(\theta_{u_s}, \phi_{u_s})$. The angle $\theta_{u_s}(m)$ is generated from the OTHG PDF in Eq. (7) as [61]

$$ u_{\theta} = \int_{0}^{\theta_{u_s}(m)} p_{\theta_s}(\theta_s, m) \sin(\theta_s) d\theta_s, \quad (18) $$

where $u_\theta \sim U[0, 1]$. Also, the azimuthal scattering angle $\phi_{u_s}$ is typically described by a uniform PDF, and it is generated as [57]

$$ p_{\phi_s}(\phi_{u_s}) = \frac{1}{2\pi}, \quad \phi_{u_s} = 2\pi u_{\phi} \quad (19) $$

where $u_\phi \sim U[0, 1]$. After scattering, the ray travels a new distance $\mu_{u_1}$ with a new direction $\vec{e}_3$ before the next scattering occurs with likelihood $p_\mu(\mu_{u_1}, m)$. Compared to the seawater and the atmosphere, particle scattering takes place much more frequently in snow and sea ice layers. Typically, the optical ray is scattered few times in the seawater or atmosphere layer, however, hundreds of scattering events can typically take place in the sea ice or snow layers.

The MCNRT traces the optical rays until they are either absorbed, trapped in the ice layer, escape to the atmosphere, or diffuse back into the seawater. The diffused rays only contribute in the obtained diffusing pattern for the upward transmission and the remainder of the rays are considered as lost. For a given position and orientation for the AUV-Tx, $\Delta_t$, the normalized diffusing pattern is obtained with the intensity $I_d$ as a function of the space, angles and time as follows

$$ I_d(x_d, y_d, \theta_d, \phi_d, t_d | \Delta_t) = \text{MCNRT}\{\text{SSCL}, \Delta_t, I_o, \lambda_o, W_o\} \quad (20) $$

where, as shown in Fig. 4, the intensity $I_d$ is measured on the bottom of the sea ice surface at position $x_d$ and $y_d$, with polar

<footnote n="8">⁸Modelling using non-polarized light is typical case of scattered light.</footnote>
413 $\theta_d$, azimuth $\phi_d$ angles, and time $t_d$. As well, the DC gain of the upward transmission $G_u$ is computed using $I_d$ as

$$G_u = \int_{-L_x/2}^{L_x/2} \int_{-L_y/2}^{L_y/2} \int_0^{\pi/2} \int_0^{2\pi} \int_0^{\infty} I_d(x_d, y_d, \theta_d, \phi_d, t_d) \\ dt_d d\phi_d d\theta_d dy_d dx_d, \qquad (21)$$

415 where $L_x = \max\{L_x(m=2), \dots, L_x(m_i+m_s+1)\}$ and $L_y = \max\{L_y(m=2), \dots, L_y(m_i+m_s+1)\}$ are the considered lengths of the SSCL channel in x and y axes, respectively.

## B. DOWNWARD TRANSMISSION MODEL

420 Figure 4 shows a model for the downward transmission which corresponds to the link from the bottom of the sea ice to the AUV-Rx through the seawater channel. A diffused ray emitted from a position $(x_d, y_d, 0)$ in the direction of $\vec{e}_d$ is represented in the figure, where $\vec{e}_d$ is defined as

$$\vec{e}_d = \vec{x}_d \sin(\theta_d) \cos(\phi_d) + \vec{y}_d \sin(\theta_d) \sin(\phi_d) + \vec{z}_d \cos(\theta_d), \qquad (22)$$

425 where $(\vec{x}_d, \vec{y}_d, \vec{z}_d)$ are the unit vectors in the direction of $(X_d, Y_d, Z_d)$ axes. The impurities in the seawater cause absorption and scattering for the diffused ray $\vec{e}_d$. Under a single scattering assumption, which is reasonable here because the link is short, $\vec{e}_d$ arrives to the AUV-Rx either with LOS (i.e., non-scattering) or after one scattering with the direction $\vec{e}_r$. In the LOS path, the direction is maintained (i.e. $\vec{e}_d = \vec{e}_r$) and the ray arrives with arrival position $(x_r^o, y_r^o, z_r^o)$.

430 In the scattering path, let $(x_s, y_s, z_s)$ denote the position of the scattering event relative to $(X_d, Y_d, Z_d)$ axes, and with polar and azimuthal scattering angles $(\theta_{d_s}, \phi_{d_s})$ relative to the axes of the scattering $(X_s, Y_s, Z_s)$, as shown in Fig. 4. The scattering angles $\theta_{d_s}$ and $\phi_{d_s}$ are computed using Eqs. (18) and (19) by replacing angles $\theta_{u_s}$ and $\phi_{u_s}$ with angles $\theta_{d_s}$ and $\phi_{d_s}$, respectively. The scattered ray arrives to the receiver with polar and azimuthal arrival angles $(\theta_r, \phi_r)$ measured relative to the sea ice axes, $(X_d, Y_d, Z_d)$. For given scattering angles $(\theta_{d_s}, \phi_{d_s})$, the arrival angles $(\theta_r, \phi_r)$ are computed from $(\theta_d, \phi_d)$ as follows. Let $e_d$ be $(3 \times 1)$ vector, represented in $(X_d, Y_d, Z_d)$ as $e_d=[\sin(\theta_d) \cos(\phi_d); \sin(\theta_d) \sin(\phi_d); \cos(\theta_d)]$. Then, $e_d$ is rotated around $(Y_s, X_s, Z_s)$ axes by two angles: $\theta_y = \arcsin(\cos(\phi_{d_s}) \sin(\theta_{d_s}))$ and $\theta_x = \arcsin(\sin(\phi_{d_s}) \sin(\theta_{d_s})/\cos(\theta_y))$ respectively. Thus, $\theta_r$ and $\phi_r$ are computed as

$$\theta_r = \arccos([0, 0, 1] \mathbf{R}_x(\theta_x) \mathbf{R}_y(\theta_y) \mathbf{e}_d), \\ \phi_r = \arcsin\left(\frac{[0, 1, 0]}{\sin(\theta_d)} \mathbf{R}_x(\theta_x) \mathbf{R}_y(\theta_y) \mathbf{e}_d\right), \qquad (23)$$

450 where $\mathbf{R}_x(\theta_x)$ and $\mathbf{R}_y(\theta_y)$ are $(3 \times 3)$ rotation matrices around $(X_d$ and $Y_d$ axes, respectively [62]). The arrival vector $\vec{e}_r$ is expressed with respect to the axes of the sea ice as

$$\vec{e}_r = \vec{x}_d \sin(\theta_r) \cos(\phi_r) + \vec{y}_d \sin(\theta_r) \sin(\phi_r) + \vec{z}_d \cos(\theta_r). \qquad (24)$$

453 Vector $\vec{e}_r$ is also characterized by arrivals angles $(\dot{\theta}_r, \dot{\phi}_r)$ measured relative to the axes, $(X_r, Y_r, Z_r)$, as shown in the Fig. 4, and can be equivalently written as

$$\vec{e}_r = \vec{x}_r \sin(\dot{\theta}_r) \cos(\dot{\phi}_r) + \vec{y}_r \sin(\dot{\theta}_r) \sin(\dot{\phi}_r) + \vec{z}_r \cos(\dot{\theta}_r), \qquad (25)$$

455 where $(\vec{x}_r, \vec{y}_r, \vec{z}_r)$ are the unit vectors relative to the Rx axes $(X_r, Y_r, Z_r)$. For the given angles $(\theta_r, \phi_r)$, the angles $(\dot{\theta}_r, \dot{\phi}_r)$ are calculated from Eq. (23) by replacing $e_d$ with $e_r=[\sin(\theta_r) \cos(\phi_r); \sin(\theta_r) \sin(\phi_r); \cos(\theta_r)]$ and substituting $\theta_y = \arcsin(\cos(\phi_{in}) \sin(\theta_{in}))$ and $\theta_x = \arcsin(\sin(\phi_{in}) \sin(\theta_{in})/\cos(\theta_y))$. The scattered ray arrives at arrival position $(x_r^s, y_r^s, z_r^s)$ over the aperture of the AUV-Rx.

460 The arriving ray from the LOS or scattering path is detected if the position of arrival $(x_r, y_r, z_r)$ is located on the lens of the AUV-Rx with arrival angles $(\dot{\theta}_r, \dot{\phi}_r)$ less than half angle of the FOV. This can be compactly represented as the geometric loss $G_g$ and it is written as

$$G_g(\Delta_r) = \begin{cases} 1, & \text{if } (x_r, y_r, z_r) \in f_p(D_r, \Delta_r) \text{ and } \dot{\theta}_r \le \frac{\theta_{FOV}}{2} \\ 0, & \text{otherwise,} \end{cases} \qquad (26)$$

465 where $f_p(D_r, \Delta_r)$ represents the spatial extent of the AUV-Rx lens with respect to the sea ice axes $(X_d, Y_d, Z_d)$.

470 1) Case 1: Low Scattering Seawater

473 Consider the case of seawaters with small scattering coefficient (e.g., pure seawater) where the impact of scattering is negligible. In this case, only the LOS component need to be considered [63], [64]. In the LOS path, the direction is maintained (i.e. $\vec{e}_d = \vec{e}_r$), and the amplitude of the optical ray is attenuated according to the Beer-Lambert law. The LOS ray arrives with arrival position $(x_r^o, y_r^o, z_r^o)$, shown in Fig. 4, and is computed as [65]

$$\begin{aligned} x_r^o &= x_d + \Delta_z \tan(\theta_d) \cos(\phi_d), \\ y_r^o &= y_d + \Delta_z \tan(\theta_d) \sin(\phi_d), \\ z_r^o &= \Delta_z. \end{aligned} \qquad (27)$$

476 For rays diffused from a single point on the bottom of the sea ice $(x_d, y_d, 0)$, the CIR can be well approximated by a linear combination of LOS components as

$$P^o(t_r, x_d, y_d) \approx P_o \int_0^\infty \int_0^{2\pi} \int_0^{\pi/2} [\exp(-l_r^o c) G_g(\Delta_r) \\ I_d(x_d, y_d, \theta_d, \phi_d, t_d) \delta\left(t_r - \left(t_d + \frac{l_r^o}{\nu}\right)\right)] d\theta_d d\phi_d dt_d \qquad (28)$$

480 where the length of the LOS path is computed geometrically from the figure as $l_r^o = \sqrt{(x_d - x_r^o)^2 + (y_d - y_r^o)^2 + (z_d - z_r^o)^2}$ The symbols $t_r$ and $\nu$ are the arrival time and the light speed in the seawater, respectively, and $\delta(.)$ is the Dirac-delta function.
<figure><img src="image_6.png" /><figcaption>FIGURE 5: The proposed system architecture for the SDOC approach.</figcaption></figure>

power during the on slots. Consider a LD with green wavelength ($\lambda_o$ = 532 nm) given its relatively low attenuation in seawater [69]. A beam expander is the LD implemented using two lenses, one lens for beam diverging and other one for beam collimating. This collimated wide beamwidth optical beam helps in transmitting more optical power while keeping the constraint of the maximum permissible exposure (MPE) optical power on the eye<fnref n="9" /> regarding eye safety.

## 2) Receiver Optoelectronics

The proposed Rx architecture is shown in Fig. 5b. First, the AUV-Rx uses a hemispherical concentrator which is implemented using a hemispherical non-imaging lens coated by a bandpass optical filter as shown. Such a lens with a relatively large diameter, $D_r$, and a wide FOV, $\theta_{FOV}$, is desired to collect much of diffused rays to compensate SDOC high channel attenuation. As well, an optical filter with narrow bandwidth, $\Delta\lambda$, is preferred to eliminate the incoming background radiation from the sun during the daytime. The concentrator is an essential component in the SDOC approach especially with high background radiation levels at $\lambda_o$ = 532 nm [72]. The concentrator is quantified by its gain $G_c$ which depends on its refractive index, $n_c$, and the FOV as [73]

$$G_c(\dot{\theta}_b) = \begin{cases} \frac{n_c^2}{\sin(\theta_{FOV}/2)^2} & \text{if } \dot{\theta}_b \le \theta_{FOV}/2, \\ 0 & \text{if } \dot{\theta}_b > \theta_{FOV}/2, \end{cases} \quad (38)$$

where $\dot{\theta}_b$ is the incident angle of the background ray upon the concentrator and it is measured relative to the optical axis of the Rx, $Z_r$, as shown in Fig. 4. As well, the optical band pass filter is quantified by its transmission coefficient $T(\dot{\theta}_b)$ which

depends on the incident angle of the received ray. Such hemispherical concentrators are commercially available and have been used in optical diffusing communication systems for indoor applications<fnref n="10" />. The concentrator enlarges the effective area of the PD, $A_{ef}$, which means capturing solar noise. The effective area of the PD is obtained as [73]

$$A_{ef}(\dot{\theta}_b) = \begin{cases} A_{PD} T(\dot{\theta}_b) G_c(\dot{\theta}_b) \cos(\dot{\theta}_b) & \text{if } \dot{\theta}_b \le \theta_{FOV}/2 \\ 0 & \text{if } \dot{\theta}_b > \theta_{FOV}/2 \end{cases} \quad (39)$$

where $A_{PD}$ denotes the physical active area of the PD. Here, for simplicity, the dependence of the effective area on the incident angle $\dot{\theta}_b$ is represented by replacing $A_{ef}(\dot{\theta}_b)$ by its average $\overline{A_{ef}}$ over the incident angle, while making two assumptions. Firstly, we assume that the function $T(\dot{\theta}_b)$ can be replaced by its average, $\overline{T}$, over all incident angles. This assumption holds, especially, when the incident optical ray arrives within a wide range of the angles which is the typical case of diffusing communications [73]. Secondly, we assume a uniform PDF for $\dot{\theta}_b$. Then, the average Rx effective area is obtained as

$$\overline{A_{ef}} = \frac{2}{\pi} \int_{0}^{\theta_{FOV}/2} A_{ef}(\dot{\theta}_b) d\dot{\theta}_b = \frac{2 A_{PD} \overline{T} n_c^2}{\pi \sin(\theta_{FOV}/2)} \quad (40)$$

Note that, enlarging the FOV decreases the average effective area of the Rx.

After the hemispherical concentrator, a silicon PIN photodiode (PIN-PD) with a trans-impedance amplifier (TIA) is used. The PIN-PD converts the collected optical rays to an electrical current proportionally to its responsitivity $\Re$ and $A_{PD}$. Then, the TIA converts the small current to a high voltage proportionally to its load resistance $R_L$. In contrast to avalanche photodiodes, photo-multiplier tube and SiPM PDs, the silicon PIN-PD achieves a better performance when the background radiation is much high and dominates the receiver noises [76], [77].

<footnote n="9">⁹The typical optical powers used in underwater communication experiments are on the order of fraction of Watt [9], and are far below levels needed to alter the ice surface [70]. Though direct human contact with UAVs is possible, safety must also be considered to preserve wildlife which may interact with these optical emissions [71].</footnote>

<footnote n="10">¹⁰The optical concentrator and filter with the mentioned specifications can be implemented [74]. However, some customization may be required for use in underwater applications [75].</footnote>
### 3) Channel Equalization

Using the described setup, the Rx can overcome the effects of high channel attenuation and background noise. A low pass filter (LPF) is employed after the TIA to eliminate any out-of-band noise, where the filter bandwidth is adjusted according to the actual data rate. The bandwidth of the filter is adopted according to the link speed to maximize the system performance. The output signal of the LPF is sampled with sampling rate $T_s$, where $T_s < T_b/2$ to avoid aliasing [78]. The sampled signal is then processed by a discrete-time channel equalizer to reduce the impact of ISI. Among the available equalization schemes, the digital decision feedback equalizer (DFE) is chosen due to the mobility of the AUVs [79]. With proper training, the DFE can adapt itself to the changing channel conditions and the PO vector. As well, the DFE coupled with the least mean squares (LMS) algorithm has the advantage of simplicity and is a good choice for non-fading dispersive channels [79]. As shown in Fig. 5b, the DFE has two input branches namely, feed forward (FF) and feedback (FB). The input through the FF is the electrical signal from the output of the sampler $v_{ei}(k)$, where $k$ indicates to the index of the received sample. While, the input through the FB is the output of the OOK demodulator $\hat{v}_{eo}(k)$, where $\hat{v}_{eo}(k) \in \{1, 0\}$. The output of the equalizer is the summation of the weighted inputs as follows [79]

$$v_{eo}(k T_b) = \sum_{j_F=0}^{N_{FF}} \alpha_{j_F} v_{ei}(k T_b - j_F T_b) + \sum_{j_B=1}^{N_{FB}} \beta_{j_B} \times \hat{v}_{eo}(k T_b - j_B T_b), \quad (41)$$

where $\alpha_{j_F}$ and $\beta_{j_B}$ are the FF and FB weighting coefficients, respectively. Symbols $N_{FF}$ and $N_{FB}$ indicates the number of the tabs used in the FF and FB filters, respectively. The DFE has two operation modes, training and tracking modes. In the training mode, the Tx sends a training sequence which is known to the Rx. The DFE adopts the LMS algorithm to obtain the optimal values for $\alpha_{j_F}$ and $\beta_{j_B}$ [79]. In the tracking mode, the DFE uses the optimal values obtained for the gains vector to eliminate the effect of ISI in the transmitted OOK symbols. In the next subsection, we discuss the effect of non-ideal performance of the DFE and the different noise sources on the SINR of the SDOC system.

## B. SIGNAL-TO-INTERFERENCE-PLUS-NOISE RATIO ANALYSIS

During training, the filter coefficients are adapted based on output of OOK decision device and the training sequence [78]. In the tracking mode, assuming that training was successful, error propagation at the output of the decision device should be minimized. Assuming an absence of decision errors, a simple linear model of the DFE output voltage can

be approximated as

$$v_{eo}(k T_s, \mathbf{\Delta}_r) \approx \underbrace{v_n(k T_s, \mathbf{\Delta}_r)}_{\text{noise}} + \underbrace{P_{DFE}(k T_s) \otimes (R_L \Re [p(\tau) * P(\tau, \mathbf{\Delta}_r)] |_{\tau_s=k T_s})}_{\text{signal + residual ISI}} \quad (42)$$

where, $v_n$ is the sampled noise voltage, and $\mathbf{\Delta}_r$ is PO vector of the AUV-Rx as defined in the previous section, $P_{DFE}(k T_s)$ is the sampled system impulse response of the DFE, $\otimes$ is the discrete convolution operator, $*$ is continuous convolution and $p(t)$ is the instantaneous transmitted optical optical power. The signal in Eq. (42) can be decomposed as the sum of the desired signal, denoted by $v_s(\mathbf{\Delta}_r)$, and the residual ISI denoted by $v_{isi}(\mathbf{\Delta}_r)$ where,

$$v_s(\mathbf{\Delta}_r) \approx v_n(k T_s, \mathbf{\Delta}_r) + \sum_{k=\tau_d/T_s}^{(T_b+\tau_d)/T_s} P_{DFE}(k T_s) \otimes (R_L \Re \times [p(\tau) * P(\tau, \mathbf{\Delta}_r)] |_{\tau_s=k T_s}), \quad (43)$$

where $\tau_d$ is the time delay of the channel and $T_b$ and $\tau_d$ are assumed to be multiples of $T_s$. Also,

$$v_{isi}(\mathbf{\Delta}_r) \approx v_n(k T_s, \mathbf{\Delta}_r) + \sum_{k=(T_b+\tau_d)/T_s}^{\infty} P_{DFE}(k T_s) \otimes (R_L \Re [p(\tau) * P(\tau, \mathbf{\Delta}_r)] |_{\tau_s=k T_s}). \quad (44)$$

The noise contribution in (42) includes the effects of the thermal $v_{th}$ and shot $v_{sh}$ noises, i.e., $v_n = v_{th} + v_{sh}$. The thermal noise $v_{th}$ is well described by zero mean Gaussian distribution with variance $\sigma_{th}^2$ given as [80], [81]

$$\sigma_{th}^2 = 4 R_L K [T(m=1) + 273.15] B, \quad (45)$$

where, $K$ is the Boltzmann constant and $T(m=1)$ is the temperature of the seawater layer in Celsius as defined in Sec. II-B. Usually, the temperature of the seawater underneath sea ice is $T(m=1) \le 0$ °C, as shown in Fig. 2. On the other hand, the shot noise is associated with the superposition of the desired signal voltage $v_s$, the ISI distortion voltage $v_{isi}$, and the background radiation voltage $v_{sun}$. Due to the high intensity of the solar radiation, the shot noise can be modeled using Gaussian random process with variance given as [80]

$$\sigma_{sh}^2(\mathbf{\Delta}_r; \Lambda) = 2 R_L q [\Lambda v_s(\mathbf{\Delta}_r) + v_{isi}(\mathbf{\Delta}_r) + v_{sun}(\mathbf{\Delta}_r)] B, \quad (46)$$

where $q$ is the electron charge in electron-volt units and $\Lambda = 1$ and $\Lambda = 0$ with on and off of the OOK symbols, respectively. The value of $v_{sun}$ is quantified as [76]

$$v_{sun}(\mathbf{\Delta}_r) = \begin{cases} R_L \Re (1 - G_u) \Delta \lambda E_s \overline{A_{ef}} \cos(\theta_{inc}) \\ \exp(K_d \Delta_z) \\ 0 \quad \text{if} \quad \theta_{inc} > \theta_{FOV}/2 \end{cases} \quad (47)$$

where, $E_s$ is the spectral solar intensity with unit Watt/(m$^2$. nm), and $K_d$ is the light diffusion coefficient
in the seawater. The value of $E_s$ depends on the weather conditions, as well the zenith angle of the sun [72]. The zenith angle is in range 25° to 90° in Arctic and Antarctic regions where frozen oceans exist, and it records a minimum value during the summer seasons [82]. The light diffusing coefficient is related to the seawater parameters by $K_d = a(m = 1) + 0.03 b(m = 1)$. The factor of $(1 - G_u)$ represents the transmission coefficient of the sea ice sheet. This means a thicker sea ice sheet assists in raising the gain of the upward transmission and in reducing the received background radiations.

The mean $\eta_{\Lambda}$ and the variance $\sigma_{\Lambda}^2$, $\Lambda \in \{0, 1\}$, of the total signal and noise affects system performance are given as

$$ \left\{ \begin{aligned} \eta_{\Lambda}(\Delta_{\mathbf{r}}; \Lambda) &= \Lambda v_s(\Delta_{\mathbf{r}}) + v_{isi}(\Delta_{\mathbf{r}}) + v_{sun}(\Delta_{\mathbf{r}}) \\ \sigma_{\Lambda}^2(\Delta_{\mathbf{r}}; \Lambda) &= \sigma_{isi}^2(\Delta_{\mathbf{r}}) + \sigma_{sh}^2(\Delta_{\mathbf{r}}; \Lambda) + \sigma_{th}^2 \end{aligned} \right. \qquad (48) $$

where $\sigma_{isi}^2$ is the variance of ISI signals and it is equal to $R_L v_{isi}^2/4$. Thus, the instantaneous SINR, $\gamma(\Delta_{\mathbf{r}})$, is obtained as

$$ \gamma(\Delta_{\mathbf{r}}) = \frac{[\eta_1(\Delta_{\mathbf{r}}) - \eta_0(\Delta_{\mathbf{r}})]^2}{[\sigma_1(\Delta_{\mathbf{r}}) + \sigma_0(\Delta_{\mathbf{r}})]^2} = \frac{v_s^2(\Delta_{\mathbf{r}})}{[\sigma_1(\Delta_{\mathbf{r}}) + \sigma_0(\Delta_{\mathbf{r}})]^2} \qquad (49) $$

In the numerical results, we consider three configurations for the AUV-Rxs, namely, unequalized AUV-Rx (Rx-UE), AUV-Rx with DFE (Rx-E), and AUV-Rx with perfect DFE (Rx-PE). The BER of the Rx-E and Rx-UE systems are evaluated numerically using Monte Carlo simulations. However, the BER of the Rx-PE system is evaluated by eliminating the effect of ISI from (49), i.e. $v_{isi} = 0$, then using the well known AWGN channel as [83]

$$ p_e(\Delta_{\mathbf{r}}) = Q\left(\sqrt{\gamma(\Delta_{\mathbf{r}})}\right) \qquad (50) $$

where $Q(x) = 1/\sqrt{2\pi} \int_{x}^{\infty} \exp(-[z/\sqrt{2}]^2) dz$.

## V. NUMERICAL RESULTS AND DISCUSSIONS

In this section, we numerically evaluate the diffusing patterns of upward transmission, the CIRs of downward transmission, and overall system performance. We assume that the AUV-Tx is located at depth $z_o = 2$ m, and perfectly orientated to the bottom of the sea ice, i.e., $\theta_0 = \phi_0 = 0$ and $(x_o, y_o) = (0, 0)$. The AUV-Tx is equipped with a laser source which emits a collimated beam $I_o$ with uniform intensity, wavelength $\lambda_o = 532$ nm, average transmitted power $P_o \le 200$ mW, and a width of $W_o = 5$ cm to maintain eye-safety. However, we assume the AUV-Rx moves in the $x-y$ plane where the length of the downward transmission does not exceed the limit of the single scattering assumption. Note that the single scattering assumption is valid with lengths 13.5 and 6.6 m for clear and coastal seawaters, respectively [65]. Note that in the following results, the particular values for parameters of the AUV-Rx were not optimized for communication performance but were chosen to demonstrate the range of operating conditions of the SSCL channel.

We consider four SSCL channels, namely, a clear and coastal seawater with a snow-covered sea ice sheet (Cl-S and

Co-S channels) and the clear and coastal seawater with a bare sea ice sheet (Cl-B and Co-B channels). The snow-covered sea ice sheet has a thickness of 36 cm and it well described by Eqs. (1), (2). The bare sea ice sheet has thickness 12 cm and its temperature and salinity profiles are described by Eqs. (3), (4). We use Eqs. (1)-(9) and Table 1 to calculate the optical and roughness parameters associated with each SSCL layer, where the results are given in Table 2.

The bare-ice cases are divided into 6 layers while the snow-covered cases are divided into 9 layers$^{11}$. In all cases, each layer is assigned with the average values of the salinity and temperature using Eqs. (1)-(4), as shown in Tables 2a and 2b. Clear weather above the sea ice sheets is assumed, which is the typical case during sunny days. As shown in Tables 2a and 2b, the scattering coefficients of the snow-covered sea ice sheet and coastal seawater are higher than that for bare sea ice sheet and clear seawater, respectively. In addition, it is clear that the changes in the refractive indices and asymmetry parameters are small. In Table 2c, the RMS of the roughness and correlation length, are assumed in millimetre and centimetre ranges, respectively, as measured in [34], [35]. As well, we assume isotropic layers (i.e., $\sigma_x(m) = \sigma_y(m)$ and $l_x(m) = l_y(m)$), and the interfaces are Gaussian-correlated (i.e., $\xi = 2$) [32], [33]. The interfaces between the ice layers are assumed smooth due to fact that the variation in the effective refractive indices are negligible in the presented cases. To ensure an accurate realization for the SSCL model, the roughness is sampled with intervals and lengths with values $\delta_x(m) = \delta_y(m) = 0.1 l_x(m)$ and $L_x(m) = L_y(m) = 60 l_x(m)$ [84].

## A. RESULTS FOR UPWARD TRANSMISSION

Figure 6 shows the marginalized diffusing patterns$^{12}$ for the Cl-B and Cl-S channels with the orange and maroon colors, respectively. The diffusing pattern is measured at the bottom of the sea ice, i.e., $\Delta_z = 0$, with DC gains of $G_u = 0.26$ and 0.37 for Cl-B and Cl-S channels, respectively. These results were obtained by running the MCNRT using the ZeMax Opticstudio software [85] over $10^6$ iterations. Note that we have verified that increasing the number of iterations to $10^7$ resulted in almost identical results.

Figures 6a and 6b show the marginalized diffusing patterns versus the polar and the azimuthal angles, $I_{d\downarrow\theta}$ and $I_{d\downarrow\phi}$, respectively. As shown in these figures, the marginalized intensity is uniform with respect to (w.r.t.) $\phi_d$, however, it is oriented w.r.t. $\theta_d$ with a peak at $\theta_d \approx 45^{\circ}$. The orientation indicates non-specular diffusing due to the dense scattering occurred in the sea ice and snow. The value of $45^{\circ}$ is interrupted as follows; each diffusing point on the sea ice is an identical random variable described by Eq. (18), and the diffusing pattern is a summation of that diffusing points. Assuming the central limit theory, $I_{d\downarrow\theta}$ approaches

<footnote>$^{11}$This is done as a compromise between the accuracy and the computational complexity of the MCNRT method.</footnote>

<footnote>$^{12}$The marginalized diffusing pattern with $x_d$ variable, for instance, is obtained by integrating Eq. (20) over all remaining variables.</footnote>
<caption>TABLE 2: The parameters of seawater bare sea ice and seawater snow-covered sea ice cascaded models.</caption>

<table><caption>(a) The Cl-B and Co-B SSCL channels.</caption><thead><tr><th>Layer No. (m)</th><th>a(m) [1/m]</th><th>b(m) [1/m]</th><th>g(m)</th><th>n<sub>e</sub>(m)</th></tr></thead><tbody><tr><td>m = 6 (Clear Air, T = -14 C<sup>o</sup>, S ≈ 0 ppt)</td><td>0</td><td>0</td><td>0</td><td>1</td></tr><tr><td>m = 5 (Ice, T = -13.56 C<sup>o</sup>, S = 7.19 ppt)</td><td>0.563</td><td>480.473</td><td>0.9894</td><td>1.3494 - 0.0395i</td></tr><tr><td>m = 4 (Ice, T = -9.635 C<sup>o</sup>, S = 5.91 ppt)</td><td>0.492</td><td>422.4211</td><td>0.9906</td><td>1.3496 - 0.0395i</td></tr><tr><td>m = 3 (Ice, T = -6.38 C<sup>o</sup>, S = 5.53 ppt)</td><td>0.483</td><td>473.9</td><td>0.9923</td><td>1.3499 - 0.0395i</td></tr><tr><td>m = 2 (Ice, T = -2.6 C<sup>o</sup>, S = 5.91 ppt)</td><td>0.721</td><td>996.81</td><td>0.9946</td><td>1.3504 - 0.0395i</td></tr><tr><td>m = 1 (Clear Seawater, T = 0 C<sup>o</sup>, S = 5.91 ppt)</td><td>0.069</td><td>0.08</td><td>0.8708</td><td>1.333</td></tr><tr><td>m = 1 (Coastal Seawater T = 0 C<sup>o</sup>, S = 5.91 ppt)</td><td>0.088</td><td>0.216</td><td>0.9470</td><td>1.333</td></tr></tbody></table>

<table><caption>(b) The Cl-S and Co-S SSCL channels.</caption><thead><tr><th>Layer No. (m)</th><th>a(m) [1/m]</th><th>b(m) [1/m]</th><th>g(m)</th><th>n<sub>e</sub>(m)</th></tr></thead><tbody><tr><td>m = 9 (Clear Air, T = -12 C<sup>o</sup>, S ≈ 0 ppt)</td><td>0</td><td>0</td><td>0</td><td>1</td></tr><tr><td>m = 8 (Snow, T = -11.5 C<sup>o</sup>, S ≈ 0 ppt)</td><td>0.282</td><td>3.1593 × 10<sup>3</sup></td><td>0.8878</td><td>1.1620 - 0.0395i</td></tr><tr><td>m = 7 (Ice, T = -10.06 C<sup>o</sup>, S = 11.9 ppt)</td><td>0.532</td><td>845.81</td><td>0.9900</td><td>1.3445 - 0.0395i</td></tr><tr><td>m = 6 (Ice, T = -8.56 C<sup>o</sup>, S = 8.99 ppt)</td><td>0.4661</td><td>645.03</td><td>0.9903</td><td>1.3451 - 0.0395i</td></tr><tr><td>m = 5 (Ice, T = -6.57 C<sup>o</sup>, S = 7.77 ppt)</td><td>0.463</td><td>652.64</td><td>0.9913</td><td>1.3453 - 0.0395i</td></tr><tr><td>m = 4 (Ice, T = -4.65 C<sup>o</sup>, S = 6.66 ppt)</td><td>0.451</td><td>760.68</td><td>0.9926</td><td>1.3455 - 0.0395i</td></tr><tr><td>m = 3 (Ice, T = -3.63 C<sup>o</sup>, S = 7.961 ppt)</td><td>0.4532</td><td>724.321</td><td>0.9930</td><td>1.3454 - 0.0395i</td></tr><tr><td>m = 2 (Ice, T = -2.34 C<sup>o</sup>, S = 7.97 ppt)</td><td>0.684</td><td>1334.457</td><td>0.9943</td><td>1.3457 - 0.0395i</td></tr><tr><td>m = 1 (Clear Seawater, T = 0 C<sup>o</sup>, S = 7.97 ppt)</td><td>0.069</td><td>0.08</td><td>0.8708</td><td>1.333</td></tr><tr><td>m = 1 (Coastal Seawater T = 0 C<sup>o</sup>, S = 7.97 ppt)</td><td>0.088</td><td>0.216</td><td>0.9470</td><td>1.333</td></tr></tbody></table>

<table><caption>(c) The roughness parameters for the interfaces between the layers of the SSCL channel models [34], [35].</caption><thead><tr><th>The Interface</th><th>σ<sub>x<sub>m</sub></sub> = σ<sub>y<sub>m</sub></sub> [mm]</th><th>l<sub>x<sub>m</sub></sub> = l<sub>y<sub>m</sub></sub> [mm]</th><th>The Interface</th><th>σ<sub>x<sub>m</sub></sub> = σ<sub>y<sub>m</sub></sub> [mm]</th><th>l<sub>x<sub>m</sub></sub> = l<sub>y<sub>m</sub></sub> [mm]</th></tr></thead><tbody><tr><td>Snow-Air</td><td>1.2</td><td>50</td><td>Ice-Snow</td><td>2</td><td>75</td></tr><tr><td>Ice-Air</td><td>5</td><td>120</td><td>Seawater-Ice</td><td>2.2</td><td>100</td></tr></tbody></table>

the Gaussian with mean 45° which the mean of the range; 0-90 degrees. In addition, the marginalized intensity in case of Cl-S channel is relatively higher than the case of Cl-B channel. Specifically, the peaks of the marginalized intensities in Fig. 6a are 6 × 10⁻⁴ and 4 × 10⁻⁴ for Cl-S and Cl-B channels, respectively. Furthermore, the marginalized intensities in Fig. 6b are 3.5 × 10⁻⁴ and 2.5 × 10⁻⁴ for Cl-S and Cl-B channels, respectively. The pattern in these figures can be fit to two-dimensional Lambertian and uniform functions in θ<sub>d</sub> and φ<sub>d</sub> respectively, for both Cl-B and Cl-S channels<fnref n="13" />

$$I_{d\downarrow\theta,\phi}^{\text{Cl-S}}(\theta_d, \phi_d) = 4.438 \times 10^{-5} \cos^{6.6}(\theta_d - 0.248 \pi), \qquad (51a)$$

$$I_{d\downarrow\theta,\phi}^{\text{Cl-B}}(\theta_d, \phi_d) = 6.032 \times 10^{-5} \cos^{6.5}(\theta_d - 0.242 \pi). \qquad (51b)$$

Figure 6c shows the marginalized intensities $I_{d\downarrow x}$ versus the distance $x_d = [-0.5, 0.5]$ m. The intensities decay exponentially with peaks 1 × 10⁻³ and 1.8 × 10⁻³ at the center, $x_d = 0$ m, for the Cl-B and Cl-S SSCL channels, respectively, and almost zero value at $|x_d| = 0.5$ m. Due to the uniform value of the marginalized intensities w.r.t. $\phi_d$, the intensity profiles for $x_d$ and $y_d$ are similar and can be fitted

with the following two-dimensional functions for the Cl-B and Cl-S SSCL channels as<fnref n="14" />

$$I_{d\downarrow x, y}^{\text{Cl-B}}(x_d, y_d) = \frac{0.591}{10^3} \exp(-10.95 |x_d| - 11.3 |y_d|), \qquad (52a)$$

$$I_{d\downarrow x, y}^{\text{Cl-S}}(x_d, y_d) = \frac{1.466}{10^3} \exp(-15.41 |x_d| - 15.46 |y_d|). \qquad (52b)$$

Though, the diffusing pattern has a small spot on the bottom of the sea ice sheet (i.e., $|x_d|$ and $|y_d| \le 0.5$ m), due to the orientation with angle 45°, the spot expands out with the propagation in the seawater as shown in the next subsection.

Figure 6d shows the marginalized diffusing pattern $I_{d\downarrow t}$ (i.e., temporal dispersion patterns of the upward transmission) with $t_d = [2, 24]$ ns. The pattern of the Cl-S channel has a high peak with amplitude 14 × 10⁻³ and it decays slowly with a long dispersion time due to the thickness and much particle scattering occurred for the laser beam in the channel as can be seen from Table 2b (i.e., a larger thickness, and higher temperature and salinity values). In contrast to the Cl-S channel, the pattern of the Cl-B channel has two peaks with amplitudes 32 × 10⁻⁴ and 26 × 10⁻⁴. The time interval between the two peaks is nearly equal to the time taken by the optical ray to propagate from the bottom to the surface of the ice sheet. Thus, the shown dispersion pattern can reveal information about the thickness of the bare sea ice sheets whilst performing a communication function. The

<footnote n="13">¹³In this paper, the fitting is accomplished using the tool of curve fitting in Matlab [86, CFTOOL]. The goodness of the fit R-square= {0.8574, 0.8714} and RMSE= {5.95 × 10⁻⁶, 7.606 × 10⁻⁶} for Cl-B and Cl-S channels, respectively.</footnote>

<footnote n="14">¹⁴The corresponding goodness of the fit criteria are; R-square= {0.9053, 0.9261} and RMSE= {1.522 × 10⁻⁵, 2.505 × 10⁻⁵} for Cl-B and Cl-S channels, respectively.</footnote>
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 6: The marginalized diffusing patterns (a) $I_{d\downarrow\theta}(\theta_d)$, (b) $I_{d\downarrow\phi}(\phi_d)$, (c) $I_{d\downarrow x}(x_d)$, and (d) $I_{d\downarrow t}(t_d)$. The results are shown with the orange and maroon colors for the Cl-B and Cl-S channels, respectively. The fitted curves are indicated in black dashed and dotted lines for the Cl-B and Cl-S channels, respectively.</figcaption></figure>

shown time dispersion patterns can be fitted to a sum of Gaussian functions in $t_d$ as<fnref n="15" />

$$I_{d\downarrow t}(t_d) = \sum_{i=1}^{3} \alpha_i \exp \left( - \left( \frac{t_d - \beta_i}{\gamma_i} \right)^2 \right). \quad (53)$$

Equations (51)-(53) serve as a guide for a future analytic model for the upward transmission. The equations are shown with dotted lines in Fig. 6. We notice a good agreement between the equations and MCNRT results in space but less accuracy for the temporal dispersion patterns. Note that we also tested other fitting functions proposed in the literature for underwater CIRs in other scenarios (such as double gamma weighted [87], combination of exponential and arbitrary power [88], and Beta Prime distributions [89]), however, Eq. (53) provided a much better in fit for SDOC. In fact,

the fitting is challenging due to the dense scattering taken place in the channel. Thus, further investigation is required to obtain more accurate equation as a future work.

## B. RESULTS FOR DOWNWARD TRANSMISSION

In this subsection, we demonstrate numerical results for the CIR, DC gain $H_o$, and the delay spread $\tau_{RMS}$. The results are obtained using equations (34)-(37), and take into account the effects of the type of the sea ice, seawater, Rx configuration, and the position of the AUV-Rxs. The position and FOV parameters in the following were chosen to show there scope of operating characteristics for the SSCL channel. The optimization of these parameters for maximize communication performance is left as future work.

### 1) Impact of Sea Ice

Figures 7 shows the normalized received power versus arrival time (i.e., CIR) for the case of coastal seawaters and different types of ice sheet, namely, Co-S, Co-B and coastal-

<footnote n="15">¹⁵For Cl-B channel, the coefficients $\alpha_i$, $\beta_i$ and $\gamma_i$ are {2.405 × 10⁻³, 1.931 × 10⁻³, 1.025 × 10⁻³}, {9.093 × 10⁻⁹, 1.041 × 10⁻⁸, 1.121 × 10⁻⁸}, {7.041 × 10⁻¹⁰, 6.492 × 10⁻¹⁰, 7.043 × 10⁻¹⁰}, respectively. As well, For Cl-S channel, the coefficients $\alpha_i$, $\beta_i$ and $\gamma_i$ are {1.059 × 10⁻², 1.547 × 10⁻³, 5.333 × 10⁻⁴}, {9.274 × 10⁻⁹, 1.206 × 10⁻⁸, 1.492 × 10⁻⁸}, {1.465 × 10⁻⁹, 1.817 × 10⁻⁹, 2.262 × 10⁻⁹}, respectively.</footnote>
845

846

847

848

849

850

851

852

853

854

855

856

857

858

859

860

861

862

863

864

865

pure (Co-P)<fnref n="16" /> channels. The AUV-Rx has the parameters $\theta_{FOV} = 140°$ and $D_r = 15$ cm, and is located at the position ($\Delta_x = 2, \Delta_y = 0, \Delta_z = 3$ m). As shown in the figure, the snow-covered sea ice sheet records the highest CIR amplitude and the largest dispersion thanks to the dense scattering occurring through its layers, as given in Table 2b. The Co-B channel shows a lower CIR amplitude and a relatively narrow dispersion due to a lower scattering coefficient as compared to the Co-S channel, see Tables 2a and 2b. The CIR of the pure sea ice sheet channel records the smallest amplitude and dispersion because there are no particles to scatter from inside the sheet. This result is likely to arise when the sea ice is thinned, such as when a part of the sea ice sheet melts in the summer season. The channel time delay, $\tau_d$, takes its smallest value in the case of Co-S channel, which due to the fact that the second layer in the Co-S channel ($m = 2$) has a larger scattering coefficient with contrast to the second layer in the Co-B channel, see Tables 2b and 2a. Numerically, the peaks of the CIRs are $3.1 \times 10^{-6}, 2.4 \times 10^{-6}$ and $3.2 \times 10^{-8}$, and the delay spreads are $15 \times 10^{-9}, 8 \times 10^{-9}$ and $4 \times 10^{-9}$ sec for the Co-S, Co-B and Co-P SSCL channels, respectively.

866

867

868

869

870

871

872

873

874

875

876

877

878

879

880

881

882

883

884

885

886

## 2) Impact of Seawater

Figure 8 shows the CIR for an AUV-Rx with $\theta_{FOV} = 90°$ and $D_r = 15$ cm) at position ($\Delta_x = 3, y = 0, \Delta_z = 2$ m) below a bare sea ice sheet. The CIRs are shown for the Co-B, Cl-B and a Pu-B SSCL channels, where Pu-B denotes pure seawater cascaded with the bare sea ice sheet¹⁷ (i.e., $a(m = 1) = 0.053 \text{ m}^{-1}, b(m = 1) = 0.003 \text{ m}^{-1}$) [64]. Here, we used the bare-sea ice which has less scattering compared to snow-covered sea ice, this makes the effect of the seawater on the channel more significant. At a distance of $\Delta_x = 3$ m from the AUV-Tx, the FOV does not see the diffusing spot on the bottom of the sea ice. Thus, the amplitude of the CIR depends on beam scattering in the sea water. As shown in the figure, the case of coastal seawater has the highest amplitude and largest dispersion due to particle scattering. However, pure seawater provides the AUV-Rx with the less significant CIR. Numerically, the peaks of the CIRs are $5.2 \times 10^{-9}, 3.4 \times 10^{-9}$ and $2 \times 10^{-10}$, as well, and the delay spreads are $8 \times 10^{-9}, 8 \times 10^{-9}$ and $4 \times 10^{-9}$ sec for the Co-B, Cl-B, and Pu-B SSCL channels, respectively.

887

888

889

890

891

892

893

## 3) Effects of FOV

Figure 9 shows the DC gain and delay spread of the channel versus the Rx FOV for AUV-Rx located at position ($\Delta_x = 3, \Delta_y = 0, \Delta_z = 2$ m). In general, increasing the FOV leads to the collection of more rays and improves the DC gain. However, the rate of change in the DC gain with the FOV

<figure><img src="image_9.png" /><figcaption>FIGURE 7: The effects of the sea ice on CIR with ($\theta_{FOV} = 140°, D_r = 15$ cm) and position {$\Delta_x, \Delta_y, \Delta_z$} = {2, 0, 3} m.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 8: The effects of the seawater on CIR with ($\theta_{FOV} = 90°, D_r = 15$ cm) and position{$\Delta_x, \Delta_y, \Delta_z$} = {3, 0, 2} m.</figcaption></figure>

($\partial h_o / \partial FOV$) depends on the location of the AUV-Rx with respect to the diffusing surface. For the given case study in Fig. 9 and according to the geometry of the topology, the receiver aperture begins to receive a direct signal from the diffusing surface at a computed FOV = 102.7° and receives signals from the complete diffusing surface at a computed FOV $\ge 120.5°$. The computed FOVs are shown in the figure with values 93° and 122°, respectively, due to the impacts of the orientation of the diffusing beam with angle 45° and scattering occurring in the coastal seawater. This observation can help explain the results given in the figure as follows. When the FOV changes from 36° to 93.6°, the rate of change in $\partial h_o / \partial FOV$ is 0.0456 per degree. As the FOV increases further, it starts to collect rays with high energy from the diffusing surface. Thus, when the FOV changes from 93.6° to 122°, the rate of change increases to $\partial h_o / \partial FOV$ is 0.2137 per degree. Increasing the FOV further ($FOV \ge 122°$), there is no additional improvement in the DC gain since nearly all power is collected.

894

895

896

897

898

899

900

901

902

903

904

905

906

907

908

909

910

911

912

<footnote n="16">¹⁶Co-P SSCL channel is the coastal seawater cascaded with a free-impurity sea ice, i.e., a perfect transparent sea ice. This pure sea ice rarely exists on the frozen oceans, and it is considered here just as benchmark.

¹⁷The pure seawater rarely exists underneath the frozen oceans, and it is considered here for comparison.</footnote>
<figure><img src="image_11.png" /><figcaption>FIGURE 9: The effects of the FOV on $h_o$ and $\tau_{RMS}$ with $D_r$ = 10 cm at position $\{\Delta_x, \Delta_y, \Delta_z\}$ = $\{3, 0, 2\}$ m and Co-B channel.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 10: The effects of the depth on $h_o$ and $\tau_{RMS}$ with ($D_r$ = 10 cm, $\theta_{FOV}$ = 90°) at x-y position $\{\Delta_x, \Delta_y\}$ = $\{3, 0\}$ m and Co-B channel.</figcaption></figure>

On the other hand, the RMS delay spread depends on the time of diffusing $t_d$ (20) in addition to the time taken to propagate from the bottom of the sea ice to the lens of the Rx (i.e., the distance $\mu_{d_o} + \mu_{d_1}$, see Fig. 4). In general, the value of $t_d$ is a smaller for diffused rays that leave the sea ice close to the origin of the diffusing spot than for those rays that are further away. However, the propagation time from the sea ice to Rx for rays near the diffusing spot is longer than those further away. The RMS delay spread of the link is thus impacted by the balance of diffusing and propagation times. Qualitatively, when $\theta_{FOV} \le 97^\circ$, the Rx does not see the diffusing spot origin directly and the RMS delay spread is dominated by $t_d$. That is, the total time of propagation will be close to the mean value resulting in a smaller RMS delay spread. However, as FOV increases, i.e., $97^\circ \le \theta_{FOV} \le 107^\circ$, the received rays from the diffusing spot with longer propagation time dominate increasing the delay spread. Finally, for $\theta_{FOV} \ge 107^\circ$, the AUV-Rx receives diffused rays arriving from both the origin of the diffusing spot, $\{x_d, y_d\} \approx 0$, as well as diffused rays over a wider area of the ice sheet which contributes to a reduction in the delay spread.

## 4) Impact of Depth

Figure 10 shows the DC gain and delay spread of the channel versus the depth, $\Delta_Z$, for an AUV-Rx with lens diameter $D_r$ = 10 cm and $\theta_{FOV}$ = 90°. The AUV-Rx is located at a relatively long distance from the origin of the sea ice, ($\Delta_x$ = 3m, $\Delta_y$ = 0) m. The AUV-Rx captures more diffused rays with increasing depth in the range $\Delta_z$ = [1, 3] m, then, the power captured decays with range for $\Delta_z \ge 3$ m, as shown. This phenomena can be interpreted as follows. The spatial coverage of the diffusing pattern, in the x-y plane, extends with the depth due to two reasons. Firstly, the orientation of the diffusing pattern with the polar angle

as shown in Fig. 6a. Secondly, the scattering taking place in coastal seawater contributes more in extending the spatial coverage of the diffusing pattern. However, for $\Delta_z \ge 3$ m, the DC gain decays with the depth, due to the absorption taking place in the coastal seawater which dominates the impact of scattering. Numerically, the rate of change in the gain with the depth, ($\partial h_o / \partial(\Delta_z)$), is fixed in the range $\Delta_z$ = [1.5, 2.5] m with value $\partial h_o / \partial(\Delta_z)$ = 0.6860 per meter, however, it is higher in the range $\Delta_z$ = [2.5, 3] m with value $\partial h_o / \partial(\Delta_z)$ = 2.1419 per meter. On the other hand, the delay spread reaches to its minimum value at depth $\Delta_z$ = 2.7m as shown. This occurs since the lens (with $\theta_{FOV}$ = 90° and at location $\Delta_x$ = 3) captures the LOS rays diffused from points close to the diffusing spot on the bottom of the sea ice. These LOS rays arrive with high amplitude and small propagation times, resulting in the RMS delay spread attaining its minimum value.

## 5) The Spatial Distributions of $H_o$ and $\tau_{RMS}$

Figure 11 shows the spatial distributions of the DC channel gain and the RMS delay spread versus the position of the AUV-Rx in the x-y plane. The results are shown for Co-S channel within the area of $6 \times 6$ m². As well, Table 3 summarizes statistical values of the DC and RMS delay spread and contrasts the results with the case of the Co-B channel. The results are associated to an AUV-Rx located at $\Delta_z$ = 3 and equipped with a lens with $D_r$ = 15 and cm, $\theta_{FOV}$ = 140°. These settings for the AUV-Rx are used in the remainder of the numerical results.

As shown in Fig. 11a, the DC gain distribution is symmetric in the x-y plane around the center ($\Delta_x$ = 0, $\Delta_y$ = 0) and the DC gain value decreases monotonically with $\Delta_x$ and $\Delta_y$. The shown distribution matches with the average response from the results in Figs. 6b and 6c. As well, as given in the table, the DC gain values in case of the Co-S channel are
<table><caption>TABLE 3: Extracted statistics from Fig. 11</caption><thead><tr><th>Parameters</th><th>The Co-S channel</th><th>The Co-B channel</th></tr></thead><tbody><tr><td>Minimum h<sub>o</sub></td><td>5.82 × 10<sup>-6</sup></td><td>6.2 × 10<sup>-6</sup></td></tr><tr><td>Maximum h<sub>o</sub></td><td>1.239 × 10<sup>-4</sup></td><td>0.7 × 10<sup>-4</sup></td></tr><tr><td>Average h<sub>o</sub></td><td>3.587 × 10<sup>-5</sup></td><td>2.33 × 10<sup>-5</sup></td></tr><tr><td>Minimum τ<sub>RMS</sub></td><td>8.74 × 10<sup>-10</sup></td><td>5.50 × 10<sup>-10</sup></td></tr><tr><td>Maximum τ<sub>RMS</sub></td><td>1.53 × 10<sup>-9</sup></td><td>1.085 × 10<sup>-9</sup></td></tr><tr><td>Average τ<sub>RMS</sub></td><td>1.073 × 10<sup>-9</sup></td><td>0.77 × 10<sup>-9</sup></td></tr></tbody></table>

higher than that in case of the Co-B channel due to the dense scattering taking place in the snow cap. Numerically, from Table 3, the maximum values of the DC gain are $0.7 \times 10^{-4}$ and $1.239 \times 10^{-4}$ and the average values are $2.33 \times 10^{-5}$ and $3.587 \times 10^{-5}$ for the Co-B and Co-S channels, respectively.

In Fig. 11b, the RMS delay spread spatial distribution is shown. In the area under the diffusing surface, the main amount of the power arrives from the LOS rays from the diffusing spot. Thus, these rays take the shortest path and the corresponding RMS delay spread has the lowest value in this area. At the edge of the considered area, the majority of received power arises from scattered rays. Thus, the corresponding RMS delay spread has the highest value. In the intermediate area, the RMS delay spread value fluctuates with the position of the AUV-Rx depending on whether a LOS or scattering components dominate. As well, from Table 3, the RMS delay spread is on the order of nanoseconds, however, the values in case of the Co-S channel are higher than that in case of the Co-B channel. Numerically, the maximum values of the RMS delay spread are $1.085 \times 10^{-9}$ s and $1.53 \times 10^{-9}$ s and the average values are $0.77 \times 10^{-9}$ s and $1.073 \times 10^{-9}$ s for the Co-B and Co-S channels, respectively.

### C. LINK PERFORMANCE

In this subsection, we numerically investigate the BER performance and maximum achievable bit rate for proposed system. The AUV-Rx is equipped with $D_r = 15$ cm with $\theta_{FOV} = 140^{\circ}$ and located at y-z position $\{\Delta_y, \Delta_z\} = \{0, 3\}$ m, load resistance is $R_L = 200$ $\Omega$, and the electrical bandwidth of the Rx is considered as 0.7 GHz. A DFE equalizer is implemented using 15 taps $T_b$-spaced branches. The coefficients of the taps are obtained using 2024 training symbols, and the LMS algorithm runs with control value equal to 0.15. The ISI, shot, and thermal noises are evaluated using Eqs. (44)-(48). For the background radiation, clear weather is assumed with the sun at zenith angle equal to $\approx 60^{\circ}$ [72], [82]. For the thermal noise, the temperature of the seawater is assumed zero Celsius, see Fig. 2.

Figure 12 shows the average BER performance versus the distance $\Delta_x$ for a perfect equalizer receiver (Rx-PE), i.e., $v_{isi} = 0$, where performance limitation only arises from the Rx noise, dominated by background radiation. Here, we consider Co-B and Co-S channels, and the Rxs are equipped with optical filters with bandwidths $\Delta\lambda \in \{1, 5, 10\}$ nm. The average transmitted optical power is $P_o = 100$ mW and the bit rate is $R_b = 50$ Mbps. As shown, the BER

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 11: The distributions of $h_o$ and $\tau_{RMS}$ with ($D_r = 15$ cm, $\theta_{FOV} = 140^{\circ}$, $\Delta_z = 3$) for the Co-S channel.</figcaption></figure>

performance degrades with distance and improves by decreasing the bandwidth of the optical filter. As well, the BER performance in the case of the Co-S channel is better than Co-B channel for two reasons. Firstly, the Co-S channel has a higher upward transmission DC gain; secondly, the Co-S channel reduces impact of the solar radiations much more than the Co-B channel. For example, considering a BER threshold of $10^{-3}$ as indicated by the green line in the figure, the AUV-Tx can communicate with the Rx-PE at ranges $\Delta_x = \{4, 3, 2.75\}$ and $\{3.5, 2.75, 1\}$ m with the bandwidth $\Delta\lambda = \{1, 5, 10\}$ in cases of the Co-S and Co-B channels, respectively. In other words, scaling $\Delta\lambda$ down by 10 times raises the communication range by 45% and 250% in cases of Co-S and Co-B channels, respectively.

Figure 13 compares the normalized optical power penalty (NOPP) versus the normalized RMS delay spread (NRDS) defined as

$$NRDS = \frac{\tau_{RMS}}{T_b}$$

at BER$= 10^{-3}$ for receivers with equalization (Rx-E) and unequalized (Rx-UE). The NOPP is defined as the required transmitted optical power to achieve the desired FEC limit
<figure><img src="image_15.png" /><figcaption>FIGURE 12: The BER for the Rx-PE (perfect equalization) system versus the distance $x$ with $R_b = 50$ Mbps and $P_o = 100$ mW.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 13: The normalized optical power penalty (NOPP) versus the normalized RMS delay spread (NRDS) for Rx-E and Rx-UE (no equalization) cases.</figcaption></figure>

in cases of Rx-E and Rx-UE systems normalized by that required in case of the Rx-PE system. The RMS delay spread $\tau_{RMS}$ is computed for the AUV-Rx at position $\Delta_x = 2$ m, where $\tau_{RMS} = 1.1 \times 10^{-9}$ s and $8.5 \times 10^{-10}$ s for the Co-S and Co-B SSCL channels, respectively. As well, the bit duration is varied in the range $T_b \in [2, 100]$ ns, i.e., $R_b \in [10, 500]$ Mbps. The case of Rx-UE is used as a benchmark to highlight the benefit of channel equalization.

At low data rates, e.g., ($NRDS \le 0.05$), where the bit duration is much larger than the RMS delay spread, the effect of ISI on the system performance is limited and the performance of Rx-UE and Rx-E are nearly the same. As the data rate increases, the impact of ISI increases and Rx-E gradually outperforms Rx-UE. Specifically, for the Co-B channel at $NRDS= 0.075$, Rx-E and Rx-UE require $NOPP= 2.15$ dB and $NOPP= 2.52$ dB, respectively. For the Co-S channel at

<figure><img src="image_17.png" /><figcaption>FIGURE 14: The maximum achieved bit rate for the Rx-E system versus the distance with FEC limit of BER $= 10^{-3}$.</figcaption></figure>

the same NRDS, Rx-E and Rx-UE require $NOPP= 2$ dB and $NOPP= 2.5$ dB, respectively. At higher data rates of NRDS= 0.2, Rx-E and Rx-UE require $NOPP= 6.3$ dB and $NOPP= 3.5$ dB, respectively, for the Co-B channel. For the Co-S channel at the same NRDS, Rx-E and Rx-UE require $NOPP= 3.2$ dB and $NOPP= 5.8$ dB, respectively. These results indicate that the equalizer improves the power efficiency of the systems by nearly 3 dB, which means the required transmitted power is reduced roughly by a factor of two. In other words, the AUV with the equalized system enhances the power-efficiency of the AUVs which means more lifetime for the battery.

Figure 14 shows the maximum achievable bit rate under the constraint $BER \le 10^{-3}$ versus the distance $\Delta_x$ with average transmitted optical power $P_o \in \{100, 200\}$ mW. As shown in the figure, the maximum achievable bit rate ($R_b \approx 700$ Mpbs) is achieved directly under the diffusing surface ($\Delta_x \le 1$ m). However, as $\Delta_x$ increases, the maximum achievable bit rate decreases; the proposed system can achieve broadcast data rates on the order of $R_b = 1$ Mpbs over communication ranges of $\Delta_x = 6$ m. As indicated by the green dashed line, to maintain a communication rate of 10 Mbps, scaling the transmitted power by 2 increases the communication range by 18% and 10% in cases of Co-B and Co-S channels, respectively. This trade off between data rate and coverage distance should be considered during planning stage of the AUV swarms, based on the required data rate and range.

## VI. CONCLUSIONS

In this paper, for first time, we propose a broadband-broadcast approach suitable for networking AUVs under sea ice, albeit with limited range. We take advantage of existing ice sheets on the sea surface to establish a diffusing communication systems. The SSCL model was introduced in which the channel is represented in the form of cascaded layers with uniform optical characteristics. Due to the challenge of analytic modeling of optical signal scattering inside the ice sheet, MCNRT is used to evaluate the diffusing pattern of

<footnote>$^{18}$At distance $\Delta_x = 2$ m, a data rate of 500 Mpbs is considered as a maximum since the average transmitted power is limited to 200 mW for the Rx-E systems, as indicated in the simulation parameters.</footnote>
upward transmission. For downward transmission, the CIR was derived in the form of a quasi-analytic equation assuming single scattering light propagation. Due to the expected effects of ISI and relatively high background solar power noise, we propose a new transceiver architecture that helps in mitigating the effects of these factors. We also provide extensive numerical results to investigate the effects of water and ice types, Rx parameters i.e., FOV and optical filter bandwidth, and the Rx location on the system performance.

The challenges in implementing SDOC systems includes the transceiver size which must be carefully chosen depending on the size of the AUV. The transmitted power must also be determined according to battery-life and eye-safety constraints. Lastly, the SDOC approach is not appropriate bellow transparent sea ice sheets which rarely exist in practice on frozen oceans. Future work includes further investigations for quasi-analytic forms of the diffusing patterns, investigating a better fit for the temporal diffusing pattern as well as experimental validation of the obtained results.

## REFERENCES

[1] NSIDC, "All about sea ice," [Available at] https://nsidc.org/cryosphere/seaicelink/index.html, [Accessed on]: 2021-04-01.

[2] P. Wadhams, J. Wilkinson, and A. Kaletzky, "Sidescan sonar imagery of the winter marginal ice zone obtained from an AUV," *Journal of Atmospheric and Oceanic Technology*, vol. 21, no. 9, pp. 1462-1470, 2004.

[3] N. Milard, P. Stevenson, M. Brandon, C. Banks, K. Heywood, M. Price, P. Dodd, A. Jenkins, K. Nicholls, D. Hayes *et al.*, "Instruments and methods autonomous underwater vehicles (AUVs) and investigations of the ice-ocean interface in antarctic and arctic waters," *Journal of Glaciology*, vol. 54, no. 187, p. 661, 2008.

[4] T. Maksym, H. Singh, C. Bassett, A. Lavery, L. Freitag, F. Sonnichsen, and J. Wilkinson, "Oil spill detection and mapping under arctic sea ice using autonomous underwater vehicles," *Final Report BSEE Contract E12PC00053*, 2014.

[5] E. Sahin and A. F. Winfield, "Special issue on swarm robotics." *Swarm Intelligence*, vol. 2, no. 2-4, pp. 69-72, 2008.

[6] C. Osterloh, T. Pionteck, and E. Maehle, "Monsun II: A small and inexpensive auv for underwater swarms," in *ROBOTIK 2012; 7th German Conference on Robotics*, 2012, pp. 1-6.

[7] N. Vedachalam, R. Ramesh, V. B. N. Jyothi, V. Doss Prakash, and G. Ramadas, "Autonomous underwater vehicles-challenging developments and technological maturity towards strategic swarm robotics systems," *Marine Geore sources & Geotechnology*, vol. 37, no. 5, pp. 525-538, 2019.

[8] U. Behrje, A. Amory, B. Meyer, and E. Maehle, "System identification and sliding mode depth control of the micro AUV sembio," in *ISR 2018; 50th International Symposium on Robotics*. VDE, 2018, pp. 1-8.

[9] Z. Zeng, S. Fu, H. Zhang, Y. Dong, and J. Cheng, "A survey of underwater optical wireless communications," *IEEE Communications Surveys Tutorials*, vol. 19, no. 1, pp. 204-238, 2017.

[10] F. Miramirkhani and M. Uysal, "Visible light communication channel modeling for underwater environments with blocking and shadowing," *IEEE Access*, vol. 6, pp. 1082-1090, 2018.

[11] A. E.-R. A. El-Fikky, M. E. Eldin, H. A. Fayed, A. A. E. Aziz, H. M. H. Shalaby, and M. H. Aly, "NLoS underwater VLC system performance: static and dynamic channel modeling," *Appl. Opt.*, vol. 58, no. 30, pp. 8272-8281, Oct 2019.

[12] F. R. Gfeller and U. Bapst, "Wireless in-house data communication via diffuse infrared radiation," *Proceedings of the IEEE*, vol. 67, no. 11, pp. 1474-1486, 1979.

[13] G. Yun and M. Kavehrad, "Spot-diffusing and fly-eye receivers for indoor infrared wireless communications," in *IEEE International Conference on Selected Topics in Wireless Communications*, 1992, pp. 262-265.

[14] —, "Indoor infrared wireless communications using spot diffusing and fly-eye receivers," *Canadian Journal of Electrical and Computer Engineering*, vol. 18, no. 4, pp. 151-157, 1993.

[15] S. T. Jivkova and M. Kavehrad, "Multispot diffusing configuration for wireless infrared access," *IEEE Transactions on Communications*, vol. 48, no. 6, pp. 970-978, 2000.

[16] A. G. Al-Ghamdi and J. M. H. Elmirghani, "Line strip spot-diffusing transmitter configuration for optical wireless systems influenced by background noise and multipath dispersion," *IEEE Transactions on Communications*, vol. 52, no. 1, pp. 37-45, 2004.

[17] S. Arnon and D. Kedar, "Non-line-of-sight underwater optical wireless communication network," *J. Opt. Soc. Am. A*, vol. 26, no. 3, pp. 530-539, Mar 2009.

[18] W. Liu, D. Zou, Z. Xu, and J. Yu, "Non-line-of-sight scattering channel modeling for underwater optical wireless communication," in *2015 IEEE International Conference on Cyber Technology in Automation, Control, and Intelligent Systems (CYBER)*, 2015, pp. 1265-1268.

[19] N. Anous, M. Abdallah, M. Uysal, and K. Qaraqe, "Performance evaluation of LOS and NLOS vertical inhomogeneous links in underwater visible light communications," *IEEE Access*, vol. 6, pp. 22408-22420, 2018.

[20] T. C. Grenfell, "A radiative transfer model for sea ice with vertical structure variations," *Journal of Geophysical Research: Oceans*, vol. 96, no. C9, pp. 16991-17001, 1991.

[21] Z. Jin, K. Stamnes, W. Weeks, and S.-C. Tsay, "The effect of sea ice on the solar energy budget in the atmosphere-sea ice-ocean system: A model study," *Journal of Geophysical Research: Oceans*, vol. 99, no. C12, pp. 25281-25294, 1994.

[22] C. D. Mobley, G. F. Cota, T. C. Grenfell, R. A. Maffione, W. S. Pegau, and D. K. Perovich, "Modeling light propagation in sea ice," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 36, no. 5, pp. 1743-1749, 1998.

[23] R. W. Preisendorfer, "Secchi disk science: Visual optics of natural waters 1," *Limnology and oceanography*, vol. 31, no. 5, pp. 909-926, 1986.

[24] Z. Lee, K. Du, R. Arnone, S. Liew, and B. Penta, "Penetration of solar radiation in the upper ocean: A numerical model for oceanic and coastal waters," *Journal of Geophysical Research: Oceans*, vol. 110, no. C9, 2005.

[25] G. F. Cox and W. F. Weeks, "Equations for determining the gas and brine volumes in sea-ice samples," *Journal of Glaciology*, vol. 29, no. 102, pp. 306-316, 1983.

[26] B. Hamre, J.-G. Winther, S. Gerland, J. J. Stamnes, and K. Stamnes, "Modeled and measured optical transmittance of snow-covered first-year sea ice in Kongsfjorden, Svalbard," *Journal of Geophysical Research: Oceans*, vol. 109, no. C10, 2004.

[27] G. Carnat, T. Papakyriakou, N. Geilfus, F. Brabant, B. Delille, M. Vancoppenolle, G. Gilson, J. Zhou, and J.-L. Tison, "Investigations on physical and textural properties of arctic first-year sea ice in the amundsen gulf, canada, november 2007-june 2008 (ipy-cfl system study)," *Journal of Glaciology*, vol. 59, no. 217, p. 819-837, 2013.

[28] M. Shokr and N. K. Sinha, *Sea ice: Physics and Remote Sensing*. John Wiley & Sons, 2015.

[29] D. K. Perovich and A. J. Gow, "A statistical description of the microstructure of young sea ice," *Journal of Geophysical Research: Oceans*, vol. 96, no. C9, pp. 16943-16953, 1991.

[30] M. Uysal, C. Capsoni, Z. Ghassemlooy, A. Boucouvalas, and E. Udvary, *Optical wireless communications: an emerging technology*. Springer, 2016.

[31] A. P. Worby, C. A. Geiger, M. J. Paget, M. L. Van Woert, S. F. Ackley, and T. L. DeLiberty, "Thickness distribution of antarctic sea ice," *Journal of Geophysical Research: Oceans*, vol. 113, no. C5, 2008.

[32] M. R. Drinkwater, "Limex '87 ice surface characteristics: implications for c-band sar backscatter signatures," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 27, no. 5, pp. 501-513, 1989.

[33] A. Manninen, "Surface roughness of baltic sea ice," *Journal of Geophysical Research: Oceans*, vol. 102, no. C1, pp. 1119-1139, 1997.

[34] J. C. Landy, D. Isleifson, A. S. Komarov, and D. G. Barber, "Parameterization of centimeter-scale sea ice surface roughness using terrestrial lidar," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 53, no. 3, pp. 1271-1286, 2015.

[35] J. C. Landy, A. S. Komarov, and D. G. Barber, "Numerical and experimental evaluation of terrestrial lidar for parameterizing centimeter-scale sea ice surface roughness," *IEEE Transactions on Geoscience and Remote Sensing*, vol. 53, no. 9, pp. 4887-4898, 2015.

[36] W. Hibler Iii, "Modeling a variable thickness sea ice cover," *Monthly weather review*, vol. 108, no. 12, pp. 1943-1973, 1980.
[37] R. Pope and E. Fry, "Absorption spectrum (380-700 nm) of pure water. ii. integrating cavity measurements," *Appl. Opt.*, vol. 36, pp. 8710-8723, (1997).

[38] A. Morel, "Light and marine photosynthesis: a spectral model with geochemical and climatological implications," *Progress in oceanography*, vol. 26, no. 3, pp. 263-306, 1991.

[39] Y. I. Kopilevich, M. Kononenko, and E. Zadorozhnaya, "The effect of the forward-scattering index on the characteristics of a light beam in sea water," *Journal of Optical Technology*, vol. 77, no. 10, pp. 598-601, 2010.

[40] W. M. Irvine and J. B. Pollack, "Infrared optical properties of water and ice spheres," *Icarus*, vol. 8, no. 1-3, pp. 324-360, 1968.

[41] A. Morel and S. Maritorena, "Bio-optical properties of oceanic waters: A reappraisal," *Journal of Geophysical Research: Oceans*, vol. 106, no. C4, pp. 7163-7180, 2001.

[42] B. Light, H. Eicken, G. Maykut, and T. Grenfell, "The effect of included participates on the spectral albedo of sea ice," *Journal of Geophysical Research: Oceans*, vol. 103, no. C12, pp. 27 739-27 752, 1998.

[43] A. Bricaud, A. Morel, and L. Prieur, "Absorption by dissolved organic matter of the sea (yellow substance) in the uv and visible domains 1," *Limnology and oceanography*, vol. 26, no. 1, pp. 43-53, 1981.

[44] T. C. Grenfell and D. K. Perovich, "Radiation absorption coefficients of polycrystalline ice from 400-1400 nm," *Journal of Geophysical Research: Oceans*, vol. 86, no. C8, pp. 7447-7450, 1981.

[45] B. Light, G. Maykut, and T. Grenfell, "A temperature-dependent, structural-optical model of first-year sea ice," *Journal of Geophysical Research: Oceans*, vol. 109, no. C6, 2004.

[46] —, "Effects of temperature on the microstructure of first-year arctic sea ice," *Journal of Geophysical Research: Oceans*, vol. 108, no. C2, 2003.

[47] J. Felske, T. Charalampopoulos, and H. Hura, "Determination of the refractive indices of soot particles from the reflectivities of compressed soot pellets," *Combustion Science and Technology*, vol. 37, no. 5-6, pp. 263-283, 1984.

[48] W. J. Wiscombe and S. G. Warren, "A model for the spectral albedo of snow: i: Pure snow," *Journal of the Atmospheric Sciences*, vol. 37, no. 12, pp. 2712-2733, 1980.

[49] I. I. Kim, B. McArthur, and E. J. Korevaar, "Comparison of laser beam propagation at 785 nm and 1550 nm in fog and haze for optical wireless communications," in *Optical Wireless Communications III*, E. J. Korevaar, Ed., vol. 4214, International Society for Optics and Photonics. SPIE, 2001, pp. 26-37. [Online]. Available: https://doi.org/10.1117/12.417512

[50] M. Mellor, "Light scattering and particle aggregation in snow-storms," *Journal of Glaciology*, vol. 6, no. 44, pp. 237-248, 1966.

[51] A. Macke, J. Mueller, and E. Raschke, "Single scattering properties of atmospheric ice crystals," *Journal of Atmospheric Sciences*, vol. 53, no. 19, pp. 2813-2825, 1996.

[52] T. Nousiainen, "Scattering of light by raindrops with single-mode oscillations," *Journal of the atmospheric sciences*, vol. 57, no. 6, pp. 789-802, 2000.

[53] M. Grabner and V. Kvicera, "Multiple scattering in rain and fog on free-space optical links," *Journal of Lightwave Technology*, vol. 32, no. 3, pp. 513-520, 2014.

[54] O. Boucher, "On aerosol direct shortwave forcing and the heney-greenstein phase function," *Journal of the atmospheric sciences*, vol. 55, no. 1, pp. 128-134, 1998.

[55] M. C. Al Naboulsi, H. Sizun, and F. de Fornel, "Fog attenuation prediction for optical and infrared waves," *Optical Engineering*, vol. 43, no. 2, pp. 319-329, 2004.

[56] T. C. Grenfell, "A theoretical model of the optical properties of sea ice in the visible and near infrared," *Journal of Geophysical Research: Oceans*, vol. 88, no. C14, pp. 9723-9735, 1983.

[57] J. L. Henyey and L. G. Greenstein, "Diffuse radiation in the galaxy," in *Astrophysical Journal*, Januray 1841, pp. 70-83.

[58] L. Gate, "Light-scattering cross sections in dense colloidal suspensions of spherical particles," *Journal of OSA*, vol. 63, no. 3, pp. 312-317, 1973.

[59] Y. Cocheril and R. Vauzelle, "A new ray-tracing based wave propagation model including rough surfaces scattering," *Progress In Electromagnetics Research*, vol. 75, pp. 357-381, 2007.

[60] P. Yeh, *Optical Waves in Layered Media*, Hoboken, Ed. NJ, USA: Wiley, 2005.

[61] S. Jacques and L. Wang, *Monte Carlo Modeling of Light Transport in Tissues*. Boston, MA: Springer US, 1995, pp. 73-100.

[62] G. F. Simmons, "Calculus with Analytic Geometry," in *McGraw-Hill*, Oct 1996, pp. 249-519.

[63] H. Zhang and Y. Dong, "General stochastic channel model and performance evaluation for underwater wireless optical links," *IEEE Transactions on Wireless Communications*, vol. 15, no. 2, pp. 1162-1173, 2016.

[64] C. Gabriel, M.-A. Khalighi, S. Bourennane, P. Leon, and V. Rigaud, "Monte-carlo-based channel characterization for underwater optical communication systems," *IEEE/OSA Journal of Optical Communications and Networking*, vol. 5, no. 1, pp. 1-12, 2013.

[65] A. S. Ghazy, S. Hranilovic, and M.-A. Khalighi, "Angular MIMO for underwater wireless optical communications: Link modeling and tracking," *IEEE Journal of Oceanic Engineering*, pp. 1-17, 2021.

[66] H. Zhang, J. Cheng, Z. Wang, and Y. Dong, "On the capacity of buoy-based MIMO systems for underwater optical wireless links with turbulence," in *2018 IEEE International Conference on Communications (ICC)*, 2018, pp. 1-6.

[67] J. M. Kahn, W. J. Krause, and J. B. Carruthers, "Experimental characterization of non-directed indoor infrared channels," *IEEE Transactions on Communications*, vol. 43, no. 2/3/4, pp. 1613-1623, 1995.

[68] G. W. Marsh and J. M. Kahn, "Performance evaluation of experimental 50-mb/s diffuse infrared wireless link using on-off keying with decision-feedback equalization," *IEEE Transactions on Communications*, vol. 44, no. 11, pp. 1496-1504, 1996.

[69] L. J. Johnson, R. J. Green, and M. S. Leeson, "The impact of link orientation in underwater optical wireless communication systems," in *2014 Oceans - St. John's*, 2014, pp. 1-8.

[70] T. Sakurai, H. Chosrowjan, T. Somekawa, M. Fujita, H. Motoyama, O. Watanabe, and Y. Izawa, "Studies of melting ice using co2 laser for ice drilling," *Cold Regions Science and Technology*, vol. 121, pp. 11-15, 2016.

[71] R. T. Mihran, "Interaction of laser radiation with structures of the eye," *IEEE Transactions on Education*, vol. 34, no. 3, pp. 250-259, 1991.

[72] C. Mobley, E. Boss, and C. Roesler, "Ocean Optics Web Book (2016), [Available at] https://www.oceanopticsbook.info, [Accessed on]: 2021-07-06.

[73] J. Kahn and J. Barry, "Wireless infrared communications," *Proceedings of the IEEE*, vol. 85, no. 2, pp. 265-298, 1997.

[74] J. R. Barry and J. M. Kahn, "Link design for nondirected wireless infrared communications," *Appl. Opt.*, vol. 34, no. 19, pp. 3764-3776, Jul 1995.

[75] ThorLab, "Thorlab," [Available at] http://www.thorlabs.com/, [Accessed on]: 2021-04-01.

[76] T. Hamza, M.-A. Khalighi, S. Bourennane, P. Léon, and J. Opderbecke, "Investigation of solar noise impact on the performance of underwater wireless optical communication links," *Opt. Express*, vol. 24, no. 22, pp. 25 832-25 845, Oct 2016. [Online]. Available: http://www.opticsexpress.org/abstract.cfm?URI=oe-24-22-25832

[77] T. Hamza and M. Ali Khalighi, "On limitations of using silicon photomultipliers for underwater wireless optical communications," in *2019 2nd West Asian Colloquium on Optical Wireless Communications (WACOWC)*, 2019, pp. 74-79.

[78] J. Proakis and M. Salehi, *Digital Communications*. McGraw-Hill, 2008.

[79] T. Komine, J. H. Lee, S. Haruyama, and M. Nakagawa, "Adaptive equalization system for visible light wireless communication utilizing multiple white led lighting equipment," *IEEE Transactions on Wireless Communications*, vol. 8, no. 6, pp. 2892-2900, 2009.

[80] R. M. G. S. Karp, "Optical communications," in *John Wiley and Sons*, feburary, First 1995, pp. 249-519.

[81] F. Xu, M.-A. Khalighi, and S. Bourennane, "Impact of different noise sources on the performance of pin- and apd-based fso receivers," in *Proceedings of the 11th International Conference on Telecommunications*, 2011, pp. 211-218.

[82] I. Juszak, W. Eugster, M. M. Heijmans, and G. Schaepman-Strub, "Contrasting radiation and soil heat fluxes in arctic shrub and wet sedge tundra," *Biogeosciences*, vol. 13, no. 13, pp. 4049-4064, 2016.

[83] S. Haykin, *Communication systems*. John Wiley & Sons, 2008.

[84] J. A. Ogilvy and J. R. Foster, "Rough surfaces: gaussian or exponential statistics?" *Journal of Physics D: Applied Physics*, vol. 22, no. 9, pp. 1243-1251, sep 1989.

[85] Zemax, "Opticstudio," [Available at] https://www.zemax.com, [Accessed on]: 2021-04-01.

[86] Mathwork, "Curve Fitting Tool," [Available at] https://www.mathworks.com, [Accessed on]: 2021-04-01.

[87] S. Tang, Y. Dong, and X. Zhang, "Impulse response modeling for underwater wireless optical communication links," *IEEE Transactions on Communications*, vol. 62, no. 1, pp. 226-234, 2014.
[88] Y. Li, M. S. Leeson, and X. Li, “Impulse response modeling for underwater optical wireless channels,” *Appl. Opt.*, vol. 57, no. 17, pp. 4815–4823, Jun 2018.

[89] R. Boluda-Ruiz, P. Rico-Pinazo, B. Castillo-Vázquez, A. García-Zambrana, and K. Qaraqe, “Impulse response modeling of underwater optical scattering channels for wireless communication,” *IEEE Photonics Journal*, vol. 12, no. 4, pp. 1–14, 2020.

<figure><img src="image_18.png" /></figure>

STEVE HRANILOVIC (S'94-M'03-SM'07) received the B.A.Sc. degree with honours in electrical engineering from the University of Waterloo, Canada in 1997 and M.A.Sc. and Ph.D. degrees in electrical engineering from the University of Toronto, Canada in 1999 and 2003 respectively. He is a Professor in the Department of Electrical and Computer Engineering, McMaster University (Hamilton, Ontario, Canada) and currently serves as the Associate Dean (Academic). During 2010-

2011 he spent his research leave as Senior Member, Technical Staff in Advanced Technology for Research in Motion, Waterloo, Canada. His research interests are in the areas of free-space and optical wireless communications, digital communication algorithms, and electronic and photonic implementation of coding and communication algorithms. He is the author of the book *Wireless Optical Communication Systems* (New York:Springer, 2004).

Dr. Hranilovic is a licensed Professional Engineer in the Province of Ontario and was awarded the Government of Ontario Early Researcher Award in 2006. In 2016 the title of University Scholar was conferred upon him by McMaster University. He has served as an Associate Editor for the Journal of Optical Communications and Networking and an Editor for the IEEE Transactions on Communications in the area of Optical Wireless Communications.

<figure><img src="image_19.png" /></figure>

ABDALLAH S. GHAZY was born in Giza, Egypt, in 1983. He received the B.S. degree in electrical engineering from Al-Azher University, Egypt, in 2007. From 2008 to 2016, he has worked in communication systems in Telecome-Egypt and Avaya business partner companies, Cairo, Egypt. In 2012, he joined the electrical engineering school, Azher University, Egypt. He received the M.S degree in electrical engineering from Egypt-Japan University of Science and Technology, Alexandria, Egypt, in 2016. He is currently on leave from the Al-Azher University and is working toward the Ph.D. degree, electrical and computer engineering school, McMaster University, Ontario, Canada. His main research interests include optical communications, wireless communications and digital signal processing.

<figure><img src="image_20.png" /></figure>

HAITHAM S. KHALLAF (SM' 21) received the B.Sc. degree (Hons.) from the Faculty of Electronic Engineering, Menoufia University, Menouf, Egypt, in 2005, and the M.Sc. (Hons.) and Ph.D. degrees from the Egypt-Japan University of Science and Technology (E-JUST), Alexandria, Egypt, in 2013 and 2016, respectively. He is currently a Postdoctoral Researcher at the Department of Electrical and Computer Engineering, McMaster University (Hamilton, Ontario, Canada), on

leave from his position as an Assistant Professor at the Nuclear Research Center (NRC), Egyptian Atomic Energy Authority (EAEA), Egypt. From 2016 to 2020, he was a Postdoctoral Researcher with Zewail City of Science and Technology (Egypt), Ozyegin University (Istanbul, Turkey), and Kyushu University (Japan). His research interests include free space optical communication channel modeling, unmanned aerial vehicle network backhauling, satellite communications, and application of wireless communication systems in nuclear facilities.

MOHAMMAD-ALI KHALIGHI (Senior Member, IEEE) is Associate Professor with École Centrale Marseille, Marseille, France, and head of “Optical Communications for IoT” group at Fresnel Institute research lab. He is currently serving as Project Coordinator for the H2020 ITN MSCA VisIoN project (Visible-light-based Interoperability and Networking) and Action Chair for the COST Action CA19111 NEWFOCUS (European Network on Future Generation Optical Wireless

<figure><img src="image_21.png" /></figure>

Communication Technologies). He has co-edited the book “Visible Light Communications: Theory and Applications” (CRC Press, 2017) and was the co-recipient of the 2019 Best Survey Paper Award of the IEEE Communications Society. He is also serving as Editor-at-Large for the IEEE Transactions on Communications, and served as Associate Editor for the IET Electronics Letters as well as Guest Editor for the IEEE Open Journal of the Communications Society and Elsevier Optik journal. His main research interests include signal processing for wireless communication systems with an emphasis on the physical layer aspects of free-space, underwater, and indoor visible-light optical communications.

<figure><img src="image_22.png" /></figure>

1471