Under-Sea Ice Diffusing Optical Communications

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

II, we introduce the SDOC approach and the SSCL chan-94 nel model. In Section III, we use the MCNRT method to 95 model the upward transmission, then derive a quasi-analytic 96 equation for the CIR. We introduce and model the proposed 97 transceiver architecture in Section IV. In Section V, we nu-98 merically investigate the channel characteristics and system 99 performance. Finally, conclusions are given in Section VI. In this section we introduce the SDOC link as a new approach 103 to establish communication between AUVs operating under 104 sea ice. We discuss the temperature and salinity profile of the 105 sea ice. Then, we introduce a new approach to model optical 106 characteristics of the sea ice. 107 A. SDOC ARCHITECTURE 108 As shown in Fig. 1, we consider a group of AUVs, for exam-109 ple an AUVs-swarm 1 , navigating several meters beneath a 110 sheet of sea ice. The AUVs move together in the coordinated 111 fashion with a separation of a few meters. In the proposed 112 approach a broadcast communication link between the AUV 113 transmitter (AUV-Tx) and the AUV receivers (AUV-Rxs) is 114 accomplished in two steps: upward and downward trans-115 missions. In the upward transmission, the AUV-Tx sends a 116 narrow collimated laser beam toward the sea ice. Due to im-117 purities (particles) 2 , the transmitted beam is subject to intense 118 scattering at the surface and during propagation in the interior 119 of the ice sheet. Inside the sea ice, a portion of the power will 120 be transmitted through the sheet and lost to the atmosphere. 121 Alternatively, the transmitted light may be trapped in the 122 interior of the sheet where it is absorbed. Finally, a portion 123 of the incident light will be diffused back from the ice sheet 124 into the water. This diffused light which escapes the ice sheet 125 is the useful signal which is used to establish the broadcast 126 communication link. Given that the light is diffused inside the 127 sheet, as shown in the green ellipse in Fig. 1, a wide coverage 128 area is possible. The AUV-Tx can control the position of the 129 diffusing spot by adjusting the direction of the laser beam, 130 i.e., polar and azimuthal launching angles. For instance, if the 131 AUV-Rxs are distributed symmetrically around the AUV-Tx, 132 the beam should be vertically oriented toward the ice sheet 133 to offer a fair coverage for all AUV-Rxs, as shown for the 134 case in Fig. 1. However, if the AUV-Rxs are biased to one 135 side, the AUV-Tx can orient its beam toward the direction 136 of the AUV-Rxs to improve link quality. In the downward 137 transmission from the ice sheet, the diffused beam propagates 138 in the seawater and covers the AUV-Rxs with a large spot. 139 Regardless of the position and orientation of the AUV-Rxs, 140 each AUV-Rx receives a portion of this diffused beam, and 141 the AUV-Tx establishes a broadcast communication with the 142 AUV-Rxs. (1) S(z) = −3.24 × 10 −7 z 6 + 3.58 × 10 −5 z 5 − 1.47 × 10 −3 × z 4 + 2.74 × 10 −2 z 3 − 0.205 z 2 − 0.095 z + 13.63 , where T is the temperature in Celsius ( o C), S is the 167 salinity in part per thousand (ppt), and 0 ≤ z ≤ 36 cm. 168 The equations are shown in the figure, and there is good 169 agreement between the measured and the fitted profiles 4 . 170 Another example is a 12 cm bare-sea ice sheet whose 171 temperature and salinity profiles are shown in [29,Fig. 3]. 172 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]. 4 The corresponding goodness of the fit criteria are; R-square= {0.9916, 0.9931} for the temperature and salinity curves, respectively. Measured Salinity FIGURE 2: The temperature and salinity profiles versus the sea ice depth for a snow-covered sea ice sheet as measured by [27].
The sheet is young laboratory-grown saline sea ice. The two 173 profiles of the sheet can be well fitted in T (z) and S(z) 174 functions as 5 175 T (z) = 1.176 z − 15.61 , 0 ≤ z ≤ 12cm (3) S(z) = 0.05003 z 2 − 0.7432 z + 8.203 . (4) These two ice sheet examples will be used later in the 176 numerical results as case studies. 177 As shown in Fig. 2, the top surface of the sea ice is lower 178 in the temperature than the bottom due to a cooling of the 179 atmosphere and a warming of the seawater. As well, the 180 salinity at the top and bottom is much higher than at the 181 middle of the sea ice sheet. The vertical variations in the 182 temperature and salinity with the thickness of the ice sheet 183 result in changes in particle densities, which impact the chan-184 nel optical characteristics. Given that the scattering inside the 185 ice sheet is extensive and varies through the thickness of the 186 ice sheet, channel modeling is challenging. In the following 187 we introduce a simplified channel model. and soot). However, snow layers are composed of air as a 222 hosting medium with a fewer numbers of mixture particles 223 (e.g., snow grains, algal and non-algal particles and soot). 224 Due to these particles, the optical ray propagating inside the 225 m th layer of the SSCL model suffers from absorption and 226 scattering effects. The absorption coefficient, a(m), is the 227 weighted summation of the contribution from the mixture 228 components as [56] where a o and f vo are the absorption coefficient and the 230 volume fraction of the hosting medium, respectively. As 231 well, a j and f vj are the absorption coefficient and the vol-232 ume fraction associated with the j th particle, respectively, 233 where f vo + Jm j=1 f vj = 1. Symbol J m is the number of 234 mixture particles in layer m, and the value of J m depends 235 on the hosting medium of the layers and its surrounding 236 environment. The hosting medium does not contribute to the 237 scattering effect, thus, the scattering coefficients for each 238 layer, b(m), are weighted summations of the contribution 239 from the impurity components only as [56] where b j is the scattering coefficient associated with the j th 241 particle.

242
Based on the assumptions given in [26] and [54], the 243 one term Henyey-Greensteen (OTHG) function is a good 244 approximation to the phase scattering function [57] 245 p θs (θ s , m) = 1 4 π 1 − g(m) 2 (1 + g(m) 2 − 2 g(m) cos(θ s )) 3/2 , (7) where g(m) is the asymmetry factor and θ s is a scattering 246 angle. The asymmetry factor is obtained using the weighted 247 sum as [56] 248 249 where g j is the asymmetry factor of the j th particle. The 250 effective refractive index of the layer is computed using the 251 volume fraction f vj as [58] where n o is the refractive index of the hosting medium, and 253 n j is the refractive index of the j th particle.

254
The interfaces between the adjacent layers are assumed to 255 be rough surfaces which leads to optical surface scattering at 256  [26], [25], [26] [45, Fig. 7 Fig. 7] Air bubble ≈ Zero [26], [25], [20] [45, Fig. 7  the entrance of each layer. The surface roughness of the inter-257 face is presented with the random height in the z direction for 258 each point (x, y), which can be well described in the x and 259 y directions using the two-dimensional Gaussian distribution 260 as measured in [32], [33]. To generate a realization of the 261 ice surface, a two dimensional Gaussian random variable is 262 generated with independent components in x and y according 263 to [59] 264 where z is the height at (x, y) point, and σ x (m) and σ y (m) 265 are the RMS values in x and y directions 7 , respectively. As 266 measured in [32], [33], the correlation between heights over 267 the surface is well approximated using the two-dimensional 268 generalized power-law function. Thus, to represent the corre-269 lation in space of the surface, the Gaussian realization can be 270 filtered by a generalized power-law function. This function is 271 given with one dimension in [34] and can be generalized to 272 7 The experimental measurements in the Arctic and Antarctic regions revealed that the roughness parameters, RMS and correlation length, are in the millimetre and the centimetre ranges, respectively [34], [35]. two dimensions p ρm (ρ x , ρ y ) as where ρ x and ρ y are the distances between correlated points 274 in x and y directions, respectively, l x (m) and l y (m) are the 275 correlation lengths in x and y directions, respectively. The 276 value of ξ depends on the geographical location of the sea 277 ice sheet, and is equal to 1 and 2 in cases of exponential-278 correlated and Gaussian-correlated surfaces, respectively. 279 Note that, the surface roughness includes parts of the ice 280 suspended in seawater. Due to the low density of these parts, 281 they typically float up toward the ice sheet and settle on its 282 bottom surface [1].

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

290
In this section, we obtain an expression for the CIR of 291 links between the AUV-Tx and the AUV-Rxs considering 292 the effects of scattering, attenuation, as well as AUV-Rxs 293 VOLUME 4, 2016   While the AUV-Rx is located at (∆ x , ∆ y , ∆ z ) position with 314 aperture polar and azimuthal inclination angles (θ in ,φ in ). 315 Thus, the AUV-Rx position can be described using the 316 position and orientation (PO) vector (5 × 1) as ∆ r := 317 [∆ x ; ∆ y ; ∆ z ; θ in ; φ in ]. The AUV-Rx is equipped with a lens 318 with diameter D r and field of view (FOV) of θ F OV .
and the random distance is generated as [61] where u µ is a uniform random variable, u µ ∼ U [0, 1], and 379 c(m) is the extinction coefficient of the m th layer represent-380 ing the loss in the power of the ray. The value of the extinction 381 coefficient c(m) is computed as When a scattering event occurs, the weight of the photon 383 packet is dropped to [61] Upon scattering, the optical ray arriving from the direction e 1 386 will have a new direction e 3 determined randomly according 387 to polar and azimuthal scattering angles (θ us , φ us ). The angle 388 θ us (m) is generated from the OTHG PDF in Eq. (7) as [61] where u θ ∼ U [0, 1]. Also, the azimuthal scattering angle φ us 390 is typically described by a uniform PDF, and it is generated 391 as [57] 392 393 p φs (φ us ) = 1 2 π , φ us = 2π u φ (19) where u φ ∼ U [0, 1]. After scattering, the ray travels a new 394 distance µ u1 with a new direction e 3 before the next scat-395 tering occurs with likelihood p µ (µ u1 , m). Compared to the 396 seawater and the atmosphere, particle scattering takes place 397 much more frequently in snow and sea ice layers. Typically, 398 the optical ray is scattered few times in the seawater or 399 atmosphere layer, however, hundreds of scattering events can 400 typically take place in the sea ice or snow layers.

401
The MCNRT traces the optical rays until they are either 402 absorbed, trapped in the ice layer, escape to the atmosphere, 403 or diffuse back into the seawater. The diffused rays only 404 contribute in the obtained diffusing pattern for the upward 405 transmission and the remainder of the rays are considered 406 as lost. For a given position and orientation for the AUV-407 Tx, ∆ t , the normalized diffusing pattern is obtained with the 408 intensity I d as a function of the space, angles and time as 409 follows 20) where, as shown in Fig. 4, the intensity I d is measured on the 411 bottom of the sea ice surface at position x d and y d , with polar 412 VOLUME 4, 2016 θ d , azimuth φ d angles, and time t d . As well, the DC gain of 413 the upward transmission G u is computed using I d as where where rival angles (θ r ,φ r ) measured relative to the sea ice axes, For given scattering angles (θ ds , φ ds ), the Then, 445 e d is rotated around (Y s , X s , Z s ) axes by two an-446 gles: θ y = arcsin (cos(φ ds ) sin(θ ds )) and θ x = 447 arcsin (sin(φ ds ) sin(θ ds )/ cos(θ y )) respectively. Thus, θ r 448 and φ r are computed as where R x (θ x ) and R y (θ y ) are (3 × 3) rotation matrices Vector e r is also characterized by arrivals angles (θ r ,φ r ) 453 measured relative to the axes, (X r , Y r , Z r ), as shown in the 454 Fig. 4, and can be equivalently written as 455 e r = x r sin(θ r ) cos(φ r ) + y r sin(θ r ) sin(φ r ) + z r cos(θ r ) , where ( x r , y r , z r ) are the unit vectors relative to the Rx 456 axes (X r , Y r , Z r ). For the given angles (θ r , φ r ), the an-457 gles (θ r ,φ r ) are calculated from Eq. (23) by replacing 458 e d with e r =[sin(θ r ) cos(φ r ); sin(θ r ) sin(φ r ); cos(θ r )] and 459 substituting θ y = arcsin (cos(φ in ) sin(θ in )) and θ x = 460 arcsin (sin(φ in ) sin(θ in )/ cos(θ y )). The scattered ray ar-461 rives at arrival position (x s r , y s r , z s r ) over the aperture of the 462 AUV-Rx. 463 The arriving ray from the LOS or scattering path is de-464 tected if the position of arrival (x r , y r , z r ) is located on the 465 lens of the AUV-Rx with arrival angles (θ r ,φ r ) less than half 466 angle of the FOV. This can be compactly represented as the 467 geometric loss G g and it is written as where f p (D r , ∆ r ) represents the spatial extent of the AUV-469 Rx lens with respect to the sea ice axes (X d , Y d , Z d ).
For rays diffused from a single point on the bottom of the sea 481 ice (x d , y d , 0), the CIR can be well approximated by a linear 482 combination of LOS components as where the length of the LOS path is computed geometrically 484 from the figure as The symbols t r and ν are the arrival time and the light 486 speed in the seawater, respectively, and δ(.) is the Dirac-delta 487 function.  Figure 4 shows the diffused ray traveling in the direction e d 495 for a distance µ do then is scattered in the direction e r and 496 travel a distance µ d1 before arriving the lens. The scattering 497 position (x s , y s , z s ) and angle θ ds are given by [66] 498 This scattering results in a reduction in the photon packet 499 weight of the ray e r by a factor of b/c relative to the packet 500 of the ray e d . After scattering and traveling a distance µ d1 , 501 the ray arrives to a position (x s r , y s r , z s r ) which is obtained as Using Eqs. (22) is given as where the length of the single scattering path is computed as 509 l s r = µ do + µ d1 , and µ d1 is computed using Eqs. (30) and The overall CIR is the summation of the LOS and scattering 512 components, and it is computed using Eqs. (28) and (31) as The CIR for the link between the AUV-Tx and an AUV-Rx with PO vector ∆ r is computed by integration over all the points on the bottom of the sea ice (x d , y d ) as Equation (34) can be used to determine the link budget and 514 the induced pulse dispersion. The DC gain of a downward 515 transmission (i.e., AUV-Tx to an AUV-Rx link) is obtained 516 from CIR as [67] 517 where P o is the transmitted power as defined in the link 518 model. As well, RMS of the pulse spreading is computed as 519 where, τ o is the mean excess delay given by [67] 521 The system of Equations, (22)-(37), are used to quantify 522 the link performance between the AUV-Tx and the AUV-Rxs 523 as shown in Section V.

525
Though the proposed SDOC approach provides a broadcast 526 communication link without requirement for alignment, its 527 performance is limited by the high channel attenuation and 528 inter-symbol interference (ISI) due to multipath propagation. 529 The ISI is induced mainly by the sea ice sheet in the upward 530 transmission, but also, in the downward transmission due 531 to the scattering occurring in the seawater. In addition, the 532 performance can be degraded by background radiations due 533 to the fact that the AUVs navigate near the bottom of the 534 sea ice and the orientation of the receivers are aligned up-535 wards, as shown in Fig. 4. In this section, inspired by indoor 536 OWC systems [12]- [16], we propose appropriate Tx and Rx 537 architectures to tackle these limitations. This communication 538 architecture can be considered as a first prototype step in 539 the development of such links. We also discuss practical 540 implementation considerations of SDOC links.  Figure 5 shows the overall block diagram of the proposed 543 SDOC system, as described in the following. The proposed architecture is shown in Fig. 5a. For sim-547 plicity, the transmitted data are encoded using intensity 548 modulation direct detection (IM/DD) with non-return-to-zero 549 OOK (NRZ-OOK) modulation scheme [68]. As well, for 550 simplicity, we consider the LD to be switched fully on and 551 off corresponding to ones and zeros of the OOK symbols, 552 respectively, i.e., zero extinction ratio. The OOK symbol 553 duration is T b , the transmitted data rate is  concentrator is quantified by its gain G c which depends on 578 its refractive index, n c , and the FOV as [73] 579 whereθ b is the incident angle of the background ray upon the 580 concentrator and it is measured relative to the optical axis of 581 the Rx, Z r , as shown in Fig. 4. As well, the optical band pass . Though direct human contact with UAVs is possible, safety must also be considered to preserve wildlife which may interact with these optical emissions [71]. depends on the incident angle of the received ray. Such hemi-584 spheric concentrators are commercially available and have 585 been used in optical diffusing communication systems for 586 indoor applications 10 . The concentrator enlarges the effective 587 area of the PD, A ef , which means capturing solar noise. The 588 effective area of the PD is obtained as [73] 589 where A P D denotes the physical active area of the PD. Here, 590 for simplicity, the dependence of the effective area on the 591 incident angleθ b is represented by replacing A ef (θ b ) by 592 its average A ef over the incident angle, while making two 593 assumptions. Firstly, we assume that the function T (θ b ) can 594 be replaced by its average, T , over all incident angles. This 595 assumption holds, especially, when the incident optical ray 596 arrives within a wide range of the angles which is the typical 597 case of diffusing communications [73]. Secondly, we assume 598 a uniform PDF forθ b . Then, the average Rx effective area is 599 obtained as Note that, enlarging the FOV decreases the the average 601 effective area of the Rx.

602
After the hemispherical concentrator, a silicon PIN pho-603 todiode (PIN-PD) with a trans-impedance amplifier (TIA) is 604 used. The PIN-PD converts the collected optical rays to an 605 electrical current proportionally to its responsitivity and 606 A P D . Then, the TIA converts the small current to a high 607 voltage proportionally to its load resistance R L . In contrast 608 to avalanche photodiodes, photo-multiplier tube and SiPM 609 PDs, the silicon PIN-PD achieves a better performance when 610 the background radiation is much high and dominants the 611 receiver noises [76], [77]. 612 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]. sampling rate T s , where T s < T b /2 to avoid aliasing [78].

622
The sampled signal is then processed by a discrete-time chan- of the weighted inputs as follows [79] 638 where α j F and β j B are the FF and FB weighting coefficients, where, v n is the sampled noise voltage, and ∆ r is PO 660 vector of the AUV-Rx as defined in the previous section, 661 P DF E (k T s ) is the sampled system impulse response of the 662 DFE, ⊗ is the discrete convolution operator, * is continuous 663 convolution and p(t) is the instantaneous transmitted optical 664 optical power. The signal in Eq. (42) can be decomposed as 665 the sum of the desired signal, denoted by v s (∆ r ), and the 666 residual ISI denoted by v isi (∆ r ) where, where τ d is the time delay of the channel and T b and τ d are 668 assumed to be multiples of T s . Also, The noise contribution in (42) includes the effects of the 670 thermal v th and shot v sh noises, i.e., v n = v th + v sh . The 671 thermal noise v th is well described by zero mean Gaussian 672 distribution with variance σ 2 th given as [80], where, K is the Boltzmann constant and T (m = 1) is the 674 temperature of the seawater layer in Celsius as defined in Sec. 675 II-B. Usually, the temperature of the seawater underneath sea 676 ice is T (m = 1) ≤ 0 o C, as shown in Fig. 2. On the other 677 hand, the shot noise is associated with the superposition of 678 the desired signal voltage v s , the ISI distortion voltage v isi , 679 and the background radiation voltage v sun . Due to the high 680 intensity of the solar radiation, the shot noise can be modeled 681 using Gaussian random process with variance given as [80] 682   The mean η Λ and the variance σ 2 Λ , Λ ∈ {0, 1}, of the total 699 signal and noise affects system performance are given as isi is the variance of ISI signals and it is equal to where coastal seawaters, respectively [65]. Note that in the follow- Co-S channels) and the clear and coastal seawater with a bare 734 sea ice sheet (Cl-B and Co-B channels). The snow-covered 735 sea ice sheet has a thickness of 36 cm and it well described by 736 Eqs. (1), (2). The bare sea ice sheet has thickness 12 cm and 737 its temperature and salinity profiles are described by Eqs.

741
The bare-ice cases are divided into 6 layers while the 742 snow-covered cases are divided into 9 layers 11 . In all cases, 743 each layer is assigned with the average values of the salinity 744 and temperature using Eqs. (1)-(4), as shown in Tables 2a 745 and 2b. Clear weather above the sea ice sheets is assumed, 746 which is the typical case during sunny days. As shown in 747 Tables 2a and 2b, the scattering coefficients of the snow-748 covered sea ice sheet and coastal seawater are higher than 749 that for bare sea ice sheet and clear seawater, respectively. 750 In addition, it is clear that the changes in the refractive 751 indices and asymmetry parameters are small. In Table 2c, the 752 RMS of the roughness and correlation length, are assumed in 753 millimetre and centimetre ranges, respectively, as measured 754 in [34], [35]. As well, we assume isotropic layers (i.e., 755 σ x (m) = σ y (m) and l x (m) = l y (m)), and the interfaces are 756 Gaussian-correlated (i.e., ξ = 2) [32], [33]. The interfaces 757 between the ice layers are assumed smooth due to fact that 758 the variation in the effective refractive indices are negligible 759 in the presented cases. To ensure an accurate realization for 760 the SSCL model, the roughness is sampled with intervals 761 and lengths with values δ x (m) = δ y (m) = 0.1 l x (m) and 762 L x (m) = L y (m) = 60 l x (m) [84].  Figure 6 shows the marginalized diffusing patterns 12 for the 765 Cl-B and Cl-S channels with the orange and maroon colors, 766 respectively. The diffusing pattern is measured at the bottom 767 of the sea ice, i.e., ∆ z = 0, with DC gains of G u = 0.26 and 768 0.37 for Cl-B and Cl-S channels, respectively. These results 769 were obtained by running the MCNRT using the ZeMax 770 Opticstudio software [85] over 10 6 iterations. Note that we 771 have verified that increasing the number of iterations to 10 7 772 resulted in almost identical results. 773 Figures 6a and 6b show the marginalized diffusing patterns 774 versus the polar and the azimuthal angles, I d↓θ and I d↓φ , 775 respectively. As shown in these figures, the marginalized 776 intensity is uniform with respect to (w.r.t.) φ d , however, it 777 is oriented w.r.t. θ d with a peak at θ d ≈ 45 o . The ori-778 entation indicates non-specular diffusing due to the dense 779 scattering occurred in the sea ice and snow. The value of 780 45 o is interrupted as follows; each diffusing point on the sea 781 ice is an identical random variable described by Eq. (18), 782 and the diffusing pattern is a summation of that diffusing 783 points. Assuming the central limit theory, I d↓θ approaches 784  (c) The roughness parameters for the interfaces between the layers of the SSCL channel models [34], [35].
Though, the diffusing pattern has a small spot on the bottom 796 of the sea ice sheet (i.e., |x d | and |y d | ≤ 0.5 m), due to 797 the orientation with angle 45 o , the spot expands out with the 798 propagation in the seawater as shown in the next subsection. 799 Figure 6d shows the marginalized diffusing pattern I d↓t 800 (i.e., temporal dispersion patterns of the upward transmis-801 sion) with t d = [2, 24] ns. The pattern of the Cl-S channel 802 has a high peak with amplitude 14 × 10 −3 and it decays 803 slowly with a long dispersion time due to the thickness and 804 much particle scattering occurred for the laser beam in the 805 channel as can be seen from Table 2b (i.e., a larger thickness, 806 and higher temperature and salinity values). In contrast to 807 the Cl-S channel, the pattern of the Cl-B channel has two 808 peaks with amplitudes 32 × 10 −4 and 26 × 10 Equations (51)-(53) serve as a guide for a future analytic the fitting is challenging due to the dense scattering taken 827 place in the channel. Thus, further investigation is required 828 to obtain more accurate equation as a future work.  to the collection of more rays and improves the DC gain.

892
However, the rate of change in the DC gain with the FOV 893 16 Co-P SSCL channel is the coastal seawater cascaded with a freeimpurity 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. 17 The pure seawater rarely exists underneath the frozen oceans, and it is considered here for comparison.  Figure 11 shows the spatial distributions of the DC channel 967 gain and the RMS delay spread versus the position of the 968 AUV-Rx in the x-y plane. The results are shown for Co-S 969 channel within the area of 6 × 6 m 2 . As well, As shown in Fig. 11a, the DC gain distribution is symmet-977 ric in the x-y plane around the center (∆ x = 0, ∆ y = 0) and 978 the DC gain value decreases monotonically with ∆ x and ∆ y . 979 The shown distribution matches with the average response 980 from the results in Figs. 6b and 6c. As well, as given in the 981    performance degrades with distance and improves by de-1029 creasing the bandwidth of the optical filter. As well, the BER 1030 performance in the case of the Co-S channel is better than 1031 Co-B channel for two reasons. Firstly, the Co-S channel has 1032 a higher upward transmission DC gain; secondly, the Co-S 1033 channel reduces impact of the solar radiations much more 1034 than the Co-B channel. For example, considering a BER 1035 threshold of 10 −3 as indicated by the green line in the figure,   the same NRDS, Rx-E and Rx-UE require NOPP= 2 dB and 1065 NOPP= 2.5 dB, respectively. At higher data rates of NRDS= 1066 0.2, Rx-E and Rx-UE require NOPP= 6.3 dB and NOPP= 3.5 1067 dB, respectively, for the Co-B channel. For the Co-S channel 1068 at the same NRDS, Rx-E and Rx-UE require NOPP= 3.2 dB 1069 and NOPP= 5.8 dB, respectively. These results indicate that 1070 the equalizer improves the power efficiency of the systems by 1071 nearly 3 dB, which means the required transmitted power is 1072 reduced roughly by a factor of two. In other words, the AUV 1073 with the equalized system enhances the power-efficiency of 1074 the AUVs which means more lifetime for the battery.
1075 Figure 14 shows the maximum achievable bit rate under 1076 the constraint BER≤ 10 −3 versus the distance ∆ x with 1077 average transmitted optical power P o ∈ {100, 200} mW. 1078 As shown in the figure, the maximum achievable bit rate 1079 (R b ≈ 700 Mpbs) is achieved directly under the diffusing 1080 surface (∆ x ≤ 1 m). However, as ∆ x increases, the maxi-1081 mum achievable bit rate decreases; the proposed system can 1082 achieve broadcast data rates on the order of R b = 1 Mpbs 1083 over communication ranges of ∆ x = 6 m. As indicated by 1084 the green dashed line, to maintain a communication rate of 1085 10 Mbps, scaling the transmitted power by 2 increases the 1086 communication range by 18% and 10% in cases of Co-B and 1087 Co-S channels, respectively. This trade off between data rate 1088 and coverage distance should be considered during planning 1089 stage of the AUV swarms, based on the required data rate and 1090 range.

1092
In this paper, for first time, we propose a broadband-1093 broadcast approach suitable for networking AUVs under sea 1094 ice, albeit with limited range. We taek advantage of existing 1095 ice sheets on the sea surface to establish a diffusing commu-1096 nication systems. The SSCL model was introduced in which 1097 the channel is represented in the form of cascaded layers 1098 with uniform optical characteristics. Due to the challenge of 1099 analytic modeling of optical signal scattering inside the ice 1100 sheet, MCNRT is used to evaluate the diffusing pattern of 1101 was derived in the form of a quasi-analytic equation assum-1103 ing single scattering light propagation. Due to the expected effects of ISI and relatively high background solar power 1105 noise, we propose a new transceiver architecture that helps 1106 in mitigating the effects of these factors. We also provide 1107 extensive numerical results to investigate the effects of water and ice types, Rx parameters i.e., FOV and optical filter 1109 bandwidth, and the Rx location on the system performance.