A Flexible Fully Nonlinear Potential Flow Model for Wave Propagation over the Complex Topography of the Norwegian Coast

Coastal wave propagation and transformation are complicated due to the signiﬁcant variations of water depth and irregular coastlines, which are typically present at the Norwegian fjords. A potential ﬂow model provides phase-resolved solutions with low demands on computational resources. Many potential ﬂow models are developed for oﬀshore waves and lack of numerical treatments of coastal conditions. In the presented work, several modiﬁcations are introduced to a fully nonlinear potential ﬂow model with a σ -grid for the purpose of coastal wave modelling: Shallow water breaking criteria are included in addition to deepwater breaking algorithms to approximate breaking waves over a complete range of water depth. A new coastline algorithm is introduced to detect complex coastlines and ensure robust simulations near the coast. The algorithm is compatible with structured grid arrangement in the horizontal plane and allows for high-order discretisation schemes for the free surface boundary conditions for an accurate representation of complex free surfaces. A parallelised solver for the Laplace equation is utilised to ensure fast simulations for large domains with multi-core infrastructures. The proposed model is validated against theories and experiments for various two- and three-dimensional nonlinear wave propagation and transformation cases that represent typical coastal conditions. The simulations show a good representation of nonlinear waves, and the results compare well with experiments. Furthermore, two large-scale engineering scenarios are simulated, where the applicability of the coastline algorithm and the parallel commutation capability of the model are demonstrated.

where η is the free surface elevation, φ = φ(x, η, t) is the velocity potential at the free sur- where h = h(x) is the water depth measured from the still water level to the seabed.

212
The Laplace equation, together with the boundary conditions are solved on a σ-coordinate 213 system. The σ-coordinate system follows the water depth changes and offers flexibility for 214 irregular boundaries. The transformation from a Cartesian grid to a σ-coordinate is expressed 215 as follows: .
The Laplace equation is discretized using second-order central differences and solved using 222 a parallelized geometric multigrid preconditioned conjugate gradient solver provided by the 223 hypre library (van der Vorst, 1992).

225
The gradient terms of the free-surface boundary conditions are discretized with the 5th-226 order Hamilton-Jacobi version of the weighted essentially non-oscillatory (WENO) scheme 227 (Jiang and Shu, 1996). The implementation of the WENO scheme in the presented model is 228 described in the Appendix. 229 For time treatment, a 3rd-order accurate total variation diminishing (TVD) Runge-Kutta 230 scheme (Shu and Osher, 1988) is used. Adaptive time-stepping is used by controlling a con-231 stant time factor as an equivalence to the Courant-Friedrichs-Lewy (CFL) condition: where u max , v max are the maximum particle velocities in x and y directions at the free 234 surface, h max is the maximum water depth, g = 9.81 m/s 2 is the gravitational acceleration.

247
In the model, the vertical coordinates follow a stretching function so that the grid becomes 248 denser close to the free surface: where α is the stretching factor and i and N z stand for the index of the grid point and 250 the total number of cells in the vertical direction.

251
As an example, a general description of a progressive Airy wave can be expressed as: and function A(z) follows: which is governed only by the wave number k, which can be defined by the linear dispersion 253 relationship to the wave angular frequency: where g is the gravity acceleration.

255
A correct representation of the phase velocity depends on the correct representation of 256 the wave number. The new method is based on the assumption that a constant absolute 257 truncation error at every vertical location can preserve the correct shape of the function f (z) 258 and yield the correct wave number. Function f (z) is a Taylor expansion of free surface over 259 the depth: If the absolute error is set to a constant E for every vertical location and the function 261 f (z) and its derivatives are known, one can find a maximum cell size ∆z(η) = z − η at every 262 location (Pakozdi et al., 2021): where x is scaled to the length of the relaxation zone. The free surface velocity potential φ and the surface elevation η are increased to the analytical values in the wave generation zone: Following the same methodology, the free surface velocities potential φ and the surface 265 elevation η are reduced to zero or initial still water values in the wave energy dissipation zone 266 or numerical beach to eliminate wave reflection of the outlet boundaries.

267
Waves can also be generated at the inlet using a Neumann boundary condition where the 268 spatial derivatives of the velocity potential are defined. In this way, the velocity potential at 269 the boundary is calculated using the desired analytical horizontal velocity: where u(x, z, t) is the analytical horizontal velocity.
In addition, wavemaker motion input can also be used to generate waves in the numerical 273 wave tank (NWT Deepwater steepness-induced breaking is activated with a steepness criterion: where β = 0.3 is recommended by Smit et al. (2013) from the comparison to physical 289 tests.

290
After a wave breaking is detected, two methods are available to represent the energy dissi-291 pation during the wave breaking process. The first method is a geometric filtering algorithm 292 that smoothens the free surface for energy dissipation (Jensen et al., 1999). Here, an explicit 293 scheme is used and therefore there is no CF L constraint. Another method is to introduce 294 a viscous damping term in the free surface boundary conditions locally around the breaking 295 region (Baquet et al., 2017). When wave breaking is detected, the free surface boundary 296 conditions Eqn.
(3) then become: Handling the complex coastline has been a challenge when applying a potential flow model in 313 the coastal area. The first difficulty is efficient grid generation around the complex boundaries. curves. The accuracy of the coastline is also sensitive to the grid resolution at the coastline.

318
The second difficulty is possible numerical instability during the wave run-up process in the 319 swash zone. The vertical velocity in the free surface boundary condition Eqn. (7) cannot be 320 given directly but calculated from derivatives of velocity potential over water depth. In some 321 scenarios, there is a thin layer of water in the swash zone where the water depth can be con-322 sidered as infinitesimal. In such cases, the vertical derivative of the velocity potential at the 323 free surface tends to be ill-defined. This tends to cause unreasonably high particle velocities η is the surface elevation, h is the still water level measured from the bottom. The 332 relationship among h , h and η is illustrated in Fig. 1. If the local water depth h is smaller than a threshold h, then the local cell is identified 334 as a dry cell. When a cell is identified as a dry cell, the velocities in the cell is set to be zero: The default threshold is set to be 0.00005 m according to the practice of Zijlema et al.
where φ(x, y), η(x, y) and ν b0 are the velocity potential, the free surface elevation and  is the boundary of the coast-following relaxation zone to reduce numerical instability and customise reflection properties ofr the coastline.
The proposed coastline algorithm (see Fig. 2) is intended to have the following beneficial waves. There is no need for coastline-following (i.e. body-fitted) grid generation. Therefore, 375 the algorithm is less case-dependent.

376
2) The level-set method enables accurate capture of the shoreline positions in an implicit 377 manner. The method uses a smooth signed distance function for the coastline geometry 378 rather than representing the coastline geometry with a body-fitted grid. This is in contrast to  water is performed to prove that the model is able to represent wave transformation over 404 significant wave depth variation within a short horizontal distance. And finally, the breaking 405 algorithm is proven to be effective for a wave breaking over a mild slope.    wave with a wave steepness corresponding to 90% deepwater breaking limit.     With the chosen grid resolution, the grid at the steep ramp is shown in Fig. 13. In spite 541 of the significant bathymetry change, the σ-grid follows the topography very well.   Figure 15: Frequency spectra at wave gauge G2 in the simulation of bi-chromatic wave propagation over a steep ramp, (a) frequency spectrum near ω 3 , (b) frequency spectrum near ω 1 and ω 2 , (c) frequency spectrum nearω 4 , (d) frequency spectrum near ω 5 , (e) frequency spectrum near ω 6 .
It is seen that all theoretical frequency components are represented in the frequency spec-  One of the wave breaking events is shown in Fig. 21. Here, the wave crest increase to its maximum and the wave front becomes vertical. Since the free surface is single-valued, a 616 visualisation of overturning breaker is not within the scope of the model. However, it is seen 617 from Fig. 20 and Fig. 21 that the breaking event is correctly detected and the wave energy is 618 correctly dissipated in the numerical simulation.        The visual observation in Fig. 31 is confirmed in Fig. 33. It is seen that the significant 752 wave heights at all wave gauges inside the harbour are significantly reduced after waves pass 753 the two breakwaters BW1 and BW2 in Fig. 26a. Breakwater 2 has a much stronger effect 754 than breakwater 1 on reducing wave height at nearly all wave gauges. However, the combined 755 usage of the two breakwaters further reduces the wave height at wave gauges 6,7 and 9. In 756 fact, most small leisure boats and fishing boats are docked in the marina near wave gauge 9.

757
The combined usage of the breakwaters proves to be necessary in order to protect the small 758 boats from large motions.

759
The study of Mehamn harbour shows the model's capability to represent complicated    The spatial distribution of the H s are compared at the fish farm site as seen in Fig. 40.    The relaxation zones along the coastlines make the coastal reflection property customisable.

891
In addition, the possible numerical instability in the free surface boundary condition at the 892 very shallow water region is avoided. As a result, the model is able to include complicated to- In conclusion, the proposed fully nonlinear potential flow model and the coastline algo-913 rithm provide a working framework for coastal wave modelling in complex coastal environ-914 ments, considering accuracy, flexibility and efficiency.  vertical plane problems, International Journal for Numerical Methods in Fluids 28 (1998)