*Article* **Underwater Sparse Acoustic Sensor Array Design under Spacing Constraints Based on a Global Enhancement Whale Optimization Algorithm**

**Lening Wang 1, Hangfang Zhao 1,2,3,\* and Qide Wang <sup>4</sup>**


**Abstract:** Sparse arrays with low cost and engineering complexity are widely applied in many fields. However, the high peak sidelobe level (PSLL) of a sparse array causes the degradation of weak target detection performance. Particularly for the large size of underwater low-frequency sensors, the design problem requires a minimum spacing constraint, which further increases the difficulty of PSLL suppression. In this paper, a novel swarm-intelligence-based approach for sparse sensor array design is proposed to reduce PSLL under spacing constrains. First, a global enhancement whale optimization algorithm (GEWOA) is introduced to improve the global search capability for optimal arrays. A three-step enhanced strategy is used to enhance the ergodicity of element positions over the aperture. In order to solve the adaptation problem for discrete array design, a position decomposition method and a V-shaped transfer function are introduced into off-grid and on-grid arrays, respectively. The effectiveness and superiority of the proposed approach is validated using experiments for designing large-scale low-frequency arrays in the marine environment. The PSLL of the off-grid array obtained by GEWOA was nearly 3.8 dB lower than that of WOA. In addition, compared with other intelligent algorithms, the on-grid array designed using GEWOA had the lowest PSLL.

**Keywords:** sparse array design; underwater acoustic sensor array; peak sidelobe level; intelligence optimization algorithm

#### **1. Introduction**

Sparse sensor arrays are widely used in observation and communication systems, such as radar, sonar, satellite, etc. [1–3]. With fewer array elements, such arrays obtain better array performance, including the narrow main lobe width and the PSLL compared with equally spaced arrays. In addition, such arrays have clear advantages in terms of the installation cost and system complexity.

Generally, the design problems of sparse sensor arrays are mainly divided into two types [4]: (1) The element positions are optimized to minimize the PSLL with a fixed number of array elements and aperture [5]. (2) The array beampattern is designed to match the desired beampattern by optimizing the array element positions and weights while reducing the number of array elements [6,7]. Some arrays, such as large towed arrays and the seabed array, are prone to shift in array position and change in array configuration during operation, which causes a mismatch between the real positions and the design ones.

For example, the towed array is flexible with oil-filled polyurethane tubing. When the array is towed, it is difficult to keep the array straight; therefore, the real positions of the

**Citation:** Wang, L.; Zhao, H.; Wang, Q. Underwater Sparse Acoustic Sensor Array Design under Spacing Constraints Based on a Global Enhancement Whale Optimization Algorithm. *Appl. Sci.* **2022**, *12*, 11825. https://doi.org/10.3390/ app122211825

Academic Editors: Lin Mu, Enjin Zhao and Hao Qin

Received: 6 October 2022 Accepted: 17 November 2022 Published: 21 November 2022

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

elements deviate from the designed positions. For the seabed array mounted at the bottom of sea, the uncertainties of the element positions and the entire array configuration are induced by ship drift and subsurface currents when the array deploys to the seafloor from the layout vessel. We focus on the first type problem under the marine environment to design a random sparse array with low PSLL to enhance the detection of weak targets [8].

The existing methods to solve the first problem are mainly divided into two categories: deterministic optimization methods and stochastic optimization methods. Deterministic optimization methods, including minimum redundancy array [9], coprime array [10], and difference coarray [11], are easy to implement; however, the designed arrays are generally suboptimal. Stochastic optimization methods use intelligent algorithms to find the element positions that satisfy the design requirements. Genetic algorithm (GA) [12], particle swarm optimization (PSO) algorithm [13], and simulated annealing algorithm [14] have been successfully applied to sparse sensor array design.

However, these intelligent algorithms converge slowly and easily fall into local optima. Some improved algorithms [15] have been proposed to solve the problem for array design for example, a series of swarm intelligence optimization algorithms mimicking swarm behavior, including the wolf pack algorithm (WPA) [16], bat algorithm (BA) [17], and cuckoo search (CS) algorithm [18], which show good search accuracy and quick convergence speed. Although the stochastic optimization method suffers from global performance degradation during the design process of large-scale arrays, this method outperforms the deterministic optimization method in terms of the array performance as the calculation capacity increases.

The whale optimization algorithm (WOA) [19]—a swarm intelligence algorithm—has been applied successfully to sparse sensor array design for both fixed aperture and nonfixed aperture and has shown its advantages in sidelobe suppression and null steering for small- and medium-scale arrays [20,21]. Combining the salp swarm algorithm (SSA) and WOA, a novel swarm algorithm was proposed to improve its convergence accuracy for conformal antenna array design [22] and for dual-band aperture-coupled antenna design [23].

Nevertheless, the global search capability of WOA decreases with the expansion of the search range in large-scale array design. Researchers have attempted to improve the algorithm from several aspects in successive steps. The initialization strategy played an important role [24]. Chaos initialization improves population diversity due to its ergodic and random nature, and the chaotic particle algorithm was demonstrated to successfully reduce the PSLL of sparse sensor arrays [25].

With the development of intelligent algorithms, more chaotic initialization algorithms, including the chaotic invasive weed optimization algorithm [26], chaotic cuckoo search algorithm [27], and chaotic sparrow search algorithm [28], have been applied to sparse array design for improving global search capability. These success examples lead us to believe that WOA embedded with chaotic initialization will improve the performance of large-scale array design.

In the next step, search strategy optimization can be achieved by introducing adaptive weights [29] and search path planning [30]. WOA with adaptive weights had a good global search capability [31,32]. Search path planning is another effective method. Levy flight strategy, as a random search path strategy with global search capability, was introduced to WOA [33]. Although the method has not been widely applied to array design within WOA, it has been successful within other swarm algorithms for linear and circular array design, including the differential evolutionary algorithm [34], CS algorithm [27], InvasiveWeed Optimization [35], and biogeography-based optimization approach [36]. Unfortunately, the Levy flight strategy destroys the original whale spiral search, which causes the convergence speed to slow down; hence, a new search strategy should be introduced to maintain the speed.

In the third step to produce a new population, it was found that opposition-based learning (OBL) can also improve the search accuracy of WOA [37]. The combined strategy of adaptive weights with OBL was applied to high-dimensional global optimization

problems [38]. However, it is easy to fall into the inverse local optimum, since the search range of the inverse solution is limited. A new learning algorithm should be developed to jump out of local minima. In short, the existing methods still have some shortcomings in the large-scale array design, particularly for the insufficient global search capability. To tackle this issue, GEWOA with a three-step enhanced global search strategy is introduced in this paper.

Additionally, the sparse sensor array design is a discrete problem, while the WOA is suitable for the solution of continuous variables. Therefore, WOA needs to be customized for discrete problems. Continuous variables discretization is generally achieved by transfer functions. If the discrete array element is on an integral multiple of the half wavelength, then this is an on-grid problem, otherwise it is an off-grid problem. The former is a discrete optimization problem with popular S-shaped and V-shaped transfer functions [39,40]. The latter uses binary discrete continuous variables with the ability to achieve lower PSLL [41,42]. Considering the physical limit of the array element size, the array design needs to constrain the minimum spacing between adjacent array elements while minimizing the PSLL.

A novel approach based on GEWOA for sparse sensor array design is proposed to reduce PSLL. A three-step enhanced global search strategy is introduced in GEWOA to improve the global search capability. In the initial stage, chaotic initialization is embedded in GEWOA, which enhances the ergodicity of the algorithm. In the search stage, the conventional spiral strategy is replaced by the Archimedean spiral strategy (ASS) to avoid premature algorithm results. In the offspring selection stage, refraction learning (RL) is used to obtain the inverse solution of the offspring, which avoids falling into local optima.

Moreover, to solve the adaptation problem for the discrete array design based on GEWOA, a position decomposition method is applied to make the element positions of off-grid arrays continuous, and a V-shaped transfer function is introduced into the on-grid array design to map the search position into the discrete grid points to obtain the desired array. The effectiveness of the proposed method is validated on the design tasks for sparse sensor arrays. To recap, the main contributions of this paper can be summarized as follows:


The rest of this paper is organized as follows: Section 2 presents the problem description and the background. Section 3 introduces the proposed GEWOA. In Section 4, the effectiveness and superiority of the proposed method are demonstrated. Finally, Section 5 concludes this paper.

#### **2. Problem Description**

Let us consider that a sparse linear array has *N* isotropic radiative elements, which are distributed along the *x*-axis at **D** = [*d*0, *d*1, ... , *di*, ... , *dN*−1], *i* ∈ [0, *N* − 1]. The array aperture is *L*. The beampattern of the array is given by:

$$P = \sum\_{i=1}^{N} w\_i e^{j\frac{2\pi}{\lambda}d\_i \sin(\theta - \theta\_0)}\tag{1}$$

where j <sup>=</sup> √−1, *<sup>λ</sup>* denotes the wavelength of the operating frequency, *wi* denotes the complex excitation of the *i*th element in the array, *θ* is the steering angle, and *θ*<sup>0</sup> represents the desired beam direction. The schematic of the array elements distribution is shown in Figure 1. Sparse arrays can be divided into on-grid arrays and off-grid arrays according to the positions of the array elements. The grid points are set at positions that are integer multiples of half the wavelength.

The black solid circles in Figure 1 are the actual positions of the array elements. The red hollow circles are grid points. If all array elements fall on the grid points, it is an on-grid array. Otherwise, it is an off-grid array. The PSLL of the off-grid array could be lower than that of on-grid array. In other words, an off-grid array can use fewer elements than an on-grid array to meet the same PSLL. However, the on-grid array is easier to process with less computation charge. Therefore, different design methods can be selected with a trade-off between the PSLL and computation charge.

**Figure 1.** Schematic of the array element distribution.

In this paper, we focus on the layout problem of large-scale sparse sensor arrays, and all weights are set to 1 (*wn* = 1 for all elements). The beampattern of a uniform array is a Dirichlet function with a PSLL of about −13 dB. However, the positions of the sparse array can be distributed nonuniformly. Assuming that the weight vector of a uniform array is **w**˜ , then the weight vector **w** of the sparse array can be expressed as **w** = **w**˜ **t**, where **t** = [1, 0, 1, ... , 0, 1] can be considered a tapering vector, which consists of 1 and 0, with 1 for element and 0 for no element, and denotes the Hadamard product.

The tapering operation enables the side lobes to be lower than those of the uniform array. The optimal positions by minimizing PSLL are filled with 1 in vector **t**. The beampattern of a sparse array obtained through Equation (1) is the convolution of the Dirichlet function with the spatial spectrum of weight vector **w**, which is based on the array positions. Hence, a sparse sensor array with a low PSLL can be obtained by optimizing the array element positions. Considering the physical limit of the sensors' size, the minimum spacing of the elements must be restricted. Therefore, the array design problem can be formulated as follows:

$$\min\_{\mathbf{X}} \mathbf{PSLL}(P) \quad \text{s.t.} \begin{cases} \mathbf{x}\_i - \mathbf{x}\_j \ge d\_c > 0\\ i, j \in \mathbf{Z}, 0 \le j \le i \le N - 1\\ \mathbf{x}\_0 = \mathbf{0}, \mathbf{x}\_{N-1} = L \end{cases} \tag{2}$$

where *xi* and *xj* denote the position of the *i*th and *j*th element, respectively. *dc* is the minimum distance between the adjacent elements. **Z** denotes the set of positive integers.

#### **3. Methodology**

#### *3.1. Overview*

The flowchart of the proposed method is shown in Figure 2. For off-grid array design, the position decomposition is used to make the element positions continuous. This strategy scales the array aperture *L* to the new aperture *L*∗, and then the proposed GEWOA is

performed to obtain the optimized design result. For on-grid array design, GEWOA is directly performed to update the positions of the population, and then a V-shaped transfer function is used to map the search positions into the discrete grid points to obtain the desired array.

**Figure 2.** Flowchart of the GEWOA method.

In GEWOA, a three-step enhanced global search strategy is introduced to improve the global search capability of the array element positions as shown in the red dashed box. Chaotic initialization is used to initialize the population. Next, the position updating process is executed, where the three strategies, including the shrinking encircling mechanism, ASS, and random search, are selected according to preset rules.

After the position updating, RL is performed to seek the inverse solution to optimize the offspring. Next, the best offspring is selected by calculating the PSLL according to the objective function. Finally, the iterative optimization process is continued until the convergence condition is satisfied or the maximum number of iterations is reached.

#### *3.2. Global Enhancement Whale Optimization Algorithm*

GEWOA is a swarm intelligence optimization algorithm inspired by the foraging process of humpback whales. It is inspired by WOA and continues to follow its strategy execution process. The execution process simulates the behavior of whales attacking prey through a shrinking encircling mechanism and spiral updating position. In addition, a random search strategy is performed to avoid falling into a local solution [19]. However, the global search capability of the original strategy decreases as the search range increases. A three-step enhanced global search strategy, including chaotic initialization, ASS, and RL, is proposed to improve the search performance in GEWOA.

#### 3.2.1. Chaotic Initialization

Chaos is a phenomenon that can traverse all states without repetition in a certain range [24]. We embed the chaos model to calculate the initial value aiming at increases the diversity of the population. Conventional chaotic mapping functions include Logistic mapping, Tent mapping, Kent mapping, Henon mapping, and Sin mapping.

After setting the initial value *<sup>s</sup>*<sup>0</sup> and the parameters, the chaotic sequence **<sup>S</sup>** ( *<sup>s</sup>*1, *<sup>s</sup>*2, ... , *sN*−2) can be obtained by using the equations in Table 1, where *α* ∈ (0, 4] is the control parameter of Logistic mapping. *β* = 1.4, *ζ* = 0.3 are the control parameters of Henon mapping, which allows the equation to enter a chaotic state and generates a chaotic sequence with strong randomness.

*μ* ∈ (0, 1) is the the control parameter of Kent mapping, *η* ∈ (0, 2] is the control parameter of Tent mapping, and *γ* ∈ [0, 1] is the control parameter of Sin mapping. Then, the positions of the initial population **<sup>X</sup>** (*<sup>x</sup>* 1, *<sup>x</sup>* 2,..., *<sup>x</sup> <sup>N</sup>*−2) can be obtained by:

$$
\overline{\hat{\mathbf{x}}}\_{n} = \overline{\hat{s}}\_{n}(\mathbf{u}b - lb) + lb,\\
\mathbf{n} = 0, 1, 2, \dots, N - 1 \tag{3}
$$


**Table 1.** Chaotic mapping equations.

In this paper, different chaotic initialization strategies are simulated to determine the best one as the optimal initialization strategy that has the lowest PSLL. Compared with the randomly distributed population in WOA, chaotic initialization overcomes the shortcoming of falling into a local solution to a certain extent.

For a swarm of whales with *M* individuals, the search space of each whale position is *N*-dimensional. The position of the *j*th whale in the population after the *q*th iteration is **X**(*q*)=(*x*0,*j*, *x*1,*j*, ... , *xi*,*j*, ... , *xN*−1,*j*), *i* ∈ [0, *N* − 1], *j* ∈ [0, *M* − 1]. Here, **X** ∈ [*lb*, *ub*], *ub* is the upper boundary of the search range, and *lb* is the lower boundary of the search range.

#### 3.2.2. Shrinking Encircling Mechanism

The optimal position is regarded as the position of a target prey, and other individuals surround the optimal position. Such behavior is called the shrink encircling mechanism, and its expressions are as follows:

$$\tilde{\mathbf{D}} = |\mathbf{C}\mathbf{X}^\*(q) - \mathbf{X}(q)|\tag{4}$$

$$\mathbf{X}(q+1) = \mathbf{X}^\*(q) - \mathbf{A} \cdot \mathbf{D} \tag{5}$$

where *q* is the number of iterations, and **X**∗(*q*) is the optimal position. **X**(*q*) denotes the individual position in the population after *q*th iterations. **X**(*q* + 1) is the updated position. · denotes the multiplication of the corresponding elements. **A** and **C** are coefficient vectors:

$$\mathbf{A} = 2\mathbf{a} \cdot \mathbf{r} - \mathbf{a} \tag{6}$$

$$\mathbf{C} = 2\mathbf{r} \tag{7}$$

where **a** decreases linearly from 2 to 0, and **r** is a random vector in [0, 1]. It is worth noting that the individuals move to the optimal position by adjusting **A** and **C**.

#### 3.2.3. Random Search

To avoid falling a local optimum, an individual is randomly selected as a target prey, and the positions of the population are updated by:

$$\dot{\mathbf{D}} = |\mathbf{C} \mathbf{X}\_{rand}(q) - \mathbf{X}(q)|\tag{8}$$

$$\mathbf{X}(q+1) = \mathbf{X}\_{rand}(q) - \mathbf{A} \cdot \mathbf{D} \tag{9}$$

where **X***rand* is a randomly selected whale position.

#### 3.2.4. Archimedean Spiral Strategy

Whales swim in a spiral shape while shrinking towards their prey with the spiral position update model:

$$\mathbf{X}(q+1) = \mathbf{X}^\*(q) + \mathbf{D}\_p e^{bl} \cos(2\pi l) \tag{10}$$

where **D***<sup>p</sup>* denotes the distance between candidate whales and their prey, *b* is a constant used to control the logarithmic spiral shape, and *l* is a random number between [−1, 1].

The spiral position update method of WOA follows the logarithmic spiral model, which is characterized by equal spiral angles but unequal pitches. Such a method can find the optimal value quickly; however, it likely misses the intermediate optimal solution. Archimedean spiral curve provides a better idea to improve the original spiral search strategy, which avoids the reduced accuracy caused by long search steps. Each position dimension of a whale is generated by (11) instead of (10).

$$\mathbf{X}(q+1) = \mathbf{X}^\*(q) + (\mathbf{D}\_{\overline{p}} + b \cdot l)e^{bl}\cos(2\pi l) \tag{11}$$

Figure 3 shows the Archimedean spiral curve and the logarithmic spiral curve, where the blue solid line is the Archimedean spiral curve, and the red dashed line is the logarithmic spiral curve. It can be seen that the Archimedean spiral curve is characterized by equal pitch, and the search step can be set by setting the pitch.

In the process of WOA, the strategy selection probability is set to 50% to choose encircling prey or the spiral updating position:

$$\mathbf{X}(q+1) = \begin{cases} \mathbf{X}^\*(q) - \mathbf{A} \cdot \overline{\mathbf{D}}, p < 0.5\\ \mathbf{X}^\*(q) + (\mathbf{D}\_{\overline{p}} + b \cdot l)e^{bl} \cos(2\pi l), p \ge 0.5 \end{cases} \tag{12}$$

where *p* is a random value between (0, 1).

If *p* is less than 50%, encircling prey is chosen. Otherwise, the Archimedean spiral strategy is chosen, and the whale position is updated by (11). In encircling prey, the prey is selected to pass by **A**. If |**A**| < 1, the optimal individual in the population is selected as the prey to execute encircling prey by (4) and (5). If |**A**| ≥ 1, an individual is randomly selected as the prey, and the positions of population are updated by (8) and (9). By following the above three mechanisms, the whale swarm can gradually approach the optimal position.

**Figure 3.** The curves of an Archimedes spiral and a logarithmic spiral.

#### 3.2.5. Refraction Learning

RL is an inverse solution solving method that introduces the refraction principle of light into opposition-based learning (OBL). The superior inverse solution is screened out by OBL, and it is used to replace the existing solution for the next iteration. This strategy has a higher probability of reaching the global optimal solution [37]. However, the conventional OBL is not adjustable in the search range and cannot prevent the inverse solution from falling into the local optimal solution. To fix the defect, a scaling factor is introduced to control the search range according to the refraction principle.

Figure 4 shows the process of RL, where the center point of the search interval [*x*min, *x*max] is denoted as *O*. The light propagates from *P*(*x*, *y*) to the origin *O* in one medium and enters another medium. Refraction occurs after it enters a new medium, and the refraction point is *<sup>P</sup>*ˆ∗. The distance <sup>|</sup> −→*PO*<sup>|</sup> from the incident point to the origin *<sup>O</sup>* is *<sup>h</sup>*, and the distance | −−→ *OP*ˆ∗| between the origin *<sup>O</sup>* and the refraction point is *<sup>h</sup>*∗. According to Shell's Law, the relationship between the angle of incidence *α* and that of refraction *β* satisfies the shell formula: *n* = sin *α*/ sin *β*, where *n* is the refractive index.

**Definition 1.** *The projection x*ˆ<sup>∗</sup> *of the refraction point P*ˆ<sup>∗</sup> *on the x-axis is defined as the inverse solution of x based on the RL.*

According to the shell formula, we find

$$\begin{aligned} \sin \alpha &= \left( (\mathbf{x}\_{\min} + \mathbf{x}\_{\max}) / 2 - \mathbf{x} \right) / h \\ \sin \beta &= \left( \mathbf{\hat{x}}^{\*} - (\mathbf{x}\_{\min} + \mathbf{x}\_{\max}) / 2 \right) / h^{\*} \end{aligned} \tag{13}$$

Assuming *k* = *h*/*h*∗, (13) can be rewritten as

$$\pounds^\* = (\mathbb{x}\_{\min} + \mathbb{x}\_{\max})/2 + (\mathbb{x}\_{\min} + \mathbb{x}\_{\max})/(2kn) - \ge /kn \tag{14}$$

The inverse solution *x*ˆ∗ can be solved using (14), where *k* is the scaling factor, which is used to adjust the positions of the inverse solutions on the *x*-axis. *n* is a fine-tuning factor. It can be seen that, assuming *n* = 1, *k* = 1, (14) can be simplified as *x*ˆ<sup>∗</sup> = −*x*, which is the expression for OBL. It can be seen that OBL is a special case of RL. The refractive index of light from air to water is about 0.75. Assuming that *n* = 0.75 is fixed, the change in incident angle causes the change in the projection point *x*∗.

**Figure 4.** The process of refraction learning.

The refraction point can be moved from point *P*ˆ<sup>∗</sup> to point *P*ˆ<sup>∗</sup> <sup>1</sup> by adjusting *k*, which is closer to the optimal point. If the *P*ˆ<sup>∗</sup> <sup>1</sup> does not reach the optimal value, the optimal solution can be reached by adjusting the incidence angle and moving the refraction point to the *P*ˆ∗ <sup>2</sup> . In order to improve the projection range of the inverse solution, we set *k* in a linearly decreasing trend:

$$k = k\_{\rm min} + (k\_{\rm max} - k\_{\rm min})t/t\_{\rm max} \tag{15}$$

The *N*-dimensional positions of a whale after position updating are denoted as:

$$\widetilde{\mathbf{X}} = \left( \widetilde{\mathbf{x}}\_0, \widetilde{\mathbf{x}}\_1, \dots, \widetilde{\mathbf{x}}\_{i\prime}, \dots, \widetilde{\mathbf{x}}\_{N-1} \right) \tag{16}$$

where *xi*,*<sup>j</sup>* ∈ [*xmin*, *xmax*], *i* = 0, 1, ... , *N* − 1. *xmin* and *xmax* are the minimum and maximum values in the search range, respectively. The inverse solution **X**ˆ <sup>∗</sup> obtained by RL is

$$\hat{\mathbf{X}}^{\*} = \left( \mathfrak{x}\_{0}^{\*}, \mathfrak{x}\_{1}^{\*}, \dots, \mathfrak{x}\_{i}^{\*}, \dots, \mathfrak{x}\_{N-1}^{\*} \right) \tag{17}$$

If **PSLL**(**X**<sup>ˆ</sup> <sup>∗</sup>) <sup>&</sup>lt; **PSLL**(**<sup>X</sup>** ), **<sup>X</sup>**<sup>ˆ</sup> <sup>∗</sup> is chosen as the best candidate solution. Otherwise, **<sup>X</sup>** is chosen as the best candidate solution.

#### *3.3. Off-Grid Arrays Design Based on GEWOA*

For off-grid arrays, considering the limit of the array element size, the minimum spacing *dc* is constrained. This means that the search range is not the whole aperture *L*. To solve this problem, the positions **D**(*d*0, *d*1, ... , *dN*−1) are decomposed. Since the positions of the two array elements *x*<sup>0</sup> and *xN*−<sup>1</sup> are fixed, only the positions **D**∗(*d*1, *d*2, ... , *dN*−2) of the *N* − 2 array elements need to be designed. Thus, a position decomposition method is applied to make **D**∗ continuous:

$$\mathbf{D}^\* = \breve{\mathbf{X}} + \begin{bmatrix} d\_c \\ 2d\_c \\ \vdots \\ (N-2)d\_c \end{bmatrix} = \begin{bmatrix} \breve{x}\_1 + d\_c \\ \breve{x}\_2 + 2d\_c \\ \vdots \\ \breve{x}\_{N-2} + (N-2)d\_c \end{bmatrix} \tag{18}$$

In (19), *<sup>x</sup>* <sup>1</sup> <sup>≤</sup> *<sup>x</sup>* <sup>2</sup> <sup>≤</sup>, ... , <sup>≤</sup> *<sup>x</sup> <sup>N</sup>*−2,*d*<sup>1</sup> <sup>&</sup>gt; *dc*, *dN*−<sup>1</sup> <sup>&</sup>lt; *<sup>L</sup>* <sup>−</sup> *dc*. The solution for the array position **D**<sup>∗</sup> is converted to the solution for **X** by the the position decomposition method. Here, the search range of *N* − 2 array elements becomes *L* − 2*dc*. In addition, (*N* − 3)*dc*

needs to be subtracted because the spacing between adjacent array elements is not less than *dc*. Therefore, the search aperture becomes *L*∗:

$$L^\* = L - 2d\_c - (N - 3)d\_c = L - (N - 1)d\_c \tag{19}$$

Through the above position decomposition, the discontinuous problem can be converted into a continuous problem. Then, GEWOA can be performed to achieve off-grid array design (Algorithm 1).


The algorithmic details of the proposed method for off-grid array design based on GEWOA are summarized in Algorithm 2. First, position decomposition is performed to obtain the search range *L*∗. Secondly, GEWOA is performed to search the optimal position. In the initial stage of GEWOA, the population initialization is performed by using chaos sequences to obtain **X** , and the maximum number of iterations *tmax* is set. The optimal whale in the initialized population is selected by (2), and then the iterative optimization is performed.

The parameters **a**,**r**, *l*, *p* of each whale are set, and **A**, **C** are calculated for choosing different strategies of position updating. If *p* > 0.5, ASS for position updating is performed. Else, if |**A**| < 1, a shrinking encircling mechanism is executed. Otherwise, a random search is used. After the position updating, the inverse solution **X**ˆ <sup>∗</sup> is generated by using RL. The best solution is selected from the new offspring and its inverse solution. The iterative optimization process is continued until the convergence condition is satisfied or the maximum number of iterations is reached. Finally, the optimal array positions are output by (19).

#### *3.4. On-Grid Arrays Design Based on GEWOA*

For on-grid arrays, the positions are discrete on half wavelength grid points. Therefore, the V-shaped transfer function is used to discretize the output of GEWOA. The V-shaped transfer function maps the whale position to the [0,1] interval. When the mapped values fall in different intervals, discrete values are selected. The expression of the V-shaped transfer function is:

$$V(\mathbf{x}\_{\vec{ij}}(t)) = |\frac{e^{\mathbf{x}\_{\vec{ij}}(t)} - 1}{e^{\mathbf{x}\_{\vec{ij}}(t)} + 1}|$$

$$\mathbf{x}\_{\vec{ij}}' = \begin{cases} G & 0 \le V(\mathbf{x}\_{\vec{ij}}(t)) < r\_1 \\ G - 1 & r\_1 \le V(\mathbf{x}\_{\vec{ij}}(t)) < r\_2 \\ \dots & \\ 0 & r\_M \le V(\mathbf{x}\_{\vec{ij}}(t)) < 1 \end{cases} \tag{20}$$

where *<sup>x</sup>* is the converted discrete value. *<sup>r</sup>*1,*r*2, ... ,*rk*−<sup>1</sup> are random numbers between (0, 1), which satisfy 0 < *<sup>r</sup>*<sup>1</sup> < *<sup>r</sup>*<sup>2</sup> < ... < *rk*−<sup>1</sup> < 1. When 0 ≤ *<sup>V</sup> xij*(*t*) < *r*1, the whale position is mapped to *G*. *G* = **floor**(2 ∗ *L*/*λ*) is the largest grid point (integer multiple of half wavelength), where **floor**() is a function to find the nearest integer less than or equal to the parameter. The algorithmic details of the proposed method for on-grid array design based on GEWOA are summarized in Algorithm 2.


#### **4. Experiment**

In this section, first, two kinds of numerical experiments on sparse-sensor-array-design tasks are performed to validate the effectiveness and superiority of the proposed method. (1) For off-grid array design, the proposed three-step enhanced global search strategy in GEWOA is simulated step by step to test the effectiveness of each strategy. Furthermore, the designed results are compared with those of WOA.

(2) For on-grid array design, the proposed GEWOA is compared with some other intelligence algorithm-based design methods, including GA, PSO, and WOA, in terms of PSLL. Then, the array designed using the proposed method is verified in the marine environment, and the beampattern obtained by the array based on GEWOA is analyzed.

The basic design requirements in numerical experiments are as follows. The array aperture is 200 m, the number of array elements is 40, and the array operating frequency is 300 Hz with a wavelength of 5 m. The objective function is shown in (2), where PSLL is chosen as the evaluation criterion. The number of whales in the population *M* = 30, the dimension of search *N* = 40, the number of iterations *tmax* = 1000.

The minimum spacing of the array *dc* is set to 2 m for off-grid arrays, and the array elements are set at half-wavelength grid points for on-grid arrays. In addition, we also conduct numerical experiments for large-scale low-frequency sensor array design, and the design conditions are as follows: the array operating frequency is 300 Hz, the number of array elements is 256, and the array aperture is 1280 m.

#### *4.1. Results of Off-Grid Arrays*

#### 4.1.1. Performance of Chaotic Initialization

In this part, different chaotic initialization methods are used to design the off-grid sensor array, and the convergence speed and beampattern performance are analyzed. Five chaotic mapping functions, including Logistic, Henon, Kent, Sin, and Tent, are tested under two termination conditions.

Under the first termination condition, which depends on the maximum number of iterations *tmax*, the search accuracy of five chaotic initialization methods is compared. The parameter settings and the obtained PSLLs are shown in Table 2, where the PSLL is the minimum value among 10 experiments. The beampatterns of the obtained array are shown in Figure 5. The iterative process of five chaotic initialization methods are shown in Figure 6. It can be seen that the search accuracy is improved after chaotic initialization compared to the original WOA (denoted as Origin in Figures 5 and 6). In addition, the Kent initialization method yields the array with the lowest PSLL.

**Table 2.** PSLLs and the numbers of iterations for five chaotic initialization methods.


**Figure 5.** Beampatterns of the obtained array based on five chaotic initialization methods.

Under the second termination condition that the algorithm stops if the PSLL remains unchanged beyond 100 iterations, the convergence speed of five chaotic initialization methods is compared. The numbers of iterations are shown in Table 2. It can be seen that Origin<Kent<Sin<Tent<Logistic<Henon in terms of the convergence speed.

Here, the Kent initialization method is adopted in GEWOA because it has the best search accuracy and faster convergence speed compared with the original method. Considering that the search accuracy is the most important criterion for array design, the selection is reasonable.

**Figure 6.** The iterative process of five chaotic initialization methods.

#### 4.1.2. Performance of Archimedean Spiral Strategy

In this part, we compare three search path planning strategies, including the original spiral position update, ASS, and Levy flight strategy. The Archimedean spiral spacing is set to 10 m. The beampatterns and the iterative process are shown in Figures 7 and 8, respectively. It can be seen that the PSLL of ASS is −21.02 dB, while the PSLL of Levy flight strategy is −20.76 dB. After adding ASS, the PSLL is 2.37dB lower than the original strategy. ASS improves the original strategy, and has a lower PSLL than the Levy flight strategy.

**Figure 7.** Beampatterns under three search path planning strategies.

#### 4.1.3. Performance of Refraction Learning

On the basis of Kent chaotic initialization and ASS, RL is used to further optimize the offspring to obtain an array with the lower PSLL. The setting parameters of RL are *kmax* = 1.2, *kmin* = 0.8, *tmax* = 10. The beampatterns of 40-element off-grid arrays and the iterative process before and after RL optimization are shown in Figures 9 and 10. It can be seen that the PSLL of array without RL optimization is −21.02 dB, and the PSLL of array with RL optimization is −22.47 dB. Compared with WOA, the PSLL obtained by GEWOA is reduced by 3.81 dB.

**Figure 8.** The iterative process under three search path planning strategies.

**Figure 9.** Beampatterns with and without RL.

**Figure 10.** The iterative process with and without RL.

In order to test the influence of RL parameter on the search accuracy, experiments with different RL parameters are performed. The three sets of RL parameters are (1) *k* = 1; (2) *kmax* = 1.2, *kmin* = 0.8, *tmax* = 5; and (3) *kmax* = 1.2, *kmin* = 0.8, *tmax* = 10. The array beampatterns and iterative process are shown in Figures 11 and 12, where the PSLLs of three sets of RL parameters are −21.47, −22.30, and −22.47 dB. The PSLLs of the designed arrays improve as the search range of the inverse solution is expanded.

**Figure 11.** Beampatterns of off-grid arrays (40-element) under three sets of RL parameters.

**Figure 12.** The iterative process under three sets of RL parameters.

#### *4.2. Results of On-Grid Arrays*

For on-grid arrays, the grids are set according to the half wavelength *λ*/2 = 2.5 m, and the number of grid points is 81. The proposed GEWOA discretized by a V-shaped transfer function is used to design on-grid arrays. In order to validate the superiority of the proposed GEWOA, three intelligence algorithms, including GA, PSO, and WOA, are compared with GEWOA.

The beampatterns of 40-element on-grid arrays and iterative process are shown in Figures 13 and 14, respectively. The PSLLs of GA, PSO, WOA, and GEWOA are −14.90, −15.10, −15.66, and −17.14 dB, respectively. It can be seen that the on-grid array designed using GEWOA obtains the lowest PSLL, which proves the advantage of the GEWOA in terms of the global search capability.

The proposed array design method, including the off-grid array design method and on-grid array design method, are successfully applied to the large-scale low-frequency arrays (256-element) design, and the obtained arrays have low PSLL. The beampatterns of the designed arrays and iterative process are shown in Figures 15 and 16, respectively. The results show that the PSLLs of the off-grid array and the on-grid array obtained by the proposed methods are −27.71 and −21.47 dB, which verifies the feasibility of the algorithm in large-scale low-frequency array design.

**Figure 13.** Beampatterns of on-grid arrays (40-element) designed using four intelligence algorithms.

**Figure 14.** The iterative process of four intelligence algorithms.

**Figure 15.** Beampattern of the on-grid large-scale low-frequency senor array (256-element) obtained using the proposed method.

**Figure 16.** The iterative process of the on-grid large-scale low-frequency senor array (256-element) based on the proposed method.

In order to compare the performance of the proposed method with typical configuration arrays, such as coprime arrays and nested arrays, three different arrays were designed at the same aperture and the same number of elements. The 40-element on-grid array was designed using the proposed method. The nested array has three levels of nesting. The first level is a 20-element uniform array with half-wavelength element spacing. The second level of nesting is a 10-element uniform array with a wavelength element spacing. The third level of nesting is a 10-element uniform array with two times wavelength element spacing.

To ensure the same aperture, the 40-element coprime array was designed, and its elements were arranged on three and four times half-wavelength grids. Figure 17 shows the beampatterns of the three arrays. It can be seen that the PSLL of the array designed using GEWOA are −14.79 dB, and the PSLLs of the coprime array and the nested array are −5.82 and −6.02 dB, respectively. The main lobe widths of three arrays are approximately the same. However, the sparse array designed using GEWOA has a lower PSLL, which proves the advantage of the proposed method compared with the determined methods.

**Figure 17.** Beampatterns of the coprime array, nested array, and GEWOA array.

To summarize, the above experimental results show that the proposed GEWOA with a three-step enhanced global search strategy had a good global search capability, and the array designed using GEWOA had a lower PSLL when compared with the other algorithms and deterministic methods. In addition, the proposed array design method was successfully applied to the large-scale low-frequency array design, and the obtained array had a low PSLL.

#### *4.3. Discussion and Analysis*

To demonstrate the superiority of the proposed method, we compared the performance of the sparse array designed using GEWOA with those of arrays obtained by traditional deterministic and hybrid methods, including the GWO array [23] and the different coarray (DC) designed using the method in [11]. Array design via these methods was performed at the same aperture.

The number of array elements is 40, and the array aperture is 200 m. The beampatterns of the designed arrays are shown in Figure 18, and their PSLLs are shown in Table 3, which indicates the lowest PSLL of our method. To further illustrate the repeatability of our method, we conducted 50 Monte Carlo simulations of the algorithms under the same conditions.

The results show that our method had a 96% probability of making PSLL below −17 dB. We compared the time efficiency of these methods on a CPU Intel E3-1220 v6 @ 3.00 GHz, and −15 dB was chosen as the convergence threshold. The computation times of GEWOA, GWO, and DC were 14.39, 20.54, and 122.41 s, and it can be seen that our method achieved the best computational efficiency.

**Figure 18.** Beampatterns of sparse arrays designed using different methods.



#### 4.3.1. Mutual Coupling Analysis

Due to the limitations of the natural frequencies, the mode shapes of frequencies, and inertial effects, the low-frequency acoustic sensors, such as transducers are usually large [43]. For engineering implementation, the minimum spacing of the array elements is constrained. The distance between adjacent element of the off-grid array designed using GEWOA can be less than half the wavelength, which results in the array element falling in the radiated sound field of that adjacent element.

The mutual coupling phenomenon occurs between adjacent array elements. This increases with the decrease of array element spacing, and directly affects the setting of the minimum spacing of the array element. According to the mutual radiation impedance model [44], the mutual radiation impedance *XAB* of an array element in the radiation field of another array element is

$$X\_{AB} = \rho c Ska^2 / d\_{AB}. \tag{21}$$

where *ρ* is the medium density. *c* is the speed of sound. *S* is the surface area of the array element. *a* is the radius of the sensor element. *dAB* is the distance of two elements. *k* is the wave number assuming the resonant frequencies of the array element is 3.5 kHz. The radius of the array element is 0.4 m. The curve of the radiation impedance changing with the array element spacing is shown in Figure 19.

It can be seen that *d*/*λ* ≥ 0.4, mutual-impedance is less than its self-impedance. In the off grid array design, the operating frequency is 300 Hz, and the minimum spacing of array elements is set to 2 m due to engineering requirements. Since the minimum spacing is greater than 0.4 times the wavelength of the resonant frequency, the mutual radiation effect is less than the self-radiation, and it can be disregarded.

**Figure 19.** Curve of the radiation impedance changing with the array element spacing.

#### 4.3.2. Array Gain Analysis

Array gain is an important indicator of array performance. Assuming that the sparse arrays receive a 300 Hz sinusoidal signal with a duration of 1 s. The array gains of the 40-element off-grid array and on-grid array designed using GEWOA are calculated to verify the feasibility of the proposed method under the Gaussian white noise background. The theoretical value of the array gain for the 40-element array is 16.02 dB.

The change trend of the array gain with the input signal-to-noise ratio (*SNRin*) is shown in Figure 20. It can be seen that the array gain decreases at a lower input SNR. When the *SNRin* is −30 dB, the array gain of both arrays decreases by about 1.5 dB compared with the theoretical value. When the *SNRin* is higher than −20 dB, the array gain of both arrays decreases by, at most, 0.3 dB compared with the theoretical value, which proves the feasibility of the designed array based on the proposed method from the perspective of array gain.

#### 4.3.3. Position Uncertainty Analysis

To verify the effect of array performance when the array element position is varied, Gaussian perturbation is added to the element position of the on-grid array designed using GEWOA in Section 4.2 to compare the effect of the PSLL performance. The outer

layer of most fiber optic towing arrays is the polyurethane tubing, which is filled with oil. This connection limits the offset of the array element position. Assume that the position perturbation with Gaussian distribution <sup>N</sup> (0,( *<sup>λ</sup>* <sup>12</sup> )2), which can guarantee 99.73% probability within a quarter-wavelength error range.

**Figure 20.** Curve of the radiation impedance changing with the array element spacing.

The beampattern of the sparse array designed using GEWOA before and after adding the position perturbation is shown in Figure 21. It can be seen that the PSLL of the designed array is −17.58 dB, while the PSLL of the array is −15.87 dB in the case of perturbation of all array elements. We performed 100 experiments using Monte Carlo simulations. The average value of PSLL for the position perturbation case is 16.23 dB. A quarter-wavelength Gaussian position perturbation causes the PSLL of the array to degrade by about 1.35 dB.

**Figure 21.** The beampattern after position perturbation.

#### *4.4. Experiment in the Marine Environment*

To verify the effectiveness of the proposed method, the performance of the designed arrays in the marine environment was analyzed. Due to the long production cycle and high investment of large-scale low-frequency sensor arrays, 40-element on-grid arrays were designed for the experiments. We extracted 40 elements from an 81-element uniform array to obtain the required arrays. The positions of the arrays designed using GA, PSO, WOA, and GEWOA are shown in Figure 22.

All array elements are set on the grid points. The beampatterns of designed arrays are shown in Figure 23. The PSLLs and main lobe widths of the beampatterns are shown in Table 4. By analyzing the performance of the sparse-sensor array designed based on GEWOA in the marine environment and comparing with the performance of the arrays obtained by GA, PSO, and WOA, the effectiveness of the algorithm was verified. The experimental conditions were as follows: the experimental array was a uniform hydrophone array with a half wavelength of 0.4167 m as shown in Figure 24. The design frequency of the array was 1800 Hz, and the number of array elements was 81.

The experimental data were from the 2020 South China Sea Experiment. The designed arrays were used to receive the transmit signal. The received data of 40-element on-grid arrays obtained by GA, PSO, WOA, and GEWOA were extracted from the 81-element uniform array. The relative azimuth angle between the sound source and the array was 15◦. The frequency of the transmitted signal was 1800 Hz. Its duration was 1 s, and the transmission period was 40 s.

**Figure 22.** The positions of arrays designed using GA, PSO, WOA, and GEWOA.

**Figure 23.** Beampatterns of the on-grid arrays (40-element) designed using GA, PSO, WOA, and GEWOA.

**Figure 24.** Profile display of equipment in the marine environment.

The beampatterns were obtained by conventional beamforming processing of the direct wave signals. The PSLLs of the beampatterns were analyzed to evaluate the array performance. The beampatterns obtained by the 40-element on-grid array designed using GA, PSO, WOA, and GEWOA are shown in Figure 25. The PSLLs and main lobe widths of the beampatterns in marine environment are shown in Table 4. It can be seen that the PSLLs of the beampatterns in the marine environment are degraded, due to the deviation of the array position from the designed ideal position.

The main lobe widths vary relatively little. The PSLL of the array designed using GEWOA is −17.10 dB in the marine environment. According to Table 4, it can be seen that compared to GA, the main lobe width expanded by nearly 0.06 degrees, which was only increased by 4.5% in the marine experiments. However, the PSLL decreased by 3.9 dB, which is a decrease of nearly 29.5%.

Compared with WOA, the main lobe width expanded by 0.01 degree, which is an increase of only 0.7%, while the PSLL was reduced by 2.5 dB, which is a decrease of 17.1%. The PSLL of GEWOA was 3.9 dB lower than that of the traditional GA algorithm with a small change in the performance of the main lobe width. In comparison with PSLLs obtained using other arrays, the array designed using GEWOA had the lowest PSLL.

**Figure 25.** Beampatterns of the on-grid arrays (40-element) designed using four intelligence algorithms in the marine environment.


**Table 4.** The PSLLs of arrays designed using different algorithms.

#### **5. Conclusions**

In this paper, we proposed a novel approach based on GEWOA for sparse sensor array design to suppress PSLL under spacing constrains. A three-step enhanced global search strategy was introduced into GEWOA to improve the global search capability. In the initial stage, chaotic initialization was embedded in GEWOA to enhance the ergodicity of the algorithm.

In the search stage, the conventional spiral strategy was replaced by the ASS strategy to avoid premature algorithm results. In the offspring selection stage, RL was used to obtain the inverse solution of the offspring, which prevents falling into local optima. Moreover, in order to solve the adaptation problem for discrete array design based on GEWOA, a position decomposition method and a V-shaped transfer function were introduced into the off-grid arrays and on-grid arrays, respectively.

The effectiveness and superiority of the proposed method were validated by experiments on large-scale low-frequency sparse sensor array design tasks. The experimental results show that the proposed GEWOA with a three-step enhanced global search strategy had a good global search capability, and the array designed using GEWOA had the lowest PSLL compared with other intelligent algorithms. Additionally, the array designed using the proposed method was further verified in the marine environment, where the proposed GEWOA still achieved the lowest PSLL. In the future, we will attempt to expand our application of the proposed method on 2-D or 3-D arrays.

**Author Contributions:** Conceptualization, L.W. and Q.W.; methodology, L.W.; software, L.W.; validation, L.W., H.Z. and Q.W.; formal analysis, L.W.; investigation, L.W.; resources, Q.W.; data curation, H.Z.; writing—original draft preparation, L.W.; writing—review and editing, H.Z. and Q.W.; visualization, L.W.; supervision, H.Z.; project administration, H.Z.; funding acquisition, H.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Maurizio Soldani \* and Osvaldo Faggioni**

Istituto Nazionale di Geofisica e Vulcanologia, Via di Vigna Murata 605, 00143 Roma, Italy

**\*** Correspondence: maurizio.soldani@ingv.it

**Abstract:** Sea level changes in coastal areas significantly influence port activities (e.g., the safety of navigation). Along Italian coastlines, sea level variations are mainly due to astronomical tides (well known, due to gravitational attraction between Earth, Moon and Sun); however, during the last fifteen years, a high number of "anomalous" tides has been observed: the study of the phenomenon has allowed to attribute its cause to variations in atmospheric pressure (the so-called meteorological tides: sea level drops when atmospheric pressure increases and vice versa); the statistical analysis of acquired data made it possible to evaluate the hydrobarometric transfer factor (a local parameter which represents the correlation between atmospheric pressure changes and consequent sea level variations): it was found that it is usually much larger within gulfs or port basins than offshore areas, where a pressure change of 1 hPa results in a sea level variation of about 1 cm; the statistical analysis described in the following, and aimed at correctly estimating the hydrobarometric transfer factor in harbors, can play a fundamental role in optimizing the management of port waters: its results allow to forecast meteorological tides and therefore future sea level (and depth) variations in a given port basin. The results of the study conducted in the port of La Spezia (North Western Italy) are presented here, together with possible applications on port activities and harbor water management.

**Keywords:** meteorological tides; marine environmental monitoring; sea level forecasting; harbor water management; port navigation safety

#### **1. Introduction**

The knowledge of sea level fluctuations in coastal areas is fundamental to increasing the safety of people who work in ports or onboard ships, to best managing port activities, and to reducing economic losses and environmental damage (improving the safety of navigation, optimizing ships' cargo and mooring, managing the refloating of stranded vessels, dimensioning maritime works, planning dredging activities, checking the chemical and physical parameters of water), as well as for civil protection purposes (mitigating the risk of flooding at the mouth of rivers) [1–7].

Sea level variations along Mediterranean coastlines are mostly due to astronomical tides (up and down motions caused by the gravitational attraction between the earth, moon and sun, then periodic and predictable deterministically), in particular to the diurnal and semi-diurnal components, shown in Table 1 [8,9].

However, a large number of anomalous tides have been observed over the last fifteen years inside many Italian harbors. The study of the phenomenon allowed us to associate these events to changes in atmospheric pressure above the sea basin under examination; in particular, sea level lowers/rises following an increase/decrease in atmospheric pressure (good/bad weather) and then in the weight of the overlying air column [10–17]: these low frequency oscillations, called meteorological tides, represent the geodetic adjustment (Newtonian compensation) of sea surface (induced effect), which compensates for atmospheric pressure perturbations (inducing cause), as shown in Figure 1 [18–23].

**Citation:** Soldani, M.; Faggioni, O. Observing Meteorological Tides: Fifteen Years of Statistics in the Port of La Spezia (Italy). *Appl. Sci.* **2022**, *12*, 12202. https://doi.org/10.3390/ app122312202

Academic Editors: Enjin Zhao, Hao Qin and Lin Mu

Received: 24 October 2022 Accepted: 24 November 2022 Published: 29 November 2022

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).


**Table 1.** The diurnal and semi-diurnal tide components and their periods expressed in hours (from [9]).

**Figure 1.** An increase in atmospheric pressure inducing a low meteorological tide (modified from [19]).

The evidence of the phenomenon was observed by means of ISPRA's (Italian Institute for Environmental Protection and Research) meteo-mareographic station located in La Spezia harbor (Italy); for example, on 13 November 2020 at 10:15 UTC (Universal Time Coordinates) and on 8 December 2020 at 09:25 UTC, the same astronomical tide was present (nearly 20 cm), but a pressure decrease of nearly 27.1 hPa resulted in an increase of about 34 cm in the sea level acquired.

In many Italian ports, we have verified that the hydrobarometric transfer factor, a parameter that represents the correlation between the atmospheric pressure variation and the consequent change in sea level, assumes often much larger values (even double) compared with offshore areas, where, as is well known, 1 hPa (about equal to 1 mBar) of atmospheric pressure variation corresponds to approximately 1 cm of change in sea level (the so-called inverted barometer effect); so, within harbors, a few hPa of decrease/increase in atmospheric pressure can cause several cm of sea level rise/fall; in fact, a sea basin behaves in the same way as a semi-constrained domain: by hindering the horizontal movement of the water mass, it amplifies its vertical displacement. Therefore, meteorological tides can cause exceptional changes in sea level within a port basin if they occur in-phase with astronomical ones.

To be able to forecast sea level in harbors, it is therefore necessary to estimate the correlation between atmospheric pressure and sea level, represented by the hydrobarometric transfer factor; it depends on a series of local parameters (first of all the morphology of the basin examined and the atmospheric dynamics over it), so it is not described by a deterministic law valid everywhere but obtained port by port through a statistical analysis of local data acquired [24–31].

In this work, we present the analysis carried out starting from data acquired since 2006 by ISPRA's monitoring station located inside the port of La Spezia (Eastern Ligurian Sea, Italy); based on the measurements of atmospheric pressure and sea level, the statistical analysis described below aims at estimating the hydrobarometric transfer factor in La Spezia harbor.

We also describe a prototype application developed to provide support to local authorities and port communities, in order to improve port navigation safety (obviously, the low tide hinders the port navigation, while the high tide facilitates it): based on the sea level measured or forecasted, the application updates water depths inside a port basin ("real" port bathymetric map) and detects, by means of a simple and intuitive graphic interface, hazardous areas for a certain ship moving inside the harbor at a given moment [32–38].

#### **2. Materials and Methods**

The starting point of this study is the monitoring of environmental parameters in the port of La Spezia, performed by means of the meteo-mareografic station working in the position 09◦51 27.52 E, 44◦05 47.79 N (see Figure 2) and belonging to the National Tidegauge Network managed by ISPRA [39–45].

**Figure 2.** The position of the monitoring station (pictures from Google Earth): (**a**) inside the port of La Spezia; (**b**) within the Eastern Ligurian Sea, Italy.

The instrumentation used consists of a hydrometer and a barometer; the first measures the sea level on the basis of the round trip time taken by a sequence of radar pulses sent from the air towards the sea surface; since the speed of propagation of electromagnetic waves in the air, the round trip time of the pulses, and the position of the radar transducer are known, sea levels are calculated; a typical sampling interval to acquire tide data is at least 10 min.

The barometer measures the atmospheric pressure by evaluating the deformations undergone by a silicon capacitive transducer; when the atmospheric pressure changes, the distance between the two plates and therefore the electrical capacitance also varies; measurements are typically made hourly (atmospheric pressure is characterized by very slow variations).

The date and time are expressed in UTC, while the sea level refers to IGM's (Italian Military Geographic Institute) 0 level. The data used in this work and further information about the ISPRA's monitoring station are available on the website www.mareografico.it (data from 2006 to 2009 are not available on the website; they have been kindly provided by ISPRA) (accessed on 28 April 2022).

Atmospheric pressure and sea level measurements (one sample every hour with resolutions equal to 10−<sup>1</sup> hPa and 1 cm, respectively) were compared with each other in

order to evaluate the hydrobarometric transfer factor, a local parameter that represents the correlation between the two quantities. Figure 3 shows the signals acquired between 28 May and 7 June 2009 inside the port of La Spezia.

**Figure 3.** Measurements carried out between 28 May and 7 June 2009 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level.

The sea level is the overlapping of different contributions, among which the main ones are astronomical and meteorological tides, as shown in the power spectral densities plotted in Figure 4, as regards the measurements performed in 2009: Frequency components due to diurnal and semidiurnal components are at approximately at 1.2 and 2.3 × <sup>10</sup>−<sup>5</sup> Hz respectively, while meteorological components are characterized by lower frequencies (slower oscillations) related to atmospheric pressure spectrum (DC components are not plotted because they do not carry useful information).

**Figure 4.** Power Spectral Densities of measurements carried out in 2009 within the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level.

First of all, it was therefore necessary to filter acquired data in order to remove highfrequency components due to causes different than meteorological ones (there could also be a residual of the disturbance due to sea waves, although the hydrometer works inside a stillpipe): atmospheric pressure and sea level signals are subjected to low-pass filtering with

an appropriate cut frequency (10−<sup>5</sup> Hz), so that only the contributions of meteorological origin survive. Figure 5 shows the result of the low-pass filtering applied to the data shown in Figure 3.

**Figure 5.** Measurements carried out between 28 May and 7 June 2009 inside the port of La Spezia, after Low-Pass filtering: (**a**) atmospheric pressure; (**b**) sea level.

By examining filtered signals, it is evident that a decrease Δ*p* in atmospheric pressure equal to about 17.2 hPa causes an increase in low-frequency sea level Δ*h* (low meteorological tide) equal to about 35.4 cm. So, for the event analyzed, the hydrobarometric transfer factor *Jph* can be calculated as:

$$J\_{\rm ph} = \frac{\Delta h}{\Delta p} = \frac{35.4 \text{ cm}}{17.2 \text{ hPa}} = 2.1 \text{ cm} \ast \text{hPa}^{-1} \text{ } \tag{1}$$

which is more than double the offshore case; the gradients Δ*h* and Δ*p* are expressed in absolute values.

The analysis just described for a single case was repeated for all the events that occurred in the port of La Spezia since 14 March 2006 (the hydrometer of ISPRA's monitoring station has been working since 13 January 2006, but the barometer was not installed until two months later) in order to obtain an estimate of the hydrobarometric transfer factor for the water basin under examination, as described in the next paragraph.

#### **3. Results**

The analysis described in the previous paragraph was replicated for every significant hydrobarometric event occurred in the port of La Spezia from March 2006 (data from 24 January 2015 to 6 March 2019 are not available) to the end of 2021, to produce a fifteenyear statistics (this study was realized starting from the installation of ISPRA's monitoring station, in 2006).

The values of Δ*p*, Δ*h* and *Jph* estimated event by event are listed in the Appendix A, in Table A2. Only events with Δ*p* greater than about 5 hPa are taken into account. The highest observed value of meteorological tide in this period was 52.9 cm (between 28 February and 4 March 2020).

The acquisitions related to these events are shown, in Figures A1–A4 (in the Appendix A) and Figures S1–S45 (in the Supplementary Material); events occurring in the presence of wind were not considered in the statistical analysis, in order to exclude some phenomena due to causes such as storm surges, anyway not predominant in the site examined (the anemometer to acquire wind data has been working since 30 June 2010) and then not analyzed in this work; for example, the event shown in Figure 6 (from 28 September to

5 October 2020, Δ*p* = 20.4 hPa, Δ*h* = 43.5 cm) has been excluded from the statistics because there was a wind coming roughly from the South East (then from the mouth of the gulf), stronger than 10 m/s and persistent for more than one day.

**Figure 6.** Measurements carried out between 28 September and 5 October 2020 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level; (**c**) wind direction, clockwise from the North; (**d**) wind intensity; for (**a**,**b**) solid lines represent low-frequency components survived the Low-Pass filtering; dashed lines represent high-frequency components removed by the Low-Pass filtering.

First, the event shown in Figure S5 was removed from the statistics because its *Jph* values differ from the mean value by more than 3 times the standard deviation; for this reason, it is considered outlier (rare event).

After doing that, the mean value and the standard deviation of the statistical distribution are 1.9 cm ∗ hPa−<sup>1</sup> and 0.3 cm ∗ hPa<sup>−</sup>1, respectively.

Then, a first estimate of the hydrobarometric transfer factor for the port of La Spezia can be represented by its average: in this case the mean value is about double the offshore; after this, the mean hydrobarometric transfer factor can be used to forecast a future sea level variation Δ*h* starting from the measured atmospheric pressure change Δ*p*: it is sufficient to

multiply Δ*p* by *Jph* (taking into account that when the pressure goes up the level falls and vice versa); this corresponds to suppose a linear dependence:

$$
\Delta \hbar = J\_{\text{ph}} \ast \,\,\Delta p \,,\tag{2}
$$

as represented by the red straight line in Figure 7, whereas the black dots correspond to the measured pairs (Δ*p*, Δ*h*) listed in Table A2.

A better estimate of the relationship between atmospheric pressure and sea level gradients can be obtained by considering a non-linear law Δ*h* = *f*(Δ*p*); using the least squares method applied to the pairs (Δ*p*, Δ*h*) the fixed degree polynomial can be extrapolated that best fits the data cloud (the hydrobarometric transfer factor becomes a sequence of coefficients); for example, assuming a cubic dependence of Δ*h* on Δ*p* and imposing the passage to the origin (because Δ*p* = 0 implies Δ*h* = 0), this estimate of Δ*h* is obtained:

$$
\Delta h = J\_{\text{pli3}} \ast \,\Delta p^3 + J\_{\text{pli2}} \ast \,\Delta p^2 + J\_{\text{pli1}} \ast \Delta p \,\,\,\tag{3}
$$

where: *Jph*<sup>3</sup> = 0.001 cm ∗ hPa−3, *Jph*<sup>2</sup> <sup>=</sup> −0.05 cm ∗ hPa−2, *Jph*<sup>1</sup> = 2.3 cm ∗ hPa−1, as plotted in Figure 7 (green line).

The couples (Δ*p*, Δ*h*) measured from 2006 to 2021 in the port of La Spezia, linear and cubic approximations are shown in Figure 7, together with the comparison with the "inverted barometer effect" of the offshore case.

Once the change in atmospheric pressure is measured or predicted, the expected sea level variation can be derived from the Equation (3).

Obviously, the two different estimates (linear and cubic) lead to two different errors between the measured value and the estimated trends: starting from data in Table A2, an average error on the expected sea level (difference between forecast and measurement) of 11.3% was obtained in the case of linear approximation, while in the case of cubic function the mean error is 9.6% (against a greater computational load), which is satisfactory for our purposes.

Therefore, the fundamental role of the hydrobarometric transfer factor consists in converting a variation of atmospheric pressure into a forecasted sea level change (meteorological component) for a given port; finally, meteorological tides will have to be added (or subctracted) to the astronomical tides (predictable by means of tide charts) to obtain the forecasted sea level in that basin.

Contributions to changes in sea level due to other causes such as, for example, seiches, storm surges or wind effects, which are not predominant in the site examined, are not taken into account in this work.

It should be emphasized that the hydrobarometric transfer factor must be updated year after year and periodically recalculated, e.g., following events that modify the topography of the basin examined (dredging operations, coastal erosion, bottom subsidence, sediment deposition).

Moreover, a multi-decade statistics is necessary to examine any changes in the hydrobarometric transfer factor due to climate change, as well as the variation over the years of occurrence frequency of events observed [46,47].

#### **4. Discussion**

As seen in the previous paragraph, the hydrobarometric transfer factor is usually much greater within port basins than in offshore areas: A few hPa of atmospheric pressure variation can cause several cm of astronomical tide that, if in phase with the astronomical one, can generate anomalous sea level variations.

The knowledge of the hydrobarometric transfer factor allows for correctly estimating expected meteorological tides in harbors and, together with the joint prediction of astronomical components from tide charts, forecasting sea level within port basins, an aspect of fundamental importance to better managing port operations.

In fact, the monitoring/forecasting of sea level (and then of water depth) in coastal areas is extremely important for managing, planning, and optimizing:


as well as preventing the risk of flooding at the mouths of rivers (civil protection purposes) by providing early warning to the population involved.

For this reason, the results of this study can have important applications in coastal areas: A software tool has been developed by the research group to which the authors belong with the aim of providing useful operational support to port communities, local authorities, and decision makers; it dynamically updates the initial port bathymetric map (usually acquired through multibeam surveys and updated after changes in the harbor topography, e.g., following dredging operations) based on sea level measured or expected and, if the draught of a certain ship is known, identifies the permitted/alert/prohibited areas for that same ship at a given time. An intuitive graphical interface implements what are called "virtual traffic lights" by coloring the forbidden areas red, the alert zones yellow, and the allowed ones green, based on two thresholds that in their turn depend on the vessel draught. Red areas are those with depths less than the lower threshold (usually equal to the vessel's draught), green areas those with depths greater than the upper threshold (much greater than the vessel's draught); finally, yellow areas are the intermediate ones.

The application continuously recalculates the "real" bathymetry (water depth variable over time) of a harbor using sea level data acquired in real time (by downloading them from the monitoring station via an Internet connection), measured in the past and saved in a dataset (to analyze a posteriori past events of particular importance such as the stranding of ships) or forecasted in the future by means of the hydrobarometric transfer factor and tide charts, in order to signal in advance potentially dangerous areas and avoid critical situations induced by sea level changes for a given ship. It can be very useful, e.g., for supporting a certain ship in choosing the best route to follow or the best time to enter or leave a port or the most suitable quay to moor.

For example, in Figure 8, virtual traffic lights in the middle of the port of La Spezia on 30 June 2010 at 12:00 UTC are shown, when the tide gauge was measuring −0.54 m. Thresholds were chosen equal to 12 and 13.5 m, e.g., for a Panamax cargo ship, whose draught is about 12 m (areas with depth less than 12 m are prohibited, those with depth greater than 13.5 m are permitted, the others are warning areas).

**Figure 8.** Virtual traffic lights in the port of La Spezia on 30 June 2010 at 12:00 UTC (threshold levels 12 and 13.5 m).

Easting and northing are expressed in UTM (Universal Transverse Mercator) coordinates, zone 32T; grid spacing is 2 m; depth refers to IGM's 0 level, depth resolution 1 cm; bathymetric data are courtesy of the Port System Authority of the Eastern Ligurian Sea, La Spezia, Italy.

Instead, Figure 9 refers to 24 November 2010 at 18:00 UTC (sea level = 0.83 m), for the same ship (and then the same thresholds).

Note how, following the rise of 1.37 m in sea level, the forbidden area narrows and many warning positions become allowed, while a yellow waterway appears in the west side of the port.

In particular, the state of the position indicated by the mouse pointer (coordinates 567,956, 4,883,930 m) in the middle top of the map switches from alert to allowed, as indicated by the color of the traffic light in Figures 8 and 9, since its depth increases from 12.45 to 13.82 m.

It is worth highlighting that the increase in sea level is partly due to a 9.5 hPa fall in atmospheric pressure between 30 June 2010 at 12:00 UTC and 24 November 2010 at 18:00 UTC (from 1016.9 to 1007.4 hPa).

Obviously at the same instant, for a vessel with greater draught (for example 14 m for a container ship), thresholds would be higher, and forbidden areas would expand, as shown in Figure 10; the traffic light for the position indicated by the mouse pointer becomes red.

**Figure 9.** Virtual traffic lights in the port of La Spezia on 24 November 2010 at 18:00 UTC (threshold levels 12 and 13.5 m).

**Figure 10.** Virtual traffic lights in the port of La Spezia on 24 November 2010 at 18:00 UTC (threshold levels 14 and 15.5 m).

Therefore, this interface represents a useful tool for detecting potentially dangerous areas for a given ship at a certain moment.

Additionally other contributions due to different phenomena can be passed as input to the application, such as seiches, storm surges, or wind effects, which particularly in certain locations must be taken into account.

#### **5. Conclusions**

In Mediterranean harbors, tides are mainly due to astronomical and meteorological components; while the first ones are well known and predictable everywhere by means of a deterministic law (tide charts), the second ones (due to atmospheric disturbances) need more study: starting from monitoring environmental parameters in La Spezia harbor (Ligurian Sea, North-Western Italy), we performed a statistical analysis over the last fifteen years (starting from the installation of ISPRA's monitoring station in 2006) to evaluate the hydrobarometric transfer factor, which represents the correlation between atmospheric pressure (cause) and sea level variations (effect) and therefore can play a fundamental role in forecasting meteorological tides in harbors.

We found that the hydrobarometric transfer factor in La Spezia harbor is larger (sometimes more than double) than in offshore areas (where 1 hPa of atmospheric pressure gradient induces nearly 1 cm of sea level rise/fall); thus, some hPa of atmospheric gradient can induce several cm of sea level rise/fall.

In particular, we found that using *Jph* = 1.9 cm ∗ hPa−<sup>1</sup> results in an estimate of the average error of about 11.3%, while approximating the dependence of sea level on atmospheric pressure by means of a nonlinear (cubic) law reduces the mean error to about 9.6%. This error very rarely induces confusion between different colors along the edges of adjacent areas in the map representing virtual traffic lights described above.

Furthermore, meteorological tides can cause exceptional changes in sea level if they occur in conjunction with astronomical components: the observation of the phenomenon allowed us to highlight anomalou" tides, sea level changes that are very different from the expected astronomical tides.

The hydrobarometric transfer factor allows to forecast meteorological tides in the sea basin examined and then, simply by adding the contribution of the astronomical tides, to know in advance the sea level (and therefore water depth) expected in the near future; this can represent a useful tool for optimizing port activities and logistic operations by planning them in advance.

Since sea level changes in coastal areas affect in a relevant way the safety of port communities, a prototype application has been developed that updates the port bathymetric map based on sea level acquired in real time or forecasted, with the aim of planning and optimizing port activities or managing emergencies (it is also able to load old measurements stored in a dataset, e.g., to analyze accidents occurred in the past); the software tool classifies the port bathymetric map in forbidden (red), warning (yellow), or permitted (green) areas for a certain ship at a given moment, based on its draught (virtual traffic lights); the graphical interface is able to detect and easily signal hazardous situations in harbors in order to provide a useful support to decision makers (port authorities, coast guards), with the aim of increasing safety for people working in ports.

**Supplementary Materials:** The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/app122312202/s1, Figure S1. Measurements carried out between 21 and 27 February 2007 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 10 hPa, <sup>Δ</sup>h = 16 cm, Jph = 1.6 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S2. Measurements carried out between 25 February and 1 March 2007 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 9.8 hPa, <sup>Δ</sup>h = 16.4 cm, Jph = 1.7 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S3. Measurements carried out between 9 and 16 May 2007 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 5.9 hPa, <sup>Δ</sup>h = 16.1 cm, Jph = 2.7 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S4. Measurements carried out between 31 May and 2 June 2007 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 4.9 hPa, <sup>Δ</sup>h = 14.4 cm, Jph = 2.9 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S5. Measurements carried out between 1

and 4 June 2007 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 4.9 hPa, <sup>Δ</sup>h = 16.8 cm, Jph = 3.4 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S6. Measurements carried out between 3 and 10 August 2007 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 12.2 hPa, <sup>Δ</sup>h = 20 cm, Jph = 1.6 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S7. Measurements carried out between 20 and 27 October 2007 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 13.7 hPa, <sup>Δ</sup>h = 27.6 cm, Jph = 2 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S8. Measurements carried out between 3 and 13 April 2008 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 17.5 hPa, Δh = 40.2 cm, Jph = 2.3 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S9. Measurements carried out between 10 and 17 April 2008 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 16.9 hPa, Δh = 30.8 cm, Jph = 1.8 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S10. Measurements carried out between 26 August and 9 September 2008 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 7.8 hPa, Δh = 19 cm, Jph = 2.4 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S11. Measurements carried out between 26 November and 2 December 2008 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 28.9 hPa, Δh = 51.6 cm, Jph = 1.8 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S12. Measurements carried out between 29 January and 9 February 2009 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 28.2 hPa, Δh = 45.6 cm, Jph = 1.6 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S13. Measurements carried out between 28 May and 7 June 2009 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 17.2 hPa, Δh = 35.4 cm, Jph = 2.1 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S14. Measurements carried out between 6 and 18 September 2009 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 13.8 hPa, Δh = 24.6 cm, Jph = 1.8 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S15. Measurements carried out between 28 November and 1 December 2009 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 20.5 hPa, <sup>Δ</sup>h = 38.6 cm, Jph = 1.9 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S16. Measurements carried out between 18 and 20 February 2010 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 19.4 hPa, <sup>Δ</sup>h = 33.9 cm, Jph = 1.7 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S17. Measurements carried out between 4 and 15 June 2010 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 8.2 hPa, Δh = 16.4 cm, Jph = 2 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S18. Measurements carried out between 7 and 16 August 2010 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 7.8 hPa, Δh = 16.5 cm, Jph = 2.1 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S19. Measurements carried out between 27 October and 02 November 2010 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 17.3 hPa, Δh = 35.7 cm, Jph = 2.1 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S20. Measurements carried out between 26 January and 8 February 2011 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 18.9 hPa, Δh = 30.3 cm, Jph = 1.6 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S21. Measurements carried out between 18 and 28 June 2011 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 11.3 hPa, Δh = 24.1 cm, Jph = 2.1 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S22. Measurements carried out between 18 September and 2 October 2011 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 20.8 hPa, <sup>Δ</sup>h = 33 cm, Jph = 1.6 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S23. Measurements carried out between 21 and 27 October 2011 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 15 hPa, Δh = 33.7 cm, Jph = 2.2 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S24. Measurements carried out between 25 October and 3 November 2011 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 17.3 hPa, Δh = 28.8 cm, Jph = 1.7 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S25. Measurements carried out between 28 October and 7 November 2011 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 17.2 hPa, Δh = 29.4 cm, Jph = 1.7 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S26. Measurements carried out between 12 and 20 July 2012 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 10.5 hPa, Δh = 18.1 cm, Jph = 1.7 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S27. Measurements carried out between 8 and 16 October 2012 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 14.2 hPa, Δh = 23.4 cm, Jph = 1.6 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S28. Measurements carried out between 5 and 12 May 2013 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 5.7 hPa, <sup>Δ</sup>h = 17.3 cm, Jph = 3 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S29. Measurements carried out between 23 and 28 June 2013 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 9.3 hPa, <sup>Δ</sup>h = 17.2 cm, Jph = 1.8 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S30. Measurements carried out between 19 and 26 January 2014 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 17.2 hPa, <sup>Δ</sup>h = 35 cm, Jph = 2 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S31. Measurements carried out between 10 and 26 February 2014 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 23.6 hPa, <sup>Δ</sup>h = 41.4 cm, Jph = 1.8 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S32. Measurements carried out between 12 and 16 October 2019 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 12.3 hPa, <sup>Δ</sup>h = 20.5 cm, Jph = 1.7 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S33. Measurements carried out between 15 and 18 October 2019 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 9.3 hPa, <sup>Δ</sup>h = 15.5 cm, Jph = 1.7 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S34. Measurements

carried out between 17 and 21 October 2019 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 4.8 hPa, <sup>Δ</sup>h = 11.5 cm, Jph = 2.4 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S35. Measurements carried out between 19 and 23 October 2019 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 6 hPa, <sup>Δ</sup>h = 13.3 cm, Jph = 2.2 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S36. Measurements carried out between 19 and 25 November 2019 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 12.3 hPa, <sup>Δ</sup>h = 28.7 cm, Jph = 2.3 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S37. Measurements carried out between 27 November and 9 December 2019 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 21 hPa, <sup>Δ</sup>h = 39.2 cm, Jph = 1.9 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S38. Measurements carried out between 28 February and 4 March 2020 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 28.6 hPa, <sup>Δ</sup>h = 52.9 cm, Jph = 1.8 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S39. Measurements carried out between 31 May and 8 June 2020 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 20.3 hPa, <sup>Δ</sup>h = 38.2 cm, Jph = 1.9 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S40. Measurements carried out between 7 and 17 June 2020 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 6.3 hPa, <sup>Δ</sup>h = 11.6 cm, Jph = 1.8 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S41. Measurements carried out between 21 and 31 August 2020 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 15.3 hPa, <sup>Δ</sup>h = 26.3 cm, Jph = 1.7 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S42. Measurements carried out between 21 and 29 June 2021 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 7.4 hPa, <sup>Δ</sup>h = 1.4 cm, Jph = 1.9 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S43. Measurements carried out between 21 and 27 July 2021 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 8.5 hPa, Δh = 16.7 cm, Jph = 2 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S44. Measurements carried out between 31 July and 24 August 2021 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 14 hPa, Δh = 23.4 cm, Jph = 1.7 cm <sup>∗</sup> hPa<sup>−</sup>1). Figure S45. Measurements carried out between 18 and 25 September 2021 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 11 hPa, Δh = 18.8 cm, Jph = 1.7 cm <sup>∗</sup> hPa<sup>−</sup>1).

**Author Contributions:** Conceptualization, M.S. and O.F.; methodology, M.S. and O.F.; software, M.S. and O.F.; validation, M.S. and O.F.; formal analysis, M.S. and O.F.; investigation, M.S. and O.F.; resources, M.S. and O.F.; data curation, M.S.; writing—original draft preparation, M.S.; writing—review and editing, M.S. and O.F.; visualization, M.S. and O.F.; supervision, O.F.; project administration, M.S. and O.F.; funding acquisition, M.S. and O.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded within the framework of the MENFOR Project (MEteo-tide Newtonian FORecasting) by several Italian port authorities, in particular the Port Authority of La Spezia (now named the Port System Authority of the Eastern Ligurian Sea) for what we described in this article and from the European Union and Regional Government of Liguria (Italy) by means of the Regional Plan of Innovative Actions—European Funds for the Regional Development.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Meteo-mareographic data used in this article, acquired by means of the monitoring station in La Spezia belonging to ISPRA's National Tidegauge Network, are mostly available on the website www.mareografico.it (accessed on 28 April 2022); data from 2006 to 2009 are not available on the website; they have been kindly provided by ISPRA; additionally, archives of environmental data available on weather websites were consulted.

**Acknowledgments:** The authors wish to thank the Port Authority of La Spezia (now named the Port System Authority of the Eastern Ligurian Sea) for having kindly provided bathymetric data (a special thanks to D. Vetrala and I. Roncarolo) and D.A. Leoncini for his past contribution in the development of the software tool. Part of this research was carried out when the first author was at OGS—National Institute of Oceanography and Applied Geophysics (Trieste, Italy). Finally, the authors thank the anonymous reviewers whose comments and suggestions helped improve this work.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

**Table A1.** Estimated values for representative events occurred from 2006 to 2021 in the port of La Spezia.



**Table A2.** Estimated values for representative events occurred from 2006 to 2021 in the port of La Spezia.

Figures A1–A4 (in this Appendix A) and Figures S1–S45 (in the Supplementary Material) show data acquired during the events listed in Table A2; solid lines represent low-frequency components survived the Low-Pass filtering; dashed lines represent highfrequency components removed by the Low-Pass filtering.

**Figure A1.** Measurements carried out between 11 June and 01 July 2006 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 13.5 hPa, <sup>Δ</sup>h = 24 cm, Jph = 1.8 cm <sup>∗</sup> hPa<sup>−</sup>1).

**Figure A2.** Measurements carried out between 16 and 23 August 2006 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 10.3 hPa, <sup>Δ</sup>h = 21 cm, Jph = 2 cm <sup>∗</sup> hPa<sup>−</sup>1).

**Figure A3.** Measurements carried out between 12 and 26 October 2006 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 15.8 hPa, <sup>Δ</sup>h = 31.8 cm, Jph = 2 cm <sup>∗</sup> hPa<sup>−</sup>1).

**Figure A4.** Measurements carried out between 23 October and 12 November 2006 inside the port of La Spezia: (**a**) atmospheric pressure; (**b**) sea level (Δp = 23.8 hPa, <sup>Δ</sup>h = 45.3 cm, Jph = 1.9 cm <sup>∗</sup> hPa<sup>−</sup>1).

#### **References**

