Alternative Forms of Enhanced Boussinesq Equations with Improved Nonlinearity

Abstract

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.

1. Introduction

In the past 3 decades, great strides have brought Boussinesq-type 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–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–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 third-order 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 second-order 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.

2. Mathematical Model

2.1. Governing Equations Fully Nonlinear to 𝑂(𝜇²)

A set of governing equations with second-order full nonlinearity was presented in [6]. In nondimensional form, these equations can be expressed as 𝜂_𝑡𝐮𝐮+∇⋅ℎ+𝜀𝜂=0,_𝑡1+∇𝜂+2𝜀∇𝐮²+𝜇²Λ₂₀+𝜀Λ₂₁+𝜀²Λ₂₂+𝜀³Λ₂₃=𝑂𝜇⁴, (1)(2) where Λ₂₀=ℎΓ_𝑡,Λ₂₁=−𝜂Γ_𝑡()Γ()+∇⋅ℎ𝐮+∇𝐮⋅ℎΓ−𝜂∇⋅ℎ𝐮_𝑡+12(∇⋅(ℎ𝐮))²,Λ₂₂=16𝜂²∇∇⋅𝐮_𝑡−13()∇()Γ1𝜂∇ℎ𝐮∇⋅𝐮+∇⋅𝜂𝐮+∇−𝜂𝐮⋅Γ−2𝜂²∇⋅𝐮_𝑡()()((,Λ+𝜂∇⋅𝐮∇⋅ℎ𝐮−𝜂𝐮⋅∇∇⋅ℎ𝐮))₂₃1=−3∇()−1𝜂∇⋅𝜂𝐮∇⋅𝐮+∇⋅3𝜂²()+1𝐮⋅∇∇⋅𝐮2𝜂²(∇⋅𝐮)², (3) with 1Γ=6()−1ℎ∇∇⋅𝐮2∇(∇⋅(ℎ𝐮)). (4) In the previous equations, ℎ is still water depth, 𝜂 is surface elevation, and 𝐮 is depth-averaged velocity. ∇=(𝜕/𝜕𝑥,𝜕/𝜕𝑦) is the two-dimensional gradient operator. 𝜀=𝐴₀/ℎ₀ and 𝜇=ℎ₀/𝐿₀ are parameters that denote nonlinearity and dispersion of equations, 𝐴₀, ℎ₀, and 𝐿₀ are the characteristic wave amplitude, water depth, and wave length, respectively.

2.2. 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 𝛽₁ℎ², which yields 𝛽₁ℎ²𝜇²∇∇⋅𝐮_𝑡1+∇𝜂+2𝜀∇𝐮²=𝑂𝜇⁴,𝜀𝜇⁴. (5) The second step is to multiply (2) by ℎ and apply the operator ∇∇⋅ and multiply the result by 𝛽₂ℎ, which yields 𝛽₂ℎ𝜇²∇∇⋅ℎ𝐮_𝑡1+∇𝜂+2𝜀∇𝐮²=𝑂𝜇⁴,𝜀𝜇⁴. (6) 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 𝛽₁ and 𝛽₂ are yet to be determined.

All the terms introduced in (5)-(6) have similar expression to those in Λ₂₀ and Λ₂₁ in original equation (2). As such, the aforementioned procedure actually amounts to rearrange second-order (𝑂(𝜇²)) terms, as reviewed in [17–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.

2.3. 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 ℎ in the procedure for obtaining (5) and (6); thus, we get 𝛽₁ℎ+𝜀𝛾𝜂²𝜇²∇∇⋅𝐮_𝑡1+∇𝜂+2𝜀∇𝐮²=𝑂𝜇⁴,𝜀𝜇⁴,𝛽₂ℎ+𝜀𝛾𝜂𝜇²∇∇⋅ℎ+𝜀𝛾𝜂𝐮_𝑡1+∇𝜂+2𝜀∇𝐮²=𝑂𝜇⁴,𝜀𝜇⁴. (7) When (7) are inserted into the left of (2), the final governing momentum equations become 𝐮_𝑡1+∇𝜂+2𝜀∇𝐮²+𝜇²𝐆₂₀+𝜀𝐆₂₁+𝜀²𝐆₂₂+𝜀³𝐆₂₃=𝑂𝜇⁴, (8) with 𝐆₂₀=Λ₂₀+𝛽₁ℎ²∇∇⋅𝐮_𝑡+∇𝜂+𝛽₂ℎ∇∇⋅ℎ𝐮_𝑡,𝐆+ℎ∇𝜂₂₁=Λ₂₁+2𝛽₁𝛾ℎ𝜂∇∇⋅𝐮_𝑡+1+∇𝜂2𝛽₁ℎ²∇∇⋅∇𝐮²+𝛽₂1ℎ∇∇⋅2ℎ∇𝐮²+𝛾𝜂𝐮_𝑡+∇𝜂+𝛽₂𝛾𝜂∇∇⋅ℎ𝐮_𝑡,𝐆+∇𝜂₂₂=Λ₂₂+𝛽₁𝛾²𝜂²∇∇⋅𝐮_𝑡+∇𝜂+𝛽₁𝛾ℎ𝜂∇∇⋅∇𝐮²+12𝛽₂ℎ∇∇⋅𝛾𝜂∇𝐮²+𝛽₂1𝛾𝜂∇∇⋅2ℎ∇𝐮²+𝛾𝜂𝐮_𝑡,𝐆+∇𝜂₂₃=Λ₂₃+12𝛽₁𝛾²𝜂²∇∇⋅∇𝐮²+12𝛽₂𝛾²𝜂∇∇⋅𝜂∇𝐮² (9)Λ₂₀, Λ₂₁, Λ₂₂, and Λ₂₃ 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 (𝑂(𝜇²)) 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 𝜀𝜇², 𝜀²𝜇², 𝜀³𝜇², 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.

2.4. Theoretical Analysis

In this section, we conduct Fourier analysis of enhanced equations (1) and (8) on a horizontal bottom. The coefficients 𝛽₁, 𝛽₂, and 𝛾 will be determined thereby. The following Stokes-type expansions [6] are introduced for this purpose: 𝜂=𝑎₁cos𝜃+𝜀𝑎₂cos2𝜃+𝜀²𝑎₃₁cos𝜃+𝑎₃₃cos3𝜃,𝑢=𝑢₁cos𝜃+𝜀𝑢₂cos2𝜃+𝜀²𝑢₃₁cos𝜃+𝑢₃₃cos3𝜃,𝜔=𝜔₀+𝜀²𝜔₃, (10) where 𝜃=𝑘𝑥−𝜔𝑡. After substituting these expansions into the one-dimensional (1-D) version of the equations and collecting terms of order 𝑂(1), 𝑂(𝜀), 𝑂(𝜀²), 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 𝛽₁+𝛽₂ can be determined by matching 𝐶 to the exact solution 𝐶^Stokes=ℎtanh(𝜅)/𝜅. Setting 𝛽₁+𝛽₂=−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]: Err=∫^𝜋₀𝐶−𝐶^Stokes2𝑑𝜅∫^𝜋₀𝐶^Stokes2+∫𝑑𝜅^𝜋0𝐶_𝑔−𝐶^Stokes_𝑔²𝑑𝜅∫^𝜋₀𝐶^Stokes_𝑔²,𝑑𝜅 (11) 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, 𝛽₁+𝛽₂=−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.

(a)

(b)

Figure 1 

Comparisons of phase celerity (a) and group velocity (b) between equations (solid line) and Padé [2,2] approximation (dashed line).

Specific values of 𝛽₁ and 𝛽₂ cannot be determined yet by this process, which should be completed by analyzing shoaling property of the model. Following the same procedure in [6], 𝛽₂ can be determined by finding the minimum value of the integral ∫^𝜋₀(𝛼_𝑠−𝛼^Stokes_𝑠)²𝑑𝜅 (here, 𝛼_𝑠 is the shoaling coefficient, and the superscript “Stokes” represents the analytical solution), and 𝛽₂=−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<𝜅<𝜋.

Figure 2 

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

Nonlinear property is analyzed by solving second- and third-order problems. Parameter 𝛾 is found to contribute to second-(𝑎₂), third-order (𝑎₃₃) transfer functions, as well as amplitude dispersion 𝜔₃. In this process, the joint error (defined by (11)) of the second-order transfer function (𝑎₂) and amplitude dispersion (𝜔₃) is chosen as the criterion to determine 𝛾. 𝜔₃ is chosen as the criterion instead of 𝑎₃₃ because the amplitude dispersion is more important than 𝑎₃₃ for nonlinear wave evolution. 𝜔₃ is the factor that controls the nonlinear process of amplitude dispersion, whereas 𝑎₃₃ is a smaller nonlinear quantity compared with 𝑎₂, 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.

(a)

(b)

(c)

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).

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].

3. 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.

3.1. 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.0 m), 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].

Figure 4 

The experimental setup of Luth et al. [26]. All dimensions are in (m).

Figure 5 

Wave surface elevations of Luth’s experiment (circle) [26] and the numerical results of the enhanced model (solid) and the original model (dashed).

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 𝑎₂ 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.

(a)  First harmonic

(b) Second harmonic

(c) Third harmonic

(d) Fourth harmonic

(e) Surface elevation

Figure 6 

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).

3.2. 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: 𝑠(𝑡)=s₀|cos(Ωt)|cos𝜔₀𝑡𝜔withΩ=⁰,20 (12) where 𝑠₀ is the forcing amplitude, 𝜔₀ is the carrier frequency, and Ω is the modulation frequency. The spectrum of the signal achieves its maximum at the carrier frequency 𝜔₀, plus several smaller peaks at discrete frequencies spaced by 𝜔₀/10 with three dominant modes. The period of the carrier wave is 0.9 s and the corresponding value of 𝑘ℎ 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 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 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 high-frequency 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.

(a)

(b)

(c)

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).

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 (𝑂(𝜀³)) 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.

3.3. 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, ℎ₁, the topography is specified as ℎℎ₁=ℎ_minℎ₁+ℎ1−^minℎ₁1((,ℎcoshtan𝜋𝑥/2𝐿))𝑥≤𝐿,ℎ₁=ℎ_minℎ₁,𝐿<𝑥≤2𝐿, (13) where ℎ_min/ℎ₁=0.2 and 𝐿=50ℎ₁. Initial surface elevation is given as𝜂ℎ₁=𝑎_𝐼ℎ₁cos2𝜋𝑁_𝑤𝑥/𝐿((,𝜂coshtan𝜋𝑥/2𝐿))𝑥≤𝐿,ℎ₁=0,𝐿<𝑥≤2𝐿, (14) where 𝑎_𝐼/ℎ₁ is the initial amplitude, and the number of waves 𝑁_𝑤=10. In the simulation, ℎ₁ 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.

(a)

(b)

(c)

Figure 9 

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

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<𝑥/ℎ₁<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.

Figure 10 

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

4. 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.

Appendix

Analysis of the equations for monochromatic waves is as follows.

In one dimension, and on a horizontal bottom, (1) and (8) simplify to 𝜂_𝑡+ℎ𝑢_𝑥+𝜀𝜂𝑢_𝑥𝑢=0,_𝑡+𝜀𝑢𝑢_𝑥+𝜂_𝑥+𝜇²𝛽₁+𝛽₂−13ℎ²𝑢_𝑥𝑥𝑡+𝛽₁+𝛽₂ℎ²𝜂_𝑥𝑥𝑥+𝜀𝜇²2𝛽₁2𝛾−3ℎ𝜂𝑢_𝑥𝑥𝑡−ℎ𝜂_𝑥𝑢_𝑥𝑡+13+3𝛽₁+3𝛽₂ℎ²𝑢_𝑥𝑢_𝑥𝑥+𝛽₁+𝛽₂−13ℎ²𝑢𝑢_𝑥𝑥𝑥+𝛽₂𝛾ℎ𝜂𝑢_𝑡+𝜂𝜂_𝑥𝑥𝑥+2𝛽₁𝛾ℎ𝜂𝜂_𝑥𝑥𝑥+𝛽₂𝛾ℎ𝜂𝑢_𝑡+𝜂_𝑥𝑥𝑥+𝜀²𝜇²𝛽₁+16𝜂²𝑢_𝑥𝑥𝑡−13ℎ𝜂𝑢_𝑥𝑢_𝑥𝑥−13ℎ𝑢_𝑥𝑥𝜂𝑢_𝑥+ℎ𝜂𝑢²_𝑥𝑥−12𝜂²𝑢_𝑥𝑡𝑥−23ℎ𝜂𝑢𝑢_𝑥𝑥𝑥+𝛽₁𝛾ℎ𝜂𝑢²_𝑥𝑥𝑥+12𝛽₂𝛾𝜂𝑢²_𝑥𝑥𝑥+𝛽₂1𝛾𝜂2ℎ𝑢²_𝑥+𝛾𝜂𝑢_𝑡+𝑢𝜂_𝑥𝑥=𝑂𝜀³𝜇²,𝜇⁴. (A.1) Substituting the Stokes-type expansions (10) in the main text into (A.1) and using the following notations 𝑚^(𝑛)₁₁=𝑛𝜔₀,𝑚^(𝑛)₂₁=−𝑛𝑘1−𝛽₁+𝛽₂(𝑛𝜅)²,()𝑚𝜅=𝑘ℎ,𝑛=1,2,3_(𝑛)12=−𝑛𝜅,𝑚^(𝑛)₂₂=𝑛𝜔₀1−𝛽₁+𝛽₂−13(𝑛𝜅)² (A.2) leads to the following perturbation equations.

(i) First-order equations 𝑚⁽¹⁾₁₁𝑚⁽¹⁾₁₂𝑚⁽¹⁾₂₁𝑚⁽¹⁾₂₂𝑎₁𝑢₁=00. (A.3) From (A.3), we obtain the amplitude of the first order velocity 𝑢₁=𝜔₀𝑎₁/𝜅 and the celerity 𝜔𝐶=⁰𝑘=ℎ1−𝛽₁+𝛽₂𝜅²1−𝛽+𝛽₂−13𝜅². (A.4)

(ii) Second-order equations 𝑚⁽²⁾₁₁𝑚⁽²⁾₁₂𝑚⁽²⁾₂₁𝑚⁽²⁾₂₂𝑎₂𝑢₂=𝐹⁽²⁾𝐺⁽²⁾ (A.5) with 𝐹⁽²⁾=−𝑘𝑎₁𝑢₁,𝐺⁽²⁾=−𝛾𝛽₁+52𝛽₂𝑘³ℎ𝑎²₁+12−2𝛽₁+𝛽₂𝜅²𝑘𝑢²₁+𝛽₁+52𝛽₂5𝛾−6𝜔₀𝑘𝜅𝑢₁𝑎₁. (A.6) From the equations, we obtain the amplitudes of the second-order oscillating surface elevation and velocity 𝑎₂=𝐹⁽²⁾𝑚⁽²⁾₂₂−𝐺⁽²⁾𝑚⁽²⁾₁₂𝑚⁽²⁾₁₁𝑚⁽²⁾₂₂−𝑚⁽²⁾₁₂𝑚⁽²⁾₂₁,𝑢₂=𝐺⁽²⁾𝑚⁽²⁾₁₁−𝐹⁽²⁾𝑚⁽²⁾₂₁𝑚⁽²⁾₁₁𝑚⁽²⁾₂₂−𝑚⁽²⁾₁₂𝑚⁽²⁾₂₁. (A.7)

(iii) Third-order equations 𝑚⁽³⁾₁₁𝑚⁽³⁾₁₂𝑚⁽³⁾₂₁𝑚⁽³⁾₂₂𝑎₃₃𝑢₃₃=𝐹₃⁽³⁾𝐺₃⁽³⁾𝑚₍₁₎₁₁,𝑚⁽¹⁾₁₂𝑚⁽¹⁾₂₁𝑚⁽¹⁾₂₂𝑎₃₁𝑢₃₁=𝐹₁⁽³⁾−𝜔₃𝑎₁𝐺₁⁽³⁾−𝜔₃𝑚⁽¹⁾₂₂𝑢₁ (A.8)(A.9) with 𝐹₃⁽³⁾=32𝑘𝑎₂𝑢₁+𝑎₁𝑢₂,𝐺₃⁽³⁾=−𝛾²14𝛽₁+𝛽₂𝑘³𝑎³₁+14𝛽₁+𝛽₂𝛾²₂−13𝜔₀𝑘²𝑢₁𝑎²₁−134𝛽₂+2𝛽₁𝛾𝑘³ℎ𝑢²₁𝑎₁−9𝛾𝛽₁+2𝛽₂𝑘³ℎ𝑎₂𝑎₁+13𝛽₂+8𝛽₁𝛾−143𝜔₀𝑘²ℎ𝑢₂𝑎₁+5𝛽₂+𝛽₁4𝛾−3𝜔₀𝑘²ℎ𝑎₂+121−27𝛽₁+𝛽₂𝜅²+32𝑘𝑢₂𝑢₁,𝐹₁⁽³⁾=12𝑘𝑎₂𝑢₁+𝑎₁𝑢₂,𝐺₁⁽³⁾=−𝛽₂+14𝛽₁𝛾²𝑘³𝑎³₁+𝛽₂+14𝛽₁𝛾²−13𝜔₀𝑘²𝑢₁𝑎²₁−4𝛽₂+7𝛽₁𝛾ℎ𝑘³𝑎₂𝑎₁+54𝛽₂−2𝛽₁𝛾+1ℎ𝑘³𝑢²₁𝑎₁−8𝛽₁+5𝛽₂2𝛾−3𝜔₀𝑘²ℎ𝑢₂𝑎₁−𝛽₁+𝛽₂2𝛾−3𝜔₀𝑘²ℎ𝑎₂𝑢₁+123−𝛽₁−𝛽₂𝜅²+1𝑘𝑢₂𝑢₁. (A.10) From (A.8), we can obtain the solutions for 𝑎₃₃ and 𝑢₃₃ as 𝑎₃₃=𝐹₃⁽³⁾𝑚⁽³⁾₂₂−𝐺₃⁽³⁾𝑚⁽³⁾₁₂𝑚⁽³⁾₁₁𝑚⁽³⁾₂₂−𝑚⁽³⁾₁₂𝑚⁽³⁾₂₁,𝑢₃₃=𝐺₃⁽³⁾𝑚⁽³⁾₁₁−𝐹₃⁽³⁾𝑚⁽³⁾₂₁𝑚⁽³⁾₁₁𝑚⁽³⁾₂₂−𝑚⁽³⁾₁₂𝑚⁽³⁾₂₁. (A.11) When we eliminate 𝑢₃₁ from (A.9) to get the equation for 𝑎₃₁, the left hand side of the resulting equations will equate to zero automatically due to the first-order equation (A.3); thus, the right hand side of the resulting equation must also equate to zero from the solvable condition. This yields the equation for the amplitude dispersion 𝜔₃𝜔₃=𝐹₁⁽³⁾𝑚⁽¹⁾₂₂−𝐺₁⁽³⁾𝑚⁽¹⁾₁₂𝑚⁽¹⁾₂₂𝑎₁−𝑚⁽¹⁾₁₂𝑚⁽¹⁾₂₂𝑢₁. (A.12) All the corresponding analytical solutions from Stokes theory could be found in [6].

Acknowledgments

The authors would like to thank the finical support from the National Natural Science Foundation of China (Grant nos. 51009018 and 51079024), the State Key Laboratory of Coastal and Offshore Engineering, Dalian University of Technology (Grant no. LP1105), National Marine Environment Monitoring Center, State Oceanic Administration (Grant no. 201206), and Key Laboratory of Water-Sediment Sciences and Water Disaster Prevention of Hunan Province (Grant nos. 2012SS02 and 2013SS02). Special thanks are given to three anonymous reviewers. Their valuable suggestions greatly help improving paper’s quality.