FULLY NONLINEAR AND DISPERSIVE MODELING OF SURF ZONE WAVES : NON-BREAKING TESTS

With the objective of modeling coastal wave dynamics taking into account nonlinear and dispersive effects, an accurate nonlinear potential flow model is studied. The model is based on the time evolution of two surface quantities: the free surface position and the free surface velocity potential (Zakharov, 1968). The spectral approach of Tian and Sato (2008) is used to resolve vertically the velocity potential in the whole domain, by decomposing the potential using the orthogonal basis of Chebyshev polynomials. The model mathematical theory and numerical development are described, and the model is then validated with the application of three 1DH test cases: (1) propagation of nonlinear regular wave over a submerged bar, (2) propagation of nonlinear irregular waves over a barred beach, and (3) wave generation and propagation after an abrupt deformation of the bottom boundary. These three test cases results agree well with the reference solutions, confirming the model’s ability to simulate accurately nonlinear and dispersive waves.


INTRODUCTION
Coastal and ocean engineers have a growing need for accurate and rapid wave models that are capable of simulating the propagation of waves in the coastal zone, including interactions at the shoreline and with coastal structures.The nonlinear and dispersive characteristics of waves are particularly important in the nearshore region and thus are the focus of the present work.A model based on fully nonlinear potential flow theory is developed as a compromise between computationally expensive CFD (Computational Fluid Dynamics) approaches based on the full Navier-Stokes equations, and Boussinesq, Serre-type or Green-Naghdi models that simplify the vertical structure of the dynamics and are only partially nonlinear and/or dispersive.The model is described in part 2 of this article.The model is then tested and applied to three nonbreaking wave test cases: the propagation of (1) regular nonlinear waves over submerged bar (Dingemans, 1994) and (2) irregular nonlinear waves over a barred beach (Becq-Girard et al., 1999), and (3) wave generation from an impulsive upthrust of the bottom (Hammack, 1973).These test cases are presented and analyzed in part 3.

MODEL PRESENTATION Mathematical model
Assuming irrotational flow of an inviscid and homogeneous fluid with constant density, the velocity potential Φ(x, z, t) (where x = (x, y)) can be defined.This velocity potential must satisfy the Laplace equation in the fluid domain, which is supplemented by boundary conditions at the free surface z = η(x, t) (kinematic and dynamic conditions), at the bottom z = −h(x, t) (impermeability condition), and at the lateral boundaries (Dirichlet or Neumann conditions).By assuming a non-overturning free surface and by setting the atmospheric pressure equal to 0 at the free surface, the kinematic and dynamic free surface boundary conditions are: where is the horizontal gradient operator.By defining the velocity potential at the free surface Φ(x, t) ≡ Φ(x, z = η(x, t), t) and rewriting equations ( 1) and (2) as function of the free surface potential, the Zakharov equations (Zakharov, 1968) are obtained: where the vertical velocity at the free surface w(x, t), is defined as: w(x, t) ≡ ∂Φ ∂z z=η(x,t) . (5) This set of equations governs the temporal evolution of the free surface elevation η(x, t) and the free surface velocity potential Φ(x, t).To integrate equations (3-4) in time, it is necessary to determine the vertical velocity at the free surface w(x, t) from the two surface quantities η(x, t) and Φ(x, t), a problem typically called "Dirichlet-to-Neumann" (DtN).

Numerical model
In this work, the model (3-4) is limited to only one horizontal dimension (i.e.x = x).The numerical simulation of the mathematical model requires (a) an algorithm to compute first and second order spatial derivatives, (b) a temporal scheme to advance in time, and (c) a method to solve the DtN problem.The three components are chosen as follows: (a) High order finite difference approximations (fourth order here) are used to calculate first and second order derivatives in x following the method of Fornberg (1988).The domain is discretized in N X nodes in the x dimension that can be regularly or irregularly spaced (constant or variable ∆x).
(b) The explicit four-step fourth-order Runge-Kutta scheme with a constant time step is used to integrate (3-4) in time.
(c) The numerical resolution of the DtN problem is achieved by solving the Laplace boundary value problem for the velocity potential Φ in the entire domain, using a spectral method in the vertical dimension (z), following the work of Tian and Sato (2008).First, the fluid domain is transformed into a rectangular domain by introducing the vertical coordinate s, which varies from s = −1 at the bottom to s = 1 at the free surface: where h + (x, t) = h(x, t) + η(x, t) and h − (x, t) = h(x, t) − η(x, t).Following Tian and Sato (2008), the vertical variation of the velocity potential then can be approximated by a linear combination of Chebyshev polynomials of the first kind (T n ): where N T is the maximum order of the Chebyshev polynomials.After introducing (7) in the system of equations solving the Laplace boundary value problem (including the boundary conditions), and applying the Galerkin method to eliminate the s coordinate, a new linear system of equations is obtained for the a n (x) coefficients.Once the system is solved for the coefficients a n (x), the vertical velocity at the free surface w(x) can be estimated by: Finally, equations (3-4) are stepped forward in time.
The DtN problem is resolved at each step in the Runge-Kutta algorithm using the direct linear solver MUMPS (MUltifrontal Massively Parallel Solver) (Amestoy et al., 2001).The size of the linear system depends on the number of nodes in the x dimension (N X ) and on the maximal order of the Chebyshev polynomials (N T ), with N X (N T + 1) equations.

VALIDATION TEST CASES
The model is validated by applying it to three test cases.The first two cases simulate the propagation of respectively regular and irregular nonlinear waves over a submerged bar, and the third test case simulates waves generated by an impulsive vertical motion of the bottom to study tsunami-like wave dynamics.

Regular waves over a submerged bar
This first test case consists of simulating the propagation of regular nonlinear waves over a submerged bar, based on the wave flume experiments of Beji and Battjes (1993) and Dingemans (1994).Waves are generated at the left domain boundary and propagate over a trapezoidal bar (Figure 1).The results presented here correspond to case A (H = 2.0 cm and T = 2.02 s) of Dingemans (1994), with long incident waves (kh = 0.67 offshore) of relatively small steepness ( = kH/2 = 0.017).The model domain has a regularly spaced grid with ∆x = 0.05 m, and the waves were propagated during 25 periods (i.e.50.5 s) with a time step ∆t = T/100 = 0.0202 s.Waves were generated in a 7-m wide relaxation zone at the left boundary of the domain using linear wave theory.An additional 7-m wide relaxation zone was applied at the right boundary as an absorption zone to avoid reflexions.Simulations with N T = 3, 4, 5, 7 and 10 were completed to evaluate the impact of this parameter on the results.Eleven wave probes were used in the experiments (Figure 1), and Figure 2 shows a comparison of the simulated and observed wave time series at five of these probes (probes 4, 6, 8, 10, 11).When the water depth decreases, the wave height and steepness increase (see probe 4 in Figure 2a).The wave profiles are asymmetric due to nonlinear interactions that create higher frequency bound harmonic components.These harmonics are then released in the shallowest region and on the back slope of the bar, and then propagate at their own phase speed (see probes 6 and 8 in Figure 2b-c).After the bar, the wave profiles measured at each probe vary significantly from one location to another due to the differences in free wave celerity (see probes 10 and 11 in Figure 2d-e).At the two last probes, the model performs well in reproducing the complex wave profiles, including the dispersive (high frequency) components.
The maximum order N T of the Chebyshev polynomials also effects the results, as seen in Figure 2. Along the tank up to the end of the submerged bar (i.e. up to probe 8), only the results of the run with the lowest order N T = 3 differs significantly from those with higher N T , which are superimposed and agree well with the measured time series.At the last three probes, where dispersive effects are more important, the results with N T = 3 deteriorate, and results with N T = 4 and 5 also exhibit some differences with the measurements.Results with N T = 7 and 10 remain however superimposed and in good agreement with the measurements.This validates the ability of the model to simulate this nonlinear and dispersive case with N T values as low as 7.

Random waves over a barred beach
The second test case simulates the propagation of irregular nonlinear waves over a submerged bar, based on the wave flume experiments of Becq-Girard et al. (1999).In the experiments, waves are generated with a piston-type random wave-maker using a JONSWAP wave spectrum with a peak-enhancement factor of γ = 3.3.The bathymetric profile of these experiments (a barred beach, Figure 3) was specifically designed to study nonlinear wave effects in shallow water, in particular triad interactions.The bathymetric profile was made with smooth metal sheets, so the bottom friction dissipation was assumed negligible, and wave reflections from the beach were minimized with a beach absorber on the upper part of the beach.Resistive-type wave probes measured the free surface elevation at 16 locations in the domain (Figure 3) during the 40-minute experiment with a sampling time step ∆t = 0.07 s.Non-breaking irregular waves are simulated with a significant wave height of H m0 = 3.4 cm and a peak period of T p = 2.39 s in the deepest part of the domain (h = 0.65 m).The computational grid extends from x = −5 m to x = 25 m (with the foot of the bar at x = 0 m), with a regular mesh with ∆x = 0.05 m (601 nodes).The total simulation time is 2380 s (approximately 39.7 min), with a time step equaling the sampling time step of the free surface elevation probes, ∆t = 0.07 s.The model is forced at the left boundary and in a 5-m wide relaxation zone with a time series of the free surface position and wave potential, which was reconstructed using linear wave theory to sum the components of the wave spectrum obtained from the free surface measurements at probe 2 (located at the foot of the submerged bar).A 10-m long relaxation zone was added at the right boundary to act as an absorption zone.For these tests, the maximum order of the Chebyshev polynomial was set to N T = 7.
From the simulated free surface elevation time series, the wave variance spectrum at each probe is compared to the wave variance spectrum obtained from the experimental data (Figure 4, spectra at probes 2, 4, 7, 9, 11, 13, 15 and 16).With a decrease in water depth, the main spectral peak increases due to shoaling effects between probes 2 and 4. In addition, in this zone, energy is transferred from lower to higher frequencies, particularly from the peak frequency to its super-harmonics.The second harmonic peak (2 f p ), which was hardly visible at probe 4, is well developed at probe 7, as are the third (3 f p ) and fourth (4 f p ) harmonic peaks.In the spectra at probes 9 and 11, a peak at the fifth (5 f p ) harmonic appears, and its amplitude is well represented by the model.When the water depth increases again on the back side of the bar, one can observe a decrease in energy of the higher harmonics, as this energy is transferred back to lower harmonics (in particular to the second harmonic).At probe 13, the signature of the fifth harmonic disappears, and at probe 15, the peaks of the third and fourth harmonics also disappear, leaving only the second harmonic and main spectral peaks visible.At probe 16, as water depth decreases again, nonlinear interactions cause an increase in amplitude of the third (3 f p ) harmonic peak., 4, 7, 9, 11, 13, 15 and 16 (see Figure 3 for the probe positions).The frequency scale is normalized by the peak frequency ( f p ) to identify harmonic peaks of frequency (e.g. 2 f p , 3 f p , ...).
Nonlinearities also modify averaged integral parameters such as the significant wave height (H m0 ), the mean energetic wave period (T m−1,0 ), and the horizontal and vertical asymmetry.The model globally estimates well these four parameters (Figure 5).The significant wave height increases as the water depth decreases (shoaling effects), before decreasing in the trough, and then finally increasing again as the water depth decreases nearing the beach.Along the bathymetric profile, the simulated H m0 agrees well with the measured values, with a maximum difference of 8.8%, which appears after probe 7, where the simulations overestimate the significant wave height.The observed shoaling phenomenon is accompanied by a reduction of the energy contained in the low frequency spectrum part and an augmentation for the high frequencies resulting in a decrease of the mean period (T m−1,0 ).The subsequent release of higher harmonics in the trough leads to an increase in the mean wave period that persists along the tank.The largest differences in mean period occur at the last three probes, but only represent an error of less than 3.5%.Spatial evolutions of horizontal and vertical asymmetry along the bathymetric profile agree well with the experimental data.This test case validates ability of the model to simulate the generation, propagation, and absorption of irregular, non-breaking waves, including nonlinear processes such as wave shoaling and wave interactions transferring energy back and forth between higher and lower harmonics.

Waves generated by impulsive bottom motion
The third test case investigates waves generated by an impulsive vertical motion of the bottom to study tsunami-like wave dynamics.Hammack (1973) studied the case of a moving bottom with physical experiments, an analytical study of the linear solution, and numerical simulations with a KdV (Kortewegde Vries) model.Fuhrman and Madsen (2009) (abbreviated F&M09 hereafter) also tried to reproduce numerically these experiments with a high-order Boussinesq model.Both upthrust and downthrust of the bottom were studied using different expressions for the temporal evolution of the bottom movement.Here, impulsive exponential bed upthrust is simulated and compared to the results of F&M09.The fluid domain is initially motionless, with an undisturbed water depth h 0 = 1 m.During the simulation, a step in the bottom between x = 0 m and x = b m is rapidly upthrust (Figure 6, here b/h 0 = 12.2), and the water depth (h(x, t)) evolves following: where H is the Heaviside step function.The step is upthrust over a total vertical distance of ζ 0 /h 0 = 0.1, and the exponential decay constant and critical time are respectively α = 1.11/t c and t c = 0.148/ gh 0 .The domain is discretized on a regular grid extending from 0 to 2500h 0 = 2500 m, with a spatial step ∆x = h 0 /5 = 0.20 m (12501 nodes).The simulation length is the nondimensional time t g/h 0 = 2375 (i.e.t ≈ 75828 s), with nondimensionnal time step ∆t g/h 0 = 0.20 (i.e.∆t ≈ 0.0638 s).The resulting ICCE 2014 CFL number is CFL = gh 0 ∆t/∆x = 1.Fully reflective vertical boundaries are applied at both ends of the domain.The maximum order of the Chebyshev polynomial is N T = 7.Time series of the free surface elevation at four positions in the domain ((x − b)/h 0 = 0, 20, 180, and 400) are compared to the experimental measurements (red line) and to the analytical integral solution (dotted blue line) based on the linearized potential flow equation of Hammack (1973) in Figure 7.
Near the generation zone (Figures 7a and b), the simulated free surface position (black line) agrees well with the global trend of the experimental data, with a slight overestimation at (x − b)/h 0 = 0 and with larger oscillations at (x − b)/h 0 = 20.These results are similar to those of Figure 2 in F&M09.Farther from the generation zone (Figures 7c and d), the differences are more significant.The linear solution diverges from the experimental measurements, showing that a model with nonlinear properties is required to simulate properly this case.The simulated results (with the present model) are closer to the measurements than to the linear solution, confirming the nonlinear ability of the model.A comparison of the far-field waves shows that they have the same global shape, but that the current model predicts faster wave propagation speeds (relative rightward shift in the time series), which is consistent with the overprediction of the wave amplitude.Hammack (1973) suggested that differences in simulated and measured wave amplitudes could be due to viscous energy losses occurring in the experiments but were neglected in both numerical models.
When the wave train propagates over long distances (up to 2500 times the water depth), the leading waves separate into two solitary waves.These solitary waves then propagate with constant shape, with the first wave traveling faster due to its higher amplitude, as observed by F&M09.In Figure 8, the free surface profile at t g/h 0 = 2375 (black) is compared to the two theoretical solitary wave profiles obtained from  the Boussinesq equations, fit to the final results of F&M09 (see their Figure 3).The current simulations predict slightly faster propagation speeds of slightly higher amplitude solitary waves.
Comparisons of the simulated free surface time series with the experimental results of Hammack (1973) confirm the ability of the model to reproduce accurately the wave disturbance dynamics, including the formation of leading waves leaving the generation zone and their separation from the wave train after a sufficient length of time, as observed by F&M09.

CONCLUSIONS AND OUTLOOK
A high-order potential flow model is presented for modeling highly nonlinear and dispersive nonbreaking waves over irregular bottom profiles.Three test cases (in one horizontal dimension) confirm the model's ability to simulate accurately nonlinear and dispersive effects such as the transfer of energy to sub and super-harmonics (up to 5 f p ), as shown in the second test case, and soliton fission after long-term propagation of an impulsive wave from bed upthrust, as illustrated in the third test case.The benefit of solving the DtN problem with the proposed method is that it has exponential convergence in the vertical with a flexible level of accuracy depending of the choice of N T , with accurate and robust order 4 numerical methods in both x and t (finite difference schemes and Runge-Kutta scheme), yet with moderate computational cost.For all of the simulations performed, highly accurate results can be obtained with N T in the range of 5 to 7.
Ongoing work includes model validation with additional 1DH test cases and extension of the model to 2DH using unstructured meshes.Additional physical processes will also be parameterized or implemented in the model, such as wave breaking and run-up at the shoreline, with the goal of representing various types of coastal and harbor structures.

Figure 8 :
Figure 8: Free surface profile at t g/h 0 = 2375.Comparison with Fuhrman and Madsen (2009) analytical solitary waves fit to their final results (their figure 3) (green) and the current model simulation results (black).