On the influence of wave reflection on shoaling and breaking solitary waves

A coupled BBM system of equations is studied in the situation of water waves propagating over a decreasing fluid depth. A conservation equation for mass and also a wave breaking criterion, both valid in the Boussinesq approximation, are found. A Fourier collocation method coupled with a 4-stage Runge–Kutta time integration scheme is employed to approximate solutions of the BBM system. The mass conservation equation is used to quantify the role of reflection in the shoaling of solitary waves on a sloping bottom. Shoaling results based on an adiabatic approximation are analysed. Wave shoaling and the criterion of the breaking of solitary waves on a sloping bottom are studied. To validate the numerical model the simulation results are compared with reference results and a good agreement between them can be observed. Shoaling of solitary waves is calculated for two different types of mild slope model systems. Comparison with reference solutions shows that both of these models work well in their respective regimes of applicability.


INTRODUCTION
Model equations for free surface water waves propagating in a horizontal channel of uniform depth have been widely studied for many years.Boussinesq models incorporate the lowest-order effects of nonlinearity and frequency dispersion as corrections to the linear long wave equation.These models are widely used for describing the propagation of non-linear shallow water waves near coastal regions.In Boussinesq theory, it is important to assume that water is incompressible and inviscid and the flow is irrotational.There are two important parameters: the nonlinearity, the ratio of amplitude to depth, represented by α = a/h 0 , and the dispersion, the ratio of depth to wavelength, represented by β = h 2 0 /l 2 .As explained in detail in [5], the Boussinesq approximation is valid only when both α and β are small and have the same order of magnitude.
The more realistic situation of an uneven bottom profile is fundamental to studies of ocean wave dynamics in coastal regions.Several authors [15,19,21,28,30,40,43] have included the effect of smooth and slowly varying bottom topographies in both Boussinesq and shallow water theory.Wave shoaling is the effect by which surface waves propagating shorewards experience a decrease in the water depth.The study of shoaling waves is of importance in the nearshore areas and in the design of coastal structures.The 'classical' Boussinesq model was applied to shallow water of uneven bottom in two horizontal dimensions by Peregrine [36], who used depth-averaged velocity as a dependent variable and derived the system where ∇ = (∂ x , ∂ y ) T , η = η(x, y,t) represents the deviation of the free surface from its rest position at time t, u = u(x, y, z,t) denotes the horizontal velocity of the fluid at some height, ū denotes the depth-averaged velocity, and the bottom is z = −h(x, y).
Several improved Boussinesq-type models have been developed, starting with Madsen et al. [27], Nwogu [32], and Wei et al. [42], among others.Madsen et al. [27] achieved an improved linearized model by rearranging higher-order terms in the classical momentum equations, which are formally equivalent to zero within the accuracy of the model.Nwogu [32] demonstrated the flexibility obtained by using the velocity at an arbitrary depth as the velocity variable.Wei et al. [42] used Nwogu's approach to derive a fully nonlinear extension of Boussinesq equations, which further extended the range of validity of Boussinesq models without the weak nonlinearity restriction.It is worth mentioning that in [8,31] the Boussinesq model ( 1) is extended to the moving bottom topography, where the bottom topography depends on x, y, and t.In [31], a Benjamin-Bora-Mahony (BBM-BBM) type system (see [3]) is derived and solved numerically using a finite element method.One aspect in which the BBM system differs from Peregrine's Boussinesq' system is that it is amenable to numerical integration.Indeed, it is much easier to define a stable numerical approximation to a system of BBM type than to other Boussinesq systems, such as the Peregrine system.On the other hand, the Peregrine system features exact mass conservation while mass conservation in the BBM-BBM type systems is only approximate.Nevertheless, in the current work, we use a system of BBM type for numerical convenience.
The main contribution of the present paper is an in-depth study of wave reflection in a shoaling analysis based on Boussinesq systems such as (1).As part of our analysis, we formulate an approximate mass balance law associated with the Boussinesq scaling developed for flat bottoms in [1].We also extend the wave breaking criterion from [4] to the case of uneven beds.The mass balance equation is used in quantifying wave reflection due to the bottom slope, and the wave breaking criterion is used to determine an approximate termination point for the shoaling curves.A significant amount of literature has focused on the use of nonlinear shallow water equations to analyse long wave shoaling on a mildly sloping beach, and both experimental and numerical investigations have been carried out.However, reflection has not been quantified.
Many experimental studies, including the early studies [6,18], were aimed partly at comparison with classical shoaling laws such as the laws of Green and Boussinesq.However, most experimental work on wave shoaling has shown that actual shoaling curves vary considerably from the predictions of both Green's and Boussinesq's law.Grilli et al. [17] solved the full Euler equations by direct numerical integration, and this work compares their shoaling results with the numerical solution obtained in the present work.
Wave breaking is also important in studying nearshore area phenomena and tsunami propagation in coastal regions, because solitary waves are often used to model steep surface waves shoaling on beaches.An enormous literature also exists on breaking waves in a number of situations, including shoaling, wave breaking in open bodies of water, and breaking induced by a wave-maker (see [12,38], for instance).Chou and Quyang [9,10] and Chou et al. [11] discussed the criterion for the breaking of solitary waves on different slopes using the boundary element method to simulate the process of wave breaking.Using the fully nonlinear potential flow wave model, Grilli et al. [17] derived a criterion for wave breaking.In this paper, a different criterion of breaking solitary waves on a sloping bottom of a BBM-BBM type system is derived based on previous work in [4].Characteristics such as the breaking index, the wave height, the water depth, and the maximum particle velocity at the breaking point are studied and the breaking indices are compared with those obtained by Grilli et al. [17] and Chou et al. [11].The relation between breaking and reflection is investigated.
The present paper is organized as follows.In Section 2, the outline for the derivation of the coupled BBM-BBM type system [31] is given, and also the mass balance equations and the wave breaking criterion are derived.In Section 3, the coupled BBM-BBM type system is solved numerically using a Fourier collocation method coupled with a 4-stage Runge-Kutta time integration scheme and the convergence of the numerical scheme is validated.In Section 4, we demonstrate the effectiveness of the numerical method applied to our model system in simulations of solitary wave shoaling on a sloping bottom.Shoaling and wave breaking are studied numerically.This paper compares two models: the coupled BBM-BBM type system derived by Chen [8] and the one in Mitsotakis [31] with respect to the evolution of solitary waves.This comparison is concerned with initial wave profiles and wave shoaling on slopes that correspond to unidirectional propagation.In Section 5, the mass balance expressions are tabulated and the reflection of a small amplitude wave propagating over a slope is examined.Finally, a short conclusion is given in Section 6.

DERIVATION OF THE SYSTEM
The main model system to be used here belongs to the family of models derived in Mitsotakis [31].In order to obtain the Boussinesq system, the full water wave problem is used.A Cartesian coordinate system (x, z) is considered, with the x-axis along the still water level and the z-axis pointing vertically upwards.The fluid domain is bounded by the sea bed at z = −h(x) and the free surface z = η(x,t).Then the system of Euler equations for potential flow theory in the presence of a free surface is used.The derivation of the Boussinesq system is only briefly sketched.For a full derivation, the interested reader may consult [8] and [31].The variables are non-dimensionalized using the following scaling: where the tilde ( ˜) denotes non-dimensional variables, and h 0 , l, and a denote characteristic water depth, wave length, and wave amplitude, respectively.Consider a standard asymptotic expansion of the velocity potential ϕ and using the Laplace condition (△ϕ = 0, −h < z < η), write the velocity potential φ in the simplest form which is a series solution with only two unknown functions φ (0) and φ (1) .Then the velocity field can be expressed as where û and ŵ are the velocities at z = 0 and given by û = φ (0) x , ŵ = (1/β ) φ (1) .
In order to establish the relation between û and ŵ, use the bottom kinematic boundary condition (ϕ z + h x ϕ x = 0 at z = −h), which has the following form after substituting the above asymptotic expressions: Now inserting (4), (5), and (6) into free surface boundary conditions, one may derive the following Boussinesq system with variable bottom ηt It is emphasized that from the above system, and in terms of û, one can extend the system in terms of other velocity variables, such as the velocity at an arbitrary z location.In this work we use a trick due to [32].Namely, a new velocity variable ũθ defined at an arbitrary water level z = − h + θ (α η + h), with 0 ≤ θ ≤ 1.
Applying the standard techniques of inversion it is not difficult to derive the following expression as an asymptotic formula for û in terms of ũθ : Switching to the variable ũθ , the following expressions are obtained: Following the methodology in [5], for arbitrary µ, ν ∈ R and using ( 9), the following equations are derived Using equations ( 7)- (10) and appropriate expansions, the following system is derived: The parameters a, b, c and d are the same as in [5], where The coupled BBM-BBM type system is derived from (11) by selecting µ = 0 and ν = 0. Disregarding terms of order O(αβ , β 2 ) and dropping the superscript θ , the system takes the following form in dimensional variables Assuming the depth h is constant, the above system reduces to the original coupled BBM system in [5].

Mass balance
As mentioned in the introduction, the use of the BBM system necessitates the derivation of an approximate mass balance law.The following mass balance derivation is based on the work in [1], where mass balance theory is presented for the Boussinesq models with an even bottom profile.Since we are interested in varying bottom topography, we provide the following derivation.The integral form of the equation of mass conservation is d dt because there is no mass flux through the bottom or through the free surface.In non-dimensional variables the above relation becomes After integration with respect to z and use of asymptotic expansion of φ , we obtain Note that if we take the limit x2 → x1 , where x2 = x 2 /l and x1 = x 1 /l, then we obtain the balance equation (17 where The derivation could also be based on the differential form of the mass conservation, such as in [2].If we use the scalings M = ρh 0 M and q M = ρh o √ (gh 0 ) qM , then the dimensional forms of these quantities are Equation ( 17) is an approximate mass balance equation.The net mass transfer to or from a control volume during a time interval △t is equal to the net change (increase or decrease) in the total mass in the control volume during △t.In [1], they proved that the maximum error in the conservation of mass is smaller than O(αβ , β 2 ) in the case of an even bottom profile when a coupled BBM system is used.In Section (5) the amount of mass reflection will be computed for different cases.

Wave breaking in the BBM model system
As waves approach the shoreline the wave length and phase velocity decrease and the wave amplitude grows larger.The wave then crashes onto the shore because it becomes too steep for the bottom of the wave to carry it.The breaking of waves mostly depends on wave steepness and beach slope.As explained in [4], if the horizontal velocity near the crest of a wave exceeds the celerity of the wave, then the wave will break.
Let us denote the propagation speed by U and the horizontal velocity by u.The horizontal velocity u can be obtained from (5a) and ( 8): It is evident that once u θ (x,t) is known, (19) can be used to approximate the horizontal velocity at any depth.
After neglecting the second-order term, the dimensional form of the equation is given by Wave breaking occurs if Since the fluid domain depends on the surface profile, the value z = η is used to approximate velocities near the surface.It is clear that the solutions η(x,t) and u θ (x,t) of system (13) and propagation speed U are needed to find the breaking criterion.

NUMERICAL METHODS
System (13) has been solved numerically using a Fourier collocation method coupled with a 4-stage Runge-Kutta time integration scheme.For numerical computations, periodic boundary conditions on the domain [0, L] are used.The problem is translated to the interval [0, 2π] using the scaling u(λ x,t) = v(x,t), η(λ x,t) = ξ (x,t) and h(λ x) = h 1 (x), where λ = L 2π .Then the BBM-BBM system (13) becomes Consider the set of N evenly spaced grid points x j = 2π j N , j = 1, ..., N in the interval [0, 2π] referred to as collocation nodes.The spectral-collocation method is implemented in the physical space by seeking approximate solutions through a global periodic interpolation polynomial of the form ) and v N (x), ξ N (x) are interpolations of the function v(x) and ξ (x), respectively, i.e., v N (x j ) = v(x j ), ξ N (x j ) = ξ (x j ) (see [14,41]).Moreover, the corresponding Fourier collocation differentiation matrices D x and D xx are given by Then at the collocation points x = x j , the system becomes [ where I N is the unit N × N matrix and D N , D N are square matrices with the dimensions N × N following from (23a) and (23b), respectively, and diag(h 1 ), diag(h 2 1 ) are the diagonal matrices of h 1 and h 2 1 , respectively.This is a system of N ordinary differential equations for ξ N and also v N .The system is solved by using a fourth-order explicit Runge-Kutta scheme with time step △t.

Convergence study
It is important to verify the convergence of the numerical scheme.This is done following [37].A numerical method is convergent if the numerically computed solution approaches the exact solution as the step size approaches 0. To test the convergence of these numerical methods, the following discrete L 2 -norm is used and the corresponding relative L 2 -error is then defined to be where ξ N (x j ) is the approximated numerical solution and ξ (x j ) is the exact solution at a time T , for j = 1, 2, . . ., N.
Supposing the case of an even bottom, the coupled BBM system features solitary-wave solutions in a closed form if θ 2 = 7 9 (see [7]).Since the analysis of the solitary wave shoaling and breaking given here depend on the exact formula for the solitary wave, θ 2 = 7  9 is used in the present work.Then the exact solitary wave solutions of system of equations ( 13) take the forms and the constants W 0 , C 0 , and κ 0 are given by where h 0 is the undisturbed depth, H 0 is wave amplitude.
To check the convergence of these methods, we determine the L 2 -error each time for n steps and set the step size as △t = (t max − t min )/n for different n values n = 20, 40, 80, ... (Table 1) and different number of grid points N = 256, 512, 1024, ... (Table 2) in the case of an even bottom topography.A representative result for a wave of amplitude 0.5 m is given in Tables 1 and 2. The numerical scheme was implemented in MATLAB.In this calculation, the solution was approximated from T = 0 s to T = 5 s and the size of the domain was L = 100 m.In the computations shown in Table 1, N = 1024 Fourier modes were used.Table 1 shows fourth-order convergence of the Runge-Kutta method in terms of the time step △t.The fourth-order convergence of the scheme is apparent up to △t = 0.0039 s, when the error became dominated by the spatial discretization and the artificial periodicity.Table 2 shows the results of some computations aimed at validating the spatial convergence of the code.As expected, spectral convergence in terms of the number of spatial grid points N was achieved in these computations.Computations were also performed for other solitary waves with heights between 0.1 m and 0.6 m, and similar results were obtained for these cases.
Table 1.L 2 -error and convergence rate for Runge-Kutta method for different fixed step sizes in case of even bottom profile BBM-BBM type system.Here the convergence rate is the ratio of consecutive L 2 -errors Table 2. L 2 -error and convergence rate due to spatial discretization in case of even bottom profile BBM-BBM type system.Here the convergence rate is the ratio of consecutive L 2 -errors To indicate the significance of the improvement, Tables 3 and 4 show the results of computing approximate solutions of the inhomogeneous BBM-BBM type system where the functions η(x,t) = 0.3 cos(x −t) and u(x,t) = 0.  3 and Table 4, where the solutions were approximated from T = 0 s to T = 5 s.These tables show that the numerical implementation of a BBM-BBM type system with the periodic bottom function h(x) is correct.Similar results can be obtained for other 2π-periodic functions u, η, and h(x).

EVOLUTION OF SOLITARY WAVES ON A SLOPING BOTTOM
Shoaling of solitary waves with different wave heights for the initial undisturbed depth h 0 = 1 m to a smaller new depth up to h = 0.1 m is considered.The maximum wave heights were computed at different locations over the slope S = 1 : 35. Figure 1 shows results for a solitary wave of 0.6 m height.It shows that wave crests become steeper while shoaling on the slope.We generally see the reflection of a small amplitude wave when a solitary wave goes through a slope.After carefully measuring wave heights over the different slopes the relative maximum local wave height H/H 0 versus the relative local depth h 0 /h are plotted in Fig. 2, where h, h 0 , H, and H 0 represent the local water depth, the constant reference water depth, local solitary wave height, and initial solitary wave height, respectively.For later reference, we define the shoaling rate to be the exponent α if the relation ) α holds.The effect of a varying bottom on water waves of this class is of obvious engineering importance and numerical solutions were obtained by Peregrine [36] and Madsen and Mei [26] using a finite difference scheme to compute the deformation of a solitary wave climbing a beach.Experimental results for wave shoaling and breaking of solitary waves were obtained by Ippen and Kulin [18], Kishi and Saeki [23], Camfield and Street [6], and Synolakis [39].Note also that Pelinovsky and Talipova [34,35] studied the shoaling curves obtained by the wave height-wave energy relation for numerical solutions of the full water wave problem found by Longuet-Higgins [24] and Longuet-Higgins and Fenton [25].In the case of a periodic sequence of solitary waves, Ostrovsky and Pelinovsky [33] found that the shoaling relation reduces to a 'nonlinear' Green's law.The experimental results of Grilli et al. in [16] and numerical studies based on the potential flow theory for the Euler equations presented by Grilli et al. in [17] concentrate on shoaling studies.It is noteworthy that the studies of Grilli et al. [16,17] give a nice picture of different shoaling regimes and predict a variety of scaling relations for the local wave amplitude ahead and beyond the breaking point.
For comparison, we have considered numerical results of Grilli et al. [17].Figure 2 shows plots of data taken from [17]. Figure 2 shows that the shoaling curves of the current work are in good agreement with the numerical results of Grilli et al. [17].It can be seen that the shoaling rate increases initially more slowly than predicted by Green's law, but then increases as the water depth keeps decreasing.Although there is no breaking point in our numerical calculation, it can be noticed that the breaking points appeared in the results obtained by Grilli et al. [17].For instance, ( 21) is used to check the breaking criterion as discussed above.
It can be seen from Fig. 2 that for waves on a mild slope (1 : 100) the water depth at breaking, h b /h 0 , will be larger than that on a steep slope (1 : 35).Furthermore, under the same wave conditions, the amplitude at breaking points H b /h 0 is larger for a mild slope than for a steep slope.In particular, the agreement of the breaking criterion (21) with the results of Grilli et al. [17] is much better on a mild slope.
Figure 3 shows plots of shoaling rates for a wave of initial amplitude 0.1 m with slopes 1 : 100, 1 : 400, and 1 : 800.It is apparent that for the slope 1 : 100, the shoaling rate is lower than predicted by Green's law for small h 0 /h and higher for large h 0 /h.However, for the smaller slopes 1 : 400 and 1 : 800, the shoaling Note that the bottom topography is smoothed near the corners.Here G denotes Green's law, B denotes Boussinesq's law, the dotted curves are our numerical results, and the solid curves are numerical results from Grilli et al. [17].Rectangular and circular symbols denote the breaking points of Grilli et al. [17] and the present work, respectively.
rate is closer to the line h −1 for large h 0 /h.Apparently, the computed curves get close to those predicted by Boussinesq's law for smaller slopes.Now the breaking criterion ( 21) is applied to the solitary wave solutions.In order to find wave breaking in these solitary wave solutions, the x-location of the maximum wave height at each time step is found.The propagation speed U is then estimated using these x-locations at each time step.Finally, if the computed horizontal velocity u exceeds the mean propagation speed U, we can conclude that around this time step the wave is starting to break.
The water depth at breaking measured under the wave crest is denoted as h b and the corresponding solitary wave height at breaking is denoted as H b .In Table 5 the relative breaking wave height H b /h b at the corresponding breaking points is compared with those of Grilli et al. [17] and Chou et al. [11].For a large wave amplitude the wave height will exceed the breaking criterion very soon after it propagates on the slope and so wave breaking occurs almost instantly without too much change in height.The ratio of relative breaking wave height is larger for small amplitude waves than for large amplitude waves.Wave breaking occurs sooner for larger initial waves.
McCowan [29] theoretically defined the breaker depth index as H b /h b = 0.78 for a solitary wave travelling over a horizontal bottom using the assumption that instability is reached when the particle velocity at the crest equals the wave celerity and that the crest angle is then 120 • .To estimate the initial breaking wave height on a mild-slope beach, this value (H b /h b = 0.78) is most commonly used in engineering practice as a first estimate.Ippen and Kulin [18] showed that the upper limit of the breaking criterion should be 0.78 for a solitary wave over a very mild slope.In this article the slope 1 : 35 is used and it was found that the relative breaking wave heights H b /h b are smaller for higher amplitude waves (Table 5).It can be seen that the relative breaking wave heights H b /h b at breaking points are well above the McCowan limit 0.78.Since the relative breaking wave heights H b /h b at breaking points are smaller than those obtained by Grilli et al. [17] and Chou et al. [11], we might consider a higher order Boussinesq model for further study.

Comparison to mild slope model systems
For comparison, the work of Chen [8] is considered.Chen presented equations for bi-directional waves over an uneven bottom, which may be written in non-dimensional, unscaled variables and disregarding terms of order O(αβ , β 2 ) as ) The models of Chen [8] and Mitsotakis [31] represent the same type of coupled BBM-BBM type system, derived in the context of the Boussinesq scaling.One can derive a number of special cases of the general Boussinesq system.Since we are interested in coupled BBM-BBM type system, the model of Chen was chosen for comparisons.The above system (28) is solved using the same numerical technique as above.
The main difference between the two systems ( 13) and ( 28) is approximation of bottom motion.In (13), the bottom motion is non-dimensionalized by h = h h 0 , and in (28), it is non-dimensionalized by h = h−h 0 a 0 , which is similar to the approximation of wave amplitude η. Figure 4 shows computations for the shoaling curves with the initial amplitude 0.4 m.It can be noticed that the shoaling curve corresponding to system (28) lies below that predicted by Green's law because of the lower order approximation.
In [8], the bottom function h(x) is assumed to be O(α) and in [31], the bottom function h(x) is assumed to be O (1).The results are in line with the assumptions used in their respective derivations.
In Fig. 5, the shoaling curves close to a breaking point of a solitary wave with an initial wave height 0.2 m are snown.A comparison between the present shoaling result using system (13) derived by Mitsotakis [31] and the numerical results of the Nwogu system [32] presented in [13] and also the numerical results of system (28) derived by Chen [8] is shown in the left panel of Fig. 5.It can be seen that the Nwogu system, which is derived in terms of amplitude-velocity, quickly over-shoals.The dashed-dotted curve represents numerical results for system (13) derived by Mitsotakis [31], the solid curve shows the numerical results from Grilli et al. [17], and the dashed curve represents numerical results for system (28).Indeed, system (28) works for small-amplitude bottom variations as expected, since the bottom function h(x) is assumed to be of order O(α).
x/h 0  In the left panel, the dashed-dotted curve represents numerical results for system (13) derived by Mitsotakis [31], the solid curve shows the numerical results from the Nwogu system [32], the circular symbols denote the data of the laboratory experiments of [16], and the dashed curve represents numerical results for system (28) derived by Chen [8].Here x is the horizontal coordinate, z is the vertical coordinate, h 0 is the constant reference water depth, and H/H 0 is the relative maximum local wave height.

MASS CONSERVATION ON A SLOPING BOTTOM
The effect of depth variations on solitary waves of shallow water wave theory is examined.The BBM-BBM type system ( 13) is simulated.In all the numerical results of this section we use N = 1024, θ 2 = 7 9 .Mass conservation is used to quantify the role of reflection in the shoaling of solitary waves.Note that the piecewise smooth linear bottom topography is used.To avoid the generation of small spurious oscillations due to the discontinuity in the derivative of the bottom function, it is smoothed near the singular points.
Consider a control volume delimited by the interval [50 m, 150 m] on the x-axis.The mass per unit width contained in this interval is defined by ∫ 150 50 M(x,t) dx and the mass flux through the boundaries of the control volume is defined by q M (50,t) and q M (150,t), where M and q M are given in (18).The quantities M and q M during the passage of a solitary wave are computed.It is observed that the mass outflux is approximately equal to the addition of mass influx and the reflection of the mass.As can be seen in Fig. 6, the mass reflection has negative values.
In Table 6, the results for various amplitudes of the solitary wave are displayed for ∆h = 0.3 m on a slope 1 : 35.For the height of the topography ∆h = 0.3 m, the mass influx through the initial boundary of the control volume is defined by 'Mass outflux = ∫ 15 0 q m (50,t) dt', the mass outflux through the final boundary of the control volume is defined by 'Mass outflux = ∫ 60 15 q m (150,t) dt', and the mass reflection through the initial boundary of the control volume is defined by 'Mass reflection = ∫ 60 15 q m (50,t) dt'.Note that the time limit may vary for other ∆h's.The error is defined by 'error = mass outflux-mass reflection-mass influx'.Table 6 suggests that mass conservation has a negligible error and that the error tends to 0 as α = a/h 0 approaches 0.
In Table 7, the results for various ∆h of water level are displayed with initial amplitude a = 0.3 m on a slope 1 : 35.It is clear from Table 6 and Table 7, that the mass conservation holds approximately for the coupled BBM system and the ratio between mass reflection and mass influx is called 'mass ratio', which is smaller for larger amplitude waves and smaller ∆h.
The reflection of a small amplitude wave when a solitary wave goes through a slope is defined as 'reflection'.To find the ratio between reflection and initial solitary wave, we use the following L 2 -norm: To calculate the L 2 -norm of initial solitary waves, the value of η is integrated with respect to x on the fluid domain [0, L] at initial time t = 0. To determine the L 2 -norm of a reflected wave, we run the solitary wave on the slope for long enough time to separate the reflection of the small wave from that of the solitary wave.The end point of the reflected waves on the x-axis is denoted by x r (see Fig. 7).Then the reflected wave is integrated on the interval [0, x r ].The corresponding 'reflection coefficient' is then defined to be .
x, m Table 8.Calculation of the amount of 'reflected' waves for different slopes and amplitudes.It shows that the 'reflection coefficient' approaches zero as the slope becomes more and more gentle.Here ∆h is the height of the topography and H 0 is the initial wave amplitude The 'reflection coefficient' approaches zero as the slope becomes more and more gentle (Table 8).Moreover, the reflection coefficient for the steep slope (1 : 35) is approximately twice the value of that of the mild slope (1 : 100).For steeper slopes the reflection coefficient is large because the wave height at breaking points is smaller for a steep slope than for a mild slope.

CONCLUSION
In this article, a coupled BBM system of equations is studied in the situation of water waves propagating over a decreasing fluid depth.A conservation equation for mass and also a wave breaking criterion, both valid in the Boussinesq approximation, were found.A Fourier collocation method coupled with a 4-stage Runge-Kutta time integration scheme was employed in this work to approximate the solution of the BBM system.It is shown that the approximate mass conservation relation is reasonably accurate.Moreover, the results from the evaluation of the approximate mass conservation law show that the ratio of mass reflection to mass influx approaches zero as the difference in flow depths (∆h) becomes small.In our previous paper [20] we showed that for waves of very small amplitude, the shoaling relation approaches Boussinesq's law for Boussinesq-type systems that are valid for waves with the Stokes number S = α/β of order 1, and in this case we measured the transition of the wave only at the initial and final stage assuming the wave undergoes an adiabatic adjustment.It is confirmed (Table 8) that the L 2 -ratio between reflection and initial solitary wave approaches zero as the slope becomes more and more gentle.This lends additional credibility to shoaling results based on adiabatic approximation.In addition, the results displayed in Fig. 3 indicate that shoaling rates for small amplitude waves are closer to those predicted by Boussinesq's law for very gentle slopes.
Considering the shoaling of finite amplitude waves, we compared shoaling curves obtained with the current method to numerical results of Grilli et al. [17] for the Euler equations based on potential flow theory.The experimental results of Grilli et al. [16], and the corresponding shoaling curve of the current work were in good agreement with the numerical results of Grilli et al. [16,17].It was found that the variation in wave height of a shoaling solitary wave initially increased at a lower rate than predicted by Green's law, but then increased similar to Boussinesq's law.Indeed, the shoaling curves achieved in this paper match the shoaling curves of Grilli et al. [17] better than the similar approximation established by Khorsand and Kalisch [22].
The comparison of shoaling curves of two model systems, (13) ( [31]) and ( 28) ([8]), with the numerical results of Grilli et al. [17] showed that each of these models works well in their respective regimes of applicability.The agreement of the breaking criterion (21) with the results of Grilli et al. [17] is much better on a mild slope.
3 sin(x −t) are used as the exact solutions and the bottom h(x) = 0.5 − (0.1) cos(x) is assumed.Then the relative L 2 -error for various pairs of combinations between the time steps ∆t = 0.1/2n for n = 1, 2, 3, ...; and N = m × 64 for m = 1, 2, 3, ... is calculated.The results are shown in Table

1 ∆ h=0. 7 Fig. 1 .
Fig. 1.Transformation of a solitary wave of initial amplitude 0.6 m on the slope S = 1 : 35.Here ∆h is the height of the topography.Note that the bottom topography is smoothed near the corners.

Fig. 2 .
Fig. 2. Computations for the shoaling curves with initial amplitudes 0.6 m (upper panel), 0.4 m (middle panel), 0.2 m (lower panel) on slopes 1 : 35 (left panel) and 1 : 100 (right panel).The relative maximum local wave height H/H 0 versus the relative local depth h 0 /h are plotted.Here G denotes Green's law, B denotes Boussinesq's law, the dotted curves are our numerical results, and the solid curves are numerical results from Grilli et al.[17].Rectangular and circular symbols denote the breaking points of Grilli et al.[17] and the present work, respectively.

Fig. 3 .
Fig. 3. Computations for the shoaling curves with the initial amplitude H 0 = 0.1 m on different slopes.The relative maximum local wave height H/H 0 versus the relative local depth h 0 /h is plotted.

Table 5 .
Comparison of the relative breaking wave height H b /h b and the wave height at breaking points H b /h 0 for waves with initial amplitudes 0.2 m, 0.25 m, 0.3 m, and 0.4 m on slope 1 : 35 from Chou et al. [11], Grilli et al. [17], and the present work

Fig. 4 .
Fig.4.Computations for the relative maximum local wave height H/H 0 versus the relative local depth h 0 /h plotted with initial amplitude 0.4 m on a slope 1 : 35.The dashed-dotted curve represents numerical results for system (13) derived by Mitsotakis[31], the solid curve shows the numerical results from Grilli et al.[17], and the dashed curve represents numerical results for system(28).Indeed, system(28) works for small-amplitude bottom variations as expected, since the bottom function h(x) is assumed to be of order O(α).

Fig. 5 .
Fig.5.Computations for the shoaling curves with the initial amplitude 0.2 m on a slope 1 : 35 close to the breaking point.Right panel: computational setup.Left panel: close-up of shoaling curves near breaking points.In the left panel, the dashed-dotted curve represents numerical results for system (13) derived by Mitsotakis[31], the solid curve shows the numerical results from the Nwogu system[32], the circular symbols denote the data of the laboratory experiments of[16], and the dashed curve represents numerical results for system (28) derived by Chen[8].Here x is the horizontal coordinate, z is the vertical coordinate, h 0 is the constant reference water depth, and H/H 0 is the relative maximum local wave height.

Fig. 6 .Table 6 .
Fig. 6.The left panel shows a solitary wave solution for system (13) with the initial amplitude 0.3 m at time t = 60 s.The right panel shows plots of time series of the mass influx at x = 50 m (solid curve), the mass reflection at x = 50 m (dashed curve), and mass outflux at x = 150 m (dash-dotted curve), per unit span.The results are shown in the numerical domain.

Fig. 7 .
Fig. 7. Reflection of the transformation of solitary waves of initial amplitude 0.2 m on a slope S = 1 : 35 and the height of the topography ∆h = 0.3 m in the physical domain.

Table 3 .
Inhomogeneous BBM-BBM type system (27); L 2 -error and convergence rate due to temporal discretization.Here the convergence rate is the ratio of consecutive L 2 -errors

Table 4 .
Inhomogeneous BBM-BBM type system (27); L 2 -error and convergence rate due to spatial discretization.Here the convergence rate is the ratio of consecutive L 2 -errors