Equations for small amplitude shallow water waves over small bathymetric variations

A generalized version of the abcd-Boussinesq class of systems is derived to accommodate variable bottom topography in two-dimensional space. This extension allows for the conservation of suitable energy functionals in some cases and enables the description of water waves in closed basins with well-justified slip–wall boundary conditions. The derived systems possess a form that ensures their solutions adhere to important principles of physics and mathematics. By demonstrating their consistency with the Euler equations and estimating their approximation error, we establish the validity of these new systems. Their derivation is based on the assumption of small bathymetric variations. With practical applications in mind, we assess the effectiveness of some of these new systems through comparisons with standard benchmarks. The results indicate that the assumptions made during the derivation are not overly restrictive. The applications of the new systems encompass a wide range of scenarios, including the study of tsunamis, tidal waves and waves in ports and lakes.


Introduction
The propagation of water waves is mathematically described by partial differential equations derived by Leonard Euler (1707Euler ( -1783) ) in 1757 [25] for the study of incompressible flow of an ideal fluid.These equations are known as the (incompressible) Euler equations.To present the Euler equations, we first define the horizontal and vertical coordinates of three-dimensional space as x = (x, y) ∈ R 2 and z, respectively.The undisturbed water level at rest is considered to be at z = 0.The velocity of the fluid is denoted as (u(x, z, t), w(x, z, t)) = (u(x, z, t), v(x, z, t), w(x, z, t)), separating the vertical component from its horizontal components.The free surface elevation above its undisturbed level is represented by η(x, t), and the depth of the ocean floor is denoted as D(x, t).(Note that we consider time-dependent bottom topography to allow for the modeling of water waves generated by the motion of the ocean floor.)The horizontal gradient is denoted by ∇ = (∂ x , ∂ y ), while t represents time.The Euler equations consist of equations for the conservation of mass and momentum, which, for −D(x, t) < z < η(x, t), can be written as: Since our aim is to model water waves in oceans bounded above by a free surface and below by the ocean floor, Euler's equations require the following boundary conditions at the free surface (4) η t + u • ∇η − w = 0, on z = η(x, t) , and at the bottom (5) D t + u • ∇(z + D(x, t)) + w = 0 , for z = −D(x, t) .
Another important ingredient in the theory of water waves is the irrotationality condition which is unavoidable when potential flow is assumed.While the pressure inside the fluid volume is an unknown function, the pressure at the free surface is a given constant (7) p = p atm , for z = η(x, t) .
The study of these particular equations presents a great difficulty due to the fact that the fluid domain Ω t , which includes the free surface of the water and the impermeable bottom (and possibly solid wall boundaries), is practically unknown since it is bounded by the unknown free-surface z = η(x, t).Furthermore, the domain is not stationary, adding further complexity to the situation.As a result, scientists have derived approximate equations based on simplifying assumptions.The development of such approximations began after John Scott Russell (1808-1882) discovered solitary waves and published his report on Waves in 1844, [47], and this work has continued to the present day.Examples of these approximations are the Boussinesq systems, which describe long waves of small amplitude.The first Boussinesq system was derived by J. Boussinesq himself [12,13], and subsequent improvements, also referred to as Boussinesq systems, have been derived.It is important to note that there are only a few water wave systems that are mathematically and physically justified.For instance, the Euler equations have been proven to be well-posed (in the Hadamard's sense) only in unbounded domains [53,54,38], while among all the Boussinesq systems that have been derived, only the system proposed in [32] has been proven to be well-posed in a bounded domain with slip-wall boundary conditions.
Several Boussinesq systems have been extensively used in the literature to describe nonlinear and dispersive waves with variable bottom topography.One of the notable Boussinesq systems is the Peregrine system [46]: (8) and the Nwogu (or extended Boussinesq) system [44] (9) , and θ = 1/5.(The parameter θ determines the depth in which the horizontal velocity u is computed, cf.e.g. ( 33)).Although in some cases their Cauchy problems have been proved well-posed, [23,42,37], they do not preserve any meaningful approximation of the total energy of the Euler equations.While the particular property is not restrictive, it is rather desirable from physical as well as numerical point of view.Specifically, we can ensure physical validity for long periods of time when the solution of a system preserves an energy functional [26].Similarly, systems of BBM-BBM-type [42] or the systems with flat bottom topography of [8] fail at describing water waves in a basin with slip-wall boundary conditions [21,22,20].Some other Boussinesq systems even fail in having traveling waves of finite energy, proving that asymptotic justification does not imply physical relevance, [9,10].For all these reasons in [34,32] we derived and analyzed a new Boussinesq system of BBM-BBM-type which can be used in a straightforward manner in bounded domains with slip-wall boundary conditions.Moreover, its solutions preserve the same energy functional as its non-dispersive counterpart, i.e. the shallow water-waves (Saint-Venant) equations [52] making the system attractive from numerical and physical point of view.
In general, Boussinesq systems are extensions of equations that describe the propagation of water waves in one direction, such as the Korteweg-de Vries (KdV) equation.Various versions of the KdV equation have been derived to account for variable bottom topography [33,30,24].On the other hand, such models cannot describe wave reflections on walls or obstacles.Moreover, high-order Boussinesq equations, including the Serre-Green-Naghdi equations [49,28,39], have only been proven to be well-posed in unbounded domains [31].Note that the Serre-Green-Naghdi equations and their optimized counterpart of [16] preserve an energy functional.
In this article, we start by presenting in Section 2 the asymptotic derivation of a new class of Boussinesq systems that generalizes the so-called abcd-Boussinesq systems of [8,7] in the case of variable bottom topography.We also establish their consistency with the Euler equations along the lines of [8] in Section 3.These systems also generalize the systems of [14] in two-dimensional domains, and their derivation is based on the same assumption of small bottom variations.Although we do not provide any well-posedness results, these new systems appear to have analogous theoretical properties to the system of [8], and they can further be used with slip-wall boundary conditions.The well-posedness of the BBM-BBM system was proved in [32], while the well-posedness of other Boussinesq systems of the same family, such as the Bona-Smith systems, is very similar.The new system that corresponds to the Nwogu system can be applied with accurate slip-wall boundary conditions, as demonstrated in Section 5, and is validated against laboratory data in Section 6.For the Nwogu system as well as for its regularized counterpart, we have justified this set of boundary conditions in one-dimensional domains in [41].Taking into account the results of [32,34], we conclude that the assumption of small bottom variations is not restrictive, and the necessary nonlinear and dispersive behaviors of the Boussinesq systems are retained.This was also demonstrated in [40].For the sake of completeness, we present a variational derivation of the new energy-conserving abcd-Boussinesq systems in Section 4. The new systems have similar energy conservation properties to the abcd-Boussinesq systems of [8,7], and thus we generalize them in the case of variable bottom topography.

Derivation of the new Boussinesq equations
Consider the depth function D(x, t) = D 0 + D b (x) + ζ(x, t), where D 0 is a characteristic (mean) depth and D b is smooth and varies gently in R 2 .The displacement of the bottom's moving part, denoted by ζ, is small compared to the depth and does not exceed the amplitude of the free surface elevation.If λ 0 is a typical wavelength, a 0 a typical wave height, and d 0 a typical order of bottom topography variations, then we consider the (scaled) non-dimensional independent variables (10) x = and also the non-dimensional dependent variables where c 0 = √ gD 0 is the linear speed of propagation.Denoting as usual and assuming that 0 < ε ≈ σ 2 ≲ β ≪ 1 all the aforementioned Boussinesq systems are asymptotically equivalent, in the sense that all of them occur from the Euler equations after appropriate scaling and discarding high-order terms with the same order in ε and σ 2 .In the new variable the scaled depth will be D = 1 + β Db + ε ζ.Note that Dt = O(ε) and thus we assume a slowly varying bottom.The requirement β ≪ 1 is used to describe small bottom variations compared to the wavelength.
Using the new variables, we write Euler's equations ( 1)-(3) in non-dimensional form: for − D < z < εη.The irrotationality condition is also transformed to ∇ × ũ = 0 , (15)   20) is exact and it is the mass equation of the hydrostatic shallow water wave equations and of Peregrine's original system.
Note that Using the previous approximations we write the general abcd-Boussinesq system as The reason for considering such approximations is to ensure that certain combinations of the parameters a, b, c, d lead to systems of significant interest of [7] that preserve meaningful approximations of the total energy even with variable bottom topography.
Discarding the high-order terms in ( 50)-( 51) we write the equations ( 50)-( 51) in dimensional variables as for fixed µ, ν ∈ R and 0 ≤ θ ≤ 1.The first equation corresponds to an approximation of the mass conservation law, while the second equation corresponds to an approximation of the momentum conservation law.Note that the system (52) has been derived in such a way that the ∇×u is preserved for all times t ≥ 0 for stationary bottom topography, and thus the irrotationality of the water waves is respected exactly.Moreover, we observe that in most cases the presence of the regularization operators I − b∇•(D 2 ∇•) and I − d∇∇•(D 2 •) requires for well-posedness that b, d ≥ 0 [19].As we shall see later, conservation of energy requires additionally b = d ≥ 0. Water waves are dispersive by their nature.For modeling purposes we employ the dispersion relation of the linearized equations and with flat bottom topography as a measure of accuracy.Specifically, to study the dispersion characteristics of the system (52) we first consider the expansion of the linear dispersion relation for general periodic solutions of the Euler equations and of the form e i(kx−ωt) , which is where k is the wavenumber and ω the frequency of the wave.Recall that the period of a linear wave is defined as T = 2π/ω, the phase speed is ω/k, the linear speed of propagation is √ gD and the wavelength is λ = 2π/k.The system (52) with flat bottom D has linear dispersion relation given by the formula (55) If we choose the coefficients a, b, c, d of ( 53) such that the formula (55) coincides with (54) up to the term of order (Dk) 4 , then we obtain the optimized extended Boussinesq system for ν = 0, µ = 1 and θ 2 = 1/5.This is the system which we call Nwogu system and is practically the only system of the abcd class of systems with such an optimal dispersive behavior.Notable systems of the abcd class that have appeared in various works [6,7,11,19] are the following: Next, we will prove that all these systems are asymptotically equivalent.On the other hand, it worth mentioning that some of these systems cannot be used for the physical description of surface water waves.For example, the KdV-KdV system does not admit classical solitary wave solution (except perhaps embedded ones) but instead has generalized solitary waves that decay exponentially to small periodic orbits.Such solutions have infinite energy and thus they are not physically meaningful [19,9,10].

Consistency with the Euler equations
In this section we consider the Boussinesq system ( 50)-( 51) with stationary bottom topography and parameters a, b, c, d as in (53) (56) Note that in the notation of this section we have dropped the tilde in nondimensional variables.Consider also the Euler equations in the nondimensional and scaled Zakharov formulation: (57) where ψ is the velocity potential and Following [39] (Section 5.3) one can derive the following four-parameter family of Boussinesq systems of order where u is the horizontal velocity at the level (D + εη)z + εη and z = −1 Here, we will prove consistency results between the asymptotic model (56) and the Euler equations in case of small amplitude topographic variations β = O(σ 2 ).However, note that from [39] the Zakharov system (57) and the Boussinesq system (58) are consistent (see Remark 5.35 in [39]).In particular, for δ Next we will show in detail the consistency of the new system (56) with Nwogu system (59).But we shall first introduce the definition of consistency along the lines of [8].
2 ) of the new Nwogu system of Boussinesq type (60) with the initial condition (η 0 , u 0 ) satisfies The proof of existence is demonstrated in ([48, Theorem 1.1]).Applying the operator curl (•) to the second equation of the system (60) we see that ∇×u t = 0 , which implies ∇×u = ∇×u 0 , and the result follows.□ Lemma 3.2.Let u, η be solutions either of (59) or (60) and D = 1 + βD b be a bottom parametrization, then Proof.The first identity follows from the fact that D(x) = 1 + βD b (x).The other two follow from the system (60).□ Here, we present the main result of this section.For simplicity, we take β = ε = O(σ 2 ) for σ 2 ≪ 1.
Proof.For the sake of simplicity we drop the index ε from the notation of (η ε , u ε ) but we will still assume solutions of (59).From the first equation of (59), it suffices to prove that Using Lemma 3.2 we have where in the sense of consistency.Thus, (η ε , u ε ) is consistent with the first equation in (56) with a residual of order ).The same situation holds for the second equation of ( 59).Here we denote c = Bν, d = B(1 − ν), ν ∈ R and B = c + d.From Lemma 3.1, and again using Lemma 3.2 we see that in the sense of consistency.Therefore, (η ε , u ε ) is consistent with the second equation of (56) with a residual of order β 2 in L ∞ ([0, T ε ]; H s−4 (R 2 ) 2 ).We conclude that (η ε , u ε ) is consistent with (56) in the sense of definition 3.1 and the proof is complete.□ The Euler equations are consistent with the System (59) at order ε 2 in C([0, T ε ]; Ḣs+10 × H s+5 ) as it was proved in [39,Corollary 5.31].A direct consequence of this is that the Euler equations are consistent with the new class of systems (56).
Assume moreover that there is h min > 0 such that the non-cavitation condition h = 1 + εη ε 0 + βD b ≥ h min is satisfied, and that there exists T > 0 such that − There is a unique family of solutions ) of the Euler equations (57) with initial condition U 0 .− There is a unique family of solutions The error estimate (62) follows using the triangle inequality and the results of [8, Corollary 2] for both (η ε E , u ε E ) and (η ε B , u ε B ). Indeed, performing symmetrization like in [8] (using appropriate nonlinear change of variables) to the new Boussinesq systems we obtain a class of systems with symmetric nonlinearity.This is the class of the so-called symmetric systems (see also [8]), which is denoted by S ′ θ,µ,ν .Let (u, η) be a solution consistent with S ′ θ,µ,ν and (θ, µ, ν) fixed satisfying a = c.By applying the nonlinear change of variables, we have that (u, η) is consistent with the class S ′ θ,µ,ν .Denote by Σ = S ′ θ,µ,ν the corresponding subclass of symmetric non-linear and dispersive systems.Moreover, for all ε > 0, let (u ε Σ , η ε Σ ) be the exact solution of a system of the class Σ, then where Corollary 2] for the solution of the Boussinesq system (u ε B , η ε B ) for all ε > 0, we obtain I ≤ C ε 2 t.Now, thanks to the Proposition 3.1 and [8, Corollary 2] for the solution of the Euler equations (u ε E , η ε E ), for all ε > 0, yields II ≤ C ε 2 t, which implies (62) and completes the proof.□ Remark 3.1.We have established that the new Boussinesq systems are asymptotically equivalent.However, Boussinesq systems have different well-posedness requirements for different sets of the parameters a, b, c, d.In the previous theorem we consider p large enough so as to embrace all valid Boussinesq systems.The maximal time of existence also varies depending on the particular system but all systems follow the general results of [20,48].The situation is much different when bounded domains are considered.The results then depend on the choice of the system as well as on the boundary conditions.A discussion on justified slip-wall boundary conditions for Boussinesq systems is presented in Section 5. Furthermore, the shape of the solitary waves and their dynamics such as interactions and collisions depend also on the choice of the parameters a, b, c, d.

Energy conservation and variational derivation
In this section we consider Ω = R 2 while the derivations hold true in an bounded domain as well, under the hypothesis that the solutions along with their derivatives vanish appropriately on the boundary.For example, the boundary conditions are adequate to carry the following analysis.Once again we consider stationary bottom topography as there is no meaning to talk about energy conservation with a moving bottom.
Consider the system ( 50)-( 51) with stationary bottom D. Setting and using the boundary conditions we have Therefore, for a ≤ 0, c ≤ 0 and b = d ≥ 0 we define the energy functional ( 63) which in dimensional variables can be written as ( 64) Using the particular energy function, we can derive energy conservative systems of the class (52) using the variational method of [15].We define the kinetic energy as and the potential energy as Then we define the action integral where the mass conservation is imposed with the help of the Lagrange multiplier ϕ(x, t).Then the Euler-Lagrange equations for the action integral I are the following: Eliminating ϕ in (67) using (66) we obtain the momentum conservation equation In this way we re-derived the subclass of systems ( 52) that preserve a reasonable form of energy, and these must have necessarily b = d.From (66) we observe that the new systems imply potential flow, and the conservation of the ∇×u follows.

Valid slip-wall boundary conditions
Practical problems involving water waves are usually posed in bounded polygonal domains Ω.For example, waves in a port or waves interacting with the sides of a basin.Wall boundary conditions must then be imposed alongside the model equations.When the waves can slip on the wall without any friction, a slip-wall boundary condition of the form u • n = 0 needs to be imposed.Here, n is the outward unit normal vector on the boundary ∂Ω.While no other boundary conditions are usually necessary, it is interesting to observe that the Euler equations implicitly obey additional boundary conditions.For example, the Neumann condition ∇η • n = 0 holds true on ∂Ω when Ω is a polygonal domain, [36].To derive this particular boundary condition, first observe that ∇p • n = 0 on ∂Ω.This can be obtained by taking the inner product of equation ( 3) with n on ∂Ω.Writing the dynamic boundary condition (7) as p(x, η(x, t), t) = p atm , and taking the horizontal gradient operator, we obtain ∇p + p z ∇η = 0.After multiplying by the normal vector n and using the fact that ∇p • n = 0 on ∂Ω, we obtain the boundary condition ∇η • n = 0 on ∂Ω.If the domain is a general domain where well-posedness can be established, appropriate polygonal approximations Ω n ⊂ Ω can be constructed such that lim Ω n = Ω, and the particular Neumann condition can be used in more general situations.This particular Neumann condition is required for the well-posed Boussinesq systems of Bona-Smith and BBM-BBM type in bounded domains with slip-wall boundary conditions, and based on the previous analysis, it appears to be a physical condition [34,32].Note that for the original Nwogu system (9), in addition to the slip-wall boundary conditions u • n = 0 must be satisfied on the boundary ∂Ω, [51].This boundary condition is satisfied by the solutions of the Euler equations only when the bottom is flat, while for general bottoms, we obtain a simpler condition.
From the irrotationality condition (6), we have u z = ∇w, and thus ∇w•n = 0 on ∂Ω.Moreover, differentiating the mass equation ( 1), we obtain ∇(∇•u) = −∇w z , which implies the boundary condition ∇(∇•u) • n = 0 on ∂Ω.Similar conditions hold true for BBM-BBM-type systems.For example, the slip-wall conditions imply that ∂Ω, then this is satisfied for t ≥ 0. Due to the assumption of mild bottom variations in the derivation of the particular BBM-BBM system and also for the new Nwogu system, we have that ∇(∇•u) • n ≈ 0, which means that even if the boundary condition is not satisfied exactly, it is a reasonable approximation of the exact boundary condition.The same will be true for Nwogu-type systems of the form (52).This shows that some of the new abcd-systems (at least the BBM-BBM, Bona-Smith, Nwogu, and regularized Nwogu systems) can be used in practical situations in bounded domains with more accurate boundary conditions compared to the other Boussinesq systems.It is worth mentioning that of all the known Boussinesq systems, only the BBM-BBM system has been proved to be well-posed in bounded domains with slip-wall boundary conditions [32], while the well-posedness in unbounded domains follows the work of [20] in a similar manner.The Bona-Smith systems can also be proved to be well-posed in bounded domains, but this will be discussed in a separate work.

Experimental validation
The mild slope assumption used in the derivation of the new abcd-Boussinesq systems with variable bottom topography raises a validity question: What are the minimum restrictions on the bottom variations to ensure the accuracy of the new systems?This question requires a more detailed study than the one presented here.However, in this section, we demonstrate that the particular assumption is not a significant barrier.On the contrary, the simplicity of the equations and their accuracy make them highly attractive for the mathematical modeling of long waves.
Among the systems we study in this section is the new Nwogu system.Although this system does not preserve any form of energy, it has optimal linear dispersion relation for θ 2 = 1/5, µ = 1 and ν = 0 and outperforms other systems in experiments where the bottom variations as well as the steepness of the waves stress the validity of Boussinesq systems.For the sake of convenience we rewrite the new Nwogu system in the form (68) We also consider the new BBM-BBM system (θ 2 = 2/3, ν = µ = 0) (69) which preserves the energy functional (64) with a = c = 0. We also study the new Peregrine system (θ 2 = 1/3, ν = µ = 0) (70) All the experiments primarily examine the impacts of bottom topography on long one-dimensional waves, which serve as standard benchmarks in the literature of Boussinesq systems.However, it should be noted that while the selected experiments somewhat deviate from the scope of Boussinesq systems, in the sense that the waves are of larger amplitude compared with what Boussinesq systems have been derived to describe, the results appear to be remarkably promising.Future work could focus on conducting more advanced experiments to further enhance our understanding.For all the numerical simulations, we utilized previously developed Galerkin/Finite Element methods.For the Nwogu system, we employed the numerical method proposed in [41].The method proposed by [3] was utilized for the BBM-BBM system.Lastly, for the Peregrine system with variable bottom topography, we used the method introduced in [1], which was extended accordingly.In all experiments, we employed cubic spline elements and the classical explicit fourth-order Runge-Kutta method with four stages.6.1.Shoaling and reflection of solitary waves.In the first experiment, we investigate the capability of the new models to depict the shoaling behavior of waves on gentle slopes, as outlined in the study [18].Specifically, we focus on a one-dimensional computational domain Ω = [−100, 20] with a bottom topography characterized by the following function: For the numerical simulation of this experiment, we employ wall boundary conditions at the endpoints x = −100 and x = 20.All variables used in the simulation are expressed in SI units.
We begin by investigating the reflection of two right-traveling solitary waves with amplitudes A 1 = 0.07 m and A 2 = 0.12 m. e two solitary waves have been horizontally shifted so that their maximum at t = 0 is reached at x 0 = −30.These waves initially propagate over the sloping bottom and then collide with the wall at x = −20, where they are reflected and propagate to the left.We record the free surface elevation at three wave gauges positioned at x = 0, x = 16.25, and x = 17.75.The results of these experiments for the solitary wave with A 1 = 0.07 m are presented in Figure 1, while the results corresponding to the other solitary wave are shown in Figure 2. We observe that the solutions of the Reflection of shoaling solitary wave of amplitude A = 0.12 m by a vertical wall BBM-BBM system exhibit better agreement with the experimental data compared to the other systems.This behavior was also observed in [34] with a BBM-BBM system that does not preserve any energy functional.In both cases, the values of ε and σ are small but not exactly within the range of validity of Boussinesq systems.For the first case, the characteristic values of the initial solitary wave are ε ≈ 0.23 and σ ≈ 0.035.For the second case, the parameters are ε ≈ 0.4 and σ ≈ 0.05.Moreover, the bottom variations parameter is β = 0.028.It is important to note that these values change during the transformation of the shoaling solitary wave.Particularly, the value of ε increases as the wave passes over the slope.This leads to a significant amount of error, which explains the deviations between the numerical and experimental results in Figures 1 and 2.
We conclude the study of shoaling waves by examining an experiment conducted by [29].In this particular experiment, our focus is on the shoaling behavior of a solitary wave with an amplitude of A = 0.2 propagating over a bottom topography characterized by the following function: We recorded the free surface elevation at five wave gauges that are located at x = 20.96,22.55, 23.68, 24.68, 25.91, respectively.It is important to note that in this experiment, all variables are scaled.This experiment proves to be more demanding compared to the previous one, requiring appropriate calibration of the location of the initial data in order to align the numerical solutions with the experimental data.This calibration is necessary because solitary waves of the same amplitude travel at different speeds for different Boussinesq systems [19].As a result, solutions from different systems reach the wave gauges at different times.
To demonstrate this, let's examine the shoaling wave at the gauge located at x = 25.91, which has an approximate amplitude of a 0 ≈ 0.5, while the depth is D(x) ≈ 0.26, resulting in ε ≈ 1.92.Although this value of ε exceeds the range of validity for any Boussinesq system, the obtained results are still quite satisfactory after performing the appropriate calibration of the initial condition.The characteristic values for the initial solitary wave and the bottom in this case are ε = 0.2, σ ≈ 0.05, and β ≈ 0.028.
Surprisingly, the BBM-BBM system required the least amount of phase calibration and exhibited the most accurate amplification of the solution among the tested systems.The results are illustrated in Figure 3. Once again, the BBM-BBM system demonstrates superior approximation of shoaling waves compared to the other systems.Similar observations were reported for the classical Nwogu system in [27], emphasizing the need for initial condition calibration in all Boussinesq systems (except for the Serre-Green-Naghdi equations [43]).
In the comparison between the classical Peregrine system (8) and the new Peregrine system, we found that we needed to apply the same calibration to the initial conditions for both systems in order to align the phase of the numerical solutions with the experimental data.However, when we utilized the Serre-Green-Naghdi system for the same experiment, there was no need for calibration of the initial data's location [43].This implies that even in extreme situations involving Boussinesq systems, the assumption of small bottom variations is not a disadvantage, as any Boussinesq system will inevitably lose accuracy in such extreme experiments.It is important to note that the phase error observed in the first experiment was negligible, and therefore, no calibration was necessary.Please note that the numerical method used to solve the BBM-BBM system in our study is not conservative.However, we were able to measure the total energy of the solution (64), which was approximately 0.162 for the first experiment, 0.362 for the second experiment, and 0.135 for the third experiment.6.2.Periodic waves over a submerged bar.In this scenario, we assess the performance of the aforementioned Boussinesq systems in a highly challenging laboratory experiment designed to investigate the nonlinear and dispersive properties of Boussinesq systems with variable bathymetry [4].This specific experiment has been extensively used to validate different Boussinesq systems, including the original Nwogu system, and is regarded as a standard benchmark for numerical models [17,35,50].
In this experiment, small-amplitude periodic waves are generated by a wavemaker and propagate over a bathymetry defined by the function: x ∈ [14, 17] 0.4, elsewhere which is depicted in Figure 4.
During the interaction of the generated waves with the submerged bar's uphill portion, they undergo shoaling and steepening.Subsequently, the waves propagate at an increased speed over the downhill part of the bar, encompassing a broad spectrum of long wave numbers.In this experiment, the input wave has a period of τ = 2.02s and an amplitude of A = 0.01m.The water's free surface is recorded at several wave gauges positioned at x = 10.5, 12.5, 13.5, 14.5, 15.7, 17.3, 19.0, and 21.0.For numerical simulations, the waves are generated using a wavemaker located at x = 2.01, following the methodology described in [34].
The locations of the wavemaker and the wave gauges are depicted by vertical broken lines in Figure 4. Experimental measurements of the free surface elevation at the different gauges are presented in Figure 5.In the same figure, numerical results obtained using the standard Galerkin/Finite Element method with cubic splines, combined with the classical fourth-order Runge-Kutta method, are displayed for solving the one-dimensional initial-periodic boundary-value problem for the systems (68), (52) corresponding to the new Nwogu system, BBM-BBM system, and Peregrine system.Further details regarding the specific numerical method can be found in [2] and [41].It is noteworthy that the numerical solution exhibits an impressive level of agreement with the experimental data.

Figure 1 .
Figure 1.Reflection of shoaling solitary wave of amplitude A = 0.07 m by a vertical wall

Figure 2 .
Figure 2. Reflection of shoaling solitary wave of amplitude A = 0.12 m by a vertical wall

Figure 4 .Figure 5 .
Figure 4.The bathymetry in relation to the sponge layers (light blue regions), wavemaker and wave gauges locations (distances in metres)