Enhanced Mild-Slope Wave Model with Parallel Implementation and Artificial Neural Network Support for Simulation of Wave Disturbance and Resonance in Ports

Abstract

Ensuring sea surface tranquility within port basins is of paramount importance for safe and efficient port operations and vessels’ accommodation. The present study aims to introduce a robust numerical model based on mild-slope equations, capable of accurately simulating wave disturbance and resonance in ports. The model is further enhanced by the integration of an artificial neural network (ANN) to address partial reflection, and its efficiency is optimized by developing a parallel algorithm based on OpenMP, allowing for a reduction in the required simulation times for real port areas spanning several kilometers horizontally. Numerous numerical experiments focusing on wave reflection against a vertical wall were conducted to develop the ANN. This neural network was designed to determine the appropriate value of the eddy viscosity coefficient, a crucial parameter in the momentum equation of the mild-slope model, tailored to incident wave characteristics. The model’s validity was confirmed through rigorous validation against experimental measurements, covering wave disturbance, rectangular harbor resonance, and Bragg resonance. The model consistently demonstrated a more than satisfactory performance across all considered scenarios. In a practical application, the model was deployed in the Port of Rethymno, Crete Island, Greece, effectively capturing and describing dominant phenomena within the port area. The implementation of a parallel algorithm significantly reduced the simulation times by ~92%, compared to the serial algorithm, thereby enhancing the model’s efficiency and applicability in real-world port environments.

1. Introduction

Ports are vital in facilitating global trade and connecting nations across the world. They serve as gateways for imports and exports, constituting vital hubs for economic activity. To this end, ensuring port tranquility and safe berthing is of utmost importance since it guarantees the smooth flow of goods, minimizing delays and disruptions in supply chains. Efficient berthing operations enable the timely loading and unloading of cargo, enhancing productivity and reducing costs. Additionally, port tranquility ensures the safety of vessels, crew members, and the surrounding environment, minimizing the risk of accidents, collisions, or environmental hazards. Ultimately, ensuring low wave agitation levels and safe berthing is of paramount importance with strong implications for the economy’s and community’s safety.

Considering the above, the numerical modelling of wave agitation inside port basins is indispensable for engineers desiring to obtain predictions of wave heights at berthing positions to assess if the operability of the port is compromised. Several numerical models have been developed that are capable of carrying out this task, each one associated with specific advantages and limitations [1]. However, all models selected to carry out this task should be able to resolve a multitude of processes that drive wave transformation inside the port basin and the nearshore area, most notably wave diffraction and wave reflection.

When investigating wave agitation within port basins, researchers frequently utilize models that rely on solving the Boussinesq equations. These models leverage the advantageous features of the Boussinesq equations, such as their high level of nonlinearity and ability to address multiple significant wave phenomena like diffraction, refraction, reflection, and wave–structure interactions, among others. Building on the foundation laid by [2], numerous studies have significantly contributed to the development of Boussinesq-type wave models. These efforts aimed to enhance the models’ applicability to deeper waters, incorporate the effects of energy dissipation, and intensify the nonlinearity of the equations (e.g., [3,4,5,6,7,8,9]). While Boussinesq-type wave models offer a detailed description of wave transformation processes within port basins, their practical application is hindered by heightened computational resource requirements and frequent instability issues. The latter becomes particularly pronounced in areas characterized by rapidly changing bathymetry. Challenges stem from the higher order spatiotemporal derivatives in the governing equations, restricting the effective use of these models in real-world port and coastal environments.

An alluring alternative to Boussinesq-type models are those based on the mild-slope wave equation of [10], which describes the propagation of linear water waves in a mildly sloping bed. The authors in [11] formulated an equivalent version of the elliptic mild slope, while [12] introduced a set of three interconnected first-order hyperbolic equations. Building upon Copeland’s [12] work, [13] enhanced the model by incorporating harmonic time variations and implementing an Altering-Direction-Implicit (ADI) algorithm with variable time steps. Another approach involves the utilization of a time-dependent hyperbolic equation, as suggested by [14]. Extended mild-slope equations have been developed to account for the presence of steep or rapidly changing bottom topography. This extension was proposed by researchers such as [15,16,17]. The solution of the hyperbolic-form of the mild-slope wave equation is more efficient than the Boussinesq equations, with less strict stability restrictions. A crucial characteristic of hyperbolic approximation models is their ability to simulate all predominant phenomena nearshore, encompassing full wave reflection.

In addition to accounting for full reflection, which can be directly modeled, addressing partial reflection is crucial for accurately predicting wave disturbance, and, therefore, numerous research efforts have focused on this matter. The authors in [18] introduced a formulation utilizing the characteristic method to address partial wave reflection. Subsequently, [12,19] applied the numerical approach established by [18] to the hyperbolic mild-slope equation, specifically for addressing the partial reflection of water waves. Madsen [20] proposed a theoretical approach akin to the methodology employed by [21] in studying wave transmission through porous structures. The focus was on analyzing the reflection of linear shallow-water waves from a homogeneous porous wave absorber with a vertical face situated on a horizontal bottom. The reflection coefficient was derived as a function of various parameters, including those characterizing the incoming waves and the absorber properties—specifically, the friction factor, porosity, and the product of the wavenumber and the absorber width. Isaacson and Qu [22] presented a method based on linear diffraction theory for predicting wave agitation in harbors containing partially reflecting boundaries in the manner proposed by [23]. In their research, these boundaries are schematized as verticals, and partial reflection was introduced into the model by using a mixed-boundary condition instead of a full-reflection condition. Dalrymple et al. [24] explored the linear theory concerning waves impacting a vertically sided porous structure at an oblique angle. Considering both wave transmission and obliqueness, they derived an analytical solution for partial reflection and transmission. Their solution expressed reflection and transmission coefficients in relation to the inertia term, porosity, and friction factors. Karambas and Bowers [25] modeled the energy losses resulting from friction effects on the rough slope of a rubble mound structure by incorporating a friction area ahead of a vertically faced homogeneous porous structure within a hyperbolic linear wave model. Thereafter, [26] introduced a boundary condition for partial reflection in water waves. This equation found applications in the works of [27,28], and it is employed for addressing partial reflection in elliptic mild-slope equations (e.g., [29]). Chun et al. [30] developed a numerical method for the partial reflection and transmission of water waves in a hyperbolic mild-slope equation, by introducing a friction factor to the continuity equation for the treatment of partial reflection, and a formulation of the depth-integrated velocity for the transmitted wave was added to simulate wave transmission. A comparison of their numerical findings with the analytical solutions proposed by [31,32] showed good agreement. Son et al. [33] presented an energy-controlling boundary condition for partial wave reflection in the mild-slope equations. According to their research, energy dissipation is determined by the imaginary part of the wave number and the wave phase velocity. Partial reflection can be simulated by adjusting two variables: the width of an energy-controlling area and the imaginary wavenumber. Diaz-Hernandez et al. [34] incorporated partial reflection in their elliptic mild-slope model by implementing the equations for a complex constant, related to the complex reflection coefficient on the boundary, as proposed by [35]. The determination of wave reflection was achieved by adopting the formulation presented by [36], whilst a categorization of the waterfront types in the inner perimeter of the harbor should be accomplished. Karambas and Samaras [37] and Makris et al. [38] incorporated partial reflection in a hyperbolic mild-slope model by introducing an artificial eddy viscosity coefficient in the momentum equations to simulate energy loss in front of waterfronts. The values of the eddy viscosity coefficient were estimated from the method developed by [25], using the reflection coefficient values proposed by [39]. However, despite the time-efficient solutions provided by the hyperbolic mild-slope wave equations, the preprocessing procedures for determining the appropriate eddy viscosity coefficients often demand significant effort [40]. Moreover, when there is a requirement to simulate irregular waves approaching and penetrating a port, the necessary simulation times can be significant.

In the present paper, the main features and advancements of a hyperbolic mild-slope wave model, referred to as HMS hereinafter, designated to deal with partial wave reflection and wave disturbance and resonance in port basins are presented. A novel distinguishing quality of the model is that it is supported by an artificial neural network (ANN) that is tasked with estimating the appropriate eddy viscosity values depending on the desirable reflection coefficient to simulate, in a concise and reliable manner, partial wave reflection. Moreover, the model’s algorithm has been parallelized to substantially decrease the necessary simulation times within real coastal areas that extend across several kilometers in two horizontal dimensions. The model’s validation process involved thorough assessments against a diverse set of experimental measurements with very satisfactory results, rendering it capable of simulating irregular and regular wave propagation and transformation inside port basins. Furthermore, the model was applied in the Port of Rethymno, Crete Island, Greece, effectively simulating all the dominant phenomena within the port area.

2. Materials and Methods

2.1. Governing Equations and Model Features

The present model solves the hyperbolic mild-slope equation, describing the evolution, in time and space, of the free sea-surface elevation and the velocities. Based on the works of [12,19,37,41], the continuity and momentum equations for linear harmonic wave propagation in coastal waters of mildly sloping beds with energy dissipation terms are as follows:

$$\frac{\partial \zeta }{\partial t}+\nabla ·\left(\mathit{U}h\right)=-w_b\frac{\partial \zeta }{\partial t}$$

(1)

$$\frac{\partial \left(\mathit{U}h\right)}{\partial t}+\frac{c^2}{n}\nabla \left(n\zeta \right)=v_h{\nabla }^2·\mathit{U}-w_f\mathit{U}$$

(2)

where $\zeta$ is the sea-surface elevation; $h$ is the water depth; $\mathit{U}=(U,V)$ is the mean velocity vector with $U$ and $V$ being the mean horizontal velocities in the x and y directions, respectively; $c$ is the phase celerity; $n=(1/2+kh/\sinh kh)$; $k$ is the wave number; $v_h$ is the horizontal eddy viscosity coefficient responsible for replicating partial wave reflection; and $w_f$ denotes energy dissipation due to bottom friction, while $w_b$ denotes energy dissipation due to depth-induced breaking.

Energy dissipation $w_f=E_f/E$ ($E=H_{rms}^2/8$, and $H_{rms}$ is the root mean square wave height) due to bottom friction is given by the following [42,43]:

$$\begin{array}{c}{E}_f=\frac{1}{6\pi }\frac{c_{fw}/2}{g}{\left(\frac{\omega H}{\sinh kh}\right)}^3\text{for} \text{regular} \text{waves} \text{and} \hfill \\ E_f=\frac{1}{8\sqrt{\pi }}\frac{c_{fw}/2}{g}{\left(\frac{\omega H_{rms}}{\sinh kh}\right)}^3\text{for} \text{irregular} \text{waves}\hfill \end{array}$$

(3)

where $c_{fw}$ is a wave friction coefficient [44], $g$ is the acceleration of gravity, and $\omega$ is the angular frequency.

Energy dissipation $w_b=E_b/E$ due to depth-induced breaking is calculated as shown in [45]:

$$E_b=\frac{a}{4}{Q}_b{f}_p{H}_{max}^2$$

(4)

where $a$ is a calibration coefficient; $f_p$ is wave frequency at the spectral peak; and $Q_b$ is the fraction of breaking waves, obtained by solving the following relationship $\frac{1-Q_b}{\ln Q_b}=-{\left(\frac{H_{rms}}{H_{max}}\right)}^2$, assuming a Rayleigh distribution truncated at the breaker height $H_{max}$, calculated based on a Miche-like criterion [46].

The model also includes a nonlinear dispersion relationship by adopting the formulation of [47]:

$${\omega }^2=gktanh\left(kh\right)+F_1\frac{\left|A\right|}{h}+F_2{\left(\frac{\left|A\right|}{h}\right)}^2$$

(5)

with $F_1=g{k}^{2}h\frac{f_2\left(kh\right)}{{\cosh }^{2}kh}$, $F_2=g{k}^3{h}^{2}tanh\left(kh\right)\left(f_1\left(kh\right)D-\frac{{f_2}^2\left(kh\right)}{{\cosh }^{2}kh}\right)$, $D=\frac{\left(\cosh 4kh+8-2{\tanh }^{2}kh\right)}{8{\sinh }^{4}kh}$, $f_1\left(kh\right)={\tanh }^5\left(kh\right)$, and $f_2\left(kh\right)={\left[kh/\sinh kh\right]}^4$.

The numerical solution of the equations is accomplished by employing an explicit scheme implemented on a grid staggered between the cell values of the surface elevation and mean velocities. At a given spatial point i and time n, the partial differential equations are approximated on a staggered orthogonal grid. The calculation of the free surface elevation ($\zeta$) takes place at the center of each grid cell, while the velocities $(U,V)$ are estimated on the sides of the grid cell [41].

The harmonic wave generation is placed internally along any lateral boundary line by adding the following quantity to the water surface elevation at each time step [48]:

$${\zeta }^{*}=2\frac{H}{2}\cos \left(ky\sin \theta -\omega t\right)c\frac{dt}{dx}\cos \theta$$

(6)

where $\theta$ is the angle-of-incident wave direction.

Sponge layers are placed perimetrically on the numerical domain to absorb the wave energy travelling outwards by using an exponential damping factor as follows:

$$DF\left(x\right)=exp\left[\left(2^{-x/\Delta x}-2^{x_s/\Delta x}\right)\mathit{ln}b\right],0<x\le x_s$$

(7)

where $x_s$ is the width of the sponge layer, $b=1+r_s+{exp}^{-\left(1/r_s\right)}$, $r_s=10/N_s$, and $N_s$ represents the number of grid points in the sponge layer.

At the solid boundaries, which are treated by default as fully reflective by the model, the slip boundary condition is applied, setting the water particle velocity normal to the boundary as zero. For partial reflection, the eddy viscosity term is enabled.

The model can also simulate irregular waves by adopting the principle of linear superposition for the main variables $\zeta , U,$ and $V$. First, the input spectrum is discretized in distinct regular components with a specific wave height and period. For each of the wave components, Equations (1) and (2) are solved at each timestep, and subsequently, the variables are superposed in the entire domain. In this manner, the model treats the propagation and transformation of unidirectional irregular waves. If multidirectional waves are to be simulated, the input spectrum is further discretized in distinct regular components with a specific wave height, period, and direction, and the aforementioned methodology is applied. The root mean square wave height is calculated as $H_{rms}=2{\left(〈2{{\zeta }_{tot}}^2〉\right)}^{1/2}$ at each cell of the numerical domain, where angular brackets indicate the time-averaged quantity and ${\zeta }_{tot}$ the sum of the surface elevation of each component through the linear superposition. The model can represent the frequency spectrum by utilizing a JONSWAP or TMA [49] source function. These spectrum source functions are discretized by adopting the single summation model of [50]:

$$f_{mj}={\widehat{f}}_m-\frac{1}{2}\Delta f+\left(j-1+{RAN}_{mj}\right)\Delta f/N_{\theta }$$

(8)

$$a_{mj}=\sqrt{2\widehat{S}\left({\widehat{f}}_m,{\theta }_j\right)\Delta f\Delta \theta }$$

(9)

$$\Delta f=\frac{f_{max}-f_{min}}{M_f}, \Delta \theta =\frac{{\theta }_{max}-{\theta }_{min}}{N_{\theta }}$$

(10)

In the above equations, ${\widehat{f}}_m$ is the center frequency of the $m$th frequency band; ${(f}_{min},f_{max})$ and ${(\theta }_{min},{\theta }_{max})$ are the ranges of frequency and direction of the incident directional wave spectrum, respectively; $M_f$ and $N_{\theta }$ are the number of frequency and directional bands in the discretized directional spectrum, respectively; ${RAN}_{mj}$ is a randomly generated number within the range (0, 1), and it is integrated to introduce a random component to $f_{mj}$; $a_{mj}$ is the amplitude of each component; and $\widehat{S}\left(f,\theta \right)$ is the directional wave spectrum. It is noted that a random phase is considered for each individual component comprising the spectrum.

2.2. Proposed Methodology for Incorporating Partial Reflection

To accurately simulate the partial reflection of waves across diverse waterfront types, a methodology is outlined herein, structured around the following three steps (as illustrated in Figure 1):

(1)

An estimation of the incident wave characteristics on the waterfront and a determination of the anticipated wave reflection. The latter depends on a multitude of parameters [40], the most important of which are the incident wave height and period, the depth at the structure’s toe, and the type of the front (partially dissipative or fully reflective). An additional parameter that influences the wave reflection is if the incident sea-state is considered as regular or irregular [51]. For a preliminary estimation of the reflection coefficient, informative diagrams (e.g., [52]) may be used, while more precise determinations can be made using reflection coefficient formulas. A plethora of formulas for predicting the reflection coefficient have been developed, based on either limited datasets related to specific structures (e.g., [46,53,54,55,56]) or more extensive ones (e.g., [57,58,59]) that consider various types of structures derived from extensive experimental datasets. An even more accurate approach involves leveraging advanced tools based on machine learning techniques, such as the neural network developed by [60,61] for predicting wave reflection from coastal and harbor structures.

(2)

The anticipated reflection coefficient, determined in the previous steps, cannot be directly fed to the mild-slope model. To address this, we propose that the horizontal eddy viscosity coefficient, $v_h$, is responsible for replicating partial wave reflection by dissipating a portion of the incident wave energy. A relationship between the eddy viscosity coefficient and the anticipated reflection coefficient is essential, and it is not known a priori, as it depends on the incident wave characteristics, the depth at the structure’s toe, and the type of the waterfront. Therefore, as a second step, we suggest implementing an ANN (whose development is detailed in the following paragraphs). This ANN can provide the eddy viscosity coefficient based on the anticipated reflection coefficient and sea state characteristics.

(3)

The determined eddy viscosity coefficients are provided in a separate file, serving as input data for the mild-slope wave model, which ultimately executes the simulation. Within this file, the eddy viscosity coefficient uniformly assumes a zero value throughout the entire numerical domain, except for the zones corresponding to the waterfronts’ frontal areas where partial reflection is anticipated. In these specific regions, the values generated by the artificial neural network (ANN) are incorporated. It is noted that, in case of a sloping structure, the depth along the sloping region has to be set as constant and equal to the depth at the structure’s toe.

It is noted that the proposed methodology enables the simulation of partial reflection both from vertical fronts and sloping structures, such as breakwaters, either permeable or impermeable, rough or smooth, and with varying steepness. For the case of vertical fronts, the proposed methodology can be directly applied. However, in the presence of sloping structures, an additional step is necessary, adopting an approach similar to the one proposed by [25], further embraced by [37] and subsequently by [38]. Specifically, the bathymetric grid is modified at the front of the sloping structures. For a width of $L_e$ (a specific width of the cells seaward of the land cell), where the eddy viscosity coefficient is to be applied, the depth must remain constant to avoid simultaneous energy loss from the eddy viscosity and breaking and friction processes. Therefore, the bathymetric grid in front of the cells representing the shoreline of the sloping structure is adjusted, with a constant depth equal to the depth of the most seaward cell where the turbulent diffusion coefficient will be applied. With this approach, the energy loss due to breaking and friction as well as the reflected wave energy can be controlled, given that the desired reflection coefficient is known in advance.

2.3. Model Parallelization

Based on the above, it becomes evident that in engineering applications, particularly those involving coastal areas spanning multiple kilometers along two horizontal axes, the simulation of multidirectional seas with numerous components in both frequency and directional domains poses a significant computational challenge. This is primarily due to the necessity of solving the entire numerical domain for each component at every time step. Consequently, the motivation behind parallelizing the algorithm was to alleviate this computational burden and achieve a substantial reduction in the required simulation times.

In recent decades, there has been a notable emphasis on the development of parallel algorithms for simulating wave propagation. The aim is to accelerate the generation of results and enhance the applicability of numerical models in real coastal areas spanning several kilometers, all while mitigating the need for substantial CPU time and memory resources.

Sitanggang and Lynett [62] developed a parallel Boussinesq wave model utilizing the Domain Decomposition (DD) technique in conjunction with the Message Passing Interface (MPI). Employing the DD technique involves dividing the domain into multiple regions, each containing an overlapping area of ghost cells. Individual subdomains are assigned to separate processor cores, and MPI facilitates data exchange in the overlapping regions between neighboring processors. The parallelization of tridiagonal matrix systems, a challenging aspect of the serial scheme, was accomplished using a parallel matrix solver [63] with an efficient data transfer scheme. Similarly, [64] adopted an approach to parallelize the well-known Boussinesq-type model FUNWAVE-TVD, employing the DD technique. MPI, with non-blocking communication, was utilized for exchanging data in overlapping regions between neighboring processors. The solution of tridiagonal matrices was achieved through the implementation of the parallel pipelining tridiagonal solver, as described by [65]. Additionally, [66] provided insights into certain aspects of parallel implementations of the widely used spectral wave model SWAN [67] on distributed-memory architectures using MPI. The authors in [68,69] explored the parallel efficiency of both the OpenMP and MPI versions of SWAN. Their findings indicate that the MPI version was less efficient than the OpenMP version within a single node. Consequently, they proposed a hybrid version of SWAN that combines the efficiency of the OpenMP version on shared memory, with the capability of the MPI version to distribute memory across computational nodes. Rautenbach et al. [70] implemented and compared OpenMP- and MPI-distributing-memory architectures, focusing on the scalability of SWAN. Their results show that the MPI version outperformed the OMP version across three performance metrics: efficiency, the speed-up ratio, and the timesaving ratio. The SWASH model [71], a multi-layer nonlinear Shallow-Water-Equation solver, capable of simulating wave fields and rapidly varied flows in coastal waters, was parallelized using strip-wise grid partitioning and MPI. This involved dividing the computational domain into strips along the two horizontal axes, with each strip being assigned to a different processor. Message passing was implemented using a high-level communication library. Bova et al. [72] developed a dual-level parallel analysis for the elliptic, mild-slope wave model CGWAVE [73], using both MPI and OpenMP. The model treated each incident wave component independently, resulting in a set of independent differential equations solved in parallel and distributed to multiple processors via MPI. OpenMP was utilized to parallelize the matrix-vector product and inner-product kernels, which are the most time-consuming calculations. Hence, MPI was used to simultaneously obtain solutions for multiple incident wave components, while OpenMP was used to accelerate the solution to each component. Gerostathis et al. [74] developed a parallel implementation using the message-passing paradigm [75] in a Coupled-Mode, phase-resolving model [76,77] for the transformation of a wave spectrum over steep 3D topographies. Concurrent calculations were achieved by separating problems corresponding to different frequencies and directions and decomposing resulting systems of linear algebraic equations with complex coefficients. Calculations were performed in parallel, on the basis of decomposing the computational domain into as many parts as the number of employed processes. Each process is executed on a different computer node and communicates with the others, exchanging results through messages.

In this paper, the parallelization of the proposed model is implemented using OpenMP, primarily due to its ease of implementation. OpenMP is particularly well suited for shared-memory architectures, such as those found in common personal computers, where multiple processors share access to a unified memory pool. It offers a straightforward and efficient method for parallelizing code through compiler directives, enabling the distribution of workload among available cores [78]. These directives allow for the identification of the parallelizable sections of code, such as loops or sections, and express parallelism without the need of significantly modifying the code, decomposing the domain or delving into low-level details of thread creation and synchronization. Hence, the current approach can be readily implemented in similar mild-slope models. Furthermore, the performance evaluation of the parallel implementation (facilitated with OpenMP) reveals a notable decrease in actual simulation times, rendering the model valuable for engineering applications.

The code for the model presented here has been developed in Fortran 95, with a set of compiler directives and runtime library routines to indicate parallel regions without having to deal with the low-level details of thread creation, synchronization, and communication. The most important directive is the one defining the parallel regions, which are blocks of code that are going to be executed by multiple threads running in parallel. Out of these regions, the code is executed by only one thread, i.e., the serial regions. When a thread encounters a parallel region while executing a serial region, it establishes a team of threads and assumes the role of the team’s master thread. The master thread is not only the initiator of the team but also actively participates in the computations. Within the parallel region, each thread is assigned a distinct thread number, starting from zero for the master thread and extending up to N_th − 1, where N_th represents the total number of threads in the team (Figure 2).

The implementation of OpenMP is a straightforward process that can be easily accomplished by incorporating specific clauses at the beginning of each parallel region, outlined by a loop. These clauses play a crucial role in determining the behavior of the parallel region. Notably, the employed clauses herein include PRIVATE, SHARED, and NUM_THREADS. The PRIVATE clause is used to designate variables as local to each thread. This means that these variables will possess distinct values for each thread. Examples of such variables include quantities related to the spatial derivatives present in the governing equations as well as mean values of wave characteristics, which are calculated at each time step. Conversely, the SHARED clause encompasses variables that should be accessible to all threads, either because their values are required by all threads or because all threads need to update their values. Lastly, the NUM_THREADS clause is employed to control the number of threads utilized in each parallel region. Typically, the number of the maximum available threads is utilized.

2.4. ANN Parametrization and Training for Calculating the Corresponding Eddy Viscosity Coefficient for Partial Reflection

This section presents the parametrization and training of the developed artificial neural network, which is tasked with the prediction of the proper eddy viscosity coefficients to achieve a target wave reflection value. The wave reflection coefficient depends on a multitude of parameters [40], the most important of which are the incident wave height and period, the depth at the structure’s toe, and the type of the front (partially dissipative or fully reflective). An additional parameter that influences the wave reflection is if the incident sea-state is considered as regular or irregular [51]. Considering this, a parameter signifying if the sea-state is regular or irregular was included as an input parameter at the ANN. Consequently, the eddy viscosity values to be predicted from the ANN are a function of the following input parameters (IPs, hereafter):

$$v_h=f({sea state, H}_s, T_p, h, L_e, C_R)$$

(11)

where $H_s$ is the incident significant wave height of an irregular sea-state or the wave height of a regular sea-state, $T_p$ is the peak wave period (irregular waves) or the wave period (regular waves), $h$ is the water depth at the face of the vertical wall, and $L_e$ is the width of the cells seaward of the vertical wall, where the corresponding eddy viscosity values must be applied to achieve the desirable wave reflection coefficient $C_R$.

A methodology was conceptualized in order to train and implement an artificial neural network (ANN) that will subsequently possess the capability to predict the suitable eddy viscosity coefficients to be incorporated in the model and accurately simulate partial wave reflection. The individual steps constituting the methodology are outlined as follows:

• The first step concentrates on defining a number of combinations of the IP $\left\{{ H}_s, T_p, h, L_e, v_h\right\}$ to be inserted in the numerical modelling simulations. It is noted that the eddy viscosity coefficient has to be prescribed as an IP for the numerical modelling simulations. The goal was to define many combinations of the IP in order to cover the wide array of possible values in an arbitrary dataset of the offshore sea-state wave characteristics, the water depth at the toe of the structure, and the eddy viscosity values. Hence, the lower and upper boundaries of the IP were defined, and are shown in

  Table 1. Lower and upper limits of IP considered for the numerical simulations.

  	${\mathit{H}}_{\mathit{s}}\left[\mathbf{m}\right]$	${\mathit{T}}_{\mathit{p}}\left[\mathbf{s}\right]$	$\mathit{h}\left[\mathbf{m}\right]$	${\mathit{L}}_{\mathit{e}}\left[\mathbf{m}\right]$	${\mathit{v}}_{\mathit{h}}[{\mathbf{m}}^2/\mathbf{s}]$
  minimum	0.5	3.5	2	0	0
  maximum	8.5	13.5	25	30	25

  , which cover typical values in real life field cases, and a Saltelli sampling method, employed usually in Global-Sensitivity-Analysis applications, was utilized to obtain a set of 16,000 possible combinations of the IP. To ensure realistic pairs of $H_s$ and $T_p,$ an additional restriction was imposed on the values of the peak wave period, which was to fall within the limits of the correlation between $H_s$ and $T_p$, as presented in [79]. To reduce the sheer number of combinations, the widely used Fuzzy C-Means clustering algorithm [80] was implemented to reduce the combinations of the IP to 3800. Additionally, 300 complementary scenarios were considered to cover the upper- and lower-bounding values of the IP, raising the final number of possible combinations to 4100. It is noted that to avoid creating several bathymetries corresponding to each water depth, this parameter was assigned discrete values within the range shown in Table 1, in particular $h$=$\left\{2, 5, 7.5, 10, 15, 25\right\}$.
• Six idealized 1DH wave flumes were created (Figure 3), corresponding to the six still-water depths, as defined in step 1. The flume has constant horizontal dimensions of 875 m × 25 m (350 × 10 cells). The grid cells have dimensions of dx = dy = 2.5 m. The wave maker is placed at x = 100 m. For each combination of IP, two discrete set of simulations with the HMS model were carried out: (a) no reflection, and sponge layers are positioned on either side of the flume; and (b) full or partial reflection, and a vertical wall is considered on the east side of the flume, with or without the eddy viscosity along a distance of $L_e$ seaward of the vertical wall.
• The simulation sets for the combinations of the IP are executed considering both regular and irregular unidirectional waves. The spectrum of the irregular unidirectional waves utilized a JONSWAP (with a peak factor γ = 3.3) source function and was discretized using 6 discrete wave components in the frequency domain. The time step is set equal to dt = 0.05 s, and the model runs until steady-state conditions are achieved (>4000 s).
• The reflection coefficients ($C_R$) are then calculated, in a different manner for regular and irregular waves. Specifically, for regular waves the reflection coefficient, $C_{Rreg}$, is estimated as follows:

  $$C_{Rreg}=\frac{{\zeta }_{max}-{\zeta }_{min}}{{\zeta }_{max}+{\zeta }_{min}}$$

  (12)

  where ${\zeta }_{max}$ and ${\zeta }_{min}$ are the maximum and minimum values of the envelope amplitude. Indicatively, Figure 4 illustrates the envelope plot of the surface elevation for a regular wave scenario (d/L = 0.16, H/L = 0.01) for three test cases: a) no reflection—the sponge layers are on both sides of the flume ($v_h=0$ m²/s, $C_{Rreg}=0.0$); b) full reflection—the vertical wall is on east side of the flume without the eddy viscosity ($v_h=0$ m²/s, $C_{Rreg}=1.0$); and c) partial reflection—the vertical wall is on east side of the flume with the eddy viscosity ($v_h=20$ m²/s, $C_{Rreg}=0.77$).

In the case of irregular waves, the reflection coefficient, $C_{Rirr}$, is estimated as follows:

$$C_{Rirr}=\sqrt{{\left(H_{se}/H_{si}\right)}^2-1}$$

(13)

where $H_{se}$ is the wave height by enabling the eddy viscosity coefficient and $H_{si}$ the incident wave height without enabling the eddy viscosity coefficient, calculated as average quantities in the middle of the flume.

Proceeding to the training of the ANN, two datasets are needed: the training and the validation dataset. Together, the training and validation sets form the generalization set, which includes the required data for constructing the ANN. In the present study, the generalization dataset comprises 4100 events, as described previously, divided into 3690 cases for the training dataset and 410 randomly selected cases for the validation dataset.

Subsequently, a Multilayer Feed-forward artificial neural network (ANN) is employed to predict the association between the eddy viscosity and reflection coefficient. The effectiveness of the ANN hinges on the number of units (neurons) and hidden layers within its architecture. Commonly used approaches for determining the number of hidden neurons involve trial and error as well as experimentation.

The chosen ANN architecture was selected after evaluating various configurations, considering the number of hidden layers and the corresponding number of neurons. Implementation was carried out using the Keras library [81] in the Python programming language. The architecture comprises the following components:

• An input layer transmitting input data to the network.
• A hidden layer consisting of 24 units with a rectified linear unit (relu) activation function.
• Another hidden layer comprising 24 units, also utilizing a relu activation function.
• An output layer with an identity function.

Input and output parameters are normalized utilizing a minmax normalization procedure, ensuring that they consistently fall within a designated range [0,1]. This normalization is implemented to mitigate errors associated with the specific characteristics and magnitudes of the data. Each respective value of the input parameters:$\left\{H_s, T_p,h,L_e,v_h\right\}$ is normalized according to the specified bounding limits outlined in Table 1.

The normalized values of the input parameters propagate from the input layer (IL) through the hidden layers ($\mathrm{HL}1$, $\mathrm{HL}2$) to the output layer ($\mathrm{OL}$). The output is computed as follows:

$$x_j^{\mathrm{HL}1}=f\left[\left(\sum _{i=1}^{N1}{w}_{ij}^{\mathrm{HL}1}{x}_i^{\mathrm{IL}}\right)-{\theta }_j^{\mathrm{HL}1}\right] j=1,...,N1$$

(14)

$$x_j^{\mathrm{HL}2}=f\left[\left(\sum _{i=1}^{N2}{w}_{ij}^{\mathrm{HL}2}{x}_i^{\mathrm{HL}1}\right)-{\theta }_j^{\mathrm{HL}2}\right] j=1,...,N2$$

(15)

$${\left(x_j^{\mathrm{OL}}\right)}^n=f\left[\left(\sum _{i=1}^4{w}_{ij}^{\mathrm{OL}}{x}_i^{\mathrm{HL}2}\right)-{\theta }_j^{\mathrm{OL}}\right] j=1,1$$

(16)

where $N1$ and $N2$ represent the number of units in HL1 and HL2, respectively. $F$ is the activation or transfer function,${\theta }_j$ is the bias, and $w_{ij}$ is the numerical weight between the input and hidden layers $w_{ij}^{\mathrm{HL}1}$, $w_{ij}^{\mathrm{HL}2}$.

The initial weights and biases are assigned randomly, and the ANN undergoes training using the IP to determine the optimal weights and biases. The effectiveness of the training procedure is evaluated by computing the Mean Squared Error (MSE) between the output and the target values, as shown below:

$$\mathrm{MSE}=\frac{1}{N}\sum _{i=1}^N{\left[{\left(x_j^{\mathrm{OL}}\right)}^n-{\left(x_j^{\mathrm{target}}\right)}^n\right]}^2$$

(17)

where the subscript j denotes a specific record in the dataset, and N represents the number of generalization data points, i.e., 3200 in this particular case. The MSE is subsequently employed to adjust weights and biases through the gradient descent method:

$$w_{ij}^{\prime }=w_{ij}-\kappa \frac{\partial \mathrm{MSE}}{\partial w_{ij}}$$

(18)

$${\theta }_{ij}^{\prime }={\theta }_{ij}-\kappa \frac{\partial \mathrm{MSE}}{\partial {\theta }_{ij}}$$

(19)

where $w_{ij}^{\prime }$ and ${\theta }_{ij}^{\prime }$ represent the updated weights and biases after each iteration, and $\kappa$ is the learning rate, set to 0.001. It was determined that setting the number of iterations (epochs) to 100 is sufficient to achieve minimal values of the MSE.

Ultimately, upon the conclusion of the ANN training and validation, MSE values close to zero were achieved, accompanied by a correlation factor $R$ approaching 0.996. Figure 5 depicts the ANN-predicted values plotted against the corresponding targets in the generalization dataset for both the regular and irregular waves.

From Figure 5, it can be deduced that the predicted values of the ANN for the irregular sea-state deviates in some cases from the generalization dataset targets compared to the regular waves.

3. Model Validation

In this section, the capability of the HMS model to simulate important wave transformation processes in the nearshore and inside port basins is evaluated. For this reason, a comparison of the model results with three benchmark experiments is carried out, with the first experiment concerning the Bragg scattering of regular waves over a patch of sinusoidal bars [82], the second one dealing with wave-induced oscillations due to frequency resonance in a rectangular harbor, and the third one focusing on irregular wave diffraction from a breakwater gap [83].

3.1. Bragg Resonant Reflection of Regular Waves (Davies & Heathersaw, [84])

The effect of Bragg resonant reflection of incident waves propagating over an undulating seabed with sinusoidal bars has been studied extensively by numerous researchers for several decades (e.g., [84,85,86,87]). When the wavelength of the periodic bottom undulations is approximately about half of the wavelength of the incident waves, a majority of the incoming waves are reflected by the periodic undulating seabed, which thereafter leads to a significant decrease in the transmitted waves propagating towards the shoreline, having strong implications for sediment transport and the dampening of wave energy.

To assess the capability of the HMS model to simulate the Bragg scattering of regular waves over a bed with sinusoidal undulations, the experiment of [82] was reproduced. The experiments were conducted in a wave tank with dimensions of 45.72 m × 0.91 m × 0.91 m. A patch of sinusoidal bars with the wavelength of S = 1.0 m was built into a false bottom in the tank. Two wave gauges (G1 and G2) were placed in front of the bars to calculate the reflection coefficient, according to the two-point method proposed by [88], at respective distances of about 5.125 m and 4.875 m, measured from the starting tip of the bar. Test 1 of [40]’s experiments were reproduced numerically with the HMS model. In this particular test, the number of bars is N = 4, the water depth was kept constant at 0.156 m, while the amplitude of the bars was set to 0.05 m. The experimental test contains a series of regular incident wave conditions with various wavelengths.

The numerical grid was 42 m long and was discretized with a spatial step of 0.025 m, with the total simulation time set to 100 s. Two sponge layers were placed at both sides of the numerical domain with a length of 6 m to minimize the effect of reflected waves from the tips of the numerical tank. The wavemaker was placed at a distance of 10 m from the starting point of the sinusoidal undulations. The numerical configuration of the experiment used in the HMS model is showcased in Figure 6.

Figure 7 presents the variations of both the measured and the simulated reflection coefficients with respect to the dimensionless parameter 2S/L, where L denotes the wavelength of the incident waves, and S denotes the wavelength of the undulations.

In general, a satisfactory agreement is observed between the experimental data and model results. In particular, the model slightly overpredicts the reflection coefficients for 0.9 ≤ 2S/L ≤ 1.1, indicating that the intensity of the bottom reflection is overestimated by the model near the resonance point. Overall, the model adequately reproduces the effect of the Bragg resonant reflection of regular waves over sinusoidal bars, validating its use for this type of application.

3.2. Wave Resonance in Fully Open Rectangular Harbour (Ippen & Goda [89]; Lee [90])

Wave-induced surface oscillations in harbors due to frequency-dependent resonance have been studied extensively both experimentally [89,90] and numerically (e.g., [91,92,93,94]). Such oscillations in harbors can cause significant damage to moored vessels and the corresponding structures.

To evaluate the HMS model’s ability to simulate the dynamic response of a rectangular harbor to wave-induced oscillations, the benchmark experiments of [89,90] were reproduced. The rectangular harbor basin was 0.311 m long and 0.060 m wide, with the excitation provided by regular incident waves of varying wavelengths. The numerical grid was discretized with a spatial step of 0.01 m, with the total simulation time set to 100 s. Two sponge layers were placed at the lateral boundaries of the numerical domain and behind the wavemaker to absorb the reflected wave energy.

Figure 8 showcases the variations of both the measured and the simulated amplification factor ($R$) in relation to the parameter kl, where k denotes the wavenumber and l the length of the basin. As a supplement, Figure 8 also contains the definition of the dimensions of the rectangular harbor.

Overall, a highly satisfactory agreement is evident between the experimental data and the model results. The model slightly overpredicts the amplification factors at the resonant periods. Overall, the model is well equipped to simulate wave-induced oscillations in harbors due to the effect of resonance.

3.3. Diffraction of Irregular Waves through a Breakwater Gap (Yu et al. [83])

Yu et al. [83] conducted a series of experiments focusing on the irregular wave diffraction through a breakwater gap. The experiments were carried out in the wave basin at the Dalian University of Technology, measuring 55 m in length, 34 m in width, and 1.3 m in depth. The breakwater was positioned parallel to the wavemaker, situated 7 m away from the wavemaker’s location. The still-water depth, h, inside the basin remained constant at 0.4 m. The breakwater, featuring a thickness of 0.35 m, had a gap in the center formed by two semi-circular tips. Two gap widths were explored in the experiments: 3.92 m and 7.85 m. To minimize wave reflection, absorbing layers with a width of 0.8 m were placed waveward of the breakwater. Various angles of wave attacks relative to the breakwater’s placement (i.e., 90°, 75°, 60°, and 45°) were investigated to study the combined effects of refraction and diffraction. Wave heights were measured at two cross sections parallel to the breakwater and at a normalized distance of $\mathrm{Y}/\mathrm{L}$ = 3.0 and $\mathrm{Y}/\mathrm{L}$= 6.0 shoreward of the breakwater (with $L$ representing the wavelength corresponding to the peak wave period).

The experiments concerning irregular unidirectional and multidirectional (for narrow and broad banded directional spectrum) direct wave attacks were reproduced with HMS, for the layout, with a breakwater gap of 3.92 m. The simulated cases of [41], simulated with HMS, are shown in

Table 2. Incident wave and spectrum characteristics for the three cases in Yu et al., 2000.

Case	${\mathbf{H}}_{\mathbf{s}}$ (cm)	${\mathbf{T}}_{\mathbf{p}}$(s)	h (cm)	${\mathit{\theta }}_{\mathit{o}}$ (°)	α	γ	s
U1	5.0	1.20	40	0.0	0.0081	4.0	-
N1	5.0	1.20	40	0.0	0.0081	4.0	19
B1	5.0	1.20	40	0.0	0.0081	4.0	6

, where ${\theta }_o$ is the incident wave angle, α is the spectrum Philips constant, γ is the spectrum peak-enhancement factor, and s is the directional spreading parameter.

Absorbing sponge layers, with a width of 2.3 m, were placed at the south and north boundaries of the numerical domain to minimize wave reflection. The spatial step was set equal to 0.05 m in both directions. To discretize the wave spectrum, 25 wave components were considered, and the total simulation time was set to 150 s, which was sufficient to achieve a steady-state solution. The layout, showcasing the numerical domain configuration of the experiment, is shown in Figure 9.

For each of the three simulated scenarios, wave heights were extracted along each of the two discrete sections. The obtained model results concerning diffraction coefficients ($K_d=H/H_o$) along both sections are compared with the experimental measurements and are presented in Figure 10.

It can be observed that a very good overall agreement is obtained between the HMS simulation results and the experimental measurements. For the unidirectional wave propagation case (U1), for the section with the largest proximity to the breakwater (Y/L = 3), the model predicts very well the diffraction coefficient at the center of the breakwater, while a slight overprediction of the focusing of the diffracted wave energy is observed for the section that is located further from the breakwater gap (Y/L = 6) for the positions between −1 ≤ X/L ≤ 1. For both multidirectional irregular wave propagation cases, the model predicts, in a very satisfying manner, the extent of the wave diffraction along both sections; however, once again, the peaks of the wave energy near the center of the breakwater gap are slightly overpredicted. The model results indicate that the performance for the broad-band directional spectrum case (B1) is better than the one for the narrow-band directional spectrum (N1). Ultimately, the model reproduces the experimental measurements in a very satisfying manner, validating its capability of simulating the irregular wave diffraction through a breakwater gap, a case of significant implication for port agitation purposes.

4. Harbor Agitation Application and Implementation of the ANN

In order to thoroughly evaluate the ability of the HMS model to simulate wave penetration in harbors, the physical experiment presented in [95], within the framework of the project “Benchmark tests for harbour models” [96], was numerically reproduced with HMS. The importance of this case is twofold, since it will not only serve as a further validation of the HMS model, but it will also showcase the capability of the model to simulate partial wave reflection, aided by the ANN presented in Section 2.

The experiments were conducted with the purpose of creating a schematized version of a typical port geometry. In total three different port layouts were examined, with each layout having an increased complexity. The third case was simulated numerically with the HMS wave model, with the layout consisting of a main basin, a closable side basin connected to the main basin at a 45° angle, and a composite breakwater protecting the main basin.

The wave basin, where the physical experiments were conducted, was approximately 40 m × 40 m, with a uniform water depth, h, of 0.44 m. The main basin had dimensions of 8.66 m × 14.53 m, while the side basin was 3.07 m × 10.49 m. The composite breakwater had a length of 4.6 m and a height of 0.7 m. To simulate partial reflection, a damping slope was considered both seaward and at the lee side of the breakwater. At either side of the main basin, as well as in its northern side, gravel slopes were placed to facilitate wave energy dissipation. In order to reproduce the physical experiment with HMS, considering that the ANN was trained for real-field values, a prototype was constructed with a Froude scaling of 1:20, compared to the physical experiment. Correspondingly, the rescaled incident wave characteristics for case T079, which was simulated with the HMS model, are shown in

Table 3. Incident wave characteristics for the T079 case (Van Mierlo et al., 2014) simulated with the HMS model.

Case	${\mathbf{H}}_{\mathbf{s}}$ (m)	${\mathbf{T}}_{\mathbf{p}}$(s)	h (m)	${\mathit{\theta }}_{\mathit{o}}$ (°)	Sea State
T079	2.12	6.67	0.44	0	Irregular unidirectional

.

In total, 27 gauges were used to extract the significant wave heights, and most of them were placed near the vertical walls, at positions with significant importance for loading and unloading operations conducted in ports. Based on the material used to construct the composite breakwater and the gravel slopes, the following reflection coefficients ($C_R$) were determined [95]:

• Vertical walls: $C_R$ = 1.0
• Breakwater: $C_R$ = 0.4
• Gravel slopes: $C_R$ = 0.2

The bathymetry of the domain had a constant depth of 8.8 m and was discretized using a spatial step size in both directions equal to 1.5 m. Also, a sponge layer with a width of 105 m was placed behind the wavemaker. The numerical setup of the prototype is showcased in Figure 11.

Having created the bathymetry of the numerical domain, the determination of the eddy viscosity values to achieve the desirable reflection coefficients follows. It is noted that in order to estimate the eddy viscosity values, a preliminary estimation of the wave height at the toe of the breakwater and the gravel slopes must first be undertaken. Considering that the depth of the basin is constant and thus no refraction and diffraction takes place, it is determined that the incident wave characteristics at the toe of both the breakwater and the gravel slopes are the same as the ones generated at the wavemaker. For the case of the breakwater, it was determined to apply the eddy viscosity coefficient at the area coinciding with the slopes of the structure (with a width of about 28 m). For the gravel slopes, considering that the inner gravel slope is located at about 20 m away from the measuring gauges 9, 12, and 16, a total width of 15 m was inserted as an input to the ANN to achieve the desirable eddy viscosity values. Ultimately, the input parameters inserted at the ANN are denoted in

Table 4. Input parameters provided to the ANN, depending on the geometry configuration.

Configuration	${\mathbf{H}}_{\mathbf{s}}$ (m)	${\mathbf{T}}_{\mathbf{p}}$(s)	$\mathit{h}$ (m)	${\mathit{L}}_{\mathit{e}}$	Sea State
Breakwater	2.12	6.67	8.8	28.0	Irregular
Gravel slope	2.12	6.67	8.8	15.0	Irregular

, while Figure 12 showcases the relationship between the eddy viscosity and the reflection coefficients for the breakwater structure and the gravel slopes.

From Figure 12, it can be deduced than an eddy viscosity value of 4 has to be applied at the slopes of the breakwater to achieve a desirable reflection coefficient of 0.4, whereas for 10 cells in front of the gravel slopes, an eddy viscosity value of 7 m²/s has to be implemented. The vertical walls of the main and side harbor basins are considered fully reflective and thus have an eddy viscosity value of 0. Figure 13 illustrates the significant wave height results simulated by HMS.

Observing Figure 13, reflected waves from the breakwater structure protecting the main basin can be observed, mostly near station 1 due to the oblique placement of the breakwater in relation to the incident wave direction. Diffracted waves from the breakwater tip can also be distinguished at the main basin, while significant reflection patterns are especially noticeable at the side connected basin. The simulated wave heights at each wave gauge are plotted against the experimental measurements in Figure 14.

Concerning the complexity of the physical processes present in the experimental case, the HMS wave model can reproduce the experimental measurements in a very satisfying manner across all measuring gauges. In general, the model predicts reasonably well the wave heights at the offshore wave gauges (1–4, 10, and 23–25), with the largest discrepancies observed at gauge 25, where the model underestimates the reflected waves. Of particular importance are gauges 5, 23, and 24, which are located in the vicinity of the breakwater and are influenced by reflected wave trains from the structure. A perfect agreement between measurements and the model results is obtained at gauge 23, while the model slightly underpredicts the wave heights at gauges 5 and 24. Taking this into account, it is determined that the implemented eddy viscosity coefficients provided by the ANN lead to a satisfactory model performance in simulating partial reflection from the structures. This is further supported by the comparisons at gauges 9, 12, and 16, which are located closely to the dissipating gravel slopes (see Figure 11). It should be noted that the model slightly overpredicts the impact of diffraction in this case, as can be observed mostly in gauges 26, 27, and 11. This can be attributed to the implementation of a constant eddy viscosity value along the slopes of the breakwater; hence, inserting the slopes of the breakwater at the leeside and near its tip can potentially improve the model results, albeit this requires more effort in the model pre-processing. Finally, the model predicts, in a very satisfying manner, the increase in the wave energy at the gauges situated near the vertical walls (gauges 14–20), validating its use for wave agitation purposes. Summarily, the model can simulate, in a concise manner, all the dominant processes driving wave penetration inside harbors, rendering it an accurate and efficient solution, both of which are further enhanced by the underlying support of an artificial neural network.

5. Parallel Implementation Performance

To investigate the performance of the parallel algorithm, numerical simulations of a real case are tested with different numbers of threads on a common desktop PC with a Windows 10 Pro 64-bit Operating System, equipped with a 10th Generation Intel^® Core™ i7 Processor-10700 with a base frequency at 2.90 GHz and with 8 total cores, 16 total threads, and an 8.0 GB installed RAM. A GNU Fortran Compiler version 10.2.1 was used to compile and build the algorithm.

For the purposes of the current research, the port of Rethymno, situated on the island of Crete, Greece, in the Mediterranean Sea, was selected as a test case. The port’s infrastructure and the adjacent shores were constructed in a bathymetric grid of 2660 × 1200 cells, covering an area of 6650 m × 3000 m. The grid cells have equal dimensions along the two horizontal axes, specifically, dx = dy = 2.5 m. The time step was set equal to dt = 0.05 s, and the model ran for 600 s. The incoming waves were considered irregular, adopting the JONSWAP spectrum with a peak factor of γ = 3.3. In the frequency dimension, the spectrum was discretized into five components. The offshore significant wave height was specified as H_s = 2 m, and the peak wave period was set to T_p = 7s, coming from NNW (330° N). The windward breakwater of the port was constructed with Tetrapods, and, hence, along this front, an eddy viscosity coefficient equal to 18 m²/s was applied, as provided by the ANN, for H_s = 2 m, T_p = 7s, L_e = 10 m (4 * dx), and an anticipated reflection coefficient equal to 0.4. The inner works of the port consisted of vertical non-dissipative quay walls, rendering them fully reflective. Hence, the eddy viscosity at these boundaries, as well as towards the adjacent beach, was zero.

Initially, the simulation was implemented using a serial algorithm (N_th = 1). Subsequently, the same simulation set up was executed employing multiple threads ranging from 2 to 16 (i.e., N_th = 2, 4, 6, 8, 10, 12, 14, and 16). Figure 15 shows the speedup factor of the model in relation to the number of threads. The performance exhibited excellence for a thread number up to 8, scaling almost proportionally with the increase in the number of threads. However, inefficiencies in parallelization become apparent for thread counts exceeding 8, attributable to overhead costs. Despite these inefficiencies, the overall performance remained satisfactory. The simulation duration using a single thread was 5.5 h, but with 16 threads, the time significantly decreased to 0.46 h, representing a reduction of ~92%. This enhancement is especially valuable for engineering applications that entail the execution of numerous simulations across a wide array of scenarios, which can represent varying incident wave characteristics, different sea water level shifts, and alternative port layouts. It should be noted that the results were consistently identical both for the serial and the parallel algorithms. Figure 16 illustrates the nearshore wave field results, revealing distinct reflection and diffraction patterns in the vicinity of the port infrastructure. Additionally, noticeable shoaling, refraction, and breaking phenomena are observed along the adjacent shore.

6. Conclusions

In this study, we presented a novel hyperbolic mild-slope wave model (HMS) designed to address partial wave reflection, wave penetration, disturbance, and resonance within port basins. The key features of the HMS model include its ability to solve the hyperbolic form of the mild-slope wave equation, allowing for the simulation of both regular and irregular unidirectional and multidirectional waves. The parallelization of the algorithm using OpenMP substantially reduced simulation times, making it applicable to real port areas spanning several kilometers. The model, supported by an artificial neural network (ANN), demonstrated significant advancements, particularly in estimating eddy viscosity values, efficiently and accurately simulating partial wave reflection. The training of the ANN involved a comprehensive dataset covering a wide range of sea-state characteristics, water depths, and eddy viscosity values. The results indicated high accuracy, with Mean Squared Error values close to zero and high correlation factors.

The model’s validation against benchmark experiments demonstrates its capability to accurately simulate important wave transformation processes in nearshore areas and inside port basins. Regarding the Bragg scattering of regular waves propagating over a sinusoidal bar, the model results demonstrate satisfactory agreement with experimental measurements, offering an accurate estimation of wave energy reflection near the resonance point. This case study indicates that apart from port engineering applications, the model is also well equipped to handle wave propagation in a nearshore with complex bathymetries. In relation to wave-induced oscillations in a rectangular harbor due to frequency resonance, the numerical model very satisfactorily reproduced the amplification of wave energy at the backwall of the harbor basin, observed at the wave resonant periods. Hence, the model is deemed capable of simulating the long-wave resonant oscillations generated by fully reflecting waves. Furthermore, regarding wave penetration through a breakwater gap, diffraction coefficients obtained by the numerical model’s simulations for the cases of unidirectional and multidirectional wave fields with narrow and broad directional spreading, respectively, were, in general, in a very good agreement with the experimental measurements. This case has strong implications for the capability of the model to accurately simulate the effect of diffraction through port entrances.

Furthermore, wave penetration in a schematized harbor basin, protected by a sloping breakwater, was reproduced numerically with the HMS model. This case study was notably significant, as it provided an additional validation of the ANN’s capability to accurately predict the level of partial reflection from the sloping breakwater. A comparison of the model results with experimental data from 27 wave gauges was deemed very satisfactory, confirming the model’s capability to simulate all dominant wave-transformation processes occurring in wave-penetration applications inside ports and harbors. Future research efforts will delve more into partial reflection under oblique wave attacks. Lastly, the application of the HMS model to the Port of Rethymno, Crete, showcased its ability to simulate wave disturbance in real-case applications, while the parallel implementation of the model exhibited excellent performance, resulting in a significant reduction of ~92% in simulation times. Future research endeavors could explore the implementation of MPI or a hybrid approach to further reduce simulation times.

In summary, the presented HMS model, supported by an ANN, offers a robust and efficient solution for simulating wave disturbance and resonance in port basins, making it a valuable tool for engineers and researchers in the field of port and coastal engineering.