# Time Evolution of Underwater Sensor Networks Coverage and Connectivity Using Physically Based Mobility Model

Authors: Fatma Bouabdallah
Journal: Wireless Communications and Mobile Computing
 (2019)

---

## Abstract

Underwater mobile acoustic sensor networks (UW-ASNs) require the design of new networking protocols due to fundamental differences with terrestrial wireless sensor networks. The performance of these protocols is highly impacted by the mobility of sensors, especially when they are freely floating. In such mobile UW-ASNs, nodes move with the water currents but are constrained by the gravitational weight of the sensor along with the water resistance and the buoyant force. A realistic mobility model that can reflect the physical movement of randomly scattered and freely floating sensor nodes under ocean currents provides clearer understanding of the communication challenges and hence helps conceiving efficient communication protocols. In this paper, we first propose an exhaustive physically inspired mobility model which meticulously captures the dynamics of underwater environments. We, then, study the resulting time evolution of network coverage and connectivity. Our objective is to provide the underwater network research community with a realistic mobility model that could be exploited in conceiving networking communication protocols such as routing, localization, and medium access. Namely, we show that the network mobility effect on coverage and connectivity is more significant in intermediately dense UW-ASNs. Less effect is recorded on the coverage and connectivity for low- and high-density UW-ASNs.


---

## Body

## 1. Introduction

Underwater acoustic sensor networks (UW-ASNs) are gaining an increasing interest within the research community since they are the enabling technology for a broad range of applications. Indeed, UW-ASNs brought many applications to life, such as; deep offshore exploration, tsunami warning, and especially oil field detection [1]. Consequently, designing networking protocols is of paramount importance to enable the deployment and the proper functioning of UW-ASNs. One of the most important challenges that faces the design of networking protocols is the underwater environment mobility characteristics. Thus conceiving an underwater mobility model that takes into consideration most of the physical aspects impacting the mobility of an underwater sensor is crucial for the underwater network community as it provides clearer understanding of the underwater environment mobility challenges and hence helps conceiving efficient communication protocols.Earlier works on UW-ASNs tried to model and address the harsh challenges imposed by the underwater channel [2–7] such as high attenuation, limited bandwidth, and high and variable propagation time. Then, researches have started work on networking solutions such as synchronization [8], localization [9, 10], routing protocols [4, 11, 12], energy efficiency, and MAC issues [13–16], but mainly under limited mobility constraint. Note that the performance of all these networking protocols can dramatically decrease with a realistic mobility model that can capture the physical movement of the sensor nodes with ocean currents. Consequently, understanding and considering the mobility features of a freely floating underwater sensor are needed to design efficient networking protocols for underwater acoustic sensor networks.In this paper, we study underwater mobile acoustic sensor networks that consist of free-floating sensors with networking capability. We present an exhaustive and more realistic physically inspired mobility model for underwater environments. Indeed, as opposed to previous works on mobility [10, 17–19], which simply suppose that the sensor mobility follows exactly the current mobility, in our model we consider all the physical forces applied to a freely floating underwater sensor. This results in a different mobility pattern.In more detail, we consider freely floating sensor nodes that are randomly scattered from a boat at the water surface to perform measures of an oceanic environment. Thus, we consider a field of study that covers several kilometers. In such mobile UW-ASNs, nodes move with the force of the water currents but are constrained by the gravitational weight of the sensor along with the water resistance and the buoyant force. In the literature, some works have modeled the mobility of freely floating sensor nodes subject to the underwater environment challenges. One of the earliest works modeling underwater mobility was proposed in [20]. The authors study the trajectory of a sinking underwater sensor under the influence of water current where they assume a constant current velocity. However, the real physical problem is more complex. First, the water density depends on the depth and, second, the water current velocities depend on the sea-shore environment. A widely used mobility model was proposed in [18]. The authors succeed to mathematically model the ocean current and hence the sensors mobility is deduced by simply assuming that the sensor mobility model is the one of the water current. This neglects any physical force applied to a freely floating underwater sensor that may modify its mobility model especially compared to the current mobility model. Indeed, the water current should be considered as a force, among others, acting on the mobility of an underwater sensor rather than defining it. Moreover, due to inertia effects, there is a delay in the sensor response relatively to the current force. The sensor does not response instantaneously to it. This delay increases as the variation of the current velocity becomes important. Note that there has been a lot of localization work based on [18], such as [9, 21–23]. In [19], based on [18], the authors further introduce a ring like mobility pattern that captures the behavior of some ocean currents like the ones observed in Lake Pontchartrain [24]. Here again, the authors assume that the underwater body will adopt the same mobility pattern as the water currents and hence ignore the underwater sensor physical forces that impact the mobility of the underwater object. Authors in [17] and [10] used the kinematic model proposed in [25] to develop a mobility aware routing algorithm and a prediction based localization technique, respectively.In this paper, we are interested in deriving the trajectory of sinking freely floating underwater sensors that were initially deployed at the water surface. In such mobile UW-ASNs, nodes’ trajectories depend not only on the surface winds and subsurface currents as claimed by [18], but also on the gravitational weight of the sensor along with the water resistance and the buoyant force. This mobility response accounts for the sensor inertia, which becomes significant with nonconstant current velocity. Note that, depending on these forces’ values, each underwater sensor may sink to the ocean bottom or to a given depth resulting thus in a time variant 3D deployment. Once our mobility model is set and studied, we investigate its impact on the dynamic network connectivity and coverage. Note that the time variant network connectivity and coverage are considered as the most important performance criteria as they affect the network reliability and lifetime. Moreover, the effectiveness of any underwater communication protocol either routing or medium access is directly related to these two issues. Thus studying the dynamics of these two performance metrics under a realistic mobility model is of paramount importance in UW-ASNs.The paper is organized as follows: in Section2, we define our physically inspired mobility model. In Section 3, we present the network model, as well as the simulation outcomes of the temporal impact of the mobility model on the network connectivity and coverage. Section 4 concludes the paper and presents the main objective of our work.


## 2. Mobility Model

### 2.1. Kinetics of Underwater Sensors

We are interested in predicting the trajectory of a sinking underwater sensor. The positionP of the sinking object is given by the three time-dependent coordinates (x,y,z). The velocity and the acceleration of the sensor can be written as v→=x˙i→+y˙j→+z˙k→ and a→=x¨i→+y¨j→+z¨k→, respectively. Four forces are acting upon the sensor as introduced in [20]. The first force is the weight F→W, which reads(1)F→W=ρsVsg→whereρs is the sensor’s density, Vs is the senor’s volume, and g→ is the terrestrial gravitational acceleration. The weight accelerates the sensor downwards. It is the main sinking force. It is countered by the buoyant force F→B, which is equal to the weight of the displaced water. It is written as(2)F→B=-ρwVsg→ where ρw is the water density, which can depend on the water depth z, i.e., ρwz. If the sensor density is higher than the water density, the weight overtakes the buoyant force and the sensor sinks down. However, the water density is mostly increasing with the water depth. Thus, the buoyant force increases as the sensor sinks while the sensor’s weight is constant. At a certain depth the water density can equal the sensor density and, consequently, the buoyant force balances the weight. Hence, this water depth will be the asymptotic value of the sensor’s z-coordinate.Unlike the weight and the buoyant forces, the third and fourth forces depend on the sensor’s kinematics. In other words, the weight and buoyant forces mainly depend on the water and sensor densities and the sensor’s shape. However, the third and fourth forces depend on the sensor and water velocities. The third force is the one applied by the water current to drive the sensor in the same direction as the current. It writes(3)F→C=CσACv→C-v→// where C is a constant, σ is a shape factor, AC is the sensor’s cross-section area facing the current, v→C is the velocity of the water current, and v→// is the projection of the sensor’s velocity in the same plane as the current velocity. More precisely, any upwelling or downwelling effects are neglected as justified in [20]. Thus, no current is considered in the vertical direction. The water current is simplified to v→C=vCxi→+vCyj→. Consequently, v→//=x˙i→+y˙j→. The current’s force is the only force that governs the horizontal sensor’s motion. It vanishes as the sensor’s velocity tends to equal the current’s velocity. The last force is the water resistance’s force which is applied normally to the current’s velocity plane. It reads(4)F→R=-KρwμARv→-v→// where K is a constant, μ is a shape factor, and AR is the sensor’s cross-section area perpendicular to the current. This fourth force tends to prevent any sensor’s motion perpendicular to the water current. The weight, buoyant and water resistance forces govern the vertical motion. As the sensor sinks towards the water depth that has a density equal to the sensor’s one, the weight, and buoyant forces balance each other. At the same time, the vertical sensor’s speed decreases.The free-body diagram is shown in Figure1. Applying the second Newton’s law yields(5)F→W+vecFB+F→R+F→C=ρsVsa→ or(6)ρs-ρwVsg→+αCv→C-v→//-αRv→-v→//=ρsVsa→ where αC=CσAC and αR=KρwμAR. Projecting (5) on the three space directions leads to the following dynamic system equations:(7)x¨+αCρsVsx˙=αCρsVsvCx,(8)y¨+αCρsVsy˙=αCρsVsvCyand(9)z¨+αRρsVsz˙=ρs-ρwρsg These equations are rewritten as follows:(10)x¨+Cxyx˙=CxyvCx,(11)y¨+Cxyy˙=CxyvCxand(12)z¨+Czz˙=CzvCz.where Cxy=αC/ρsVs, Cz=αR/ρsVs, and vCz=ρs-ρwg/ρsCz.Figure 1
Free-body diagram.


### 2.2. Numerical Solution

If the current velocitiesvCx and vCy and the water density ρw are constant, the above system of dynamic equations turns to a system of three linear independent second-order differential equations that can be solved analytically as done in [20]. Namely,(13)xt=x0+vCxt-t0+x˙0-vCxCxy1-e-Cxyt-t0,(14)yt=y0+vCyt-t0+y˙0-vCyCxy1-e-Cxyt-t0and(15)zt=z0+vCzt-t0+z˙0-vCzCz1-e-Czt-t0.where t0 is a reference time, x0, y0, and z0 are the sensor coordinates at t0, and x˙0, y˙0, and z˙0 are the sensor velocity components at t0.However, the real physical problem is more complex. The water density depends on the depth. In this work, the water density is assumed to vary linearly with depth. Mainly,(16)ρwz=ρ0+βz where ρ0 and β are two constants. Moreover, the current’s velocities depend on the sea-shore environment. Unlike the work proposed in [20], the current’s velocity is assumed nonconstant. Authors in [10, 17] used the kinematic model proposed in [25]. Specifically, the current velocities vCx and vCy depend on time t and the horizontal coordinates of the current position x and y. Namely,(17)vCx=k1λvsin⁡k2xcos⁡k3y+k1λcos⁡2k1t+k4(18)vCy=-λvcos⁡k2xsin⁡k3y+k5,where k1,k2,k3,k4,k5,λ, and v are parameters that depend on the underwater environment such as tide and bathymetry as explained in [17]. Within this framework, (10) and (11) are now coupled. Moreover, (10) and (11) are nonlinear differential equations, because of the x and y dependency of the water current velocity described in (17) and (18). Though (12) is still linear, its coefficients are now variable.In this work, we tackle the problem of solving the above system of differential equations using a numerical incremental iterative approach. The position of the sensor, which is given by the solution of (10), (11), and (12), is here calculated on a discrete set of times tj=jΔt;j∈N, where Δt is a constant time step. The time step is chosen so that the water density ρw and the current velocity components vCx and vCy are almost constant. Thus, the solution written in (13), (14), and (15) can be used for only one step. Let us assume that the sensor coordinates xj, yj, and zj and its velocity components x˙j, y˙j, and z˙j are known at a given time tj. Let us also consider that, during the period of time tj,tj+1, where tj+1=tj+Δt, ρw, vCx and vCy are almost constant. Within this hypothesis the sensor coordinates xj+1, yj+1, and zj+1at tj+1 are evaluated using (13), (14), and (15). Namely,(19)xj+1=xj+vCxΔt+x˙j-vCxCxy1-e-CxyΔt,(20)yj+1=yj+vCyΔt+y˙j-vCyCxy1-e-CxyΔtand(21)zj+1=zj+vCzΔt+z˙j-vCzCz1-e-CzΔt.Moreover, a simple time differentiation of (13), (14), and (15) and substituting tj for t0 give the sensor velocity components at time tj+1. More precisely,(22)x˙j+1=vCx+x˙j-vCxe-CxyΔt,(23)y˙j+1=vCy+y˙j-vCye-CxyΔtand(24)z˙j+1=vCz+z˙j-vCze-CzΔt.Knowing the initial position and velocity of the sensor, it is possible to calculate the position and velocity at any time using a time incremental procedure. This procedure can independently be repeated for each sensor to obtain the time positions of all sensors of the network.


## 2.1. Kinetics of Underwater Sensors

We are interested in predicting the trajectory of a sinking underwater sensor. The positionP of the sinking object is given by the three time-dependent coordinates (x,y,z). The velocity and the acceleration of the sensor can be written as v→=x˙i→+y˙j→+z˙k→ and a→=x¨i→+y¨j→+z¨k→, respectively. Four forces are acting upon the sensor as introduced in [20]. The first force is the weight F→W, which reads(1)F→W=ρsVsg→whereρs is the sensor’s density, Vs is the senor’s volume, and g→ is the terrestrial gravitational acceleration. The weight accelerates the sensor downwards. It is the main sinking force. It is countered by the buoyant force F→B, which is equal to the weight of the displaced water. It is written as(2)F→B=-ρwVsg→ where ρw is the water density, which can depend on the water depth z, i.e., ρwz. If the sensor density is higher than the water density, the weight overtakes the buoyant force and the sensor sinks down. However, the water density is mostly increasing with the water depth. Thus, the buoyant force increases as the sensor sinks while the sensor’s weight is constant. At a certain depth the water density can equal the sensor density and, consequently, the buoyant force balances the weight. Hence, this water depth will be the asymptotic value of the sensor’s z-coordinate.Unlike the weight and the buoyant forces, the third and fourth forces depend on the sensor’s kinematics. In other words, the weight and buoyant forces mainly depend on the water and sensor densities and the sensor’s shape. However, the third and fourth forces depend on the sensor and water velocities. The third force is the one applied by the water current to drive the sensor in the same direction as the current. It writes(3)F→C=CσACv→C-v→// where C is a constant, σ is a shape factor, AC is the sensor’s cross-section area facing the current, v→C is the velocity of the water current, and v→// is the projection of the sensor’s velocity in the same plane as the current velocity. More precisely, any upwelling or downwelling effects are neglected as justified in [20]. Thus, no current is considered in the vertical direction. The water current is simplified to v→C=vCxi→+vCyj→. Consequently, v→//=x˙i→+y˙j→. The current’s force is the only force that governs the horizontal sensor’s motion. It vanishes as the sensor’s velocity tends to equal the current’s velocity. The last force is the water resistance’s force which is applied normally to the current’s velocity plane. It reads(4)F→R=-KρwμARv→-v→// where K is a constant, μ is a shape factor, and AR is the sensor’s cross-section area perpendicular to the current. This fourth force tends to prevent any sensor’s motion perpendicular to the water current. The weight, buoyant and water resistance forces govern the vertical motion. As the sensor sinks towards the water depth that has a density equal to the sensor’s one, the weight, and buoyant forces balance each other. At the same time, the vertical sensor’s speed decreases.The free-body diagram is shown in Figure1. Applying the second Newton’s law yields(5)F→W+vecFB+F→R+F→C=ρsVsa→ or(6)ρs-ρwVsg→+αCv→C-v→//-αRv→-v→//=ρsVsa→ where αC=CσAC and αR=KρwμAR. Projecting (5) on the three space directions leads to the following dynamic system equations:(7)x¨+αCρsVsx˙=αCρsVsvCx,(8)y¨+αCρsVsy˙=αCρsVsvCyand(9)z¨+αRρsVsz˙=ρs-ρwρsg These equations are rewritten as follows:(10)x¨+Cxyx˙=CxyvCx,(11)y¨+Cxyy˙=CxyvCxand(12)z¨+Czz˙=CzvCz.where Cxy=αC/ρsVs, Cz=αR/ρsVs, and vCz=ρs-ρwg/ρsCz.Figure 1
Free-body diagram.


## 2.2. Numerical Solution

If the current velocitiesvCx and vCy and the water density ρw are constant, the above system of dynamic equations turns to a system of three linear independent second-order differential equations that can be solved analytically as done in [20]. Namely,(13)xt=x0+vCxt-t0+x˙0-vCxCxy1-e-Cxyt-t0,(14)yt=y0+vCyt-t0+y˙0-vCyCxy1-e-Cxyt-t0and(15)zt=z0+vCzt-t0+z˙0-vCzCz1-e-Czt-t0.where t0 is a reference time, x0, y0, and z0 are the sensor coordinates at t0, and x˙0, y˙0, and z˙0 are the sensor velocity components at t0.However, the real physical problem is more complex. The water density depends on the depth. In this work, the water density is assumed to vary linearly with depth. Mainly,(16)ρwz=ρ0+βz where ρ0 and β are two constants. Moreover, the current’s velocities depend on the sea-shore environment. Unlike the work proposed in [20], the current’s velocity is assumed nonconstant. Authors in [10, 17] used the kinematic model proposed in [25]. Specifically, the current velocities vCx and vCy depend on time t and the horizontal coordinates of the current position x and y. Namely,(17)vCx=k1λvsin⁡k2xcos⁡k3y+k1λcos⁡2k1t+k4(18)vCy=-λvcos⁡k2xsin⁡k3y+k5,where k1,k2,k3,k4,k5,λ, and v are parameters that depend on the underwater environment such as tide and bathymetry as explained in [17]. Within this framework, (10) and (11) are now coupled. Moreover, (10) and (11) are nonlinear differential equations, because of the x and y dependency of the water current velocity described in (17) and (18). Though (12) is still linear, its coefficients are now variable.In this work, we tackle the problem of solving the above system of differential equations using a numerical incremental iterative approach. The position of the sensor, which is given by the solution of (10), (11), and (12), is here calculated on a discrete set of times tj=jΔt;j∈N, where Δt is a constant time step. The time step is chosen so that the water density ρw and the current velocity components vCx and vCy are almost constant. Thus, the solution written in (13), (14), and (15) can be used for only one step. Let us assume that the sensor coordinates xj, yj, and zj and its velocity components x˙j, y˙j, and z˙j are known at a given time tj. Let us also consider that, during the period of time tj,tj+1, where tj+1=tj+Δt, ρw, vCx and vCy are almost constant. Within this hypothesis the sensor coordinates xj+1, yj+1, and zj+1at tj+1 are evaluated using (13), (14), and (15). Namely,(19)xj+1=xj+vCxΔt+x˙j-vCxCxy1-e-CxyΔt,(20)yj+1=yj+vCyΔt+y˙j-vCyCxy1-e-CxyΔtand(21)zj+1=zj+vCzΔt+z˙j-vCzCz1-e-CzΔt.Moreover, a simple time differentiation of (13), (14), and (15) and substituting tj for t0 give the sensor velocity components at time tj+1. More precisely,(22)x˙j+1=vCx+x˙j-vCxe-CxyΔt,(23)y˙j+1=vCy+y˙j-vCye-CxyΔtand(24)z˙j+1=vCz+z˙j-vCze-CzΔt.Knowing the initial position and velocity of the sensor, it is possible to calculate the position and velocity at any time using a time incremental procedure. This procedure can independently be repeated for each sensor to obtain the time positions of all sensors of the network.


## 3. Performance Evaluation

### 3.1. Network Model

Our mobile network can be modeled as a time varying graphG=Vt,Et where Vt represents the set of deployed sensor nodes and Et represents the communication link between them. Note the temporal dependence of Vt and Et. In fact, since our network is mobile, sensor nodes are continuously moving leading thus to variable links establishment between them as function of time as shown in Figures 2, 3, and 4. Indeed, u,v∈Et if node v is within the transmission range Rt of u at time t and hence a data communication is possible between them at time t. Our network model consists of a large number of freely floating sensor nodes N moving in a 3D rectangular domain D. Moreover, due to nodes’ mobility, sensor may leave the considered domain D, and thus the set of nodes, Vt, in our network depends on t.Figure 2
The positions of 400 sensors right after the deployment.Figure 3
Sensor nodes dispersion after 2 days and a half.Figure 4
Sensor nodes dispersion after 5 days.In our model, we suppose that the same transmission power is adopted by all the sensors in the network. The considered transmission power allows establishing communication between nodes that are within a communication rangeRt=2000m. Moreover, since we are interested in studying the network mobility over several days, we consider a large domain D=[0,600]×[400,4] km. Besides, we suppose that all the nodes will be randomly and uniformly deployed on D0=[0,80]×[-2,2] km, all at the same time as shown in Figure 2.


### 3.2. Numerical Simulation Parameters

The numerical solution presented in Section2.2 provides the trajectory of one sinking object where the speed parameters of (17) and (18) are assumed normally distributed as used in [10, 17]. More precisely, the mean and the standard deviation of k1, k2, k3, λ, and v are 0.3π,π,2π,1,and0.3 and 0.03π,0.1π,0.2π,0.1,and0.03, respectively. The parameters k4 and k5 are random variables whose mean and standard deviation are 1 and 0.1, respectively. The sensor volume Vs is set equal to 0.5L. Now, to get the water density ρw(z) as function of the depth z, ρ0 and β are set equal to 1025 and 0.02, respectively. Moreover, ρs is randomly chosen for each sensor between [1025,1045]. Note that, in order to numerically simulate a network of sensors, N sensor nodes will be created with different deployment positions and with different ρs.Figures2, 3, and 4 provide the 3D time evolution of the positions of 400 sensor nodes randomly scattered in a domain D0=0,80×-2,2km, along x, y, and z axis, respectively, as shown in Figure 2. After two days and a half, the sensor nodes dispersion along the x and y axis is between [220,320]km and [150,200]km, respectively, which shows the nodes’ displacements with the water currents. Two days and a half later, the nodes’ dispersion along the x and y axis progressed to reach ranges of [460,560]km and [392,398]km, respectively. Note however that the nodes’ dispersion along z covers anytime the range [0,4]km.It is worth pointing out that our numerical simulations consider two study cases: 3D mobility model and 2D one. In the 3D mobility model, sensor nodes having different densities and under water currents will end up moving in a 3D domain. However, if the deployed sensor nodes have the same density, they will end up moving in a well-defined water depth and hence their mobilities will be confined in a 2D environment. In the latter case, it is more reasonable to study the network coverage within the 2D domain relative to the water depth where the deployed sensors will sink.


### 3.3. Analysis of the Mobility Model

Our mobility model implies thaty˙ and z˙ approaches zero for large y and z, respectively. Thus nodes’ mobility is mainly oriented along the x-axis. Consequently, we rather focus on studying the horizontal movements of nodes, which is called sensors’ advection. To do so, we measure the absolute dispersion along the x-axis, as introduced in [18, 26]. Accordingly, the absolute dispersion reads (25)Ds2t,t0=xit-xit02=1N∑i=1Nxit-xit02where N=V is the number of sensors in the network, ⋯ indicates average over the sensor nodes, xit is the x coordinate of the ith sensor at time t, and t0 is the deployment time. Ds2 simply measures the average dispersion of sensor nodes along the x-axis relative to their initial position as a function of time. In other words, Ds2 provides the average temporal displacement of sensors from their initial deployment positions.


### 3.4. Mobility Model and Network Connectivity

To study the temporal network connectivity, we are interested in studying the behavior of the largest connected componentLCCt as it is our active network representative component susceptible to deliver reports to the offshore station. Note that, as expected, due to nodes’ movement, the largest connected component is time dependent as nodes may leave or incorporate the largest connected component over time. Our study of the LCCt includes its size variation over time, the first time its size drops below 90% of N, and finally the time variation of its bounding box. The bounding box functions of LCCt are defined as follows: xminLCCt=mini∈LCCtxit; xmaxLCCt=maxi∈LCCtxit. We aim at comparing the bounding box of LCCt with the bounding box of the whole network, i.e., xminGtmini∈Vtxit and xmaxGt=maxi∈Vtxit for both study cases, namely, 3D and 2D mobility models.Figures5 and 6 show the time evolution of xminGt, xmaxGt; the bounding box of the whole underwater sensor network and xminLCCt, xmaxLCCt; the bounding box of the largest connected component, for both study cases, namely, 3D and 2D mobility models, respectively. Note that these functions describe the motion of nodes in the whole network as well as the movement of the largest connected component. The plots show the displacement of nodes along the x-axis in a period of 5 days.Figure 5
Time evolution of theLCCt and the whole network bounding boxes. N=400. 3D mobility model.Figure 6
Time evolution of theLCCt and the whole network bounding boxes. N=400. 2D mobility model.Right after the deployment, the whole network is connectedxminGt=0=xminLCCt=0, and xmaxGt=0=xmaxLCCt=0, for both study cases. It is worth pointing out that functions xminGt and xmaxGt have the same slope, which is the maximum speed of the water current vm=1m/s, meaning that some sensors are in the heart of the water current. Moreover, for both study cases, the bounding box of the largest connected component is contained in the bounding box of the whole network where xminLCCt and xmaxLCCt move between xminGt, xmaxGt with an average velocity of vm which depicts that the largest connected component is covering the bounding box of the whole network with respect to time and that nodes in theLCC are in the middle of the current. Note that, for both study cases, the square root of the absolute dispersion Ds has also an average speed equal to vm, which point out that all the nodes are in the center of the current and move with the maximum speed. Moreover, it is worth pointing out that Ds scales linearly with time Ds∝t and thus the absolute dispersion Ds2∝t2 which assesses that our transport process is ballistic as defined in [26].Figures7 and 8 plot the size of the largest connected component with respect to time for a different number of total deployed sensor nodes for both study cases, namely; 3D and 2D mobility models, respectively. Most importantly, note that when the total number of sensor nodes is too high (N=800) or too low N=100 for the 3D case, the size of the largest connected component is rather constant. Indeed, for a very small number of sensor nodes, nodes are far from each other. Thus the average distance separating each pair of nodes is high. Consequently, nodes’ movement will have almost no impact on the average mutual distance between nodes. Similarly, when the total number of nodes is high and so is the network density, the average distance separating any pair of nodes is so small such that nodes displacement will not lead to disconnection and hence nodes’ movement will have no impact on the size of the largest connected component, which is so high. However, for medium values of the total number of deployed sensors, the size of the largest connected component may fluctuate. A very fast fluctuation is observed for N=400 for both mobility models. Indeed, for this range of total number of sensor nodes, the distance separating any pair of nodes is neither too high nor too low, such that any node movement will have an important impact on the size of the largest connected component. This is because nodes can get closer or farther rapidly and hence cause a fast variation in the size of the largest connected component. Similar behaviors are observed for the 2D case.Figure 7
Time evolution of theLCCt size. 3D mobility model.Figure 8
Time evolution of theLCCt size. 2D mobility model.In Tables1 and 2, we present the first time when the size of the largest connected component drops below 90% for different values of N for both study cases, namely, 3D and 2D mobility models, respectively. As expected, this time is increasing as a function of N as observed in Tables 1 and 2. Indeed, according to Figures 7 and 8 a faster drop in the size of the largest connected component is observed as the total number of nodes is decreasing. Note that, for N=800, the size of the largest connected component is very high and constant. Actually, for N=800, the size of the largest connected component will never drop below 90% for the 3D mobility model. Similarly, for the 2D deployment case, starting from N=600, the largest connected component will never drop below 90%.Table 1
Connectivity: first time when the size of the largest connected component drops below90%N. 3D mobility model.
N Time (s) 100 1800  200 1820  400 9000  600 23400  800 SimTime = 433800Table 2
Connectivity: first time when the size of the largest connected component drops below90%N. 2D mobility model.
N Time (s) 100 1800  200 1820  400 100800  600 SimTime = 433800  800 SimTime = 433800


### 3.5. Mobility Model and Network Coverage

Now, let us get more insights into the time evolution of the network coverage under both mobility models. A node coverage area is the area where a given node can perform some dedicated measurements or detect a specified event. In a 3D domain, a node coverage is modeled as a sphere of radiusRs centered at the sensor position. In our work, a point of the domain is considered covered, if it belongs to the coverage area of a node belonging to the largest connected component. Accordingly, a domain region that is covered by a node not belonging to the largest connected component is considered uncovered. Indeed, measurements or reports coming from isolated nodes are most likely undeliverable to the offshore station as it is hard to find a multihop path to reach the sink. Accordingly, inspired from [18], we define two measures of coverage, namely, static and dynamic coverage.Definition 1 (static area coverage).
The area coverage of a sensor network at timet is the fraction of the geographical area covered by one or more sensors in LCCt at time t.Definition 2 (dynamic area coverage over a time interval).
The area coverage of a mobile sensor network during the time interval[0,t) is the fraction of the geographical area covered by at least one sensor in LCC at some point of time within [0,t).Using either static or dynamic coverage depends mainly on the application requirements. For instance, the static area coverage is more appropriate for applications that require continuous coverage of a given underwater field under exploration. However, the dynamic area coverage is more suitable for applications that do not require continuous coverage of all points, but rather cover a given field within some time interval like event detection in a specific field.Figures9 and 10 depict the static coverage as a function of time for different values of N for both study cases, namely, 3D and 2D mobility models, respectively. As expected, for both study cases, the static coverage increases with a growing number of nodes. Recall that, in our work, we study the network coverage of the largest connected component which is able to reach the sink. Considering the coverage of each sensor node is not really meaningful since isolated sensor nodes cannot reach the sink and hence cannot be considered as part of the network which explains the low achieved static coverage rate.Figure 9
Time evolution of the static coverage. 3D mobility model.Figure 10
Time evolution of the static coverage. 2D mobility model.Similar to Figures7 and 8 the static coverage, for both study cases, exhibits a rather smooth variation when the total number of nodes is too high N=800 or too low N=100 and a fast fluctuation when the number of nodes is in between especially for N=400. Indeed, for small values of N, nodes are already far from each other and hence nodes’ displacement will not impact these distances and so the coverage will not be subject to significant variation. Likewise, for high values of N, the network density is high and hence nodes are rather close to each other. Thus nodes’ movement will not radically change the distances separating pairs of nodes and hence a smooth fluctuation is observed in the coverage rate. However, for medium sized networks, distances between pair of nodes are neither high nor small, such that the least nodes’ movement may cause a fast fluctuation in the coverage rate, which is the case especially for N=400.Now, let us get more insights into the time evolution of the dynamic coverage rate for both study cases, namely, 3D and 2D mobility models shown in Figures11 and 12, respectively. Note that the dynamic coverage rate is increasing with time. Indeed, the dynamic coverage rate is rather cumulative coverage. Thanks to the node mobility, new places will be visited over time and hence inevitably the cumulative dynamic coverage rate will be increasing.Figure 11
Time evolution of the dynamic coverage. 3D mobility model.Figure 12
Time evolution of the dynamic coverage. 2D mobility model.


## 3.1. Network Model

Our mobile network can be modeled as a time varying graphG=Vt,Et where Vt represents the set of deployed sensor nodes and Et represents the communication link between them. Note the temporal dependence of Vt and Et. In fact, since our network is mobile, sensor nodes are continuously moving leading thus to variable links establishment between them as function of time as shown in Figures 2, 3, and 4. Indeed, u,v∈Et if node v is within the transmission range Rt of u at time t and hence a data communication is possible between them at time t. Our network model consists of a large number of freely floating sensor nodes N moving in a 3D rectangular domain D. Moreover, due to nodes’ mobility, sensor may leave the considered domain D, and thus the set of nodes, Vt, in our network depends on t.Figure 2
The positions of 400 sensors right after the deployment.Figure 3
Sensor nodes dispersion after 2 days and a half.Figure 4
Sensor nodes dispersion after 5 days.In our model, we suppose that the same transmission power is adopted by all the sensors in the network. The considered transmission power allows establishing communication between nodes that are within a communication rangeRt=2000m. Moreover, since we are interested in studying the network mobility over several days, we consider a large domain D=[0,600]×[400,4] km. Besides, we suppose that all the nodes will be randomly and uniformly deployed on D0=[0,80]×[-2,2] km, all at the same time as shown in Figure 2.


## 3.2. Numerical Simulation Parameters

The numerical solution presented in Section2.2 provides the trajectory of one sinking object where the speed parameters of (17) and (18) are assumed normally distributed as used in [10, 17]. More precisely, the mean and the standard deviation of k1, k2, k3, λ, and v are 0.3π,π,2π,1,and0.3 and 0.03π,0.1π,0.2π,0.1,and0.03, respectively. The parameters k4 and k5 are random variables whose mean and standard deviation are 1 and 0.1, respectively. The sensor volume Vs is set equal to 0.5L. Now, to get the water density ρw(z) as function of the depth z, ρ0 and β are set equal to 1025 and 0.02, respectively. Moreover, ρs is randomly chosen for each sensor between [1025,1045]. Note that, in order to numerically simulate a network of sensors, N sensor nodes will be created with different deployment positions and with different ρs.Figures2, 3, and 4 provide the 3D time evolution of the positions of 400 sensor nodes randomly scattered in a domain D0=0,80×-2,2km, along x, y, and z axis, respectively, as shown in Figure 2. After two days and a half, the sensor nodes dispersion along the x and y axis is between [220,320]km and [150,200]km, respectively, which shows the nodes’ displacements with the water currents. Two days and a half later, the nodes’ dispersion along the x and y axis progressed to reach ranges of [460,560]km and [392,398]km, respectively. Note however that the nodes’ dispersion along z covers anytime the range [0,4]km.It is worth pointing out that our numerical simulations consider two study cases: 3D mobility model and 2D one. In the 3D mobility model, sensor nodes having different densities and under water currents will end up moving in a 3D domain. However, if the deployed sensor nodes have the same density, they will end up moving in a well-defined water depth and hence their mobilities will be confined in a 2D environment. In the latter case, it is more reasonable to study the network coverage within the 2D domain relative to the water depth where the deployed sensors will sink.


## 3.3. Analysis of the Mobility Model

Our mobility model implies thaty˙ and z˙ approaches zero for large y and z, respectively. Thus nodes’ mobility is mainly oriented along the x-axis. Consequently, we rather focus on studying the horizontal movements of nodes, which is called sensors’ advection. To do so, we measure the absolute dispersion along the x-axis, as introduced in [18, 26]. Accordingly, the absolute dispersion reads (25)Ds2t,t0=xit-xit02=1N∑i=1Nxit-xit02where N=V is the number of sensors in the network, ⋯ indicates average over the sensor nodes, xit is the x coordinate of the ith sensor at time t, and t0 is the deployment time. Ds2 simply measures the average dispersion of sensor nodes along the x-axis relative to their initial position as a function of time. In other words, Ds2 provides the average temporal displacement of sensors from their initial deployment positions.


## 3.4. Mobility Model and Network Connectivity

To study the temporal network connectivity, we are interested in studying the behavior of the largest connected componentLCCt as it is our active network representative component susceptible to deliver reports to the offshore station. Note that, as expected, due to nodes’ movement, the largest connected component is time dependent as nodes may leave or incorporate the largest connected component over time. Our study of the LCCt includes its size variation over time, the first time its size drops below 90% of N, and finally the time variation of its bounding box. The bounding box functions of LCCt are defined as follows: xminLCCt=mini∈LCCtxit; xmaxLCCt=maxi∈LCCtxit. We aim at comparing the bounding box of LCCt with the bounding box of the whole network, i.e., xminGtmini∈Vtxit and xmaxGt=maxi∈Vtxit for both study cases, namely, 3D and 2D mobility models.Figures5 and 6 show the time evolution of xminGt, xmaxGt; the bounding box of the whole underwater sensor network and xminLCCt, xmaxLCCt; the bounding box of the largest connected component, for both study cases, namely, 3D and 2D mobility models, respectively. Note that these functions describe the motion of nodes in the whole network as well as the movement of the largest connected component. The plots show the displacement of nodes along the x-axis in a period of 5 days.Figure 5
Time evolution of theLCCt and the whole network bounding boxes. N=400. 3D mobility model.Figure 6
Time evolution of theLCCt and the whole network bounding boxes. N=400. 2D mobility model.Right after the deployment, the whole network is connectedxminGt=0=xminLCCt=0, and xmaxGt=0=xmaxLCCt=0, for both study cases. It is worth pointing out that functions xminGt and xmaxGt have the same slope, which is the maximum speed of the water current vm=1m/s, meaning that some sensors are in the heart of the water current. Moreover, for both study cases, the bounding box of the largest connected component is contained in the bounding box of the whole network where xminLCCt and xmaxLCCt move between xminGt, xmaxGt with an average velocity of vm which depicts that the largest connected component is covering the bounding box of the whole network with respect to time and that nodes in theLCC are in the middle of the current. Note that, for both study cases, the square root of the absolute dispersion Ds has also an average speed equal to vm, which point out that all the nodes are in the center of the current and move with the maximum speed. Moreover, it is worth pointing out that Ds scales linearly with time Ds∝t and thus the absolute dispersion Ds2∝t2 which assesses that our transport process is ballistic as defined in [26].Figures7 and 8 plot the size of the largest connected component with respect to time for a different number of total deployed sensor nodes for both study cases, namely; 3D and 2D mobility models, respectively. Most importantly, note that when the total number of sensor nodes is too high (N=800) or too low N=100 for the 3D case, the size of the largest connected component is rather constant. Indeed, for a very small number of sensor nodes, nodes are far from each other. Thus the average distance separating each pair of nodes is high. Consequently, nodes’ movement will have almost no impact on the average mutual distance between nodes. Similarly, when the total number of nodes is high and so is the network density, the average distance separating any pair of nodes is so small such that nodes displacement will not lead to disconnection and hence nodes’ movement will have no impact on the size of the largest connected component, which is so high. However, for medium values of the total number of deployed sensors, the size of the largest connected component may fluctuate. A very fast fluctuation is observed for N=400 for both mobility models. Indeed, for this range of total number of sensor nodes, the distance separating any pair of nodes is neither too high nor too low, such that any node movement will have an important impact on the size of the largest connected component. This is because nodes can get closer or farther rapidly and hence cause a fast variation in the size of the largest connected component. Similar behaviors are observed for the 2D case.Figure 7
Time evolution of theLCCt size. 3D mobility model.Figure 8
Time evolution of theLCCt size. 2D mobility model.In Tables1 and 2, we present the first time when the size of the largest connected component drops below 90% for different values of N for both study cases, namely, 3D and 2D mobility models, respectively. As expected, this time is increasing as a function of N as observed in Tables 1 and 2. Indeed, according to Figures 7 and 8 a faster drop in the size of the largest connected component is observed as the total number of nodes is decreasing. Note that, for N=800, the size of the largest connected component is very high and constant. Actually, for N=800, the size of the largest connected component will never drop below 90% for the 3D mobility model. Similarly, for the 2D deployment case, starting from N=600, the largest connected component will never drop below 90%.Table 1
Connectivity: first time when the size of the largest connected component drops below90%N. 3D mobility model.
N Time (s) 100 1800  200 1820  400 9000  600 23400  800 SimTime = 433800Table 2
Connectivity: first time when the size of the largest connected component drops below90%N. 2D mobility model.
N Time (s) 100 1800  200 1820  400 100800  600 SimTime = 433800  800 SimTime = 433800


## 3.5. Mobility Model and Network Coverage

Now, let us get more insights into the time evolution of the network coverage under both mobility models. A node coverage area is the area where a given node can perform some dedicated measurements or detect a specified event. In a 3D domain, a node coverage is modeled as a sphere of radiusRs centered at the sensor position. In our work, a point of the domain is considered covered, if it belongs to the coverage area of a node belonging to the largest connected component. Accordingly, a domain region that is covered by a node not belonging to the largest connected component is considered uncovered. Indeed, measurements or reports coming from isolated nodes are most likely undeliverable to the offshore station as it is hard to find a multihop path to reach the sink. Accordingly, inspired from [18], we define two measures of coverage, namely, static and dynamic coverage.Definition 1 (static area coverage).
The area coverage of a sensor network at timet is the fraction of the geographical area covered by one or more sensors in LCCt at time t.Definition 2 (dynamic area coverage over a time interval).
The area coverage of a mobile sensor network during the time interval[0,t) is the fraction of the geographical area covered by at least one sensor in LCC at some point of time within [0,t).Using either static or dynamic coverage depends mainly on the application requirements. For instance, the static area coverage is more appropriate for applications that require continuous coverage of a given underwater field under exploration. However, the dynamic area coverage is more suitable for applications that do not require continuous coverage of all points, but rather cover a given field within some time interval like event detection in a specific field.Figures9 and 10 depict the static coverage as a function of time for different values of N for both study cases, namely, 3D and 2D mobility models, respectively. As expected, for both study cases, the static coverage increases with a growing number of nodes. Recall that, in our work, we study the network coverage of the largest connected component which is able to reach the sink. Considering the coverage of each sensor node is not really meaningful since isolated sensor nodes cannot reach the sink and hence cannot be considered as part of the network which explains the low achieved static coverage rate.Figure 9
Time evolution of the static coverage. 3D mobility model.Figure 10
Time evolution of the static coverage. 2D mobility model.Similar to Figures7 and 8 the static coverage, for both study cases, exhibits a rather smooth variation when the total number of nodes is too high N=800 or too low N=100 and a fast fluctuation when the number of nodes is in between especially for N=400. Indeed, for small values of N, nodes are already far from each other and hence nodes’ displacement will not impact these distances and so the coverage will not be subject to significant variation. Likewise, for high values of N, the network density is high and hence nodes are rather close to each other. Thus nodes’ movement will not radically change the distances separating pairs of nodes and hence a smooth fluctuation is observed in the coverage rate. However, for medium sized networks, distances between pair of nodes are neither high nor small, such that the least nodes’ movement may cause a fast fluctuation in the coverage rate, which is the case especially for N=400.Now, let us get more insights into the time evolution of the dynamic coverage rate for both study cases, namely, 3D and 2D mobility models shown in Figures11 and 12, respectively. Note that the dynamic coverage rate is increasing with time. Indeed, the dynamic coverage rate is rather cumulative coverage. Thanks to the node mobility, new places will be visited over time and hence inevitably the cumulative dynamic coverage rate will be increasing.Figure 11
Time evolution of the dynamic coverage. 3D mobility model.Figure 12
Time evolution of the dynamic coverage. 2D mobility model.


## 4. Conclusion

In this paper, we proposed a physically inspired mobility model especially tailored for underwater freely floating acoustic sensor networks. To the best of our knowledge, this work is one of the most accurate mobility model that captures most of the physical aspects impacting a freely floating sensor node mobility. After studying the proposed mobility model, we analyzed its temporal impact on the network connectivity and coverage. We showed that the network mobility effect on the coverage and connectivity is more significant on intermediately dense UW-ASNs, while a lower effect was recorded on the coverage and connectivity for low- and high-density UW-ASNs. Our work objective was to provide the underwater network research community with a realistic mobility model that could be exploited in conceiving networking communication protocols such as routing, localization, and medium access.



---

*Source: 9818931-2019-02-12.xml*