Alternative Forms of Enhanced Boussinesq Equations with Improved Nonlinearity

We propose alternative forms of the Boussinesq equations which extend the equations of Madsen and Schäffer by introducing extra nonlinear terms during enhancement.Theoretical analysis shows that nonlinear characteristics are considerably improved. A numerical implementation of one-dimensional equations is described. Three tests involving strongly nonlinear evolution, namely, regular waves propagating over an elevated bar feature in a tank with an otherwise constant depth, wave group transformation over constant water depth, and nonlinear shoaling of unsteady waves over a sloping beach, are simulated by the model. The model is found to be effective.


Introduction
In the past 3 decades, great strides have brought Boussinesqtype equations into the family of operational coastal wave prediction models because they possess both dispersion and nonlinear properties and require relatively little computation effort.The development of variants of the theory was trigged by improving linear and nonlinear properties of the model.Extensive reviews have been provided in [1,2].
Earlier works have concentrated on improving linear dispersion of classical Boussinesq equations (which are accurate only for weakly dispersive and mildly nonlinear water waves), thus allowing the model to treat a wider range of water depths.Extensive research results (e.g., [3][4][5][6][7][8]) have been published, and applicability range of the equations has been greatly improved.In [9], nonlinear shallow water equations are extended to include dispersion using an original approach based on hyperbolic approximation.The resulting equations have attractive hyperbolic characteristics and can be solved efficiently using the finite volume method [10].Among these studies, one of the notable procedures was provided by Madsen et al. [3], who established a procedure for optimizing model performance through rearrangement of dispersive terms.This procedure has been extensively utilized to improve dispersion of Boussinesq equations, for example, [6,11].Other linear properties, such as shoaling property and Bragg reflection, cannot be ignored considering the applicable range of water depth to be enlarged for a large number of new forms of Boussinesq-type equations.In numerous improved models, attention has also been given to the aforementioned properties [4,6,11,12].
Subsequent to the initial work on improved linear properties, the next topic to draw attention is the problem of improving model's nonlinear accuracy.One of the critical steps in this process is to develop equations with fully nonlinear characteristics [9,[13][14][15].Wei et al. first presented such equations with second-order accuracy [13].The fourth-order Boussinesq-type equations with fully nonlinear characteristic were lately presented [2,14,15].All the aforementioned studies show that fully nonlinear models can present better numerical results than weak nonlinear models.However, a huge inaccuracy in second-and thirdorder transfer functions and amplitude dispersion still exists for deep water.This limitation has inspired further effort to develop other formulas, such as that of Agnon et al. [16] and Madsen et al. [17].These models exhibit excellent linear and nonlinear performances even for extremely deep water.However, fourth-order nonlinear equations, as well as highly dispersive and fully nonlinear models, usually involve complicated formulas, thus creating a challenge for numerical implementation.Hence, developing a model with improved nonlinear performance and a relatively simple formula is necessary.
Efforts have been made to address this issue within second-order Boussinesq equations.Zou [7] added a constant coefficient in a second nonlinear term in the derived Boussinesq equation to improve the accuracy of the super harmonic transfer function.A similar technique was adopted in [18]; however, a variable coefficient was used instead of a constant coefficient.Yao et al. [19] further extended the aforementioned method by developing a parameter-adaptive version of Zou's equations, wherein the introduced parameter is a function of space and time which locally optimizes second-order solutions throughout the simulation domain.In [20], the reference elevation concept in the equation of Wei et al. [13] (a fully nonlinear version of Nwogu's equations [5]) is generalized by considering the contribution of instantaneous-free surface elevation, and thus, considerably improving theoretical and practical nonlinear performances.These previous studies are very constructive; however, they provide no explanation on the accuracy of adding terms theoretically.In this study, we attempt to improve the nonlinear properties of a set of second-order Boussinesq equations with full nonlinearity in alternative way.
Alternative forms of Boussinesq equations with secondorder full nonlinearity are presented.The model is based on the equations of Madsen and Schäffer [6].Nonlinear performance is improved by introducing extra nonlinear terms during the enhancement.The improved model is presented in Section 2, including its governing equations, theoretical analysis, and numerical implementation.In Section 3, three numerical simulations for nonlinear dispersive waves are conducted.Finally, conclusions are drawn in Section 4.

Governing Equations Fully Nonlinear to O(𝜇 2
).A set of governing equations with second-order full nonlinearity was presented in [6].In nondimensional form, these equations can be expressed as where with In the previous equations, h is still water depth,  is surface elevation, and u is depth-averaged velocity.∇ = (/, /) is the two-dimensional gradient operator. =  0 /ℎ 0 and  = ℎ 0 / 0 are parameters that denote nonlinearity and dispersion of equations,  0 , ℎ 0 , and  0 are the characteristic wave amplitude, water depth, and wave length, respectively.

Technique of Madsen and Schäffer for Improving Linear
Dispersion and Shoaling.The dispersion contained in (1)-( 2) is the same as that in classical Boussinesq equations, and thus these are only applicable in shallow water region.To improve dispersion characteristic, an enhancement technique is introduced and applied to (1)-( 2) in [6].The rationale of the enhancement technique is summarized as the following three steps.The first step in the procedure is to apply the operator ∇∇⋅ to (2) and multiply the result by  1 ℎ 2 , which yields The second step is to multiply (2) by h and apply the operator ∇∇⋅ and multiply the result by  2 ℎ, which yields Equations ( 5) and ( 6) are correct up to the same order as the original equation (2), which consequently can be consistently modified by the use of these equations.Finally, by adding ( 5) and ( 6) to (2), the enhanced versions of the governing equations are obtained [6].The two coefficients  1 and  2 are yet to be determined.All the terms introduced in (5)-( 6) have similar expression to those in Λ 20 and Λ 21 in original equation (2).As such, the aforementioned procedure actually amounts to rearrange second-order (O( 2 )) terms, as reviewed in [17][18][19][20].This finding implies that the manipulation is not unique, and that we can rearrange terms in an alternative manner, as shown in the next section.

Enhanced Equations with Improved Nonlinearity.
Only linear properties, that is, dispersion and shoaling, are considered for optimization in the previously discussed method.In this section, we focus instead on improving model's (( 1)-( 6)) nonlinearity using a similar technique.This procedure is akin to the manipulation approach described in Section 2.2, but in a nonlinear sense.To carry out this technique, a reference variable ℎ + ℎ ( is free parameter) is used instead of h in the procedure for obtaining ( 5) and ( 6); thus, we get When ( 7) are inserted into the left of ( 2), the final governing momentum equations become with Λ 20 , Λ 21 , Λ 22 , and Λ 23 are defined in (3).Hereafter, (1) and (2) (with ( 5) and ( 6) added to the left side) are referred to as the original equations, whereas (1) and ( 8) are called enhanced equations.Note that both original and enhanced equations remain as second-order formulas (O( 2 )) with full nonlinearity; however, their nonlinear characteristics will be different, as will be shown in Section 2.4.
Compared with ( 5)-( 6), (7) introduce extra nonlinear terms at orders  2 ,  2  2 ,  3  2 , which contribute to nonlinear performance of the model.Equations (7) are still correct up to the same order as the original equation ( 2), hence, the rationale behind the method retained.

Theoretical Analysis.
In this section, we conduct Fourier analysis of enhanced equations ( 1) and ( 8) on a horizontal bottom.The coefficients  1 ,  2 , and  will be determined thereby.The following Stokes-type expansions [6] are introduced for this purpose: where  =  − .After substituting these expansions into the one-dimensional (1-D) version of the equations and collecting terms of order (1), O(), O( 2 ), governing equations for the first-, second-, and third-order problems can be obtained (see Appendix for details).
Solving the first-order problem yields the dispersion relation of the model ((A.4) in the Appendix), and the value of  1 +  2 can be determined by matching C to the exact solution  Stokes = ℎ tanh()/.Setting  1 +  2 = −1/15 leads to a Padé [2,2] approximation of the exact dispersion, thus creating a match with Stokes solution up to approximately  = ℎ =  [6].As stated in [21], the parameter can be alternatively optimized to minimize phase and group velocity errors over the entire dispersive range of 0 <  < .This process can be achieved by introducing the following form for the joint error (joint root mean square error) [21]: where   is group velocity and superscript "Stokes" represents analytical solutions from Stokes theory.The optimization process aims to find the minimum value of error function (11).Finally,  1 +  2 = −0.0533 is found to be the optimum value with a minimum error of 1.66% over 0 <  < .
Phase celerity and group velocity of the model with this value are shown in Figure 1.The agreement between the modeled dispersion relation and the theoretical one is relatively good over the range 0 <  < .Larger discrepancies are found for the group velocity; however, maximum error remains below 10% for all the physical wavelengths that the model is able to describe.Nevertheless, the error in group velocities starts to increase rapidly when  > 2.5, thus showing that this property will be overestimated in this range.Corresponding results from Padé [2,2] expansion are also shown in Figure 1 for comparison.Overall accuracy of phase celerity and group velocity is improved over 0 <  <  by using (11) as the control criterion.Specific values of  1 and  2 cannot be determined yet by this process, which should be completed by analyzing shoaling property of the model.Following the same procedure in [6],  2 can be determined by finding the minimum  ) 2  (here,   is the shoaling coefficient, and the superscript "Stokes" represents the analytical solution), and  2 = −0.0645 is found to give optimum shoaling performance.Figure 2 presents the comparisons between the shoaling coefficient for the present model and the analytical reference.They are generally in good agreements over the range 0 <  < .
Nonlinear property is analyzed by solving second-and third-order problems.Parameter  is found to contribute to second-( 2 ), third-order ( 33 ) transfer functions, as well as amplitude dispersion  3 .In this process, the joint error (defined by ( 11)) of the second-order transfer function ( 2 ) and amplitude dispersion ( 3 ) is chosen as the criterion to determine . 3 is chosen as the criterion instead of  33 because the amplitude dispersion is more important than  33 for nonlinear wave evolution. 3 is the factor that controls the nonlinear process of amplitude dispersion, whereas  33 is a smaller nonlinear quantity compared with  2 , as argued in [22]. = 2.521 is determined as the final value.Corresponding results with this value are given in Figure 3.This figure shows that the accuracy of the transfer functions and amplitude dispersion are significantly improved, which is apparently caused by retaining extra nonlinear terms during enhancement.
2.5.Numerical Implementation.The numerical scheme mainly follows procedures presented in FUNWAVE [23]; however, several modifications are made.The main procedures are summarized as follows.
Governing equations are discretized using a staggered grid (by contrast, a uniform grid is adopted in [23]) because recent studies suggest that discretizing Boussinesq equations on a staggered grid can achieve better stability [24,25] than those on a uniform grid.Derivatives are approximated by a higher-order finite difference formula, and a composite fourth-order predictor-corrector Adams-Bashforth-Moulton integration scheme is adopted in the model to perform time marching.The independent variable wave surface  can be directly determined by solving the continuity equation, whereas the independent variable  is obtained by solving a tridiagonal linear system from the momentum equation.
The entire computation domain is enclosed by impermeable walls, wherein the horizontal velocity is set to zero.Sponger layers are placed in front of the solid walls to absorb wave energy propagating out of the computational domain.The incident waves are internally generated in the computational domain by adding a source function into the mass equation.The implicit filter [21] is used sparsely to remove higher frequency oscillations.The details of numerical implementations are referred to [23].

Numerical Results and Discussions
To validate the proposed model, three tests involving strongly nonlinear evolution of waves are chosen for the simulation, including (a) 1-D regular wave trains propagating over trapezoidal bar feature in an otherwise constant depth tank, (b) nonlinear evolution of 1-D wave groups over constantly deep water, and (c) nonlinear shoaling of unsteady waves up to the breaking point over a mildly sloping beach.Both enhanced and original models will be used in the simulations, and the computed results will be compared.Available experimental data or analytical solutions will be adopted as reference.

Propagation of 1-D Regular Wave Trains over Submerged Bars.
Wave propagation over a submerged bar is an extremely complicated process.Nonlinear interactions transfer energy from the leading wave component to higher harmonics as waves propagate onto the front slope of the bar and begin to become steep.At the leeside of the bar, water becomes deep and nonlinear coupling of higher harmonics with the fundamental wave becomes progressively weaker.Therefore, each of the Fourier components is released as free waves with their own bound higher harmonics.Each component travels with a different speed, thus resulting in a fairly complicated process.Numerical prediction of this process requires higher-order linear and nonlinear accuracies of the model.Such prediction has been widely used as a benchmark test for validating Boussinesq-type models [14,15].
Here, Case (a) of Luth et al. experiments [26] is simulated, and the corresponding experimental setup is shown in Figure 4. Regular wave trains with period  = 2.02 s and height  = 0.02 m (ℎ = 0.67) are internally generated in the computation domain.Time and spatial size for the simulation are 0.01 s and 0.02 m, respectively.Numerical results from the original and enhanced models are shown in Figure 5 and compared with the experimental data.At the first gauge ( = 10.5 m), both models exhibit similar performances and present good agreements with the experimental data.Differences begin to appear on top of the bar ( = 13.5) and downstream of the bar slope ( = 15.7 m).At these locations, the enhanced model obtains slightly better numerical results than the original model.A large discrepancy between numerical results and measurements is expected at the rear part of the wave tank ( = 19.0m), where higher harmonics already exceeds the application range of the model.Such result is also found and argued in the previous literature for second-order Boussinesq-type models [14,15].
The enhancement only improves the nonlinearity of the model, whereas dispersion and shoaling remain unchanged.As such, checking the spatial distribution of higher harmonics for this case is more reasonable.The corresponding numerical results for the two models are plotted in Figure 6 and compared with the experimental data.For the first and second harmonics, numerical results from the two models exhibit negligible difference and agree well with the measurements.This outcome is expected because the second transfer functions  2 for the two models are primarily the same (Figure 3(a); ℎ = 0.67).A discrepancy is expected for third and fourth harmonics according to the analysis shown in Figure 3(b) (ℎ = 0.67).Figure 6 shows that the enhanced model indeed presents better results than the original model, and the numerical results are in good agreement with the measurements.

1-D Nonlinear Evolution of Deep Waves over a Flat Bottom.
The experiment for nonlinear evolution of wave groups over a constant water depth [27] provides a good basis for checking third-order nonlinearity of Boussinesq models.According to Stokes wave theory, third-order nonlinearity enables the amplitude to contribute to the dispersion relation, and thus, large waves in an envelope travel faster than those with lower amplitudes.As a result, wave envelopes become skew; that is, the front part of a wave group is pressed and the back part is stretched.New free-wave components will also be generated during the process because of nonlinear resonant interaction, and the wave spectrum is expected to become wider during wave group evolution.Simulations for these nonlinear properties require high accuracy in nonlinearity and dispersion of a mathematical wave model [15,22].
In the experiment, a wave maker is driven by the following signal in a flume with a constant water depth of 0.6 m: where  0 is the forcing amplitude,  0 is the carrier frequency, and Ω is the modulation frequency.The spectrum of the signal achieves its maximum at the carrier frequency  0 , plus several smaller peaks at discrete frequencies spaced by  0 /10 with three dominant modes.The period of the carrier wave is 0.9 s and the corresponding value of kh is 3.0, thus denoting that wave motions are within the deep water regime.In the corresponding simulations, we set the spatial and temporal steps to 0.02 m and 0.01 s, respectively, and collected data until a steady state is reached.Numerical results from the enhanced and original models are presented in Figure 7 and compared with the measurements.It can be seen that close to the wave maker ( = 0.24 m), the two numerical results are almost identical and agree with the measurements.As the wave propagates along the flume, the results of the two models diverge.The original model obtains a big phase error as the waves propagate along the flume, whereas the enhanced model predicts the correct phase, thus denoting that the accuracy of amplitude dispersion has been improved greatly.However, the enhanced model still fails to reproduce the details of the wave shape.This finding may be ascribed to the fact that the motion of carrier waves and higher harmonics almost exceeds the application range of the equations.
Figure 8 shows Fourier amplitude spectra of the surface elevation of the wave group at  = 0.24 m and 9.47 m.In the experimental spectra, amplitudes increase on the highfrequency side when waves reach  = 9.47 m and the spectrum becomes asymmetric.Computed results from the enhanced and original models are plotted in the same figure.The enhanced model produces amplitude spectra similar to those of the experiments, whereas the original model produces spectra which change insignificantly by the time they reach the end of the flume.This finding shows that wave components produced by the model remain unchanged along the flume.
The asymmetric growth of the spectrum shown in the figures presents another nonlinear characteristic of nonlinear wave group evolution apart from the amplitude dispersion, that is, the generation of new free-wave components due to resonant wave-wave interaction.This is a third-order (O( 3 )) nonlinear effect which causes the energy transfer from one frequency to another in a wave spectrum and is referred to as four-wave resonant interaction [28,29].For the present case, although there are only three free-wave modes for the signal (12) near the wave maker ( = 0.24), the four-wave resonant interaction will inspire other new free-wave modes to develop as waves propagate along the tank.

Nonlinear Shoaling of Unsteady Waves.
The unsteady shoaling test is chosen to investigate the effect of the improvements made to the original equations.In this procedure,   the most nonlinear case from [20], wherein unsteady waves are specified as the initial surface elevation and the subsequent shoaling over a given topography, is simulated.By normalizing all quantities with respect to the greatest depth, ℎ 1 , the topography is specified as where ℎ min /ℎ 1 = 0.2 and  = 50ℎ 1 .Initial surface elevation is given as where   /ℎ 1 is the initial amplitude, and the number of waves   = 10.In the simulation, ℎ 1 is set as 1.0 m and Δ = 0.05 m, while Δ = 0.01 s. Figure 9 shows initial and final surface and the topography.Using the initial value problem also helps eliminate comparison difficulties that can arise from differing wave generation schemes [20], as also found in our numerical experiments.
Both the original and enhanced equations are used for simulation, and the numerical results are compared with the numerical solution from the potential flow results [30] (which could be regarded as the analytical solution to   the problem).We discover that in deeper water (where waves are weakly nonlinear) and in shallow water, Boussinesq envelopes are significantly similar to each other, and to the potential flow results.However, in the highly nonlinear shoaling zone near the front of the shelf, results differ significantly, thus providing a good comparison of nonlinear shoaling properties.Figure 10 shows the computed crest and trough envelopes compared to the potential flow results for 40 < /ℎ 1 < 60.The original equations underpredict crest elevations significantly when wave heights are high.By contrast, the enhanced equations obtain better numerical results.The same trend is also exhibited for trough elevations; however, the difference between the original and enhanced equations is not very obvious.

Conclusions
An alternative form of the Boussinesq-type model, which extends the equations in [6], is developed.By introducing the contribution of time-varying component during the enhancement, we develop a possible method for improving nonlinear performance of the model.Theoretical analysis shows that nonlinear performance of the enhanced equations is considerably improved.To demonstrate the effect of the improvement, three tests involving strong nonlinear characteristics, that is, nonbreaking regular wave propagation over a submerged bar, nonlinear evolution of deep water waves over a flat bottom, and nonlinear shoaling of unsteady waves over a sloping beach, are simulated.Numerical results show that the enhanced model demonstrates considerable improvements over the original model, at least with regard to the tests considered in this study, thus demonstrating the potential of the proposed model in describing nonlinear processes up to the breaking point.

Figure 2 :
Figure 2: Comparisons of shoaling gradient between the equations (solid line) and the analytical solution (dashed line).

Figure 3 :
Figure 3: Comparisons of second-order (a), third-order (b) transfer functions, and amplitude dispersion (c) between the enhanced equations (solid line) and the original equations (dashed line).

Figure 5 :
Figure5: Wave surface elevations of Luth's experiment (circle)[26] and the numerical results of the enhanced model (solid) and the original model (dashed).

Figure 6 :
Figure6: The spatial variation of the harmonics of Luth's experiment (circle)[26] and the numerical results of the enhanced model (solid) and the original model (dashed).

Figure 7 :
Figure 7: Experimental data (circle) of wave elevations of wave groups and the numerical results of the enhanced model (solid) and the original model (dashed).

Figure 8 :
Figure 8: Comparisons of amplitude spectra at  = 0.24 m and  = 9.47 m between the experimental data (a) and the numerical results from the enhanced model (b) and the original model (c).

Figure 9 :
Figure 9: Initial (a) and final (b) surface elevations, with crest and trough envelopes and underlying bathymetry.

Figure 10 :
Figure 10: Comparisons of crest and trough envelopes between the enhanced model, the original model, and the potential flow results.