Solitary wave solutions and their interactions for fully nonlinear water waves with surface tension in the generalized Serre equations

TL;DR

This paper investigates the effects of surface tension on fully nonlinear water waves, presenting new analytical solutions and numerical analysis of solitary wave interactions across different surface tension regimes.

Contribution

It introduces analytical peaked traveling wave solutions at critical surface tension and compares solitary wave interactions in subcritical and supercritical regimes.

Findings

Differences in solitary wave interactions across surface tension regimes

New peaked traveling wave solutions at critical surface tension

Numerical validation using high-accuracy finite element methods

Abstract

Some effects of surface tension on fully-nonlinear, long, surface water waves are studied by numerical means. The differences between various solitary waves and their interactions in subcritical and supercritical surface tension regimes are presented. Analytical expressions for new peaked travelling wave solutions are presented in the case of critical surface tension. The numerical experiments were performed using a high-accurate finite element method based on smooth cubic splines and the four-stage, classical, explicit Runge-Kutta method of order four.

Tables

Table 1. Table 1. Lax categories for the Euler , Serre and g Serre equations with B = 0.2 𝐵 0.2 B\ =\ 0.2\, .

Lax categories	(a)	(b)	(c)
Euler	$r⩽2.941$		$r>3.536$
Serre	$r⩽3.096$	$3.109⩽r⩽3.978$	$r⩾3.979$
gSerre	$r⩽3.453$	$3.469⩽r⩽5.129$	$r⩾5.130$

Equations

$$h_t+[hu{]}_x=0,$$

$$\displaystyle u_{\,t}\ +\ g\,h_{\,x}\ +\ u\,u_{\,x}\ =\ \frac{1}{3\,h}\;\bigl{[}\,h^{\,3}\,(u_{\,x\,t}\ +\ u\,u_{\,x\,x}\ -\ u_{\,x}^{\,2})\,\bigr{]}_{\,x}\ +\ \tau\,h_{\,x\,x\,x}\,.$$

$$h_t+\epsilon [hu{]}_x=0,$$

$$u_t+h_x+\epsilon u{u}_x-\frac{{\delta }^2}{3h}{\left[h^3(u_{xt}+\epsilon u{u}_{xx}-\epsilon u_x^2)\right]}_x-{\delta }^{2}B{h}_{xxx}=0,$$

$$\displaystyle\biggl{[}\,u\ -\ \frac{(h^{\,3}\,u_{\,x})_{\,x}}{3\,h}\,\biggr{]}_{\,t}\ +\ \biggl{[}\,\frac{u^{\,2}}{2}\ +\ g\,h\ -\ \frac{h^{\,2}\,u_{\,x}^{\,2}}{2}\ -\ \frac{u\,(h^{\,3}\,u_{\,x})_{\,x}}{3\,h}\ -\ \tau\,h_{\,x\,x}\,\biggr{]}_{\,x}\ =\ 0\,,$$

$$\displaystyle\biggl{[}\,h\,u\ -\ \frac{(h^{\,3}\,u_{\,x})_{\,x}}{3}\,\biggr{]}_{\,t}\ +\ \biggl{[}\,h\,u^{\,2}\ +\ \frac{g\,h^{\,2}}{2}\ -\ \frac{2\,h^{\,3}\,u_{\,x}^{\,2}}{3}$$

$$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\ \frac{h^{\,3}\,u\,u_{\,x\,x}}{3}\ -\ h^{\,2}\,h_{\,x}\,u\,u_{\,x}\ -\ \tau\mathcal{R}\,\biggr{]}_{x}=0\ ,$$

$$\displaystyle[\,h\,u\,]_{\,t}\ +\ \bigl{[}\,h\,u^{\,2}\ +\ {\textstyle{1\over 2}}\;g\,h^{\,2}\ +\ {\textstyle{1\over 3}}\;h^{\,2}\,\gamma\ -\ \tau\,\mathcal{R}\,\bigr{]}_{\,x}\ =\ 0\,,$$

$$\displaystyle\biggl{[}\,\frac{h\,u^{\,2}}{2}\ +\ \frac{h^{\,3}\,u_{\,x}^{\,2}}{6}\ +\ \frac{g\,h^{\,2}}{2}\ +\ \frac{\tau}{2}\;h_{\,x}^{\,2}\,\biggr{]}_{\,t}\ +\$$

$$\displaystyle\qquad\qquad\biggl{[}\,\Bigl{(}\frac{u^{\,2}}{2}\ +\ \frac{h^{\,2}\,u_{\,x}^{\,2}}{6}\ +\ g\,h\ +\ \frac{h\,\gamma}{3}\ -\ \tau\,h_{\,x\,x}\Bigr{)}\,h\,u\ +\ \tau\,h_{\,x}\,(h\,u)_{\,x}\,\biggr{]}_{\,x}\ =\ 0\,,$$

$$({\eta }^{\prime }{)}^2={\eta }^2\frac{c_s^2-1-\eta }{c_s^2/3-B(1+\eta )}.$$

$${\eta }^{\prime }=\pm \eta \sqrt{\frac{c_s^2-1-\eta }{c_s^2/3-B(1+\eta )}}.$$

$$-{\int }_{\eta (0)}^{\eta (x_i)}\frac{1}{\eta }\sqrt{\frac{c_s^2/3-B(1+\eta )}{c_s^2-1-\eta }}\mathrm{d}\eta ={\int }_0^{x_i}\mathrm{d}x.$$

$${\int }_{\log \eta (x_i)}^{\log \eta (0)}\sqrt{\frac{c_s^2/3-B(1+\exp (z))}{c_s^2-1-\exp (z)}}\mathrm{d}z-x_i=0,$$

$${\eta }_p(\xi )=A\exp (-\sqrt{3}|\xi |),$$

$$\eta (\xi ;B)=A\exp \left(\pm \int \sqrt{\frac{c_s^2-1-\eta }{c_s^2/3-B(1+\eta )}}\mathrm{d}\xi \right),$$

$$\underset{B\to 1/3}{\lim }\eta (\xi ;B)={\eta }_p(\xi )=A\exp (-\sqrt{3}|\xi |),$$

$$\eta (\xi ;B)\sim A\exp \left(-\sqrt{\frac{c_s^2-1}{c_s^2/3-B}}|\xi |\right),|\xi |\to \infty .$$

$$S\ =\ \{\phi\ \in\ \left.C^{\,2}_{\,\mathrm{per}}\,[\,a,\,b\,]\ \right|\ \phi\bigr{|}_{\,[x_{\,i},\,x_{\,i\,+\,1}\,]}\ \in\ \mathds{P}^{\,3}\,,\quad 0\ \leqslant\ i\ \leqslant\ N\ -\ 1\}\,,$$

$$C_{\mathrm{per}}^2=\{f\in C^2[a,b]|f^{(k)}(a)=f^{(k)}(b),0⩽k⩽2\},$$

$$({\tilde{h}}_t,\varphi )+((\tilde{h}\tilde{u}{)}_x,\varphi )=0,$$

$$\mathcal{B}({\tilde{u}}_t,\varphi ;\tilde{h})+(\tilde{h}({\tilde{h}}_x+\tilde{u}{\tilde{u}}_x),\varphi )+\frac{1}{3}\left({\tilde{h}}^3(\tilde{u}{\tilde{u}}_{xx}-({\tilde{u}}_x{)}^2),{\varphi }_x\right)+$$

$$B\left[(\tilde{h}{\tilde{h}}_{xx},{\varphi }_x)+({\tilde{h}}_x{\tilde{h}}_{xx},\varphi )\right]=0,$$

$$\mathcal{B}(\psi ,\varphi ;\tilde{h})=(\tilde{h}\psi ,\varphi )+\frac{1}{3}({\tilde{h}}^3{\psi }_x,{\varphi }_x),$$

$$\tilde{h}(x,0)=P{h}_0(x),\tilde{u}(x,0)=P{u}_0(x),$$

$$E_{\mathrm{amp}}=\frac{|H(x^{\ast }(t),t)-H_0|}{|H_0|},$$

$${\tilde{c}}_s=\frac{x^{\ast }(t)-x^{\ast }(t-\tau )}{\tau },$$

$$E_{\mathrm{speed}}=\frac{|{\tilde{c}}_s-c_s|}{|c_s|}.$$

$$E_{\mathrm{shape}}=\underset{\tau }{\min }\zeta (\tau ),\zeta (\tau )=\frac{|H(x,t^n)-h(x,\tau )|}{|h(x,0)|}.$$

$$E_{\mathrm{phase}}={\tau }^{\ast }-t^n.$$

$$E_{\mathrm{shape}}=\zeta ({\tau }^{\ast }).$$

$$\displaystyle E_{\,1}\,(t)\$$

$$\displaystyle E_{\,2}\,(t)\$$

$$\displaystyle E_{\,3}\,(t)\$$

$$\displaystyle E_{\,4}\,(t)\$$

$$({\partial }_t-c_s{\partial }_y)\eta +{\partial }_{y}u=0,$$

Full text

Key words and phrases: Serre equations; solitary waves; surface tension; peakons

MSC:

PACS:

∗ Corresponding author

Last modified:

Last modified:

1. 1 Introduction

2. 2 Numerical methods

3. 2.1 Numerical computation of solitary waves

4. 2.2 Numerical time integration of the gSerre equations

5. 2.3 Method validation

6. 3 Numerical experiments

7. 3.1 Head-on collisions

8. 3.2 Overtaking collisions

9. 3.3 Critical and transcritical cases

10. 4 Conclusions

11. Acknowledgments

Introduction

Waves propagating on the surface of a thin fluid layer can be influenced by the effects of surface tension. Theoretical and experimental studies have shown that interesting phenomena can emerge during the propagation of waves in the presence of surface tension on various fluids including the two extreme cases of water and liquid mercury, [42, 41, 24]. For example, in laboratory experiments with thin layers of liquid mercury such that its depth is small compared to the capillary length, elevation and depression solitary waves have been observed around the critical surface tension regime, [24]. For these reasons the mathematical modeling of waves in the presence of surface tension has attracted the interest of scientists that have derived several asymptotic model equations and explored the dynamics of their solutions, cf. [28, 1, 26, 27, 14, 18] and the references therein.

In this paper, we numerically study the effects of surface tension on fully-nonlinear shallow-water waves. Surface water waves are usually described by the full Euler equations of water wave theory, [48]. Due to their complexity and the difficulties arising in their theoretical and numerical study, simpler model equations have been derived as approximations to the Euler equations in the shallow water regime. There are three often studied regimes within shallow water waves: (i) the weakly nonlinear – weakly dispersive, (ii) the weakly nonlinear – fully dispersive and (iii) the fully nonlinear – weakly dispersive regime. Model equations such as weakly-nonlinear Boussinesq systems and Whitham type equations modelling capillary-gravity waves in the regimes (i) and (ii) were derived in [14, 19] (see also [15, 29]). These models extend the Boussinesq systems derived for surface waves with no surface tension in [6] and, although they incorporate surface tension effects, they are limited to small amplitude waves. Because of this approximation, some effects of surface tension cannot be observed due to the absence of higher-order nonlinear terms. For this reason, the study of higher-order models should be considered. Mathematical models appropriate for water waves with surface tension in the regime (iii) were derived in [18]. These equations extend the Serre equations [43, 44], incorporating surface tension effects, and are referred to as the generalized Serre (gSerre) equations, cf. [18]. For the derivation, justification and generalisations of the model equations for surface water waves in both regimes we refer to [18, 31]. In this paper, we focus on the study of solitary wave solutions of the gSerre equations.

We note that another higher order effect, dispersion, can play a fundamental role in the near critical surface tension regime, as in the fifth order Kawahara equation [30, 28] that describes weakly nonlinear, unidirectional shallow water waves. While higher order dispersion without surface tension for the fully nonlinear regime (ii) was recently presented in [36], its generalized, surface tension counterpart was only noted in [45] and has not been studied. Although additional model equations have been derived recently that incorporate more surface tension effects, [10, 11], we will restrict this study to the effects of surface tension on the solitary wave solutions of the gSerre equations.

The gSerre system of equations is a very accurate mathematical model of shallow water waves that is derived as an approximation to the Euler equations with surface tension and a flat bottom, [18, 31]. Here we consider a two-dimensional coordinate system $O\,x\,y$ with the horizontal axis coinciding with the still water level $y\ =\ 0\,$. A layer of a perfect, incompressible fluid (of constant density $\rho\ >\ 0$) is assumed to be bounded between a flat, impermeable bottom at $y\ =\ -\,d$ and its free surface $y\ =\ \eta\,(x,\,t)$ with air. (The variable $\eta$ measures the deviation of the free surface from the still water level.) The gSerre equations can be obtained by considering the depth-averages horizontal velocity of the fluid and approximating the pressure $p$ jump across the interface using the small slope approximation $\llbracket\,p\,\rrbracket\approx\ -\,\sigma\,\eta_{\,x\,x}\,$. Following the general lines of [12, 39], the gSerre equations for the depth averaged velocity $u\,(x,\,t)$ and the total depth $h\,(x,\,t)\ =\ d\ +\ \eta\,(x,\,t)$ can then be derived using variational methods (first derived in [18] along with other gSerre-type equations with different short-wave dispersion properties) in the form:

[TABLE]

We non-dimensionalise the surface height by the undisturbed water depth $h^{\,\prime}\ =\ h/d\,$, the deviation of the free-surface elevation by a typical wave amplitude $\eta^{\,\prime}\ =\ \eta/a$ the mean velocity by the long wave speed $u^{\,\prime}\ =\ d\,u\,/\,(a\,\sqrt{g\,d})\,$, the horizontal length by a typical wavelegth $x^{\,\prime}\ =\ x/\lambda\,$, and time by the amount of time a long wave takes to traverse a horizontal length $t^{\,\prime}\ =\ \sqrt{g\,d}\,t\,/\,\lambda\,$. Dropping the primes for brevity we obtain the non-dimensional equations

[TABLE]

where $\varepsilon\ =\ a\,/\,d\,$, $\delta\ =\ d\,/\,\lambda\,$. In this setting the total depth from the dimensionless and scaled bottom $z\ =\ -1$ is $h\ =\ 1\ +\ \varepsilon\,\eta$ while both variables $h$ and $u$ are functions of the spatial $x$ and temporal $t$ variables. The coefficient $\tau$ measures the ratio of gravity to capillary forces and is defined as $\tau\ =\ \sigma\,/\,\rho\,$, where $\sigma$ is the surface tension coefficient, $\rho$ the constant density of the fluid. We will utilize the dimensionless Bond number in the form $B\ =\ \tau\,/\,g\,d^{\,2}\ \geqslant\ 0\,$. Note that this definition of the Bond number, while common in water waves [18], is the reciprocal of the Bond number (also called the Weber or Eötvös number) in other areas [25]. When $B\ =\ 0$ (i.e. no surface tension is considered), then the gSerre equations (1.3), (1.4) reduce to the Serre equations. In the rest of this paper, we focus on the dynamical properties of the model equations where the scaling is not important and so we consider the gSerre equations in non-dimensional and unscaled form (1.3) – (1.4), i.e. we take $\varepsilon\ =\ \delta^{\,2}\ =\ 1\,$, which is the non-dimensional gSerre equations (1.3) – (1.4) with $a\ =\ \lambda\ =\ d\,$.

An important, structural property of the gSerre equations is their linear dispersion relation $\omega^{\,2}\ =\ k^{\,2}\,(1\ +\ B\,k^{\,2})\,/\,(1\ +\ k^{\,2}\,/\,3)\,$. For $B\ =\ 0\,$, the dispersion curvature $\omega^{\,\prime\prime}\,(k)$ is single-signed for positive wavenumbers $k\,$. The dispersion is convex (or concave). For $B\ >\ 0\,$, the dispersion relation is non-convex, exhibiting an inflection point at $k_{\,\ast}\ =\ \sqrt{3}\,\sqrt{1\ +\ \sqrt{(1\ +\ B)\,/\,B}}\,$. This loss of dispersion convexity can have important physical implications, for example on solitary waves, [28, 10] and undular bores [45]. This feature makes the gSerre equations particularly interesting for further study. When $B\ =\ 1/3\,$, the dispersion relation is degenerate $\omega^{\,2}\ =\ k^{\,2}$ and the gSerre equations are no longer dispersive, designating $B\ =\ 1/3$ as a critical value of surface tension.

Several conservation laws can also be derived for the gSerre system in a similar way as in the case of the Serre equations. In additions to (1.3), here we mention four of these conservation laws:

[TABLE]

where we introduced two quantities: \gamma\ =\ h\,\bigl{[}\,u_{\,x}^{\,2}\ -\ u_{\,x\,t}\ -\ u\,u_{\,x\,x}\,\bigr{]}\,, $\mathcal{R}\ =\ h\,h_{\,x\,x}\ -\ {\textstyle{1\over 2}}\;h_{\,x}^{\,2}\,$. The quantity $\gamma$ can be considered physically as the vertical acceleration of fluid particles computed at the free surface. It is noted that equation (1.5) represents the conservation of surface tangential momentum, equations (1.6) and (1.7) represent the conservation of horizontal momentum (and are equivalent), and equation (1.8) represents the conservation of energy.

It is known that the Serre equations admit solitary wave solutions of elevation type (i.e. of ${\mathrm{sech}}^{2}-$type), [4]. Capillary effects have been shown to be important for the shape of the solitary waves, [22, 18]. Specifically, the presence of surface tension makes an elevation solitary wave narrower than a solitary wave of the Serre equations with the same speed. Moreover, depending on the value of the Bond number $B\,$, the solitary waves can be either of elevation or of depression type. Depression solitary waves have negative excursion relative to the fluid background. The critical Bond number where the nature of the solitary wave changes, has been found to be $B\ =\ 1/3\,$, [18]. For $B\ <\ 1/3$ the gSerre equations admit elevation solitary waves while for $B\ >\ 1/3$ they admit depression solitary waves. It is noted that there are no known analytical formulas for solitary waves of the gSerre equations when $B\ >\ 0$ and for this reason, numerical computations are required. The change in solitary waves polarity with Bond number near the critical value $B\ =\ 1/3$ has been observed in the full water wave equations too, cf. e.g. [7, 16, 17]. As it was also suggested in [45, 7], we expect some interesting phenomena around the critical Bond number $B\ =\ 1/3\,$.

The Serre equations ($B\ =\ 0$) were first derived in [43, 44] and rederived various times since then, cf. e.g. [46]. Although there exist several studies, both theoretical and numerical, for the Serre equations without surface tension [38, 23], the behaviour of solitary waves under the effects of surface tension remain unknown [18]. This paper is focused on the properties of solitary wave solutions of the gSerre equations with surface tension and their interactions and, in general, how surface tension influences their dynamics. Special attention has been given to the critical value of the Bond number $B\ =\ 1/3\,$, where it is shown that the gSerre equations admit peaked solitary waves of elevation (peakons) and depression (antipeakons) type. The existence of these weakly singular traveling wave solutions along with their stability have been reported in [39] while here we provide a detailed description of their construction and their properties. It is noted that the gSerre equations are the first bi-directional propagation model in water wave theory that admits stable peakons and therefore we can study their head-on collision. The equations (1.3) – (1.4) are solved numerically in a periodic domain using an extension of the standard Galerkin/Finite element method of [40] to include the surface tension term. This numerical method has been proven very efficient with favourable convergence properties to smooth and non-smooth solutions [2, 40, 3].

The paper is organized as follows. In Section 2, we present the numerical methods used in this paper along with their numerical validation. Section 3 contains numerical results which analyse the various interactions between elevation and depression solitary waves for multiple values of the Bond number $B\,$. Finally, the main conclusions are discussed in Section 4.

Numerical methods

Consider the gSerre equations (1.3) – (1.4) in their dimensionless but unscaled form. Lacking analytical expressions for gSerre solitary waves, we present an efficient numerical method for their numerical generation. Then, the standard Galerkin/finite element method is presented for the numerical integration of the gSerre equations.

Numerical computation of solitary waves

It is known that the Serre equations possess solitary wave solutions traveling at constant speed $c_{\,s}$ of the form $h\,(x,\,t)\ =\ 1\ +\ \eta_{\,s}\,(x\ -\ c_{\,s}\,t)\,$, $u\,(x,\,t)\ =\ u_{\,s}\,(x\ -\ c_{\,s}\,t)$ with $\eta_{\,s}\,(\xi)\ =\ A\,{\mathrm{sech}}^{\,2}\,[\,\lambda\,\xi\,]\,$, $u_{\,s}\,(\xi)\ =\ 1\ -\ 1\,/\,(1\ +\ \eta_{\,s}\,(\xi))\,$, $\lambda\ =\ \sqrt{3\,A\,/(4\,(1\ +\ A))}\,$, $c_{\,s}\ =\ \sqrt{1\ +\ A}$ and $\xi\ =\ x\ -\ c_{\,s}\,t\,$.

On the other hand, it was shown in [18] that the gSerre equations possess classical solitary wave solutions (symmetric solitary waves that decay monotonically to zero, [21]) that satisfy the ordinary differential equation

[TABLE]

It is noted that the gSerre equations share the same relation between $u$ and $\eta\,$, i.e. $u\ =\ 1\ -\ 1\,/\,(1\ +\ \eta)$ with the Serre equations. Smooth solitary waves also share the same speed-amplitude relation $A\ =\ c_{\,s}^{\,2}\ -\ 1\,$, independent of the Bond number. This is verified by evaluating (2.1) at the point of the peak amplitude of the solitary wave. In what follows, we describe the numerical method for the computation of the symmetric traveling waves.

Without loss of generality, we search for solitary waves that are positive and symmetric about $x\ =\ 0$ for $0\ \leqslant\ B\ <\ 1/3\,$. We consider a large enough interval $x\ \in\ [\,-\,L,\,L\,]$ so that the solitary wave has decayed sufficiently close to the background value $\eta\ =\ 0$ within this interval. Under the assumption that $1\ +\ \eta\ \leqslant\ c_{\,s}^{\,2}\,\min\{1,\,1\,/\,(3\,B)\}$ or $1\ +\ \eta\ \geqslant\ c_{\,s}^{\,2}\,\max\{1,\,1\,/\,(3\,B)\}$ we can solve the equation (2.1) for $\eta^{\,\prime}$ and obtain the equation:

[TABLE]

Because of the symmetry of the solution, a nonuniform grid on $[\,0,\,L\,]$ is used, $0\ =\ x_{\,0}\ <\ x_{\,1}\ <\ \cdots\ <\ x_{\,N}\ =\ L\,$, where we assume that the function $\eta$ is decreasing on this grid. Integration of (2.2) yields

[TABLE]

Making the change of variable $\eta\ =\ \exp\,(z)\,$, then (2.3) can be simplified to the equation

[TABLE]

which defines the values $\eta_{\,i}\ =\ \eta\,(x_{\,i})$ implicitly, given $\log\,\eta\,(0)\ =\ \log A\,$, $c_{\,s}^{\,2}\ =\ 1\ +\ A\,$, for some $A\,$. Gauß–Legendre numerical quadrature is used for the approximation of the integral in (2.4) while the resulting nonlinear equations are solved with the secant method for values $\eta_{\,i}$ in the interval $(0,\,c_{\,s}^{\,2}\ -\ 1\,]\,$. Usually, the secant method converges in several iterations for a relative error tolerance of $\mathcal{O}\,(10^{\,-\,10})\,$. For the nodes $x_{\,i}\,$, we used the quadrature nodes of the Gauß–Legendre quadrature rule with $5$ nodes in a uniform grid of the computational domain. It is noted that the discretization of the inner products in the Finite Element method is based on Gauß–Legendre quadrature and therefore the numerically generated solitary waves can be used directly without using interpolation.

It is noted that in the case $B\ =\ 1/3$ the gSerre equations are degenerate with no linear dispersion. This implies that the asymptotic approximation of full water waves has broken down and higher order dispersion should be taken into account, [28]. In this case, we observe that the equation (2.2) can be simplified and one can verify that the solution

[TABLE]

satisfies (2.2) with $B\ =\ 1/3$ in a trivial way. Therefore, the critical gSerre equations with $B\ =\ 1/3$ possess peakons of elevation and depression type since the formula is valid for any $A\,$, [8]. Even if solutions of the form (2.5) satisfy the equation (2.2) for any value $A\ \in\ \mathds{R}$ it has been found numerically that only the peakons that satisfy the usual speed-amplitude relation $c_{\,s}\ =\ \sqrt{1\ +\ A}$ are stable, [39]. Such stable solutions are the limits of smooth solitary waves for $B\ \not=\ 1/3$ as $B\ \rightarrow\ 1/3\,$. This can be seen by making the dependence of the traveling-wave solutions on the parameter $B$ explicit and assuming solitary waves of the form $\eta\ =\ \eta\,(\xi;\,B)$ satisfying (2.2) with $c_{\,s}\ =\ \sqrt{1\ +\ A}\,$. Then, from (2.2) after integration we obtain:

[TABLE]

which is simplified to the peakon solution

[TABLE]

after taking the limit $B\ \rightarrow\ 1/3\,$. The existence of these weakly singular solutions has been demonstrated also using phase-plane analysis of the relevant ordinary differential equation in [39].

Remark 1**.**

Formula (2.6) can provide an asymptotic estimate for the decay rate of the solitary waves of the gSerre equations. Since solitary waves are assumed to decay at infinity, i.e. $\eta\ \ll\ 1$ as $\lvert\,\xi\,\rvert\ \to\ \infty\,$, we have from (2.6) that

[TABLE]

This exponential decay of the solitary waves is also illustrated in Figure 2.

Remark 2**.**

The asymptotic relation (2.7) can also been obtained from (2.1) written in the form $(\eta^{\,\prime})^{\,2}\ =\ F\,(\eta)$ and using the second order Taylor expansion of $F\,(\eta)\ \approx\ F^{\,\prime\prime}\,(0)\,\eta^{\,2}\,/\,2\ +\ \mathcal{O}\,(\eta^{\,3})$ where $F^{\,\prime\prime}\,(0)\ =\ 2\,(c_{\,s}^{\,2}\ -\ 1)\,/\,(c_{\,s}^{\,2}\,/\,3\ -\ B)\,$.

Numerical time integration of the gSerre equations

For the numerical approximation of the initial value problem of the gSerre equations subject to periodic boundary conditions, we implement a standard Galerkin / finite element method for the spatial discretization and the fourth-order, four-stage classical Runge–Kutta method for the discretization in time (cf. [40] for the Serre equations lacking surface tension). Consider the system (1.3) – (1.4) posed in a finite interval $[\,a,\,b\,]$ and for time $t\ \in\ [\,0,\,T\,]$ with some $T\ >\ 0\,$, and with periodic boundary conditions, i.e. $h^{\,(k)}\,(a,\,t)\ =\ h^{\,(k)}\,(b,\,t)$ and $u^{\,(k)}\,(a,\,t)\ =\ u^{\,(k)}\,(b,\,t)$ for $k\ =\ 0,\,1,\,2,\,\ldots$. We consider a uniform subdivision of the interval $[\,a,\,b\,]$ consisting of the nodes $x_{\,i}\ =\ a\ +\ i\,\Delta\,x\,$, where $i\ =\ 0,\,1,\,\cdots,\,N\,$, and $N\ \in\ \mathds{N}\,$, such that the grid size is defined as $\Delta\,x\ =\ (b\ -\ a)/N\,$. We shall consider numerical solutions of the gSerre equations in the space of cubic, periodic splines

[TABLE]

where

[TABLE]

and $\mathds{P}^{\,3}$ is the space of cubic polynomials.

The numerical solution will be denoted by $\tilde{h}$ and $\tilde{u}\,$. To state the spatial Galerkin semi-discretization, we first multiply equations (1.3) – (1.4) by $\phi\ \in\ S\,$. Integration by parts leads to the weak formulation:

[TABLE]

where $(f,\,g)\ =\ \int_{\,a}^{\,b}\;f\,g\;\mathrm{d}x\,$, the bilinear form $\mathcal{B}$ is defined for a fixed $\tilde{h}$ (and $\phi,\,\psi\ \in\ S$) as

[TABLE]

and the initial conditions are

[TABLE]

where $P$ is the $L^{\,2}$ projection onto $S$ defined by $(P\,v,\,\phi)\ =\ (v,\,\phi)\,$, for all $\phi\ \in\ S\,$.

Using the standard basis functions with B$-$splines for the space $S\,$, the equations (2.8) – (2.9) form a system of ordinary differential equations. It has been shown that the classical, explicit, four-stage, fourth-order Runge–Kutta method performs very well for the surface tensionless Serre equations, [40]. Denoting by $\Delta\,t$ the uniform time-step, we consider the temporal grid $t^{\,n}\ =\ n\,\Delta\,t\,$, for $n\ =\ 0,\,1,\,\cdots,\,K\,$, with $\Delta\,t\ =\ T/K\,$. For more information about the formulation and properties of this fully-discrete scheme, we refer the reader to [40, 2]. We only note here that the convergence rates for the present system are analogous to the convergence rates observed for the Serre system while in the case of peakons the convergence properties are similar to the case of the Camassa–Holm equation, [3]. Specifically, the convergence rate of the standard Galerkin method depends on the smoothness of the solutions and therefore the convergence rate for non-smooth peakon solutions can be reduced to order $1$ while for smooth solutions it is of order $4\,$. This is reflected in the conservation properties of the method and so it requires fine grids when dealing with non-smooth solutions. For the numerical computation of the integrals appearing in the numerical method, we use the Gauß–Legendre quadrature rule with five nodes in each spatial mesh interval. The nodes of the quadrature rule form a nonuniform grid which we used in the numerical method presented in Section 2.1.

Method validation

In order to validate the presented numerical methods for solitary wave computation and time integration, we generate and numerically evolve solitary wave solutions for various values of the Bond number $B\,$. As it was noted in [18], solitary waves for $B\ <\ 1/3$ are of elevation type while for $B\ >\ 1/3$ are of depression type, cf. Figure 2. Solitary waves are computed on the domain $[\,-40,\,40\,]\,$. Figure 2 presents a magnification of the numerically generated solitary waves for various values of the Bond number $B$ along with an exact solitary wave of the Serre system with $c_{\,s}\ =\ 1.5\,$.

The numerical method for the solution of (2.2) typically converges quickly and the convergence appears to depend slightly on the amplitude of the solitary waves and the value of the Bond number. As an indication of the speed of the numerical method described in Section 2.1, we present here the elapsed time for the computations of the solitary waves on the Gauß–Legendre quadrature nodes with $\Delta\,x\ =\ 0.1$ ($N\ =\ 800$) in $[\,-40,\,40\,]\,$. For $B\ =\ 0.1$ and $c_{\,s}\ =\ 1.5\,$, the CPU time is $10.4265$ seconds, while, for $B\ =\ 0.33$ and $c_{\,s}\ =\ 1.5$ is $11.1082$ seconds, very similar to the case $B\ =\ 0.35$ and $c_{\,s}\ =\ 0.5\,$, which is $11.5507$ seconds. As the Bond number approaches the value $B\ =\ 0$ the elapsed time is again close $11$ seconds. For example, for $B\ =\ 0.1$ and $c_{\,s}\ =\ 3$ the elapsed time is $10.9633\,$. In general, we observe that the convergence time of the solitary wave computation is comparable as $B$ and $A$ are varied.

In order to study the accuracy of the approximation of the solitary waves, we use the numerically generated solitary waves for $B\ =\ 0.1$ and $B\ =\ 0.5$ as initial conditions to the fully discrete numerical time integration scheme and study several error indicators relevant to the propagation of traveling waves. Specifically, we monitor the amplitude, speed, shape and phase errors for solitary wave propagation up to time $T\ =\ 100\,$. It is noted that for elevation solitary waves, we use the discretization parameters $\Delta\,x\ =\ 0.1$ and $\Delta\,t\ =\ 0.05\,$, while for depression solitary waves we use smaller mesh lengths $\Delta\,x\ =\ 0.01$ and $\Delta\,t\ =\ 0.005$ in order to produce stable and accurate computations. For the specific numerical method it has been shown that mild restrictions on $\Delta\,t\ \leqslant\ C\,\Delta\,x$ are adequate.

We define the normalized amplitude error as

[TABLE]

where $x^{\,\ast}\,(t)$ is the curve along which the computed approximate solution $H\,(x,\,t)$ achieves its maximum and $H_{\,0}\ \equiv\ H\,(0)$ is the initial peak amplitude of the numerically generated solitary wave. We observe that $E_{\,\mathrm{amp}}$ remains very small and practically constant during the propagation of the solitary waves, cf. Figure 3.

Additionally, we approximate the solitary wave speed $c_{\,s}$ by $\tilde{c}_{\,s}$ as

[TABLE]

where $\tau$ is a constant. The corresponding speed error is defined as

[TABLE]

The results for $\tau\ =\ 10$ in Figure 3 show that the error between the numerical values $\tilde{c}_{s}$ and the exact speed $c_{\,s}$ remained practically constant.

Two other error norms that are pertinent to solitary waves are the shape and phase errors. We define the normalized shape error as the distance in $L^{\,2}$ between the computed solution at time $t\ =\ t^{\,n}$ and the family of translated, exact solitary waves with the same amplitude/speed, i.e.,

[TABLE]

The minimum in (2.14) is attained at some critical $\tau\ =\ \tau^{\,\ast}\,(t^{\,n})\,$. This, in turn, is used to define the (signed) phase error as

[TABLE]

In order to find $\tau^{\,\ast}\,$, we solve the equation $\zeta^{\,\prime}\,(\tau)\ =\ 0$ using Newton’s method. The initial guess for Newton’s method is chosen as $\tau^{\,0}\ =\ t^{\,n}\ -\ \Delta t\,$. Having computed $\tau^{\,\ast}\,$, the shape error (2.14) is then

[TABLE]

These error norms are closely related to the orbit of the solitary wave (family of translated solitary waves) and measure properties of the solitary waves, which are often not well conserved using dissipative numerical methods.

The computed numerical errors are presented on a logarithmic scale in Figure 3. It can be observed that the errors in the propagation of both solitary waves are very small. Especially when we use grids with small $\Delta x\,$, the numerical speed of propagation of the solitary wave is almost equal to the exact speed and similarly the rest of the error indicators remain very small. This study shows that both numerical methods are very accurate and they conserve the properties of a propagating solitary wave very well. We note that the errors are analogous for solitary waves with different propagation speeds since they primarily depend on the discretization $\Delta x$ and $\Delta t\,$. We also mention that the fully-discrete scheme appears to be stable with no restrictive bounds on the ratio $\Delta t/\Delta x$ except when the solution is not very smooth. Mild restrictions, empirically found to be on the order $\Delta t/\Delta x\ \lessapprox\ 10^{\,-\,2}\,$, could be necessary for the numerical stability even of the non-smooth solutions.

As the Bond number $B$ approaches the value $1/3\,$, solitary wave solutions become more cusp-shaped, approaching the peakon solution (2.5), cf. Figure 2. For example, the solitary waves for $B\ =\ 0.33$ are very close to peaked solitary waves, [8, 34]. This phenomenon has already been explained in Section 2.1 and we verify it here numerically. As the solitary waves lose smoothness, then the mesh length $\Delta x$ must be reduced in order to maintain high resolution. We discuss further the transcritical case where $B\ \approx\ 1/3$ in Section 3.

In order to study the conservation properties of the numerical scheme, we consider the following conserved quantities:

[TABLE]

These quantities occur after integration of the conservation laws (1.5) – (1.8) over a periodic interval $[\,a,\,b\,]$ while taking into account the periodic boundary conditions. We also consider the evolution of initial conditions of the form $\eta\,(x,\,0)\ =\ C\,\mathrm{e}^{\,-\,0.1\,x^{\,2}}$ and $u\,(x,\,0)\ =\ 0$ in the interval $[\,-100,\,100\,]$ up to $T\ =\ 100$ for various values of the parameter $B\,$. In general the quantities $E_{\,i}\,(t)\,$, $i\ =\ 1,\,2,\,3$ are conserved in all of our simulations with at least $11$ digits correct for all values of $B$ and $C$ we used. For example, when we take $\Delta x\ =\ 0.1$ and $\Delta t\ =\ 0.01\,$, the quantities $\lvert\,E_{\,i}\,(t)\ -\ E_{\,i}\,(0)\,\rvert$ for $i\ =\ 1,\,2,\,3$ are always of $\mathcal{O}\,(10^{\,-\,10})\,$. The quantity $E_{\,4}$ is well conserved for small and large values of $B$ but is poorly conserved for values of $B$ close to the critical value $1/3$ due to the singularity in the first derivative of the solution. For example, for $B\ =\ 0.1$ and $C\ =\ 1$ the conserved value of $E_{\,4}$ is $214.9750606245\,$. For $B\ =\ 0.5$ and $C\ =\ -0.8$ the conserved value of $E_{\,4}$ is $195.4778205\,$. For $B\ =\ 1/3$ although the quantities $E_{\,i}$ for $i\ =\ 1,\,2,\,3$ are conserved again; the value $E_{\,4}$ is conserved only with $5$ digits and is $215.05\,$. Taking smaller values for $\Delta x$ and $\Delta t\,$, the conserved quantities can be improved. This behaviour is expected since the convergence of the Galerkin method depends on the regularity properties of the solution. The results indicate that the Galerkin method preserves the quantities $E_{\,i}$ for $i\ =\ 1,\,2,\,3$ while the presence of peaked solitary waves affect the conservation of the quantity $E_{\,4}\,$.

Numerical experiments

In this section, we study the effects of surface tension on various solitary wave interactions. Specifically, we study the head-on and overtaking collisions of elevation and depression solitary waves. We also study the generation and interaction of solitary waves when the Bond number $B\ \approx\ 1/3\,$. It is noted that in what follows, we will report the deviation of the free surface $\eta$ rather than the total depth $h\,$. Thus, we consider initial conditions that are less than $10^{\,-\,10}$ on the boundaries of the computational domain.

Head-on collisions

We first study the symmetric head-on collision for the gSerre equations with $B\ =\ 0.1$ for two identical elevation solitary waves that propagate in opposite directions. In these numerical experiments, we consider the interval $[\,-200,\,200\,]$ and take $\Delta x\ =\ 0.1\,$, $\Delta t\ =\ 0.01\,$. Here we present the solitary waves with $c_{\,s}\ =\ 1.2$ and amplitude $A\ =\ 0.4472\,$. The solitary waves are initially translated so that their maximum values are achieved at $x\ =\ -\,50$ and $x\ =\ 50\,$, respectively, exhibiting essentially no overlap in their exponentially small tails. The interaction begins at approximately $t\ =\ 40$ and the peak of the interaction occurs at about $t\ =\ 42\,$. The interaction is presented in Figure 4. We observe that after the collision, the solitary waves propagate in different directions followed by small amplitude dispersive tails.

In addition to the generation of dispersive tails, the inelastic interaction causes a phase change in the propagation of the solitary waves. During the collision there is a temporal interval in which the solution has only one peak while the maximum value of the solution during the interaction recorded is $0.9838\,$, which is greater than the sum of the amplitudes, $0.88\,$, of the solitary waves. The solitary waves after the interaction are separated and stabilized to different amplitudes $A\ \approx\ 0.4377$ compared with the initial amplitudes $A\ =\ 0.4472\,$.

Figure 5 shows the amplitude (a) and the location of the maximum values of the solution (b) recorded during the interaction. We observe that, due to nonlinear interaction, the amplitude fluctuates before it is stabilized to its new value. The dotted lines in these diagrams represent the predicted solitary wave maximum for no interaction.

We perform several symmetric head-on collisions and record the maximum value at $x\ =\ 0$ during the interaction. This value is also known as the maximum runup since, by reflection symmetry $h\,(x,\,t)\ =\ h\,(-x,\,t)\,$, $u\,(x,\,t)\ =\ -\,u\,(x,\,t)$ of the equations and the initial data, it coincides with the maximum runup of a solitary wave on a vertical wall located at $x\ =\ 0$ subject to appropriate boundary conditions [47, 38]. Figure 6 shows the values recorded for $B\ =\ 0.2$ and $0.3$ compared with the asymptotic solution of [38] for $B\ =\ 0\,$. We observe that surface tension decreases the maximum runup value. This effect is stronger for larger amplitude solitary waves for the same values of $B\,$. Moreover, the maximum runup value is decreasing as $B$ is increasing. Therefore, for larger values of $B\,$, we observe smaller maximum runup values.

We also study the unsymmetric collision of two solitary waves with amplitudes $A\ =\ 0.96$ and $0.44\,$, respectively (equivalently, speeds $1.4$ and $1.2\,$, respectively). The interaction is very similar to the unsymmetric interaction of the Serre equations and is not presented here. We only mention that the solitary waves after the interaction have amplitudes $A\ \approx\ 0.9515$ and $0.4325$ respectively. In conclusion, the head-on collision of elevation solitary waves with small surface tension is qualitatively the same as the collision of solitary waves with no surface tension. It is noted that a jet formation can be observed during the head-on collision of large amplitude solitary waves, [9, 40]. This has been observed in numerical experiments with the Euler equations while with the Serre equations can be observed a weak jet formation during the head-on collision of large amplitude solitary waves [40] but the exact jet formation is not possible for the Serre equations since the water height must be a single-valued function.

We continue with the head-on collision of depression solitary waves for Bond number greater than $1/3\,$. The interaction of depression solitary waves has never been studied before and there are no previous results to compare with. We study the symmetric head-on collision of two solitary waves with speed $c_{\,s}\ =\ 0.8$ and amplitude $A\ =\ -\,0.36$ (where we keep the minus sign to emphasise that the solution is negative) for Bond number $B\ =\ 0.5\,$. Initially, the solitary waves are translated to $x\ =\ -\,40$ and $40$ respectively and we solve the gSerre equations in the interval $[\,-100,\,100\,]$ with $\Delta x\ =\ 0.01$ and $\Delta t\ =\ 0.001\,$. Figure 7 shows the inelastic interaction between two depression solitary waves. It is interesting that the dispersive tails generated after the interaction propagate faster than the solitary waves and therefore lead in the propagation while there are no deviations of the free surface between the two new solitary pulses. This is due to the linear dispersion curvature sign change as $B$ crosses $1/3\,$.

The minimum peak amplitude of the solution $-\,0.66$ was recorded at about $t\ =\ 49.78\,$. This value is larger than the sum of the peak amplitudes of the two solitary waves. Because of the nonlinear interaction, a phase shift can be observed in both pulses, but here the resulting solitary waves have larger amplitude $A\ \approx\ -\,0.3598$ (less negative) and therefore propagate faster than the initial solitary waves. Also, the phase change is different than the subcritical case $B\ <\ 1/3\,$. In the supercritical case $B\ >\ 1/3\,$, the waves travel faster during the interaction and are separated earlier compared to the case where small or no surface tension is considered. The minimum of the solution as a function of time and the phase diagram with the location of the peak amplitudes of the solitary waves are presented in Figure 8. The effect of strong surface tension $B\ >\ 1/3$ on the head-on collision of solitary waves is to invert the dynamics of the collision relative to the weak surface tension case $B\ <\ 1/3\,$, resulting in faster solitary waves and dispersive tails propagating faster than the solitary waves.

The collision of two solitary waves of unequal size is very similar to the one with solitary waves of equal amplitudes. For example, we consider solitary waves with speeds $c_{\,s}\ =\ 0.6$ and $0.8$ corresponding to amplitudes $A\ =\ -\,0.64$ and $-\,0.36\,$, respectively. The results of the inelastic collision are presented in Figure 9. The dispersive tails are traveling faster than the solitary pulses again but now their shape is different, having deformed to solitary waves with different amplitudes post interaction: $A\ \approx\ -\,0.35974$ and $A\ \approx\ -\,0.63958\,$. The evolution of the minimum value of the solution is very similar with the one presented in Figure 8 and so it is omitted.

Overtaking collisions

We now consider a different type of interaction, the overtaking collision of two solitary waves traveling in the same direction but with different speeds. In contrast to the head-on collision case, the overtaking interaction is often referred to as the strong interaction of solitary waves due to the relative importance of nonlinearity [37]. For the Korteweg–de Vries (KdV) equation, a standard weakly nonlinear, long wave model of unidirectional waves, there are three categories of overtaking collisions classified by Lax for the two-soliton solution [32]. Each category corresponds to a distinct geometry of solitary wave interaction. These categories have been observed in experiment and computation of fully nonlinear water wave models (Euler and Serre equations) [13, 2] and in experiment and a model of viscous core-annular flows [35]. Here, we will consider the Lax categories within the context of the gSerre equations for various Bond numbers.

We label the three categories (a), (b), and (c). Lax category (a) corresponds to solitary waves of similar size whose interaction remains bimodal throughout, resulting in a small exchange of mass from the larger to the smaller solitary wave. In Lax category (c), the small solitary wave is absorbed completely by the large solitary wave, resulting in a symmetric, unimodal conglomerate at the peak of interaction. Following this, the smaller solitary wave is ejected behind the larger wave and each propagates independently. Lax category (b) is a combination of categories (a) and (c). In this case, the two solitary waves initially form an asymmetric, unimodal mass. At the peak of interaction, however, the conglomerate is bimodal. This process is undone and the smaller solitary wave is emitted behind the larger wave.

In what follows, we denote the amplitudes $A_{\,1}\,$, $A_{\,2}$ for the larger and smaller solitary waves, respectively. The amplitude ratio is denoted $r\ =\ A_{\,1}/A_{\,2}\ >\ 1\,$. For the KdV equation, interactions with $r\ \leqslant\ (3\ +\ \sqrt{5})/2$ are category (a), with $(3\ +\ \sqrt{5})/2\ <\ r\ <\ 3\ \approx\ 2.6180$ category (b) and with $r\ \geqslant\ 3$ category (c).

It is important to note that the Lax categories for the KdV equation are completely determined by the amplitude ratio $r$ of the two solitons. This is because the KdV equation admits Galilean and scaling invariances that enable one to fix the leading, slower soliton to have amplitude $1$ and trailing, faster soliton to amplitude $r\,$, both on a zero background. While the gSerre equations admit Galilean and scaling symmetries, the Lax categories for two soliton interactions functionally depend on both soliton amplitudes $A_{\,1}$ and $A_{\,2}$ separately. This is because the scaling symmetry is used to fix the total water height, here normalised to unity. In this section, we fix the amplitude of the faster gSerre soliton to unity $A_{\,1}\ =\ 1$ and vary the slower soliton’s amplitude $A_{\,2}\ <\ 1$ in order to identify the Lax categories in this restricted regime. Therefore, we consider the ratio $r\ =\ 1/A_{\,2}\,$, [13, 35].

In the first case with Bond number $B\ =\ 0.2$ similar types of interactions are observed. We use $3$ decimal digits in the calculation of the parameter $r$ and we observe that interactions with $r\ \leqslant\ 3.453$ are in the category (a), for $3.469\ \leqslant\ r\ \leqslant\ 5.129$ they are in the category (b) and for $r\ \geqslant\ 5.130$ they are in the category (c). For values of $r$ in the interval $[\,3.453,\,3.468\,]\,$, we observe a transition zone between categories (a) and (b) where the small solitary wave is absorbed and re-emitted only after the exchange of the masses towards the end of the interaction. This phenomenon has been observed also in the case of overtaking collisions of the Serre equations with $3.097\ \leqslant\ r\ \leqslant\ 3.108$ [2]. The limits for the three Lax categories for the Serre equations as reported in [2], along with the limits for the gSerre (with $B\ =\ 0.2$) and Euler equations, [13], are summarised in Table 1 (for $A_{\,1}\ =\ 1$). It is noted that the values presented in Table 1 are correct to the digits shown and are approximate values, so they can be used as an indication of where the transition occurs.

Phase diagrams for the different interactions of solitary waves in each of the Lax categories are presented in Figure 10. In Figure 10(a) it is observed that the two solitary waves maintain a distance while they exchange masses. It is also possible to observe the absorption of the small solitary wave in Figures 10(b) and (c). The dotted lines represent the paths of the solitary waves as if there were no interaction.

The interaction with $r\ =\ 2.5$ in the Lax category (a) is depicted in Figure 11, where it is shown that during the interaction there are two peaks. In Figure 12 the interaction with $r\ =\ 4$ in the Lax category (b) is presented where the small solitary wave is absorbed by the large solitary wave initially and then is re-emitted and two peaks are present during the interaction, while in the end is absorbed again and finally ejected and separated from the large solitary wave. Figure 13 shows the interaction with $r\ =\ 10\,$, which belongs to the Lax category (c) where during the interaction only one peak can be observed as the small solitary wave is absorbed by the large one until it is ejected and separated from the large one at the end of the interaction.

The maximum value of the solution as a function of time for several values of $r$ is presented in Figure 14. We observe that the maximum of the solution during the interaction does not behave monotonically with $r\,$. So we can achieve the same maximum for different values of $r\,$. For large values of $r$ and in the category (c) of Lax, the amplitude as a function of time is a smooth function. In the categories (a) and (b) the maximum appeared to have a singularity at the time $t$ where the minimum occurs.

Although the limits for the Lax categories of overtaking collisions are different when surface tension is included (cf., Table 1), the interactions are very similar and no new phenomena are observed when compared to the results reported in [38, 13, 33, 2, 40].

Like in the case of the Serre equations, small amplitude dispersive tails are generated during and after the interaction of two solitary waves. The dispersive tails are propagating mainly to the right but a small $N-$shaped wavelet is generated and propagates to the left as shown in Figure 15.

We draw the conclusion that the overtaking collision of two solitary waves for the gSerre equations with $B\ <\ 1/3$ retains the qualitative characteristics of the analogous interaction for the Serre equations. The situation for the gSerre equations with $B\ >\ 1/3$ is somewhat different, as we now demonstrate. We consider the gSerre equations with $B\ =\ 0.5$ and we test overtaking collisions for different solitary waves. In this case, we were able to observe the analogous three categories of Lax but here, the values of $r$ are totally different and cannot be compared with the analogous results obtained for small values of Bond number $B\,$. Specifically, we define here $r\ =\ \lvert\,a_{\,2}\,\rvert/\lvert\,a_{\,1}\,\rvert$ with $\lvert\,a_{\,1}\,\rvert\ >\ \lvert\,a_{\,2}\,\rvert$ i.e. the amplitude of the small solitary wave over the amplitude of the large solitary wave.

In order to compare numerical results with the predictions of Lax, we restrict our attention to solitary waves of small amplitude where the unidirectional solitary waves of the gSerre equations can be asymptotically approximated by KdV solitary waves. All the experiments were performed in $[\,-200,\,200\,]$ and the solitary waves are translated initially so as to attain their maximum values at $x\ =\ -\,50$ and $x\ =\ 50$ respectively. When $B\ <\ 1/3$ we take $\Delta x\ =\ 0.1$ and $\Delta t\ =\ 0.01$ while, for $B\ >\ 1/3$ we take $\Delta x\ =\ 0.02$ and $\Delta t\ =\ 0.01\,$.

In these experiments, the fast solitary wave has amplitude $a_{\,1}\ =\ -0.1\,$. Then, for a small solitary wave of amplitude $a_{\,2}\ =\ -\,0.2\,$, the interaction falls into Lax category (a). As can be observed in Figure 16, the solitary waves exchange masses and there are always two pulses present in the domain. On the other hand, when we take $a_{\,2}\ =\ -\,0.3\,$, the small solitary wave is initially absorbed by the large solitary wave, and then re-emitted, causing the existence of two local minima during the interaction, as described by Lax category (b). Finally, for $a_{\,2}\ =\ -\,0.4$ the small solitary wave is absorbed by the large one during the interaction and after the interaction is ejected and separated from the large one, Lax category (c).

Again, the minimum value of the solution as a function of time behaves similarly to the case of small Bond number but here the amplitude is the negative minimum of the solution. Figure 19 shows the evolution of the amplitude as a function of time for different values of $r\,$. In this figure, the amplitudes have been translated so as to be all at the same level. Again, we observe that as the interaction changes categories, the maximum amplitude is not a monotone function with respect to $r\,$.

The interaction again is inelastic and the generation of small amplitude dispersive tails is observed. Figure 20 shows the dispersive tails generated by the interaction of two solitary waves with $a_{\,1}\ =\ -\,0.1$ and $a_{\,2}\ =\ -\,0.3\,$. Again, a small amplitude oscillatory dispersive tail is generated in front of the two pulses that travels faster than the pulses. On the other hand, a very fast $N-$shaped wavelet is also generated moving to the left. The tails generated in the other cases where $B\ >\ 1/3$ are always very similar, thus we do not present them here.

Critical and transcritical cases

In this section, we study the critical case of Bond number $B\ =\ 1/3$ and the cases where $B$ is close to this critical value. As it was mentioned in Section 2, and in [18], there are no smooth traveling wave solutions known for the critical Bond number $B\ =\ 1/3$ but only the peaked solitary waves given by the formula (2.5). Additionally, the solitary waves corresponding to Bond numbers close to the critical value are very close to peaked solitary waves. Moreover, the absence of dispersive effects for the critical value $B\ =\ 1/3$ is expected to influence the behaviour of the solutions for values of $B$ close to $1/3$ is expected to be affected by. We first study the interactions of solitary waves in the cases where $B\ =\ 0.32\,,0.33\,,0.34$ and $0.35\,$.

We start with the description of the head-on collision of two equal solitary waves. For subcritical values of the Bond number, we consider solitary waves of amplitude $A\ =\ 1\,$, while we take $A\ =\ -\,0.3$ for the supercritical values. Although the two waves interact in a similar manner with the cases described in Section 3, the generated dispersive tails are not oscillatory; rather, they are similar to $N-$shaped waves. For example, when $B\ =\ 0.34\,$, the dispersive tail is ahead of the solitary wave, while in the case $B\ =\ 0.33$ the tails are behind the solitary waves. A close look at the dispersive tails shows an intermediate state between dispersive and non-dispersive waves. Figure 21 presents the results after the interaction of two equal solitary waves for $B\ =\ 0.32$ and $B\ =\ 0.34\,$.

The interactions of solitary waves in all transcritical cases we tested are very similar to the analogous subcritical and supercritical cases. The main difference observed between the two cases is the shape of the dispersive tails. Analogous behavior is observed for the overtaking collision in the transcritical cases. The overtaking collisions generated small $N-$shaped wavelets and dispersive tails, which consist of a few oscillations and are very similar to the tails generated after the head-on collision. For this reason, we do not present pictures of the overtaking collisions in the transcritical cases in this paper. We proceed with the critical case $B\ =\ 1/3\,$.

Although the critical case is numerically difficult in the sense that numerical errors might lead to false conclusions, we explore the existence of stable solitary waves when $B\ =\ 1/3\,$. It is known that localized initial conditions evolve into a series of solitary waves and dispersive tails in many nonlinear and dispersive wave models, [18, 5]. In order to examine if the gSerre equations possess stable solitary wave solutions when $B\ =\ 1/3\,$, we test the evolution of a general initial condition of the form $\eta\,(x,\,0)\ =\ a\,\mathrm{e}^{\,-\,b\,x^{\,2}}$ with zero initial velocity $u\,(x,\,0)\ =\ 0\,$. We present the results for $a\ =\ 1$ and $b\ =\ 0.1\,$, i.e. when we consider the evolution of a heap of water under gravity. The initial waveform is split into two symmetric waves that eventually evolve into a series of solitary waves. In Figures 22(a) and (b), we present the resolution of a Gaussian initial condition into a series of solitary waves where we set $\Delta x\ =\ 0.02$ and $\Delta t\ =\ 0.002\,$. Although the computation is performed in the interval $[\,-200,\,200\,]\,$, we present only the solution in the interval $[\,0,\,200\,]$ because it is symmetric. In Figure 22(b), we observe that the solitary waves coincide with peaked solitary waves when we compare the shape of the numerical solution with the analogous analytical peakon (2.5). This example serves as an indication of the existence of stable peaked solitary waves for the gSerre equations in the critical case $B\ =\ 1/3\,$.

The situation is similar when we take a negative initial condition given by the same formula with $a\ =\ -0.5\,$, $b\ =\ 0.1$ for the critical value $B\ =\ 1/3\,$. It is observed that although the initial condition is negative, it is resolved into a series of depression, peaked solitary waves (also known as antipeakons), indicating that the gSerre equations possess both stable depression and elevation peaked solitary waves when $B\ =\ 1/3\,$. In this experiment, we use $\Delta x\ =\ 10^{\,-3}$ and $\Delta t\ =\ 10^{\,-4}$ and the interval of integration is $[\,-200,\,200\,]\,$. Although further theoretical studies are required to ensure the accuracy of our conclusions, some confidence in the numerical results can be gained by the fact that both elevation and depression traveling waves appear to exist at the same time when $B\ =\ 1/3\,$, contrary to what we have experienced for large and small values of Bond number $B\ \neq\ 1/3\,$, where we are able to compute either elevation or depression solitary waves in each case.

In order to study the solitary waves generated by the evolution of a Gaussian and to ensure that they are bona fide peakons in this critical case, we isolated the solitary waves using the cleaning procedure suggested in [20, 5]. After one cleaning iteration, the solitary wave propagates without change in shape, amplitude or speed. Specifically, the elevation solitary wave of Figure 22(a) propagates with constant speed $c_{\,s}\ \approx\ 1.317$ and amplitude $A\ \approx\ 0.7337$ and satisfies the usual speed-amplitude relationship. The amplitude of the emerging depression wave is approximately $A\ \approx\ -0.3368\,$. The interaction of the cleaned solutions are studied and found to be very similar to the transcritical cases. The symmetric head-on collision of two elevation peakons is presented in Figure 23 where we observe similarities with the head-on collision of two solitary waves in the case $B\ =\ 0.32\,$. Similarly, the head-on collision of two depression peakons leads to very similar results when compared to the case $B\ =\ 0.34\,$. This suggests that at least some dynamics of the gSerre equations with $B$ near-critical are reasonably approximated by dynamics of the critical case.

We also explore the head-on collision between a peakon of elevation and a peakon of depression with amplitudes $A\ =\ 0.5$ and $A\ =\ -0.3$ (with $c_{\,s}\ =\ \sqrt{1\ +\ A}$ for both peakons), respectively. This collision appeared to be a combination of the two previous collisions. After the interaction, the depression peakon sheds a small wavelet in front of the pulse, leading its propagation while the elevation peakon sheds an analogous wavelet behind, following the propagation of the main pulse. Figure 24 shows the head-on collision between elevation and depression peakons.

Furthermore, we study the overtaking collision of two elevation and two depression peakons. Figure 26 shows the overtaking collision of two right-traveling elevation peakons of peak amplitudes $A\ =\ 2.5$ and $A\ =\ 1.1$ respectively. In this experiment, we take $\Delta x\ =\ 0.005$ and $\Delta t\ =\ 0.0005$ in the interval $[\,-100,\,100\,]\,$. The two peakons interact in an inelastic way and tails are generated during and after the interaction. It is noted that during the overtaking collisions of elevation peakons the peakons maintain some distance and exchange masses as in Lax category (a). Figure 27 shows the overtaking collision of two right-traveling depression peakons of peak amplitudes $A\ =\ -0.1$ and $A\ =\ -0.5$ respectively. In this experiment, we take $\Delta x\ =\ 0.0025$ and $\Delta t\ =\ 0.00025\,$. Once again, counter-propagating tails have develop during the overtaking collision of these depression peakons while two peaks are observed for almost the whole interaction as in the Lax category (b).

The last overtaking collision between an elevation and depression peakon is presented in Figure 25. In this experiment, we use the same peakons as in the previous experiment but in this case both of them propagate to the right. The interaction is strong, especially for the depression peakon which evolve into a much smaller peakon while two more peakons are generated after the interaction. Also, other oscillatory structures appear. The initial elevation peakon evolve into a new elevation peakon of smaller amplitude.

It is noted that the tails generated during the interactions of traveling wave solutions for $B\ \approx\ 1/3$ all include a $N-$wave propagating in one direction. This indicates that as the effects of dispersion diminish, the dispersive tails behave similar to the non-dispersive tales of the critical case.

Remark 3**.**

In order to explain the differences between the tails generated during the interactions discussed above, we study the properties of the gSerre equations linearised around the trivial solution and in a reference frame $y\ =\ x\ -\ c_{\,s}\,t$ moving with the speed of a solitary wave, which is traveling to the right with speed $c_{\,s}\ >\ 0\,$. The equations then become

[TABLE]

Plane wave solutions of the form $\eta\,(y,\,t)\ =\ a\,\mathrm{e}^{\,{\rm i}\,(k\,y\ -\ \omega\,(k)\,t)}\,$, $u\,(y,\,t)\ =\ b\,\mathrm{e}^{\,{\rm i}\,(k\,y\ -\ \omega\,(k)\,t)}$ of (3.1) – (3.2) satisfy

[TABLE]

The existence of nontrivial solution of (3.3) – (3.4) implies that

[TABLE]

where $v\,(k)\ =\ \omega\,(k)/k$ and $\phi\,(x)\ =\ (1\ +\ B\,x)\,/\,(1\ +\ x/3)\,$, $x\ \geqslant\ 0\,$. It can be seen that $\phi\,(x)\ \leqslant\ \phi\,(0)$ when $B\ <\ 1/3$ and $\phi\,(x)\ \geqslant\ \phi\,(0)$ when $B\ >\ 1/3\,$. Therefore, we find

[TABLE]

These relationships imply that in the case $B\ <\ 1/3\,$, where $c_{\,s}\ >\ 1\,$, the dispersive tails can propagate with speed $v\,(k)\ <\ 0$ and therefore follow the solitary waves, while in the case $B\ >\ 1/3\,$, where $c_{\,s}\ <\ 1\,$, the dispersive tails propagate with speed $v\,(k)\ >\ 0$ and therefore lead the solitary waves. In the critical case $B\ =\ 1/3$ there is no linear dispersion ($\phi\,(k^{\,2})\ =\ 1$) and therefore the tails are expected to propagate with speed $v\,(k)\ =\ -c_{\,s}\ \pm\ 1\,$. So if $c_{\,s}\ >\ 1$ then the generated tails should follow the traveling wave since $v\,(k)\ <\ 0$ but if $c_{\,s}\ <\ 1$ then tails can be generated in front and behind the main traveling wave. For example in Figure 23 where $c_{\,s}\ \approx\ 1.32\,$. Since $v\,(k)\ <\ 0$ the tails can propagate behind the solitary wave. Analogous indications can be extracted for the case where $c_{\,s}\ <\ 1$ for depression peakons.

Conclusions

Some effects of surface tension on gravity-capillary solitary waves of the gSerre equations were presented. Head-on and overtaking collisions were studied for subcritical, critical, and supercritical values of the Bond number. The qualitative dynamical picture of the interactions for values of the Bond number $B\ <\ 1/3$ appeared to be very similar to the analogous interactions of solitary waves of the Serre equations that neglect surface tension. The maximum runup was smaller when surface tension was taken into account, and the effect of surface tension was stronger for larger amplitude solitary waves. On the other hand, the primary observed differences for solitary wave interactions with large values of Bond number $B\ >\ 1/3$ is the relative propagation of resultant dispersive tails, which propagate faster than the solitary waves as opposed to the slower when $B\ <\ 1/3\,$. We also studied the existence of solitary waves in the critical case $B\ =\ 1/3$ where the numerical experiments indicated the existence of solitary waves of elevation and depression that are reminiscent of peaked solitary waves of the Camassa–Holm equation. Analytical formulas for peakon solutions of the gSerre equations were presented while their stability was explored. The numerical experiments related to the various nonlinear interactions of the traveling wave solutions of the gSerre equations indicate that the model at hand is not integrable for any value of the Bond number $B\,$. Because the gSerre equations are a fully nonlinear water wave model with weak or strong surface tension, the results obtained in this paper go beyond existing weakly nonlinear Boussinesq models and could be useful for understanding a variety of large-amplitude surface wave dynamics in thin fluid layers.

Acknowledgments

M. A. Hoefer was partially supported by NSF CAREER DMS-1255422. D. Mitsotakis was supported by the Marsden Fund administered by the Royal Society of New Zealand. The authors would also like to thank the anonymous referees for their valuable comments and suggestions that helped to improve the original manuscript.