A Set of Accurate Dispersive Nonlinear Wave Equations

: In this study, a set of accurate dispersive nonlinear wave equations is established , using the wave velocity and free surface elevation as variables. These equations improve upon previous equations in which the velocity potential is used as a variable by considering the rotational wave motion and by adding a second-order bottom slope term that applies to general situations, allowing the equations to consider the inﬂuence of rapidly changing, horizontal, two - dimensional bottom topographies. The problem of the inaccuracy of the integral calculations used in previous equations in nearshore areas is solved by approximating the integral calculations into diﬀerential calculations, and a set of coupled wave equations is established by keeping the free surface elevation and the horizontal velocity constant, thus allowing the calculation of nearshore wave -generated currents. The beneﬁts of the current model are conﬁrmed through comparisons with corresponding laboratory experimental ﬁndings and are illustrated throu gh a comparison with the numerical outcomes of other pertinent models.


Introduction
Waves are an important hydrodynamic factor in marine environments.In coastal engineering and certain practical problems, it is necessary and crucial to accurately determine the wave field (the spatial and temporal distribution of wave elements such as wave height, period, and wave direction) of wave propagation and deformation.Wave equation calculation models are a basic research topic in ocean engineering, and the results obtained from this research not only have important academic value for further understanding of this natural law, but also provide technical support for research on coastal engineering, shipping, offshore marine environmental protection, and marine resource development.
The commonly used wave equations in ocean engineering are the mild-slope equation [1][2][3][4][5] and the Boussinesq equation [6][7][8][9][10][11][12].However, the mild-slope equation only considers regular waves or narrow-spectrum irregular waves.The dispersion of the Boussinesq equation is approximate, but it can be improved by adding derivatives of an order of five or higher or via water depth stratification models [13,14]; however, this also increases the computational load and introduces difficulties into the numerical calculations.The stability and computational accuracy of certain higher-order derivative terms in the high-order Boussinesq equation can also be problematic in areas with complex terrain.Improving previous equations improves performance but still does not fundamentally achieve complete accuracy in equation dispersion.Therefore, it is necessary to continue to develop accurate and dispersive nonlinear wave equations so that they can be applied to a wide range of frequency domains with broad spectra.
There are two methods of developing the classical gradual-slope equation into a method suitable for broad-spectrum analysis in the current research.One is the discretespectrum multi-mode method, and the other is the continuous-spectrum method.The multi-mode method uses multiple sine wave modes to construct the velocity potential or velocity, and it establishes equations based on this.Therefore, the resulting equations are coupled equations of multiple modes, which require large amounts of computation.Nadaoka et al. [15] established a set of multi-mode mild-slope equations with weak nonlinearity, which can provide accurate results for irregular situations when the number of modes is 2 N ≥ .Tang and Ouellet [16] established a multi-mode, weakly nonlinear, mildslope equation expressed via the velocity potential.Belibassakis and Athanassoulis [17] established a multi-mode, slow-slope equation expressed via the velocity potential by applying the Luke variational principle.In recent research, Kim and Kaihatu [18] also proposed a multi-mode, parabolic, mild-slope equation.Compared to the above multi-mode models, the continuous-spectrum model has the advantage of a relatively small computational complexity, as the consideration of continuous spectra is implemented through the integration of the kernel functions, and the corresponding model only has a set of equations.Schäffer [19] established a set of linear wave models with accurate dispersion for the case of horizontal, one-dimensional conditions, based on the increase in the wave surface and velocity at the still-water surface.The model was established by transforming the relationship between the vertical velocity component and the horizontal velocity component at the still water surface from a wave-number space to a physical space.Von Groesen [20] used consistent variational modelling to obtain and verify an accurate model for uni-directional surface water waves.Mei [21] derived a modified Zakharov equation for a simple nonlinear dispersion system that takes into account the random scattering of broadband waves.Natanael [22] provides an alternative methodology for analysis of three-wave interactions under the exact dispersion relation associated with gravity waves in fluid of intermediate depth.Karambas and Memos [4] established a set of continuous-spectrum models of the same type but with weak nonlinearity.The application of these models can simulate irregular wave propagation on submerged breakwaters.Zou et al. [5] established a set of horizontal, two-dimensional, nonlinear wave equations with accurate dispersion.The nonlinearity is approximated to the fourth order, and the second-order water-slope term in the one-dimensional case is also considered.The second-order water-slope term is established by extending the method of adding a rapidly varying water-depth term in the classical mild-slope equation proposed by Kirby [23] to a wide-spectrum model.The model also provides the corresponding wave dispersion and second-and third-order nonlinear response functions through perturbation analysis.Numerical simulations of Bragg reflection experiments on sandbar terrain were performed, and the effectiveness of the second-order water bottom slope term was verified.Except for Zou et al. [5], who considered the second-order bottom slope term, other models only consider the first-order bottom slope term.For the establishment of the first-order bottom slope term, see Bingham and Agnon [24].
In this study, we establish a new set of nonlinear wave equations with accurate dispersion by using the wave velocity and free surface elevation as variables.This equation only approximates nonlinearity in the continuity equation, while the momentum equation is completely accurate without any approximation.The equation does not make any assumptions about water depth during the derivation process, so it applies to irregular waves under all water depth conditions.This equation overcomes the shortcomings of previous equations that are unable to consider rotational wave motion and adds a secondorder water bottom slope term that applies to general situations.The establishment of a coupled model through the continuity of free surface elevation and velocity can simulate physical phenomena such as nearshore circulation.The organization of this study is as follows: Section 2 presents the derivation processes of the equations; Section 3 extends the equations in the nearshore region and establishes the coupled model; Section 4 presents the numerical method; Section 5 validates the model with relevant laboratory experiments; and the conclusions are drawn in Section 6.

Derivations of Government Equations
The government equations established in this section include the momentum equation and the continuity equation.The derivation of the momentum equation in Section 2.1 is based on expressing the Euler equation at the free surface.The derivation of the continuity equation in Section 2.2 is based on the free surface boundary conditions, and the equation is closed by approximating the vertical velocity with horizontal velocity and free surface elevation.The continuity equation includes the derivation of a new second-order bottom slope term.

Derivation of Momentum Equation
A Cartesian coordinate system with the origin at the still-water level and the z-axis pointing vertically upward is adopted.The fluid is assumed to be uniform, ideal, and incompressible, and the surface tension is neglected.To make the equation suitable for rotational flow, it is derived from Euler's equation: where is the horizontal velocity vector, w is the vertical velocity, p is the fluid pressure, ρ is the fluid density, g is the acceleration due to gravity, t is the time, is the horizontal gradient operator, and η is the wave surface eleva- tion.In order to express (1) and ( 2) at the free surface, the derivative relation is introduced, as follows: where q  is the free surface variable.Using the above formula, Equations ( 1) and ( 2) can be expressed as follows: The pressure term on the right side of ( 5) can be eliminated by substituting (6) into (5) and using Then, using the boundary conditions at the free surface, and eliminating w  from (7) yields the following: Equation ( 9) is the momentum equation to be established.

Derivation of Continuity Equation
The establishment of the continuity equation is based on (8), which requires expressing 0 w in terms of  u .This can be performed by first expressing w  in terms of the vertical velocity ( 0 w ) and horizontal velocity ( 0 u ) at the still-water level, and then expressing 0 u in terms of  u .The specific procedure is given below.
The Taylor expansions of w and u on the still-water level are as follows: Using the continuity equation and the equation of zero horizontal vorticity, Let w and u be represented by 0 w and 0 u , as follows: In the formula, when 2n ∇ acts on a vector, the divergence is calculated first, and then the gradient is calculated.When 2n ∇ acts on a scalar, the gradient is calculated first and then the divergence is calculated.Taking z η = in ( 14) completes the first step; w  is expressed in terms of 0 w and 0 u .
To complete the second step, it is necessary to establish the relationship between 0 w and 0 u .The following vertical integration of Equation ( 12) is performed, and using the Leib- niz formula and the water bottom boundary condition ( h w h − = ⋅∇ u ) yields the following: In order to obtain the expression of 0 w , which contains the second-order term of the water bottom slope required for this study, it is necessary to determine the expression of u in (16).Here, the expression of u needs to be accurate to the first-order term of the water bottom slope.The expression used, given by Xu et al. [25], is as follows: where . Substituting ( 17) into ( 16) yields the following: where To separate the water-depth gradient term included in (18) by order, 0 w is expanded with the water-depth slope ( h ∇ ) as a small parameter, as follows: where (2) 0 w are the zero-order, first-order, and second-order bottom slope terms, respectively; W are their Fourier transforms; ( , ) is the wave-number vector; and  22) and ( 23) into (21) yields the following: ], In the above formula, when ∇ acts on 0 U , it is necessary to consider the change in 0 U with the water depth.In order to obtain the expression of 0 U , it is necessary to express (8) at the still water level and retain the linear term, substitute , and take the divergence of . The expressions of 0 U and 0 ∇ ⋅U are substituted into (10), and the zero-order, first-order, and second-order bottom slope terms are separated.The zero-order term is as follows: The first-order term is as follows: The second-order term is as follows: where . The expression for s α is given by (A5) in Appendix A. The first-order term above can be written in the following simple form (see Bingham and Agnon [24]): Zou et al. [5] performed the inverse Fourier transform on the zero-order and firstorder terms, obtaining the Fourier integral expressions for (0) 0 w and (1) 0 w in the xy plane: where is the integral point coordinates, and 0 ( ) are the zero-order and first-order kernel functions, the expressions of which are as follows: where is the zero-order Bessel function of the first kind.The second-order term can also be expressed in a form similar to the one above.To achieve this, the results of (27) need to be organized in the following steps: (1) Separate ∇ from A , so that A is not affected by ∇ , such as where A ∇ can be expressed in terms of A using Equation (A4) in Appendix A. (2) Combine like terms (that is, combine all the terms into three major terms: ∇ , and h ∇ ⋅k ).In order to facilitate the establishment of a second-order kernel function, the same small-angle assumption as the one in Appendix A is used for h ∇ ⋅k , which can be approximately expressed as cos sin (because α and h ∇ are of the same order under the small-angle assumption, and cos 1 α ≈ and sin 0 α ≈ can be taken in the above expression when approximating to the second-order bottom slope term).(3) Express A in terms of 0 U , so that the resulting expression for (2)   0 w is expressed in terms of 0 u .The approach is to ignore the term A ∇ proportional to h ∇ in the expression for 0 U described above, resulting in 4) Perform an inverse Fourier trans- form on (27), expressing the result in the wave-number domain as the result in the xy domain.The following substitutions can be made during the Fourier transformation: After step (4), the expression of (2)   0 w can be obtained as follows: (2) 2 2 0 0 21 0 22 (35) where 21 G , 22 G , 23 G , and 24 G are the second-order kernel functions, and the expressions are as follows: Substituting the above results into Equation ( 21), 0 w can be expressed in terms of the integral-differential operator E (i.e., 0 0 E w = − u ), which is given as follows: where ϕ represents the variable being integrated.Based on the relationship between 0 w and 0 u established above, 0 u is expressed as  u .Substituting 0 0 E w = − u into (14) and taking z η = , we can obtain w  at the free surface.This allows w  to be expressed in terms of 0 u .Then, taking z η = in (15), the right-hand side expression of the equation is  u .Through continuous iteration, 0 u is expressed in terms of  u .Substituting the result- ing expression for 0 u into (14) gives an expression for w  , which allows ( 14) to be com- bined with the nonlinearity to the fifth order (

( )
O η ).Substituting the resulting expression for w  into (8) gives a continuity equation, which can be written together with the momentum equation to obtain the wave equation to be established, as follows: where where 2 F , 3 F ,and 4 F , are the second-, third-, and fourth-order nonlinear terms, respectively, and 1 F is the linear term, which is introduced to maintain the similarity between the left-hand side of the equation and the continuous equation in integral form.To further illustrate the characteristics of the equation, the following transformation is made by expressing  u in terms of the velocity potential ( η φ ) at the free surface, which is . Equations ( 41) and ( 42) are transformed into accurate dispersive nonlinear wave equations expressed by the velocity potential at the free surface established in the literature (Zou et al. [5]), as follows: where where E =E ′ ∇ .The accurate dispersion, nonlinearity (approximated to the fourth order), and shallowing effect of (47) and (48) have been proven in the literature (Zou et al. [5]), that is, the dispersion, nonlinearity (wave amplitude dispersion, second-order and thirdorder transfer functions, etc.), and shallowing effect of the equation are consistent with theoretical Stokes waves (the shallowing performance satisfies the wave energy conservation equation).The transformation from ( 41) and ( 42) to ( 47) and (48) indicates that the equations ( 41) and ( 42) established in this study under potential flow conditions are consistent with the equations (47) and ( 48), with the velocity potential as a variable.Therefore, the dispersion, nonlinearity (approximated to the fourth order), and shallowing effect performance of Equations ( 41) and (42) established in this study are completely accurate and consistent with theoretical Stokes waves.Unlike the latter set of equations, the newly established equations are derived based on Euler equation, assuming only that the horizontal vorticity is zero.Therefore, the flow considered by the equations can have vertical vorticity, that is, the motion in the horizontal plane can be rotational, which is one of the advantages of the new equations.

Nearshore Coupling Equations
We established the accurate dispersive wave equation above.This equation is of the integral-differential type, the integral calculation of which requires a symmetric integral domain relative to the control point.The calculation of the integral value at the control point requires the values of the integrand function at each point in the integral domain.Therefore, this equation cannot be used to calculate problems such as waves climbing on the coast with dynamic boundaries.In order to solve this problem, it is necessary to convert the integral kernel function into a differential form in the nearshore region so that the integral-differential calculation contained in the equation only contains the differential calculation.Based on Equations ( 41) and (42), in this section, we present a set of weakly dispersive coastal wave equations that make up for the lack of accurate dispersive wave equations in coastal boundary calculations.By coupling these two equations, a set of coupled wave models is established, where the accurate dispersive wave equation is used for the calculation of areas with relatively large water depths, and the coastal wave equation is applied to the calculation of nearshore areas.

Equation Derivation
In order to convert the integral kernel function into a differential form in the nearshore region, the continuous Equation (42) needs to be degenerated in the nearshore region into the coastal wave equation that applies to the nearshore region.Because the coastal wave equation is used in the coupled model for the calculation of the nearshore shallow water region, here, the dispersion of Equation ( 41) is refined to 2 ( ) O µ , and all the nonlinear terms are retained, resulting in the following: In order to simplify the kernel function, the hyperbolic function contained in the kernel function is Taylor-expanded, which obtains cosh 1 2( ) ( ) / 24 kh kh kh = + + , using the following formula: Transforming ( 31) and (32) into the following form and retaining terms up to

( )
O µ , we obtain the following: Substituting ( 55) and ( 56) into (40), we obtain the following: where q is the integrand.Through the above derivation, the operator E is converted from the integral-differential form to the differential form.Substituting Equation (57) into Equation (53) and retaining the terms up to

( )
O µ , we obtain the following: According to Xu et al. [25], this equation is degenerated into a weakly dispersive equation.Similarly, after expanding and approximating (56), we can obtain the following: Equation ( 59) is the continuity equation in the coastal equation.Because the momentum (43) contains a quadratic derivative term of the free surface elevation with respect to time, which is prone to numerical oscillations at the coast, Equation ( 44) is converted at the coast by separating the vorticity contained in the equation, which yields the following: where × is the cross-product.By substituting Expression (59) into the expression above, only up to

( )
O µ , we obtain the following: The weakly dispersive coastal wave equation established in this section consists of Equations ( 59) and (61).This equation is based on the accurate dispersion equation established in Section 2, converting the integral kernel function into a differential form.

Matching of Coupling Equations
In this section, the wave Equations ((41) and ( 42)) with accurate dispersion established in Section 2 for any water depth are coupled with the coastal Equations ((59) and (61)) presented in Section 3.1 for shallow coastal waters to form a new coupled model.The planar layout of the coupled computational domain is shown in Figure 1, with the origin of the Cartesian coordinate system located at one end of the incident boundary and the positive direction of the x-axis consistent with the wave propagation direction.In this section, the computational domain is divided into three regions: region 1 is the computational domain for the accurate dispersion wave Equations ((41) and (42)), which has no limitation on the depth range and is responsible for the main computational work; region 3 is the computational domain for the coastal Equations ((59) and (61)), which is responsible for the computation of the shallow-water areas near the coast; and region 2, composed of both regions, is the coupling region, which is responsible for the exchange of variables.The interface between region 1 and region 3 is the coupling interface, as marked in Figure 1.
Ω is responsible for providing the variables required for the integration and differentiation of the accurate dispersive wave equation for region 1, while region 1 Ω is responsible for providing the variables required for the differentiation of region 2. In addition, boundary layers are set at the virtual grid cells on both sides to provide variables for the kernel function integration.L1 is the grid number required for the coastal equation differential, and L2 is the grid number required for the accurate dispersion equation integration.to region 3, as follows: In order to satisfy the requirements of the integral calculation for region 1, the 2  u and 2 η variables, calculated via the coastal wave model, pass from region L1 to region 1, as follows: The calculation steps of the coupled model are as follows.Ω to time layer n .Repeat the above steps until the solution is com- plete.

Numerical Methods
The coupled equations are solved using finite difference methods.From the expressions of the two sets of equations, it can be seen that the continuity Equation (41) in the accurate dispersion wave equation is in an explicit form, while the momentum Equation (42) contains t u  and t v  simultaneously in a component equation, which cannot be di- rectly solved explicitly.It is necessary to combine the two component equations, taking them as variables to be solved, and to obtain their expressions separately before deriving explicit solutions.The continuity Equation (41) and momentum Equation (42) in the coastal wave equation only contain a linear time term; thus, they can also be solved explicitly.The coupled equations have the following common form: where * The specific expressions for these terms can be found in the FUNWAVE model [11].In terms of the numerical format, it is necessary to ensure both high numerical accuracy for discretizing the equations and stability during computation.Therefore, the equations are discretized using a fifth-order Adams-Bashforth predictor-corrector Adams-Moulton (ABM) scheme in time and a five-point finite difference scheme in space [6][7][8][9][10][11][12].Due to the complexity of the equations and numerical implementation, it is difficult to analyze the stability and convergence of the numerical model.In order to maintain computational stability in numerical calculations, the numerical filtering method proposed by Shapiro [26] is used here, with the expression where f represents the original value of each calculation point, and f * represents the new value after filtering.Numerical filtering can ensure the stability of numerical calculations.This article uses the Fortran programming language's numerical simulation.
The boundary conditions used in the numerical solution process include a fully reflective boundary, periodic boundary, and sponge layer absorption boundary.For the treatment of a moving coastal boundary, here, we use the narrow-slit method [11].The numerical wave making method is divided into boundary wave making and internal wave making.Boundary wave making is based on the theoretical solution of waves, providing the velocity and free surface elevation at the incident boundary to achieve wave making.The advantage of this method is that the waves generated through wave making are nonlinear waves, which can be used to study nonlinear problems more easily.Their disadvantage is that when encountering structures, the incident waves will be reflected, forming secondary reflections at the wave making boundary.Internal wave making adds a mass source term to the right side of the fluid continuity equation.The mass source terms of various waveforms can be obtained by integrating the wave surface equation.Since incident waves and reflected waves flow out from the left and right boundaries, respectively, the secondary reflection problem can be effectively eliminated.However, there is the problem of waveform instability when the relative wave height is large.This study adopts the internal wave making method adopted by Kim et al. [27], which is expressed as .When 1 exp( 5) S < − , the effect of the source func- tion can be ignored.Adding s f to the right side of Equation ( 62) can achieve internal wave generation.

Results
The wave equations established in this study are characterized by their accurate dispersion.In order to reflect this characteristic in the numerical results, in this section, we compare the simulation results of the high-frequency case with those of the high-order Boussinesq equations to illustrate that more accurate numerical results can be obtained with the former model than with the latter.The high-order Boussinesq equations selected in this section are a set of high-order Boussinesq equations derived from Zou and Fang [12] based on Euler's equation.The dispersion of this equation is similar to that of Padé [4,4] and retains the fully nonlinear equation of In this section, we compare the results of the N4D4 model with those of the higherorder Boussinesq equation by considering three situations: (1) the theoretical solution of the fifth-order Stokes wave; (2) a one-dimensional submerged breakwater experiment; and (3) a two-dimensional shoal experiment.In the comparison, the same difference schemes, spatial step sizes, and time step sizes are used for the numerical calculations.The accuracy of the dispersion in the calculated results of this model also depends on the accuracy of the difference scheme.Therefore, in this section, we present the different numerical schemes used in this study to investigate the impact of their accuracies on the calculated results.In addition, we also present the application of the second-order bottom slope term added to the equation to simulate the propagation (Bragg reflection) of waves in rapidly varying water depths.Then, the model is applied to simulate the propagation and breakup of waves on a coast with a flat slope verify its climbing performance.Finally, the model is applied to simulate coastal currents and rip currents under regular and irregular wave conditions to investigate its applicability to this situation.

Comparison with the Higher-Order Boussinesq Equation
(1) The theoretical solution of the fifth-order Stokes wave: in order to prove that the accurate dispersion characteristics of this model can provide better computational results at high frequencies than the high-order Boussinesq equation, the following comparison was made between the wave surface elevations of the two models and the theoretical solution of the fifth-order Stokes wave for a water depth of 0.4 m, an incident wave height of 0.04 m, and periods of 0.5 s and 1.0 s, respectively.The theoretical solution of the fifthorder Stokes [28] wave was used as the incident wave at the incident boundary.The spatial step size was / 40 is the wavelength), and the time step size was / 80 t T ∆ = .
As can be seen from Figure 2, there is a phase difference between the free surface elevation of the high-order Boussinesq equation and the theoretical solution when T = 0.5 s and kh = 5.89, while the present model is in full agreement with the theoretical solution.When T = 1.0 s and kh = 1.72, both models' wave surface elevations are in agreement with the theoretical solution, indicating that the high-order Boussinesq equation is also accurate for this larger wave period.The following further comparison of the wavelength (the distance between two adjacent peaks) in the calculation with the theoretical solution further demonstrates that the present model is superior to the high-order Boussinesq equation in high-frequency situations.Table 1 gives the wavelengths corresponding to the calculation results in Figure 1, as well as the theoretical Stokes wavelengths (with the theoretical wavelengths of the present model) and the wavelengths obtained from the dispersion relation (

) of the high-order
Boussinesq equation (N4D4), which agree with those of Padé [4,4].Comparing these results, it can be seen that, compared to the theoretical wavelength of the fifth-order Stokes wave, both the present model and the high-order Boussinesq equation have an error of approximately 0.4% when T = 0.5 s and T = 1.0 s, respectively.This indicates that the error of the high-order Boussinesq equation is also small at T = 1.0 s, which is close to that of the present model.However, at T = 0.5 s, the error of the high-order Boussinesq equation is larger, far exceeding that of the present model.The reason for this larger error can be directly seen in the table, as the wavelength given by the dispersion relation of the high-order Boussinesq equation in the table also has an error of the same magnitude as the Stokes theoretical wavelength, indicating that this larger error is due to errors in the dispersion relation of the high-order Boussinesq equation.The dispersion property of this model is completely accurate, so there is no such error in calculating wavelengths, even in the highfrequency situation when T = 0.5 s.   (2) One-dimensional submerged breakwater experiment: In order to prove that the model is more accurate than the high-order Boussinesq equation for simulating wave propagation in complex topographies, the following comparison with the experimental results [29] of the wave propagation on a submerged breakwater is presented (in the State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology, China).The experimental arrangement is shown in Figure 2, with a front slope of 1:10 and a rear slope of 1:5.The experimental wave conditions were a water depth of 0.4 m, an incident wave height of 0.02 m, and a period of 1.3 s.In the calculation, the spatial step length was 0.02 m x ∆ = , the time step length was 0.01s t ∆ = , and kA ε = .
Figure 3 shows the experimental free surface elevation results at six wave gauges from the top of the submerged breakwater to 15 m behind it alongside the calculated results of the two models.Due to the close agreement between the two calculated results at the front slope and the flat bottom in front of the breakwater, they are not shown here.From the results shown in the figure, for the two wave gauges on the top of the breakwater (x = 8.5 m and 9.5 m), the calculated results of the two models are closer to the experimental results, indicating that this model does not exhibit a better computational accuracy than the high-order Boussinesq equation.This is due to the fact that the higher-order harmonics in the time series of the free surface elevation at the breakwater top mainly exist in the form of bound waves, which correspond to special solutions of the wave equation and have the same propagation velocity as first-order waves.Therefore, the error in the propagation velocity caused by higher-order harmonics is masked by first-order waves and cannot be revealed.However, this is not the case for the wave conditions behind the submerged breakwater slope and behind the slope.The calculated results of the higherorder Boussinesq equation show larger errors than those of this model, as is clearly shown in the enlarged portion of the figure.The reason for these larger errors is that the higherorder harmonics in the region behind the submerged breakwater slope and behind the slope not only contain bound waves but also generate large-amplitude free waves, which correspond to homogeneous wave equation solutions and have different propagation velocities than first-order waves.Therefore, their wave surfaces cannot be covered by firstorder waves and instead emerge from the main wave surface elevation.The higher-order Boussinesq equation has larger errors in calculating the propagation velocity of these free waves at high frequencies due to its large dispersion error.This error is precisely why there is a large deviation between the results of the higher-order Boussinesq equation and the experimental results, namely, the large error in calculating the propagation velocity of free waves (especially for third-order (kh = 4.23) and fourthorder (kh = 7.62) harmonics).This error causes a large phase difference between the free waves and actual free waves in experiments, resulting in a large deviation in the overall shape of the time series of the calculated wave surfaces from the experimental results, as shown in the results of the last four wave gauges in Figure 3.However, this dispersion is completely accurate in this model, and there is no dispersion error problem; thus, its calculated results are closer to the experimental results.
In order to quantitatively compare the numerical results with the experimental results, the error index proposed by Wilmott [30] is used for description: where ( ) x j is the experimental result, ( ) y j is the numerical result, x is the average value of the experimental result, and WIA is the error index between the numerical result and the experimental result ( 1 WIA = means that the numerical result is in complete agreement with the experimental result, and 0 WIA = means that the numerical result is in complete disagreement with the experimental result).According to the above formula, Table 2 gives the error index between the numerical result and the experimental result of the two models at different locations.From Table 2, it can be seen that the error index between the numerical result and the experimental result of this model at different locations is greater than that of the high-order Boussinesq equation, and the locations with significant differences are at x = 11.8 m (difference value: 0.03) and x = 13.5 m (difference value: 0.02), which also quantitatively illustrates that more accurate results can be obtained with this model than with the high-order Boussinesq equation in the back slope and back area of the submerged breakwater.The accuracy of the dispersion in the above calculation results of this model also depends on the accuracy of the difference scheme.Because this model adopts a fifth-order accuracy difference scheme for the time derivative and a fourth-order accuracy difference scheme for the spatial derivative, it can better ensure the realization of accurate dispersion in the calculation results.To illustrate this point, Figure 3 shows another comparison between the wave surface time series of the wave gauges at different locations using a lowerprecision, second-order accuracy difference scheme (leap-frog scheme) (using the same spatial and temporal steps as above) and the calculation results using the higher-precision difference scheme in Figure 4.As can be seen, the former difference scheme leads to significant errors in the calculation results compared to the experimental results, indicating that the calculation accuracy of this model in Figure 4 is also guaranteed by using the high-precision difference scheme mentioned above.(3) Two-dimensional shoal experiment: The above examples are one-dimensional horizontal cases.For two-dimensional horizontal cases, Whalin [31] conducted wave propagation experiments on shoals for computational simulation.In the experiment, a three-dimensional focused terrain was placed in a 25.603 m 6.1m × water tank (see Figure 5).The water depth on the deep water side of the tank was 0.457 m, and the shallow water depth at the end was 0.152 m.Three wave conditions were simulated, using incident wave heights of 0.0212 m, 0.0212 m, and 0.0298 m, with a period of 2.0 s.The width of the computational domain was consistent with that of the experiment, and the length was larger than that of the experiment, with an additional portion for a wave-absorbing sponge layer.The spatial step size was 0.078 m x y ∆ =∆ = , and the time step size was 0.03 s t ∆ = .Figure 5 shows the comparison of the calculation results of this model with and without the second-order bottom slope term with different wave heights and also shows the high-order Boussinesq equation results and the experimental results.These comparisons not only illustrate that this model has higher-accuracy calculation results than the highorder Boussinesq equation, but they also illustrate the influence of the second-order bottom slope term on wave deformation.From the results shown in the figure, it can be seen that, in the region of x = 15-22 m, the calculation results using the second-order bottom slope term are more consistent with those of the experimental results, indicating that these terms can also improve the accuracy of the calculation results for wave propagation on this shoal.From the results shown in the figure, it can also be seen that the accuracy of the amplitude of each harmonic from the first to third order in this model with the secondorder bottom slope term is better than that of the high-order Boussinesq equation, that is, the former is closer to the experimental results than the latter, especially for the thirdorder (kh = 4.32) amplitude.This indicates that for the calculation of high-frequency wave components, the accurate dispersion characteristics of this model make it superior to the high-order Boussinesq equation, which is approximate in the high-frequency range and has a certain impact on the accuracy of the results of calculating high-frequency wave components.In order to quantitatively illustrate this issue and quantitatively compare it with the results of the high-order Boussinesq equation, Table 3 presents the error index (see (64)) for each harmonic amplitude.From the table, it can be seen that after adding the second-order bottom slope term to this model, the A values of each harmonic amplitude increased, with the largest difference appearing in the third harmonic (a difference of 0.055) for H = 0.015 m.The A values of this model are all greater than those calculated via the high-order Boussinesq equation, and the largest difference appears in the second harmonic (a difference of 0.074) for H = 0.0212 m.Experimental setup (upper) and comparison of the harmonic amplitudes of the two models (below).Experiment: solid circles, first harmonic; solid triangles, second harmonic; solid squars,third harmonic.Computation: solid lines, includes second-order slope term; dotted lines, does not include second-order slope term; dash lines, N4D4 model.

Calculation Results of Bragg Reflection
The inclusion of the new second-order bottom slope term in the continuity equation allows for a consideration of the impact of rapidly changing bottom topographies, such as Bragg reflection.This phenomenon refers to the backscattering of waves from an undulating seabed due to a resonant interaction between the waves and the sea bottom, where the wavelength of the incident wave is half that of the bed's undulation wavelength (assuming an infinite number of ripples).The capability of the current model to account for this feature was confirmed through simulations replicating the Bragg reflection experiments conducted by Davies and Heathershaw [32] for single sinusoidal ripples and Guazzelli et al. [33] for double sinusoidal ripples.The expression for the corresponding water depth remains unchanged: (where T is the incident wave period).Six experimental cases were selected for the simulations, as outlined in Table 4, with reference to the parameters detailed in (73).

Propagation of Regular Waves on Flat-Sloped Coasts
In this section, we verify the ability of the coupled models to simulate wave breaking and climbing.As the wave approaches the breaking point, the wave height gradually increases.When the critical breaking point is reached, the wave exhibits strong nonlinearity, with the wave crest becoming steeper and the wave trough becoming flatter.Therefore, by simulating the change in the elevation of the free surface, the ability of the coupled models to simulate wave breaking and climbing can be effectively verified.In this section, we describe a physical model experiment based on Ting and Kirby's [34] planar beach experiment.The experiment took place in a two-dimensional wave tank with a bottom slope of 1:35 and a consistent water depth of 0.4 m in the quiet water zone.The experimental layout is shown in Figure 8, in which the origin is located at the starting point of the topographic slope and the x-axis points in the direction of the shore.The numerical simulation used the same computational domain as that of the experiment, with the following experimental wave conditions: an incident wave height of 0.127 m and a period of 2.0 s.The spatial step size selected for numerical simulation was 0.05 m x ∆ = , and the time step size was 0.01s t ∆ = .Figure 9 shows the free surface elevation time series at wave height gauge points 1, 4, 5, and 7.The selected measuring points are located at x = −1.25 m in front of the slope, at x = 3.5 m from the non-breaking point, at x = 6.4 m from the breaking point, and at x = 9.1 m in the wave breaking zone.From the results shown in the figure, it can be seen that the waveforms of the two models are similar at the four wave height gauge points.At wave height gauge points 1 and 4, the numerical results of the coupled model at the wave trough position are better than those of the N4D4 model, which is in good agreement with the experimental results.At wave height gauge point 5, there is a slight difference between the numerical results of the two models and the experimental results at the wave trough position; however, relatively speaking, the results of the coupled model are closer to the experimental results.At wave height gauge point 7 in the wave breaking zone, there is a certain deviation between the two models and the experimental results at the wave trough position.The reason for this may be that bubbles appeared during the experiment in the wave breaking zone, which cannot be simulated by numerical models due to the conservation of flow.These bubbles will be superimposed onto the free surface elevation as water bodies; thus, the numerical simulation results will show a higher wave surface elevation.The above results show that the coupled model can accurately describe the strong nonlinear characteristics of wave propagation and breakup on slopes.Although there is a slight difference in the wave surface at the trough, at the wave crest, where the nonlinearity is strong, the coupled model can well simulate its changes.Due to the precise dispersion characteristics of the coupled model, it can better describe the wave propagation process during and after wave breaking compared to the higher-order Boussinesq equation.

Numerical Simulation of Gently Sloping Coastal Rip Flow
The results of the experiment were taken from Peng and Zou [35] and were collected at the wave basin of the State Key Laboratory of Coastal and Offshore Engineering at Dalian University of Technology.The harbor basin is 55 m long, 34 m wide, and 0.7 m deep.In order to generate oblique incident waves and increase the length of the coastal model, we arranged the coastal model at a 30-degree angle with the wave making plate.The slope of the coast was 1:40, and the flat bottom still water depth in front of the slope was 0.45 m.The slope section above the still water line retained a distance of 3 m to ensure the wave climbing length.The spur dike was arranged on the edge of the upstream coastal model and perpendicular to the shoreline, extending to the top of the slope.Because the measurement of the wave field and velocity field in the experiment was based on the spur dike as a reference, the origin of the coordinate system was the intersection point of the spur dike and the still water line.The x direction was perpendicular to the still water line and pointed in the direction away from the shore, while the y direction was perpendicular to the spur dike and pointed in the upstream direction.The sandbar on the coast was a Gaussian sandbar with a width of 2 m.The sandbar located 1-2 m away from the spur dike was removed to form a split-flow trench.The wave field in the experiment was formed by superimposing oblique incident waves and reflected waves after reflection from the spur dike.The superimposed wave field consisted of a standing wave system composed of many rhombuses, with the vertex and midpoint of each rhombus being the nodes of the standing waves along the coast.The wave field simulated in this study had a period of 1.5 s; thus, the standing wave wavelength was 2.71 m, and the positions of the nodes were 1.36 m and 4.07 m away from the spur dike.In the numerical simulation, the positions of the nodes were 0 m (0.1 m in the numerical simulation) and 2.71 m away from the spur dike.The wave field was measured using a wave meter, located at both nodes, while the split-flow field was measured using acoustic Doppler velocity measurement (ADV) and float video tracking methods, as shown in Figure 10.
In the numerical calculation, the same calculation domain setting as that of the experiment was used, and the total reflection boundary was used at one side of the groin.The simulated wave case was a regular wave with a period of 1.5 s and an incident wave height of 4.1 cm.The coupling surface was 16 m away from the wave making position.The space step was 0.05 m x y ∆ =∆ = , and the time step was 0.01 s t ∆ = .The simulation duration was 720 s, and the 120-720 s time series data were selected to count the wave height, water increase and decrease, and time-averaged flow field.× ).In order to compare with the experimental results, the numerical results also use the same grid resolution.Figure 11 was generated by MATLAB software (MathWorks USA, Inc. version: MATLAB 9.14).The vorticity results in the figure reveal the presence of two prominent pairs of large-scale vortices.The first pair emerges at the channel's edges, located precisely at the sides of the initial node line, where it aligns with the channel's central line.Similarly, the second pair forms adjacent to the second node line.These distinct vortex pairs are attributed to the swift decrease in the wave breaking height from the channel's edges towards the node line or the channel's centerline, which is caused by longshore standing waves or reduced water depths within the channel.The formation of large-scale vortices around these reduction points, as indicated by Peregrine [36], is effectively captured by both the experimental observations and numerical simulations in the figure.Further validation is presented in Figure 12, where a comparison of the experimental and numerical cross-shore and alongshore velocities along the bar crest showcases the model's accuracy in simulating vortex motion induced by wave breaking in the surf zone.

Conclusions
In this study, we developed a set of wave equations with accurate dispersion in terms of free surface velocity and wave surface elevation.The nonlinearity of the equations is accurate to the fourth order, and the dispersion is completely accurate.The continuity equations of the model maintain a spin in the vertical vortex field during derivation, and the momentum equations were derived directly from the Euler equations without omission, so that the model can simulate the nearshore circulation generated by wave breaking, such as eddy motion in the rift current.To overcome the inability of the equation to simulate wave creep in the nearshore region, the equation was expanded in the nearshore region, and the dispersion of the equation was preserved to order, thereby converting the integral calculations contained in the equation into differential calculations.
The accurate dispersion equation was coupled to the expanded equation, and a coupled model was developed that applies the accurate dispersion model in deep water and the nearshore model in shallow water.The coupled model has several advantages over other models: (1) The present equation can simulate wave-current interaction phenomena, such as coastal currents and nearshore circulation, compared with velocity potential equation.( 2) Compared with the Boussinesq equation, the accurate dispersion model applied to deep-water regions is derived without any assumptions about the water depth and without any coefficients related to the wave number; thus, it can be applied to arbitrary water depths.(3) The present equation is a nonlinear model compared to the gentle slope equation, and it can therefore simulate the nonlinear effects of waves.
The present model was validated against the relevant physical model experiments.The numerical results were compared with the experimental results of the submerged dike, and the results show that the accurate dispersion model has high accuracy in terms of nonlinearity and dispersion and is able to describe the wave propagation deformation well.A comparison of the numerical results with the amplitude results for shallow beach terrain further demonstrates that the accurate dispersion model can be applied to wave propagation in complex terrain.The numerical results were compared with the experimental results for flat beaches, where the coupled model described the free surface elevation well in the breaking zone, indicating that the coupled model also has highly accurate nonlinearity.The numerical results were compared with the experimental results of the rip current flow, and the results show that the coupled model can effectively simulate the vortices generated by the rip current flow field.
In future studies of this model, its computational accuracy will be improved using more accurate numerical methods, and its reliability can be verified using more experimental results.Sensitivity analyses of each of the model's parameters will be carried out to determine the parameters and their effects on the simulation results.In the future, the model will also be used to study nearshore circulation, pollutant transport, and wavestructure interaction and will be applied to different physical environments.

w
are also expressed in the wave-number space:

Figure 1 . 2  u and 2 η 1  u and 1 η
Figure 1.Layout plan of the coupling model.In the coupled model, the wave velocity and free surface elevation at the free surface are used as the variables; thus, the matching condition for the coupling region maintains the continuity of the wave velocity and free surface elevation.The variables calculated in region 1 are 1  u and 1 η , and the variables calculated in region 3 are 2  u and 2 η .The pro- cess of information transfer depends on the requirements of the integral and differential calculations.The variables in regions 1 and 3 are transferred to region 2. In order to meet the requirements of the differential calculations in region 3, the 1  u and 1 η variables, cal- culated via the accurate dispersion equation in region 1, are transferred from interval L2 to region 3, as follows:

1 E , * 1 F , and * 1 G
are the remaining terms in the two equations except for the lin- ear time term, and ( , ) bf bf F G , ( , ) wb wb F G , and ( , ) bs bs F G are the internal wave source term, bottom friction term, wave breaking term, and mixed subgrid effect term, respectively.

0 H
is the incident wave height, s x is the location of the wave source, and s β is the wave width parameter, µ , referred to as the N4D4 model.

Figure 2 .
Figure 2. Comparison of wave surface elevations of two computation models with theoretical solution.Solid lines, theoretical solution; dash-dotted lines, this model; dash lines, N4D4 model.

Figure 3 .
Figure 3. Experiment setup (upper) and comparison of wave surface elevations of two models (below).Solid circles, measurements; solid lines, this model; dash lines, N4D4 model. 0

Figure 5 .
Figure5.Experimental setup (upper) and comparison of the harmonic amplitudes of the two models (below).Experiment: solid circles, first harmonic; solid triangles, second harmonic; solid squars,third harmonic.Computation: solid lines, includes second-order slope term; dotted lines, does not include second-order slope term; dash lines, N4D4 model.

73 )
For singly sinusoidal ripples, N = 1, and for doubly sinusoidal ripples, N = 2. 0 h represents the water depth in the region of constant depth, while b A stands for the ripple amplitude.w L and bj L (j = 1, 2) correspond to the incident wavelength and the wave- length of the jth sinusoidal ripples, respectively.1 b n refers to the number of initial sinusoidal ripples, and sx represents the x coordinate of the starting point of the ripple patch.The numerical simulations take place in the computational domain shown in Figure6, using a space step denoted by /

Figure 6 .
Figure 6.Computation domain for Bragg reflection tests.

Figure 7
Figure7shows a comparison of the reflection coefficients between the numerical and experimental results for singly and doubly sinusoidal ripples.As can be seen from the figure, for the single sandbar case, when 1 2 / 1 b k k = , the model's calculation results and the experimental reflection coefficients both reached their maximums.For the double sandbar case, the model's calculation results and the experimental reflection coefficients both reached their maximums at 1 2 / 1 b k k = and

Figure 7 .
Figure 7. Reflection coefficients ( R C ) for single and double sinusoidal ripples.Solid circles, Experi- mental results; solid lines, numerical results.

Figure 9 .
Figure 9.Time series of free surface elevation at different locations.Circles, measured results; solid lines, the current model's results; dash lines, N4D4 model results.

Figure 10 .
Figure 10.Plan view of experimental layout for the planar beach.

Figure 11
Figure 11 illustrates the spatial arrangement of both the experimental and numerical average velocity vectors over time, along with the respective vorticity fields.The red arrows shown in the figure are the flow velocities measured by the float in the experiment, the black arrows are the flow velocities calculated by numerical calculation, and the vorticity field is calculated from the flow velocities using the expression U Ω = ∇ × .The resolution of the experimental results in the figure is ( 28 9× ).In order to compare with the experimental results, the numerical results also use the same grid resolution.Figure11was generated by MATLAB software (MathWorks USA, Inc. version: MATLAB 9.14).The vorticity results in the figure reveal the presence of two prominent pairs of large-scale vortices.The first pair emerges at the channel's edges, located precisely at the sides of the initial node line, where it aligns with the channel's central line.Similarly, the second pair forms adjacent to the second node line.These distinct vortex pairs are attributed to the swift decrease in the wave breaking height from the channel's edges towards the node line or the channel's centerline, which is caused by longshore standing waves or reduced water depths within the channel.The formation of large-scale vortices around these reduction points, as indicated by Peregrine[36], is effectively captured by both the experimental observations and numerical simulations in the figure.Further validation is presented in Figure12, where a comparison of the experimental and numerical cross-shore and alongshore velocities along the bar crest showcases the model's accuracy in simulating vortex motion induced by wave breaking in the surf zone.

Figure 12 .
Figure 12.Cross-shore (a) and longshore velocities (b) at sand bar centerline locations.Solid circles, experimental results; solid lines, numerical results.

Table 1 .
Theoretical wavelengths and calculated wavelengths.

Table 2 .
Comparison of WIA values between this model and N4D4 model for submerged breakwater terrain.
Comparison of wave surface elevations of two difference schemes.Solid circles, Measurements; solid lines, ABM scheme of fifth-order prediction and sixth-order correction; dash lines, leapfrog scheme.

Table 3 .
Comparison of WIA values between this model and the N4D4 model for shallow terrain.

Table 4 .
Experimental cases of single and double sinusoidal ripples.