An improved depth‐averaged nonhydrostatic shallow water model with quadratic pressure approximation

Phase‐resolved information is necessary for many coastal wave problems, for example, for the wave conditions in the vicinity of harbor structures. Two‐dimensional (2D) depth‐averaging shallow water models are commonly used to obtain a phase‐resolved solution near the coast. These models are in general more computationally effective compared with computational fluid dynamics software and will be even more capable if equipped with a parallelized code. In the current article, a 2D wave model solving the depth‐averaged continuity equation and the Euler equations is implemented in the open‐source hydrodynamic code REEF3D. The model is based on a nonhydrostatic extension and a quadratic vertical pressure profile assumption, which provides a better approximation of the frequency dispersion. It is the first model of its kind to employ high‐order discretization schemes and to be fully parallelized following the domain decomposition strategy. Wave generation and absorption are achieved with a relaxation method. The simulations of nonlinear long wave propagations and transformations over nonconstant bathymetries are presented. The results are compared with benchmark wave propagation cases. A large‐scale wave propagation simulation over realistic irregular topography is shown to demonstrate the model's capability of solving operational large‐scale problems.

forces on various configurations of multiple vertical circular cylinders. Furthermore, simulations of marine fluid-structure interaction were performed for semisubmerged horizontal circular cylinders in tandem, 10 and nonlinear marine hydrodynamics were investigated in detail Reference 11. Broader applications of the model are also seen on the sediment transport analysis 12 and the coastal infrastructure design. 13 Typically, these simulations require relatively fine three-dimensional (3D) grids and are, therefore, more computationally demanding.
Phase-resolved modeling of the far-field wave field is important for delivering a realistic wave generation boundary condition for higher resolution near-field wave modeling. However, the far-field wave propagation toward the coast is a large-scale phenomenon, which puts a limitation on the application of the Navier-Stokes approach in spite of the increasing computational capacities. Less computationally demanding models are required to model the far-field large-scale phase-resolved wave propagation efficiently. As most coastal areas share relatively shallower water conditions, depth-averaged shallow water models have been favored for the coastal wave modeling. These models are essentially two-dimensional (2D) and, thus, require less cells. The advances of such models have been focused on developing numerical methods to accurately capture the frequency dispersion relation and the nonlinearity when the water depth increases or a rapidly varying bathymetry is involved. A common representative of shallow water models is the Boussinesq-type wave model. 14,15 Herein, the lack of vertical flow information is compensated through the Boussinesq terms, which help to calculate the correct frequency dispersion of the waves. This approach is valid from shallow to deep water, depending on the order of the Boussinesq terms. 16 However, higher order mixed time-space derivatives in the Boussinesq equations tend to cause numerical instabilities. More recently, the possibility of using nonhydrostatic shallow equations with a single layer or multiple layers in the vertical direction has been explored by  With an increasing number of vertical layers, the flow information in the vertical direction is better resolved. However, it has been shown previously that the increase of vertical layers leads to a significant increase in computational costs. For example, Monteban 21 observed that the simulation time using two layers is nearly 10 times compared with that using a single layer. Cui et al 22 improved the two-layer approach such that it has similar computational efficiency as a one-layer counterpart and, yet, maintaining a high linear dispersion accuracy. Although the commonly used vertical pressure profile is linear, a quadratic pressure approach has been presented by Jeschke et al. 23 It is stated that, with an approximation of a proposed quadratic vertical pressure profile, the model can achieve at least a good equivalence to existing fully nonlinear weakly dispersive Boussinesq models. 23 This method presents itself as an attractive alternative for modeling shallow water waves, while potentially avoiding the numerical instabilities due to higher order terms in a Boussinesq-type model and the increased computational costs from a larger number of vertical layers in a multilayer nonhydrostatic model. However, only simple scenarios such as one-dimensional (1D) standing waves and progressive solitary waves over a flat bottom have been investigated previously. 23 Herein, several terms of the derived equations are neglected, which leaves the final question of reliability open. It is reported by Jeschke 24 that it is challenging to incorporate the vital term involving the varying bathymetry into her numerical model. As a result, the model's accuracy is seen to be less ideal than the theoretical expectations when changing bottom is present. Therefore, this article includes a numerical procedure to discretize this term appropriately. This enables the authors to emphasize the accuracy gain from the quadratic pressure approximation for nonconstant bathymetries.
The accuracy of shallow water models has been improved over the last years. High-order numerical schemes are employed in the development of Boussinesq-types models. Wei and Kirby 25 applied a fourth-order accurate Adams-Bashforth-Moulton scheme for the time discretization and a mixed fourth-order and second-order scheme for the spatial discretization. Shi et al 26 employed a mixed finite volume and finite difference method (FDM) using a fourth-order accurate monotone upstream-centered schemes for conservation laws reconstruction technique for the advection term and a third-order Runge-Kutta scheme for temporal discretization. However, few high-order implementations are presented for nonhydrostatic models. Zijlema et al 20 present their model using a second-order discretization scheme in space and a second-order leapfrog algorithm in time. Jeschke et al 23 implement the quadratic pressure model with the second-order P 1 NC − P 1 finite element method 27,28 for the advection terms and a Leapfrog method for the time stepping. In a recent development, Jeschke 24 also implemented a second-order discontinuous Galerkin scheme in the model. Thus, high-order numerical implementations are left to be fulfilled in order to advance the development of nonhydrostatic models.
In addition, parallel computations are incorporated in many shallow water models in case of computationally demanding simulations. Shi et al 26 presents a parallelized Boussinesq model following the domain decomposition strategy with a message passing interface (MPI). Good scaling characteristic is observed up to 48 cores. Zijlema et al 20 also uses the same parallelization technique and achieve linear scalability up to eight cores. However, the newly proposed quadratic pressure approximation 23 has not been incorporated into any parallel code. A good scalability up to hundreds of processors is also not presented in the literature regarding shallow water models in general. For large-scale operational engineering applications, such scalability is in great demand.
Ensuring high-quality input waves is another important aspect in the development of a shallow water model. The typical practice is to impose the surface elevation and the depth-averaged velocities to the boundary. 14,15,20,26,29,30 Periodic boundary conditions are also widely used, for example, a spatial periodic boundary condition is applied by Madsen et al, 31 and a double periodic boundary condition is implemented in Reference 23. Another popular wave generation method is the relaxation method, 1,32 which has high flexibility and tends to result in less reflected waves. 33 This method has been widely implemented in Navier-Stokes solvers, 34 but remains absent in the development of shallow water models. The feasibility of using a relaxation method for the wave generation and absorption in a nonhydrostatic shallow water model remains to be explored.
In the presented article, REEF3D::SFLOW is introduced as a novel nonhydrostatic shallow water model following the quadratic pressure approximation. 23 Developed as a part of the REEF3D framework, the proposed model has direct access to all the existing numerical schemes and parallelization algorithms in REEF3D. Thus, the model presents itself as the first nonhydrostatic shallow water model with high-order discretization schemes, for example, a fifth-order weighted essentially nonoscillatory (WENO) scheme in spatial discretization and a third-to fourth-order Runge-Kutta scheme for the temporal discretization. The model also innovatively employs the relaxation method 1 for the wave generation and absorption. With a model equipped with high-order numerical methods, this article presents for the first time the simulations of nonlinear long wave propagations over varying bathymetries using the quadratic pressure approximation. In these simulations, the equations with the depth-related terms are solved and the overall performance gain from the quadratic pressure approximation is investigated comprehensively. Computational scalability up to multihundred cores is demonstrated with the proposed model. An expanded validation process is then presented, including several well-known benchmark cases incorporating wave propagation over changing topographies and wave-structure interactions. In addition, a large-scale coastal wave propagation over a natural topography is presented to demonstrate the model's capability for engineering applications.

NUMERICAL THEORY
The mass and momentum conservation for an incompressible inviscid flow leads to the continuity and Euler equations in three dimensions: where U, V, and W are velocities in x, y, and z directions, is the constant density, P T represents the total pressure, and g is the gravitational acceleration. Additional source terms such as bottom friction and turbulent stresses are omitted here, but are straightforward to include if needed. The water depth h = d + consists of two parts: the still water depth d and the free-surface elevation , as displayed in Figure 1. Defining the horizontal velocity vector as U = (U, V), the kinematic boundary conditions at the free-surface and the bottom are: F I G U R E 1 Basic definitions in the shallow water model: the water depth h, the still water depth d, the free-surface elevation , the coordinates system, and the schematics of the assumed linear pressure profile and quadratic pressure approximation The shallow water assumption, that is, the horizontal acceleration is much greater than the vertical acceleration, implies a hydrostatic pressure. In order to get a hydrodynamic pressure correction, the total pressure P T is assumed to consist of a hydrostatic part P and a hydrodynamic part Q. The pressure and its boundary condition at the free-surface is given by: The velocities and the dynamic pressure are depth-averaged by integrating over the water depth: By contrast to previous models, 20 where the pressure is solved at the bottom, the proposed model consists of only depth-averaged quantities. A relation between the depth-averaged pressure q and the pressure at the bottom Q| −d needs to be defined in order to close the system. If the linear pressure profile 17,20 is assumed, the pressure at the bottom is simply twice the depth-averaged pressure, or: Consequently, the governing equations with only depth-averaged variables are: Jeschke et al 23 replaces the linear assumption with a quadratic vertical pressure profile as shown in Equation (15).
Following the quadratic assumption, the governing equations with depth-averaged variables become: The governing equations with the boundary conditions are solved on a structured staggered grid using a FDM. Chorin's projection method 35 is applied for the solution of the velocities. The fifth-order conservative finite difference WENO scheme proposed by Jiang and Shu 36 is used for the discretization of convective terms for the velocities u, v, and w. The total variation diminishing third-order Runge-Kutta explicit time scheme developed by Shu and Osher 37 is employed for time discretization. It involves the calculation of the spatial derivatives and the dynamics pressure three times per time step. The information containing pressure is solved using the Poisson equation: Herein, the parameter h p denotes the water level in the center of the cell. In a staggered grid arrangement, this is where the dynamic pressure q, the vertical velocities w, and the free-surface location are solved. The horizontal velocities are solved at the faces of the cells. The high-performance solver library HYPRE 38 is employed to solve the Poisson pressure equation using the PFMG-preconditioned BiCGStab algorithm. 39 The dynamic pressure q is then used to correct the velocities in a correction step. Hence, the corrections of the velocities with the quadratic pressure approximation are where u * , v * , w * are intermediate-step velocities with only hydrostatic pressure. The term Φ on the right-hand side of Equations (18) to (20) is treated with a procedure following the principles of the fractional step method of Le and Moin. 40 Assuming the dynamic pressure does not change significantly within one Runge-Kutta substep, the intermediate velocities u * , v * , w * are corrected with the dynamic pressure gradients of the previous substep: where q n,rk is the dynamic pressure from the previous Runge-Kutta substep. The spatial derivatives of Φ are updated with the corrected velocities u * * , v * * and w * * in Equation (16), which is then inserted into Equations (22) to (24) to obtain the velocities at the new time step. The time derivative term inside Φ is then calculated with simple finite differences: where is the increment factor of the corresponding Runge-Kutta substep and u n,rk , v n,rk , w n,rk are the velocities from the previous Runge-Kutta substep. Parallel computation is enabled by decomposing the simulation domain into smaller subdomains. The communication between these domains is achieved through a ghost cell approach. The MPI is then used for the communication at the subdomain boundaries.
The location of the free-surface is determined based on the divergence of the depth-integrated horizontal velocities as given in Equation (17). The free-surface is reconstructed using the fifth-order WENO scheme. 36 The solutions of the stencils are weighted, that is, a coefficient or weight is assigned to the solution of each stencil. The scheme assigns the largest weight to the smoothest solution and can therefore handle large-gradient free-surface changes caused by the varying bathymetry. As an example, the discretized form of Equation (17) whereĥ i+1∕2 is the water level at the cell face i + 1∕2.ĥ i+1∕2 is reconstructed with the WENO procedure: The ± sign indicates the upwind direction. The nonlinear weights ± n are calculated for each ENO stencil based on the smoothness indicators. 36 For the upwind direction in the positive i-direction, the three possible ENO stencilsĥ 1 ,ĥ 2 , andĥ 3 are:ĥ Wetting and drying are handled by setting the velocities in cells below a certain user-defined threshold of the water level to zero: The default threshold is set to be 0.00005 m, which is used throughout the presented work. The approach tracks the variation of the shoreline accurately and avoids numerical instabilities by ensuring nonnegative water depth. 19,41 Wave generation and absorption are carried out with the relaxation method as described in Bihs et al. 4 The relaxation function formulated by Jacobsen 1 is used in the proposed model: wherex is scaled to the length of the relaxation zone. The velocities u, v, the surface elevation , and the pressure p are increased to the analytical values in the wave generation zone and reduced to zero or initial still wave values in the wave energy dissipation zone: All types of wave theories, type of wavemakers and wave signal input available in the existing code are applicable to the proposed shallow water model as well.
A breaking wave criterion is introduced Reference 42 to represent the wave breaking process. The wave breaking is initialized when the vertical velocity of the free-surface exceeds a fraction of the shallow water celerity: At the same time, the dynamic pressure is neglected and remains so at the front of the breaker. For the persistence of the wave breaking, the coefficient (0 < < ) is introduced in Equation (42) instead of to stop the wave breaking process. The computations become nonhydrostatic again when the vertical velocity of the free-surface falls out of the range of the criterium. = 0.6 and = 0.3 are recommended as they work well with most of the waves. 42 By introducing the wave-breaking criterion and removing the dynamic pressure during breaking, the momentum is well conserved, the energy dissipation is well represented, and the asymmetry and skewness of nonlinearity are respected. 42

VERIFICATION
The proposed numerical model REEF3D::SFLOW is first verified for the wave propagation in a 28 m long 1D flume as shown in Figure 2. The wave generation zone of one wavelength is at the inlet of the flume, and a wave energy dissipation zone of two wavelengths is located at the outlet. Four different wave cases are simulated with the proposed model.  Table 1. A monotonic reduction of the error can be observed with the refinement of the grids.  Table 2. CFL number of 0.2 matches both the trough and crest well and errors approach to the ones with CFL number 0.1. As a result, CFL = 0.2 will be used in all the following simulations of this article. Figure 4A shows that the linear progressive wave is well represented by the solver at an intermediate water depth. Both, the wave height and phase are matching satisfactorily. It is also noticeable that the relaxation method dissipates the wave energy well at the wave energy dissipation zone where the surface elevation remains constant at the still water level and no artificial reflection is observed.

Linear progressive wave propagation over constant bathymetry
The advantage of the quadratic pressure approximation is demonstrated by comparing the surface elevation with quadratic pressure approximation with the linear pressure profile in References 17,20 (see Figure 4B). It is observed that, with a linear pressure assumption, the wave phase starts to shift shortly after the waves propagate outside the generation zone. By contrast, the quadratic pressure approximation improves the phase accuracy significantly and approximates the theoretical value more precisely due to a better representation of dispersion.

Second-order Stokes wave propagation over constant bathymetry
Next, a second-order Stokes wave 43 of H = 0.1 m and L = 4 m is simulated in the same 1D numerical flume. The grid convergences study is shown in Figure 5A. Similar to the previous study, the cell size dx = 0.02 m is found to be suitable for   Table 3. With the quadratic pressure approximation, the asymmetry due to the high-order approximation is well presented, and both, the wave height and phase match well with the theory. It shows that the model provides a good representation of the nonlinearity of progressive waves. In comparison, the simulation with linear pressure profile shows an increasing difference in phase overtime compared with the theoretical result.

Cnoidal wave propagation over constant bathymetry
A fifth-order cnoidal wave 43,44 of H = 0.21 m and L = 4 m is investigated in the 1D numerical flume to test steep periodic wave propagation in shallow water. The steepness of the wave is H∕L = 0.0525, the wavelength to depth ratio is H∕d = 0.42, which is about 65% of the breaking limit suggested by Laitone. 45 As shown in Figure 6A, dx = 0.02 m is still a suitable cell size to capture the wave surface elevation accurately despite the increased wave steepness. Following the same methodology as in Section 3.1, the relative error and L2 norms are computed and shown in Table 4. The wave profiles obtained with the quadratic pressure approximation and the linear pressure assumption are also compared in Figure 6B. The wave troughs start to show slight deformation, while the crests are still well preserved with the wave height to depth ratio closer to the breaking limit. The geometry of the steep cnoidal wave is kept constant during the propagation. It is also observed that the phase misalignment from the linear pressure assumption amplifies with the increase of wave steepness because the linear pressure profile assumption deviates further from the physical pressure distribution.  Figure 7A. The relative errors and L2 norms are also computed and shown in Table 5. Furthermore, simulations with the quadratic pressure approximation and the linear pressure assumption are simulated with dx = 0.02 m. The numerical computations are compared with the analytical values at propagation time 10, 20, 30, and 40 seconds, shown in Figure 7B. It is seen that the numerical results with the quadratic pressure remain in good agreement during the entire wave propagation process. Small amplitude waves propagate in opposite direction and trailing waves start to form during the simulation with the linear pressure. Simultaneously, the wave height increases during the process due to weaker dispersion from the linear assumption. These findings are in agreement with the investigations of Jeschke et al. 23 The It is empirically assumed that the scaling is linear within 16 processors, that is, one physical node on the cluster. Therefore, the computation time with one processor is linearly extrapolated from the 16-processor simulation. The computational speed of the one-processor simulation is considered as the base reference. The simulation time on one processor divided by the simulation time on multiple processors is defined as a speed-up factor. The relation between the speed-up factor and the number of processors as well as the number of cells per processor are plotted in Figure 8. It shows that the performance increases almost linearly with the number of processors within the chosen range.

VALIDATIONS AND APPLICATIONS
The evolution of waves over a nonconstant bathymetry is complicated, and the performance gain from the quadratic pressure approximation in a general setting was recommended as future work by Jeschke et al. 23 To fill the research gap, wave propagations over nonconstant bathymetries of various configurations are simulated and validated with the available experimental data. A wave-structure interaction study is also validated against the benchmark. Jeschke et al 23 suggest the quadratic pressure approximation has the best performance when the water depth to wavelength ratio is below 0.25. The selected benchmark cases all share the water depth condition within the suggested range. In addition, a large-scale wave propagation over a natural topography is presented based on an engineering scenario.

Wave propagation over a submerged bar
First, the well-known benchmark case of wave propagation over a submerged bar 47 is tested. The configuration of the numerical setup based on the experiment is shown in Figure 9. A 2D wave tank of 38 m is equipped with a wave generation zone of 5 m to the left end and a wave energy dissipation zone of 9.5 m to the right end. The beginning of the submerged bar is located 6 m downstream from the wave generation zone. Eight wave gauges are located above the submerged bar with the x-coordinates being 11,16,17,18,19,20,21, and 22 m, as shown in Figure 9. The incident wave height is H = 0.021 m, and the wave period is T = 2.525 seconds. A grid convergence study is performed at gauge 2 and 6, before and after the crest of the submerged bar, as shown in Figure 10I,J. A cell size of dx = 0.02 m is found to sufficiently represent the phenomena and shows good agreement with the experimental data. A simulation time of 60 seconds is used. The numerically predicted time series of the surface elevations at gauge 1 to 8 are compared with the experimental data in Figure 10. The results match well with the experimental measurements before the waves reach the submerged bar and during the shoaling process, for example, at gauges 1 and 2. It demonstrates that the model can represent the dispersion relations well with changing bathymetry. At the crest of the bar, no wave breaking happens, but the wave decomposition takes place and results in higher harmonic wave components. The wave decomposition phenomenon is observed at wave gauges 3 to 5, where the numerical results show accurate agreement with the experimental measurements as well. On top of the relatively steep downslope, the waves undergo a deshoaling process as the water depth increases. During this process, it is observed that the numerical results start to show differences in phase from the experimental data. The discrepancies accumulate from wave gauge 6 to 7. When the waves reach wave gauge 8, a significant difference is observed. This shows a less discussed limitation of existing shallow water approximations for deshoaling processes. Furthermore, the results are also compared between the quadratic and the linear pressure profile assumptions. As an example, the comparisons of the surface elevations at gauge 3 and 5 are shown in Figure 11. At both gauges, the quadratic assumption shows good alignment in phase with the experiment, while the linear assumption tends to predict a faster moving wave front. The observation is consistent with the investigation in Section 3.

Solitary wave interaction with a rectangular abutment
In this benchmark study, the solitary wave interaction with a surface-piercing rectangular abutment is investigated. Based on the experiments, 48,49 the numerical wave tank is defined as shown in Figure Figure 13. The first peak in the distributions is the result of the incoming solitary wave impact on the abutment. After the incident solitary wave passes the abutment, it is reflected from the wall at the end of the tank and interact with the abutment again, resulting in the second peak. The grid convergence study shown in Figure 13J is performed at gauge 7, which is located at the downstream side of the abutment. At this location, both, the interaction between the structure and the incoming waves and the properties of the reflected waves can be well observed. It indicates that the cell size dx = 0.05 m sufficiently captures the details of the wave pattern and gives good results compared with the experiments. At gauge 1 and 2, the first peaks show the solitary wave propagates without much interruption and, therefore, preserves its wave height. A second minor peak is noticed right after the peak, which is due to the partially reflected waves from the abutment. Gauge 3 shows an increase of the wave height due to the narrowing of the channel, while gauge 4 presents a further increase of the peak because of the interaction with the abutment. The peaks increase to about 0.11 and 0.13 m at gauge 3 and 4, respectively. Since gauge 5 is located in the constricted part of the channel, the flow velocity increases, and the pressure decreases. As a consequence, the wave surface drops. At gauge 6, the first peak occurs right after the wave crest passes the abutment, while the depth-averaged solution tends to smooth out the results in the sheltered region behind the abutment. At gauge 8 and 9, two peaks of equal heights are observed, indicating that the reflected wave shares the same wave height as the incoming wave. This shows that there is no damping of the soliton and the model provides an accurate representation of the solitary wave propagation. Similarly, the two peaks also share similar height at gauge 7, where no wave transformations occur before and after the wave reflects from the vertical wall. When the reflected wave reaches the abutment, a second peak occurs at gauge 6. After the reflected wave passes the abutment, gauge 4 also witnesses the second peak. In general, the wave patterns from gauge 6 to 4 mirror each other.
Finally, the second peak at wave gauge 5 and the first peak at wave gauge 7 are compared with the quadratic and the linear pressure approximation in Figure 14. Similar to the previous observations, the linear approximation predicts an increased phase velocity, while the quadratic approximation matches the experiment well in phase.
The details of the free-surface during this process is also visualized in Figure 15. Figure 15A shows the free-surface at simulation time t = 7 seconds, right before the solitary wave reaches the abutment. The solitary wave preserves its waveform. After the wave passes the abutment, a vortex is observed at the downstream behind the abutment, as can be seen in Figure 15B. When the reflected wave reaches back toward the abutment from the right-hand side, the wave crest meets the vortex from the last interaction before a second interaction, as seen in Figure 15C. After the reflected wave passes the abutment, two vortices are observed on both sides of the abutment. Figure 13 reveals that the resolution of the vortex is smoothed out at gauge 4 and 6, while the other wave gauges are well represented.
It might be interesting to notice that the 2D shallow water model is as accurate as the computational fluid dynamics study in Reference 4 except for the vortices representation in the wakes of the abutment. Herein, the results of simulations based on the 3D Navier-Stokes equations show a slightly better match with the experiments. The cost of

F I G U R E 16
The numerical wave tank setup of the wave breaking over a sloping bed, view from the side. The water depth is constant at 0.5 m, the gray-shaded object is the sloping bed with a slope of 1:35. Four wave gauges are arranged near the breaking point the computational resource, however, is significantly lower using the proposed shallow water model. This benchmark case is simulated with 16 processors on the Vilje supercomputer about 56 times faster than the 3D simulation with the same configuration.

Plunging breaking waves over a sloping bed
In Section 4.1, nonbreaking waves over a submerged bar are modeled. In a more extreme situation, where the shoaling is so strong that the wave steepness increases over a certain threshold, the wavefront becomes unstable and breaking takes place. The numerical wave tank is initialized based on the experiments in References 51,52 to model a breaking wave scenario. The wave tank has a total length of 40 m and a height of 1 m. A wave generation zone of 9.8 m is located at the inlet of the tank; a wave energy dissipation zone of the same length is arranged at the outlet. An inclined bed with a slope of 1:35 is located 4 m away from the wave generation zone. The obstacle increases to 0.748 m at the right end of the tank. The water depth is constant at 0.4 m. Wave gauges 1 to 4 are located on the slope, 10, 11, 12, and 12.3 m away from the wave generation zone, respectively. A fifth-order cnoidal wave with wave height H = 0.128 m and wave period T = 5 seconds is propagated in this simulation, which is supposed to result in a plunging breaker on the slope according to the experiment. A simulation time of 40 seconds is used Figure 16.
The sensitivity to the grid resolution is investigated with different cell sizes of dx = 0.0025, 0.005, 0.01, 0.02, and 0.05 m. The wave surface elevation at wave gauge 4 is chosen for comparing the results from different cell sizes. As can be seen in Figure 17E, the simulations capture very steep wavefronts as well as instabilities at the wave crest with all cell sizes. It is not possible to observe the overturning process because the shallow water model represents the free-surface with The view on the wave crest is shown in more detail in Figure 17F, where it is visible that dx = 0.005 m captures the peak values most accurately. The simulated wave elevations at different wave gauges with dx = 0.00 5m are compared with the experimental data in Figure 17 in order to assess the model's capacity to resolve the surf-zone wave transformations. The wave crests increase significantly when the waves propagate from gauge 1 to 2, showing an increasing shoaling process. As the waves evolve on the slope, an unstable wave crest is seen at gauge 3 and the wave height decreases slightly compared with that at gauge 2. The instability at the crest remains as the waves approach gauge 4 and a further decrease of the wave crest is noticed. These time series suggest that the breaking happens between gauge 2 and 3. To identify the breaking point, the wave elevation profile at different time are compared in the same plot ( Figure 18). It is seen that at x = 21.580 m, the wave crest is the highest while the wavefront becomes vertical for the first time indicating the location of the breaking point. Correspondingly, a breaking height of h b = 0.208 m is measured at x = 21.580 m. In the experiment, the breaking point is detected at x = 21.595 m and a breaking height of h b = 0.196 m is measured. Both, the predicted breaking point and are very close to that in the experiment. The wave surface elevation profile is shown in Figure 19. As can be seen in Figure 19A, the wave height increases significantly, the wave shape becomes narrower, the crest becomes unstable, and the wavefront becomes vertical, indicating a breaking process. At a later time, the wave energy dissipates and the wave height decreases dramatically. An attempt to simulate the breaking wave using the linear pressure approximation leads to a numerical failure. It indicates that the quadratic pressure approximation is superior for the simulation of breaking waves.  Figure 21B. Strongly reflected waves can be seen at the tips of the peninsulas that reach out northward into the ocean. Stripes of submerged reefs in the north-south directions create strong shoaling, as higher waves are shown to be following the same pattern of the submerged reefs. When the waves propagate southward, refraction occurs and bend the wave rays toward the shore. When the waves start to reach the harbor, the narrowing entry causes diffraction. A fraction of the diffracted waves manages to bypass the curve-shaped peninsulas and enter the inner harbor. The complicated wave transformations and their interactions are well demonstrated in the simulation results. Finally, the model's computational performance including a complicated bathymetry with wetting and drying and the breaking algorithm is determined in a similar manner as described in Section 2. The simulations are conducted for 500 iterations with the number of processors fixed to 16,32,64,128,256, and 512 on the supercomputer Vilje. The computational time with one processor is linearly extrapolated from the 16-processor simulation and is used as a base reference for the speed-up factor. The relation between the speed-up factor and the number of processors as well as the number of cells per processor are then plotted in Figure 22. It shows that with the presence of a complex topography and the wetting-drying scheme, the model is as computationally efficient as with a constant bottom within 200 processors, while it slows down compared with the ideal scaling characteristics afterward.

CONCLUSION
The shallow water model REEF3D::SFLOW has been presented in this article. The model solves the depth-averaged shallow water equations with nonhydrostatic extensions and a quadratic vertical pressure profile approximation. 23 In comparison to well-known Boussinesq-type models, the proposed model treats the pressure terms differently. A typical Boussinesq model adds higher order terms to express the hydrodynamic pressure. The proposed model adds nonhydrostatic extensions to the shallow water equations and solves for the hydrodynamic pressure explicitly from a Poisson equation. This equation is solved iteratively using an implicit scheme. Thus, the proposed model offers simpler numerics and indicates higher numerical stability by avoiding the high-order pressure terms of a Boussinesq model. The current model assumes a quadratic pressure approximation for a better representation of dispersion and always solves the depth-averaged pressure. This is in contrast to the multilayer approach that uses vertical layers to represent dispersion and solves the pressure at the lower layer interface. Thereby, the presented approach saves the additional computational costs from the increasing number of layers. High-order numerical methods are incorporated into the new model. Consequently, it is the first model with the quadratic pressure approximation that combines high-order schemes and fully parallelized computation. The wave generation and absorption are achieved using a relaxation method, which is absent in the current literature. The approach proves to generate various wave types with correct amplitude and dispersion, and no artificial reflections are observed in the numerical wave tank. The accuracy of the high-order scheme is confirmed for 1D and 2D wave propagation cases with a constant bathymetry. The 2D large-scale simulation of a wave propagation over constant bathymetry presents a near-linear scaling of the computational speed with an increasing number of processors up to 512. Furthermore, the model shows an almost linear scaling up to 128 processors if a natural topography is included in the numerical wave tank. The speed-up is reduced with a further increase of computational units due to the complex boundary treatment from the topography.
Overall, the study confirms the advantage of the quadratic pressure approximation over the linear pressure assumption for multiple validation cases. The linear pressure assumption leads to an overshooting phase velocity for all the regular wave tests in the articles. It also causes a secondary wave during the solitary wave propagation. The quadratic pressure approximation improves the phase information for progressive waves significantly and removes the unrealistic free-surface disturbances.
A key advancement presented in the current work is the inclusion of the varying bathymetry and structures in a nonhydrostatic shallow water model with the quadratic pressure approximation. A fractional step method is applied in the proposed numerical model in order to meet the challenge of incorporating the term Φ that appears in the bottom pressure calculation. Thus, the simulations of the nonlinear long wave propagation over varying topographies using a nonhydrostatic model with the quadratic pressure assumption are possible for the first time. The wave transformations over varying topography are well represented and in good agreement with the experimental data. The model can represent the complex free-surface during wave-structure interactions and predicts the breaking wave height and locations accurately. The quadratic pressure approximation again provides a better representation of the free-surface than the linear pressure assumption for the wave propagation over varying bathymetries. The challenges of representing the deshoaling process using a nonhydrostatic shallow water model is also discussed, and the study confirms the findings from previous research. 53 It can be concluded that, within the applicable range of the quadratic assumption, 23 the quadratic pressure approximation presents better results both with a constant and a varying bathymetry. The large-scale engineering application shows a good computational scaling character with the wetting and drying of complex topography included. In general, the model presents itself as a good alternative to shallow water modeling with robust and efficient numerical methods. The model also serves as an additional option within the hydrodynamics code REEF3D. As a consequence, an integrated wave-modeling cascade is more easily adaptable because different submodels are developed on a single platform and the information exchange can be made more convenient.