The Role of Non-Hydrostatic E ﬀ ects in Nonlinear Dispersive Wave Modeling

: Surface water waves is an important research topic in coastal and ocean engineering due to its inﬂuences on various human activities. In this study, our purpose is to gain a deeper insight on the e ﬀ ects of non-hydrostatic (NHS) pressure on surface wave motions and its role in numerical modeling, based upon the high-order NHS model and optional vertical accelerations. The relative contribution of non-hydrostatic e ﬀ ects ( P nhs / P hs ) and its sensitivity on phase celerity and amplitude of dispersive waves are quantiﬁed. The vertical structure of P nhs / P hs clearly indicates stronger NHS e ﬀ ects in deeper waters and its signiﬁcance near the surface. The NHS e ﬀ ects mainly slow down wave celerity and maintain incident amplitude for linear dispersive waves. The NHS e ﬀ ects are also responsible for increased amplitude and phase speed under strong non-linearity. The inter-relation between (un)realistic physical responses and model errors is discussed. Further, four experimental conditions for waves with complicated interactions are examined. Overall, the NHS e ﬀ ects play a critical role in side-band generation of bi-chromatic waves, and increased celerity and amplitude during nonlinear shoaling, as well as velocity moderation under co-existence of depth-varying currents. Possibly owing to weaker wave–wave interactions, however, wave directionality does not strongly interfere with FNHS / QNHS (Fully / Quasi Non-HydroStatic) e ﬀ ects on a fast-modulated nonlinear evolution of spatial focusing or di ﬀ raction waves.


Introduction
Surface water waves with significant influences upon human activities is an important research topic in the field of coastal and ocean engineering [1][2][3]. In the past several decades, analytical and experimental efforts [4][5][6] have been made to reveal basic characteristics of surface wave motions. Driven by the horizontal gradient of non-hydrostatic pressure (associated with vertical acceleration of fluid), dispersion is the most essential feature in wave propagation (i.e., longer waves traveling faster than shorter waves in a wave group) [4]. Besides, nonlinear transformation can arise owing to (mutual or self) interactions among larger wave components, such as the formation of side-band or higher harmonics in deep or shallow waters [1,4,[6][7][8][9][10]. A Doppler shift in dispersion can also be found if there are currents co-existing with waves [11,12]. In the meantime, along with theoretical works, a great amount of effort has also been devoted to surface wave modeling [13][14][15].
Accordingly, complicated wave motions in more realistic conditions (e.g., combination of shoaling, refraction, and diffraction over irregular bottom) can be resolved. However, proper representations of dispersion and nonlinearity required for accurate and efficient wave modeling are challenging [14]. With a solid background from the literature, this study aims at gaining a further and deeper insight into the role of non-hydrostatic effects in complicated wave transformation. In particular, a better understanding of the detailed contribution of non-hydrostatic pressure effects would assist continuous development and applications of nonlinear dispersive wave modeling.
To date, according to the vertical dimension treatment in governing equations, wave modeling frameworks can be classified into two major groups: (i) depth-integrated models, including Boussineq equation (e.g., [16][17][18][19][20][21]), mild-slop equation (e.g., [22][23][24]), and their variants, as well as (ii) the multi-layer models based on Laplace's equation derived from the potential flow theory with irrotational assumption (e.g., [25][26][27][28][29][30]), and the complete Navier-Stokes equations (e.g., [31][32][33][34][35][36]). In the former group, subject to limited computer capacity in the 1970s, the classical Boussinesq-type equation [16] and mild-slope equation [22] utilize a polynomial and depth integration factor to describe the velocity profiles for weakly nonlinear dispersive waves and linear waves over an entire range of depth, respectively. The constrains have been relaxed by several approaches, such as the optimal expansion levels, higher-degree polynomials, higher-order derivatives, and multiple layers [17][18][19][20][21]23,24]. In contrast to depth integration, with the vertical structure of pressure and velocity (potential) fields, accurate description of nonlinear dispersive waves can be achieved naturally, since the non-hydrostatic pressure has been directly included. Without any approximations or assumptions, non-hydrostatic modeling is the most general approach, especially for flows involving shear, vortex, and turbulence [37][38][39]. A theoretical framework for the interpretation of wave dispersion/nonlinearity in non-hydrostatic models can refer to Bai and Cheung [40]. However, a demanding computation was typically required in non-hydrostatic modeling.
In earlier works, the non-hydrostatic effects on surface wave motions were demonstrated simply by using a qualitative comparison (i.e., a noticeable discrepancy) between the prediction and analytical solutions (e.g., [34,35]). Chen [54] further quantified hydrostatic and non-hydrostatic components for linear seiche oscillations. The scale analysis indicated that the non-hydrostatic pressure gradient was not necessarily small, although the non-hydrostatic pressure was about 2-3 orders of magnitude smaller than the hydrostatic component. For wave transformation with various dispersion/nonlinearity and complicated interactions (among waves, bottoms, and currents, as well as directionality); however, the role of non-hydrostatic effects has not yet been characterized, as far as the authors are aware.
The objective of this paper is to clearly identify the effects of non-hydrostatic pressure on surface wave motions and its important role in numerical modeling. Based upon the high-order NHS model [49,51], the non-hydrostatic contributions are quantified by optional vertical accelerations, that is, hydrostatic (HS, no acceleration), quasi non-hydrostatic (quasi non-hydrostatic defined in this paper is different from "quasi hydrostatic", which represents the surface calculation step of the pressure-split algorithm in a fully non-hydrostatic model [55]), (QNHS, local acceleration only), and fully non-hydrostatic (FNHS, local and convective acceleration). Specifically, the normalized vertical profile of non-hydrostatic pressure (P nhs /P hs ) and its effects on wave celerity and amplitude are examined against the linear and non-linear analytical solutions. Further, the non-hydrostatic effects on complicated wave problems are elucidated using the high-order NHS model and available experiments, that is, non-linear bi-chromatic waves in deep water, non-linear wave shoaling, non-linear waves interacting with a depth-varying current, and 3D focusing/diffracting waves.

Model Description
The non-hydrostatic model [49,51] solves the unsteady, incompressible, Navier-Stokes equations to simulate free surface wave motions based on the staggered finite difference Crank-Nicholson scheme in the transformed σ domain (see Figure 1). The implicit numerical algorithm relaxes the stability restriction of computational time-steps in comparison to the explicit ones. Besides, the algorithm treats a 3D problem as a series of 2D vertical planes to yield block hepta-diagonal system matrixes that can be solved directly without iterations. The boundary-fitted co-ordinate σ = (z * − η)/(h + η) ensures an accurate representation of the free surface elevation (i.e., z * = η or σ = 0) and the irregular bottom topography (i.e., z * = −h or σ = −1). Most importantly, the higher-order top-layer pressure treatment based upon a cubic polynomial interpolation can accurately resolve phase velocity. The developed NHS model has been carefully validated against either analytical solutions or experimental data for various wave problems from shallow to deep water depth [49,51], clearly demonstrating the model's efficiency (requirement of only 2-5 vertical layers) and accuracy (dispersive degree up to Kh = 3-15). In the following, the model description including the governing equations and numerical algorithm is briefly presented here. Details can refer to Young et al. [49] or Young and Wu [51]. The transformed 3D governing equations can be written as ∂v ∂t with where the superscript asterisk indicates the physical domain; u, v, and w are the velocity components; P is the normalized pressure divided by a constant reference density; g is the gravitational acceleration; ν is the kinematic viscosity; and δ l and δ c are the optional flags for the vertical accelerations to determine the HS, QNHS, and FNHS experimental conditions. Several boundary conditions are required to solve the Equations (1)- (4). At the bottom boundary, a kinematic condition for a fixed bed is applied. For the free surface, a conservative form of the free-surface equation that results from the integration of the continuity equation over the water depth with associated kinematic conditions can be expressed as where H = h + η is the total water depth. In addition, the continuity of normal stress is enforced, that is, P| σ = 0 = P a , where P a is the atmospheric pressure and taken as zero here. Notice that the staggered pressure is arranged at the center of the top layer rather than the free surface. Using the integral method [15,[43][44][45]49,51], the top-layer pressure boundary condition can be derived as consisting of the hydrostatic component caused by the free-surface elevation (i.e., the first term on the right-hand side) and the non-hydrostatic one attributed to the vertical acceleration (i.e., the rest terms). At the inflow boundary, the horizontal velocity component with a ramp function is specified to generate waves without the undesired noises due to an impulse motion [41,45]. For the outflow boundary, a combination of the sponge layer technique and the Sommerfeld radiation boundary condition is used to minimize wave reflection [44]. The model domain can be discretized by a set of computational cells with non-uniform (or uniform) grid-spacing along each direction ( Figure 1). The velocity and pressure variables are defined under a staggered grid system. The overall numerical algorithm consists of three major steps. First, the discretization step applies the Crank-Nicholson scheme to the continuity equation, momentum equations, and free-surface equation, expressing the vertical velocity, and free-surface elevation as functions of horizontal velocities and providing the relation of the lower-layer pressure to its upper one and adjacent velocities. Second, pressure treatment step uses Equation (6) to specify the top-layer pressure. Particularly, a cubic function with fourth-order truncation error is employed to address the vertical profile of dispersive waves for better accuracy. With the top-layer pressure and the relation between each layer, the pressure over the entire depth can be represented as a function of horizontal velocities. Third, the solving step solves the horizontal velocities in each vertical plane from the block hepta-diagonal system matrix, which is obtained by eliminating the vertical velocity and pressure in the discretized horizontal momentum equations and is then solved using an efficient direct matrix solver. The rest variables can be simultaneously obtained by simple back-substitution.

Non-Hydrostatic Effects on Linear/Nonlinear Waves
In the following, the relative contribution (P nhs /P hs ) and sensitivity of non-hydrostatic effects on the phase celerity and amplitude in linear/nonlinear dispersive waves are examined.

Linear Wave Dispersion
Linear (frequency) dispersion means that water waves with different wavelengths λ (or wave periods T) propagate at different phase speeds c in a given depth h [1]. Incident wave conditions with infinitesimal amplitude (a/h = 0.001) and four dimensionless water depth Kh (ranging from shallow to extremely deep-water depth, i.e., Kh = π/25, π/5, π, and 5 π) are considered here. We conduct numerical experiments in a vertical 2D tank ( Figure 2) with a length L = 10λ and an undisturbed water depth h = 1 m. A sponge layer with a length L SP = 2λ is employed to avoid reflection. The simulation is carried out with 40 grids per wave length along the x direction and "stretched" five vertical layers [31]. The time-step ∆t is determined by setting the Courant number Cr = c∆t/∆x = 0.5. The viscosity is set to zero.

Contribution and Sensitivity Analysis
The vertical profile and relative contribution of non-hydrostatic pressure are analyzed here, instead of a qualitative comparison of pressure contours [34,35,44] or pressure time-series [56] in most previous works. As shown in Figure 3a, the HS condition reveals no non-hydrostatic pressure and the FNHS/QNHS distribution under the crest perfectly agree with the analytical solutions for all four different conditions of Kh. In comparison to the horizontal velocity that has a positive profile with a maximum at z * = η (not shown here), the magnitude of the negative non-hydrostatic pressure remains zero at the free surface and then increases along the depth. Sharper variations appear in the deeper water conditions (e.g., Kh = 5π). Figure 3b compares the non-hydrostatic pressure (P nhs ) with the hydrostatic component (P hs ). The orders of a relative ratio P nhs /P hs are O(10 −5 -10 −2 ) from Kh = π/25 to Kh = 5π, showing that the non-hydrostatic pressure is influential, especially in the dispersive deep-water waves. In addition, the location of the maximum contribution of P nhs occurs near the surface. Subsequently, the role of non-hydrostatic pressure in wave propagation speed is revealed in Figure 4. As the dispersive degree increases, the stronger non-hydrostatic effects lead to a decreased phase speed (see black lines). In terms of sensitivity, we consider the variation of phase speed in response to the change of non-hydrostatic pressure. Apart from the hydrostatic condition, the phase speed of intermediate-water waves (π/10 < Kh < π) reduces from (gh) 1/2 to 0.56 (gh) 1/2 , while the relative contribution P nhs /P hs increases from O(10 −5 ) to O(0.003). In the cases of deep-water waves (up to Kh = 20), the wave celerity keeps decreasing down to 0.22 (gh) 1/2 as the P nhs /P hs ratio grows to O(0.012). In addition to the physical response, model errors (calculated from the simulated free surface elevation) may interfere in the aforementioned relationship. Two relative errors for the phase speed and top-layer non-hydrostatic pressure are also shown in Figure 4, that is, ε c = (c model − c exact )/c exact and ε P = (P nhs ) model TOP − (P nhs ) exact TOP /ρgη (see blue lines for the identical FNHS/QNHS results). It can be found that the five-layer NHS model shows great accuracy for a wide range of dispersive degree (up to Kh = 5π). Compared to earlier works (e.g., [31,45]), the high-order top-layer pressure treatment is proved to be critical for efficient and accurate modeling, although the ratio (P nhs /P hs ) is small (e.g., O(10 −5 -10 −2 )). As the dimensionless water depth increases up to 20 (or the top-layer P nhs /P hs > 0.01), the 3% error in the non-hydrostatic pressure leads to a 6% under-estimation of phase speed. Such an accuracy issue in the extremely deep-water conditions can be tackled by using more vertical layers or higher-order numerical treatment in NHS or Boussinesq-type models (e.g., see [21]).

Free-Surface Spatial Profile
For a better understanding of the non-hydrostatic effects on wave celerity and amplitude, the predicted free-surface spatial profiles are demonstrated in Figure 5. For different water depths (i.e., Kh = π/25 − 5π), both FNHS and QNHS simulations are in excellent agreement with the analytical solutions. Almost identical results indicate that the local vertical acceleration dominates the behavior of linear waves. The role of non-hydrostatic pressure is further revealed through the comparison of the HS condition. While appropriate for shallow water (Kh = π/25), neglecting the NHS effects leads to unrealistic long-wave patterns once the dimensionless water depth increases (Kh > π/5), that is, in regard to the over-predicted wavelength and underestimated amplitude. Such discrepancies can be physically explained using the concepts of phase velocity and energy flux [4]. For a deep-water condition (Kh = π), the incident wave with a theoretical phase speed 0.56 (gh) 1/2 is forced to propagate at the speed of (gh) 1/2 toward the downstream of the HS numerical tank. As a consequence, the faster celerity yields the longer wavelength (i.e., λ HS = 1.77 λ Inc ). Further, the energy flux should be conserved between the incident boundary and any given location in the tank. Substituting the theoretical and HS group velocity C g into (ρgH 2 C g ) Inc = (ρgH 2 C g ) HS , one can easily obtain a smaller wave height H HS = 0.54 H Inc (or a HS = 0.27 H Inc ). Overall, the results indicate that the NHS effects attributed to the vertical acceleration would slow down celerity and remain incident amplitude for dispersive waves.

Contribution and Sensitivity Analysis
Similarly, the vertical profile and relative ratio of non-hydrostatic pressure are compared in Figure 6. In general, the overall features are quite close to those in the linear cases, except for the larger magnitude O(10 2 ) and slightly different contribution pattern near the surface. While the local vertical acceleration is dominant, the convective one plays a non-negligible role in the non-hydrostatic pressure for the cases of strong non-linearity (see the shifted profile and under-estimated top-layer significance as aK increases to 0.20 or 0.30). Subsequently, Figure 7 compares the wave celerity and top-layer non-hydrostatic pressure versus wave steepness (rather than dimensionless water depth in those linear cases). As the degree of nonlinearity extends, the stronger non-hydrostatic effects yield an increased phase speed. For the highly non-linear case (aK = 0.30), the wave celerity and the contribution ratio P nhs /P hs increase to about 1.05 c linear and 0.30, respectively. The FNHS pressure representation fully considers the vertical acceleration, and thus yields precise phase speed for Stokes waves (i.e., ε c < 0.5%). Dropping the physics of convective effects, the QNHS mode creates 2% error in the non-hydrostatic pressure, leading to 3% under-estimation in non-linear wave celerity (computed from the simulated free surface elevation).  [45]. The full influence of non-hydrostatic pressure can be observed from the difference between FNHS and QNHS results. The QNHS mode gives smaller non-hydrostatic effects (see Figure 6), further affecting wave amplitude and phase speed, especially under stronger non-linearity. Without convective vertical acceleration, the free-surface elevation presents decreased wave steepness in the simulation results, that is, ε = aK ≈ 0.20 (or η crest /H Inc ≈ 0.55). The resulting phase speed in the QNHS prediction can be explained using the Stokes theory [5]. According to Equation (8), the nonlinear incident waves with smaller ak would propagate at only 97% of phase speed c (see 3% error in Figure 7 or the lagged time series in Figure 8c) when traveling to the downstream of the QNHS numerical tank. In other words, the fully NHS effects considering both local and convective vertical acceleration would increase the wave steepness (or amplitude), as well as propagation speed for nonlinear waves.

Non-Hydrostatic Effects for Waves with Complicated Interactions
In this section, the non-hydrostatic effects on complicated wave problems are further elucidated using the high-order NHS model and available experiments, that is, non-linear bi-chromatic waves in deep water, non-linear wave shoaling, non-linear wave interacting with a depth-varying current, and 3D focusing/diffracting waves.

Bi-Chromatic Waves
We choose the laboratory experiment of bi-chromatic waves [9] to demonstrate the FNHS/QNHS pressure effects on the evolution of a nonlinear deep-water wave group with strong wave-wave interactions. The non-breaking case B11 [57] is considered, while its wave amplitude and phase speed were not easily resolved using the multi-layer Boussinesq model [21] and fourth-order non-linear Schrödinger equation [58].
According to the experimental setting of Hwung and Chiang [9], the numerical wave tank is 300 m long with an undisturbed water depth of 3.5 m. A sponge layer is used to replace the wave-absorbing beach. The incident waves consist of the low-and high-frequency components with periods of T l = 1.7 s and T h = 1.5 s (K l h = 4.87 and K h h = 6.26). The imposed wave heights are the same, that is, H = 0.1 m (aK l = 0.070 and aK h = 0.089). The computation domain is discretized by 1500 uniform horizontal cells (∆x = 0.2 m) and five stretched vertical layers. The time step is determined by setting Cr = c∆t/∆x = 0.1, where c is defined as (gh) 1/2 . The total simulation time is 400 s.

Linear Bi-Chromatic Waves
We firstly carry out a linear benchmark case with exactly the same periods but smaller amplitude (a l K l = a h K h = 0.01). Figure 9 compares the simulated free surface time series with the second-order analytical solution [59] at 11 gauge stations from the dimensionless distance K c x = 23 to K c x = 359, where K c is the wave number obtained using the mean wave period T m = 1.6 s. As it can be seen, the two components of the bi-chromatic waves travel to the downstream in symmetric groups with a period of 25.5 s. Figure 10 shows the predicted amplitude spectrum with signals only at the base frequencies.
Overall, both FNHS and QNHS simulations (with similar dispersion capability) excellently describe the linear wave group since almost no wave-wave interactions occur and dispersive effects are dominant in such a small amplitude case.

Non-Linear Bi-Chromatic Waves
For the non-linear experimental condition, Figure 11 compares the predicted wave profiles with the observation data. The second-order solution is not adopted due to its difficulty in representation of energy redistribution [60]. The FNHS results reasonably capture the evolution of non-linear bi-chromatic waves. It can be clearly seen that the initial wave trains resemble a linear superposition of two imposed wave components (K c x = 23). As the wave trains travel to the downstream, a strong asymmetric envelope with steep fronts and gently sloping rears (K c x = 58) forms, consistent with several other experimental results [58,61,62]. Afterward, the wave energy is focused into the center of the wave groups (K c x = 114), followed by the recurrence of the initial stage (K c x = 161). The processes of modulation and de-modulation cyclically repeat until the end of the tank. Figure 12 further shows the predicted spectrum with two side-band components (2ω 1 − ω 2 , 2ω 2 − ω 1 ) evolved from the non-linear interactions between the base frequencies. The QNHS mode without the convective vertical acceleration would under-estimate both wave amplitude and phase speed of the non-linear incident wave group, as discussed in Section 3.2. Moreover, the side-band frequencies are hardly captured (see Figure 12b). Therefore, the evolution of non-linear bi-chromatic waves in the QNHS numerical tank is unrealistic (see noticeable phase lag and 20% smaller amplitude at the center of wave groups at the distance of K c x = 114 in Figure 11b). Based upon a multi-layer Boussinesq model and fourth-order non-linear Schrödinger equation, such a discrepancy due to under-estimated wave amplitude and/or phase speed has also been reported [21]. Overall, the role of FNHS pressure in wave-wave interactions is apparently revealed.

Wave Shoaling
Subsequently, we consider the experiments conducted by Nwogu [18] to observe FNHS/QNHS effects on wave transformation over a sloping bottom. The numerical tank is 20 m long with a 1:25 beach placed at x = 4.6 m, and the still water depth is 0.56 m (see Figure 13d). A sponge layer (at the last 2 m of the tank) is used to represent the experimental wave absorber. The incident conditions are T = 0.85 s (Kh = 3.13) and H 0 = 0.04 m (aK = 0.11). The computational domain is discretized by a set of horizontally uniform cells (i.e., ∆x = 0.05 m) and five vertically stretched layers. The time-step is determined by setting Cr = 0.5.

Linear Wave Shoaling
Similarly, we start from the linear cases with consideration of another two wave periods (a = 0.001 h and T = 5.0 s, 0.85 s, 0.5 s) to demonstrate different shoaling behaviors [17]. Based upon energy flux conservation [4], the theoretical wave envelope can be expressed as H 1 (C g1 ) 1/2 = H 2 (C g2 ) 1/2 , where the subscripts indicate different locations; H is the wave height; and C g is the group velocity [15]. Figure 13 shows the predicted free-surface spatial profiles and the corresponding dimensionless water depth, Kh. H/H 0 = 0.98 can be found. Overall, the FNHS or QNHS pressure representation is able to describe the dispersion-dominated shoaling process.

Non-Linear Wave Shoaling
In terms of non-linear wave shoaling, Figure 14a shows a fully developed surface profile from the FNHS results. Figure 14b,c compares the FNHS/QNHS predicted and measured free-surface time series. A larger incident wave interacting with the sloping bottom results in stronger nonlinearity at the shallower depth, showing the asymmetric profile with an enlarged height [20,63]. At the depth of 0.28 m (or Kh = 1.67), the local wave height is H/H 0 = 1.01 (or aK = 0.12) and the crest-and-trough asymmetric factor is µ = η crest /(η crest −η trough ) = 0.57. For an even shallower depth h = 0.07 m (or Kh = 0.67), the wave form can reach H/H 0 = 1.12 (or aK = 0.21) and µ = 0.70. Without the full consideration of non-hydrostatic effects, the QNHS model or the Boussinesq equation [18] give 6% or 10% underestimated wave celerity and amplitude, respectively. By contrast, excellent agreement between the FNHS results and experimental data clearly indicates the importance of fully non-hydrostatic pressure in non-linear wave shoaling over a sloping bottom.

Non-Linear Wave Interacting with Depth-Varying Current
We further show the FNHS/QNHS pressure effects on nonlinear wave interacting with a depth-varying current [11]. Based upon the experimental configuration, the numerical tank has a length of 20 m and an undisturbed water depth of 0.  Figure 15 presents the measured current profile U(z) and compares the predictions of wave-induced velocity u(z) under wave crest. In addition to the experimental data, the fifth-order wave-only and Doppler-shifted solutions (respectively for U(z) = 0 and U(z) = U s ) are also included [5]. Compared with the wave-only solution, the Doppler-shifted results have a larger surface velocity with a shaper distribution along the depth owing to the enlarged amplitude and shortened wavelength/increased dimensionless water depth Kh [64]. For the sheared current, the presence of vorticity at the upper depth (−0.25 m < z < 0 m) counteracts the Doppler shift by the surface current U s [65] and then yields the smaller velocity as observed in the experiments. The FNHS results and observation data are in excellent agreement. The QNHS results appear to capture the structure of horizontal velocity as the convective effects are diminished, possibly by the opposing current.

Interaction with a Positively Sheared Following Current
Similar comparisons in Figure 16 show the contrary results in case 5. The Doppler shift by a positively uniform current presents a smaller surface velocity with a milder vertical distribution because of the decreased amplitude and lengthened wavelength (lessened dimensionless depth Kh). The vorticity of the sheared current moderates the shift, giving the larger velocity for the upper depth (i.e., −0.25 m < z < 0 m). The FNHS simulation exactly resolves the horizontal velocity profile, similar to the five-layered model results in [11]. Close to the Doppler-shifted solution, the QNHS simulation without convective vertical acceleration fails to represent the interaction of non-linear waves with the strongly sheared following current. Apparently, the wave-induced velocity with a noticeable (i.e., 20%) under-estimation near the surface can be found. Overall, the inter-comparison clearly reveals the fully non-hydrostatic effects on complicated wave-current interactions.

Spatial Focusing and Diffracting Waves
Lastly, we explore the FNHS/QNHS pressure effects upon the propagation of directional waves, that is, the spatial focusing and diffracting waves in the deep water [66]. In these two cases, a three-dimensional version of the non-hydrostatic model [49] is employed. The numerical wave flume (see Figure 17) has a length of 11 m, a width of 4 m, and an undisturbed mean water depth of h = 0.6 m, according to the experimental configuration. A deep-water spectral wave group with 32 components ranging from f = 0.6855 to 1.4745 Hz (or Kh = 1.311 to 5.246) are considered. The incident waves are generated based upon the frequency focusing technique with an addition of directionality. For the spatial focusing waves, the motion of each wave paddle is prescribed by the horizontal velocity, that is, where the amplitude a n is chosen to produce a constant steepness spectrum; ω n and K n are the angular frequency and wave number of the nth component; the incident angle of each paddle is θ p = tan −1 (x f /|y p |); the theoretical focusing location x f =3.3 m, and time t f =12 s are set. For the spatial diffracting waves, a lateral variation from |y p | = 0.45 m to |y p | = 2 m is introduced into the horizontal velocity using a cosine window that tapers the amplitude down to 10% at the lateral boundaries, that is, 1 + 0.9 cos 2 y p π 4 ] (10) Figure 17. Sketch of the experimental setup by Wu and Nepf [66].
Here, we choose the non-breaking incipient conditions (i.e., a n k n = 0.0055). At the outgoing boundary, a 3 m long sponge layer in conjunction with the Sommerfield radiation boundary condition is applied to minimize wave reflection. The computational domain is horizontally discretized by a set of uniform cells with the spatial increment ∆x = ∆y = 0.05 m. In the vertical direction, five layers are used to cover the full spectral range evolution. The model is run for the total simulation time 40 s, with a small time-step ∆t = 0.001 s.

Spatial Focusing Waves
The spatial focusing waves have been regarded as an important formation mechanism of the so-called oceanic freak waves that can mysteriously occur and cause notorious hazards. Figure 18a shows the perspective view of the FNHS modeled wave field in the vicinity of the focusing event. As clearly found in the progress of wave evolution, a "hole in the sea" is followed by the "wall of waters" during t * = −1 s to −0.5 s, presenting an exceptionally large, steep, and asymmetric wave. At the consecutive stage (i.e., at t * = 0 s), an extreme freak wave is rapidly developed with a maximum amplitude twice over the average wave height of the wave packet, which is the main characteristic of freak waves [67]. Afterward, the freak wave disappears in a second (t * = 1 s), consistent with the short-live feature reported by field observations [68]. Note that owing to nonlinearity, the actual focusing event does not occur at the theoretical location prescribed by the linear theory [66]. Figure 19a compares the FNHS/QNHS model results of free-surface elevation time-series with the available recorded data at x = x f = 3.3 m and different lateral locations, i.e., y = 0, 0.3, 0.6, 0.9, and 1.5 m. The FNHS model well predicts the distinctive crescent shape of the wave crest. The FNHS (or QNHS) predicted focusing event yields a maximum wave crest of 8.04 cm (or 7.94 cm). The actual focusing locations from both computations give a small deviation of 0.05 m (a grid spacing size), that is, a slightly shorter focusing distance in QNHS. Interestingly, the FNHS/QNHS results reveal insignificant differences for the fast-modulated focusing waves, in contrast to the slowly modulated, nonlinear, bi-chromatic wave train. Reasonable explanation might be the weaker wave-wave interactions between each smaller steepness wave component (i.e., a n k n = 0.0055) in this nonlinear transient phenomenon.

Spatial Diffracting Waves
The spatial diffracting waves can be used to represent those waves in an actual field condition with small-scale heterogeneity in amplitude due to inhomogeneous winds. Similarly, Figure 18b shows the perspective view of the developing wave field modeled by the FNHS condition. One can clearly see that the wave front diffracts as it evolves to an extreme event through frequency focusing. Figure 19b further compares the FNHS/QNHS simulated free-surface elevation time-series with the experimental data. Due to different directionality, the spatial diffracting wave crest is crescent-shaped, but with an obvious centerline leading. The FNHS-and QNHS-predicted extreme events respectively yield maximum wave crests of 7.84 cm and 7.71 cm, with a 0.05 m difference in the actual locations. The FNHS and QNHS model results are quite close, consistent with the finding and explanation in spatial focusing waves. Here, the wave-wave interactions in the incipient focusing and diffracting waves are not as strong as those in the nonlinear bi-chromatic waves, while such a physical process is greatly affected by the FNHS/QNHS pressure conditions. Overall, wave directionality does not strongly interfere with the FNHS/QNHS pressure effects on the evolution of directional waves.

Summary and Conclusions
Surface water waves with significant influences upon human activities is an important research topic in the field of coastal and ocean engineering [1][2][3]. In this study, a high-order NHS model [49,51] was used to gain a deeper insight for the effects of non-hydrostatic pressure on surface wave motions and its role in numerical modeling. First, relative contribution of non-hydrostatic pressure P nhs /P hs and its effects (sensitivity) on wave celerity and amplitude were carefully evaluated. For linear wave dispersion [1], non-hydrostatic pressure contributed from the local vertical acceleration, confirmed by the identical FNHS and QNHS results. The vertical structure of P nhs /P hs clearly indicated stronger NHS effects in deeper waters and its significance near the surface. Overall, the NHS effects would slow-down wave celerity, but remain the incident amplitude for the deep-water waves, compared to the unrealistic long-wave patterns in the HS results. The sensitivity could also be easily described by the pressure and celerity errors (e.g., ε p = 3% and ε c = 6% (at Kh = 20)). In terms of wave non-linearity [5], non-hydrostatic pressure revealed similar features, except a larger magnitude and different surface contribution pattern. Non-hydrostatic pressure components caused by the convective vertical acceleration was responsible for the increased amplitude and phase speed under strong non-linearity.
The FNHS/QNHS model [49,51] was then applied to investigate the effects of non-hydrostatic pressure for waves involving complicated interactions (i.e., wave, bottom, and current, as well as directionality). In the bi-chromatic deep-water wave train [9], the FNHS results predict the modulation of non-linear wave groups. Without convective vertical acceleration, the QNHS simulation could not capture side-band frequencies, giving the unrealistic phase speed (or group velocity) and amplitude. For the non-linear wave shoaling example [18], the FNHS results resolve the wave transformation and the QNHS predictions under-estimated wave amplitude and celerity. In terms of wave-current interactions [11], both of the FNHS and QNHS results reasonably describe wave-induced velocity under a negatively sheared opposing current, possibly due to the diminished convective effects. However, unlike the good agreement between the FNHS results and experimental data, the QNHS prediction (similar to the Doppler-shifted solution) was incapable of handling non-linear wave interaction with a strongly sheared following current. Lastly, wave directionality did not strongly interfere with the FNHS/QNHS effects on the evolution of spatial focusing and diffracting waves [66]. The FNHS and QNHS model results are quite close, possibly due to much weaker wave-wave interactions. Overall, the non-hydrostatic effects on wave propagation/transformation and its role in nonlinear dispersive wave modeling have been clearly elucidated.