Solitary water waves created by variations in bathymetry

Abstract

We study the flow of water waves over bathymetry that varies periodically along one direction. We derive a linearized, homogenized model and show that the periodic bathymetry induces an effective dispersion, distinct from the dispersion inherently present in water waves. We relate this dispersion to the well-known effective dispersion introduced by changes in the bathymetry in non-rectangular channels. Numerical simulations using the (non-dispersive) shallow water equations reveal that a balance between this effective dispersion and nonlinearity can create solitary waves. We derive a Korteweg–de Vries-type equation that approximates the behaviour of these waves in the weakly nonlinear regime. We show that, depending on geometry, dispersion due to bathymetry can be much stronger than traditional water wave dispersion and can prevent wave breaking in strongly nonlinear regimes. Computational experiments using depth-averaged water wave models confirm the analysis and suggest that experimental observation of these solitary waves is possible.

1. Introduction

1.1. Wave propagation in media with a periodic structure

The importance of laminates and composite materials in engineering led to the study of elastic waves in periodically varying media. Long wavelength linear elastic waves experience an effective dispersion that arises due to the periodic variation in the material coefficients (Sun, Achenbach & Herrmann 1968). Similar conclusions regarding general linear wave propagation were reached by using Bloch expansions in Santosa & Symes (1991) and by using homogenization theory in Fish & Chen (2001). In one dimension, this effective dispersion is the result of reflection, and its strength is correlated with the degree of variation of the impedance in the medium.

For nonlinear elastic waves, a similar effect has been studied in LeVeque & Yong (2003). The induced dispersion due to reflections leads to the formation of solitary waves that behave similarly to those arising in nonlinear dispersive wave equations like the Korteweg–de Vries (KdV) equation (Korteweg & De Vries 1895). This behaviour also depends on the degree of variation in the impedance; if it is not strong enough then there is little effective dispersion and shocks tend to develop, as they would in a homogeneous medium (Ketcheson & LeVeque 2012; Ketcheson & Quezada de Luna 2020). In the multi-dimensional setting, effective dispersion arises not only from variation in the impedance, but also from variation in the linearized sound speed (Quezada de Luna & Ketcheson 2014). This latter source of effective dispersion can also lead to solitary wave formation for nonlinear waves (Ketcheson & Quezada de Luna 2015). These waves have been called diffractons since they appear as a consequence of diffraction due to changes in the sound speed.

The first objective of the present work is to investigate similar effects for water waves. Periodic variation of the medium is introduced through bathymetry that varies periodically in one direction. In § 2 we derive effective homogenized equations for small-amplitude waves in this setting. These equations show that such waves experience an effective dispersion, similar in nature to the effects discussed above for elastic waves. We refer to this as bathymetric dispersion. These effective equations describe waves varying in two horizontal dimensions; if restricted to plane waves propagating transverse to the variation in bathymetry (as depicted in figure 1a) they are similar to those derived in Chassagne et al. (2019). The amount of dispersion increases with the variation in the bathymetry, which is also correlated with variation in the linearized impedance and sound speed for water waves.

Figure 1. An example of the kind of bathymetry studied in this work. (a) Bathymetry that changes periodically in one direction. (b) Narrow channel with non-flat bathymetry.

In § 3, we perform numerical simulations of the shallow water equations with periodic bathymetry and obtain solitary waves. Our analysis and qualitative results apply to general bathymetric profiles that are periodic in one direction. For concrete illustration, in most of the numerical examples we consider piecewise-constant bathymetry, as depicted in figure 1(a).

We refer to these solitary waves as bathymetric solitary waves, since they appear only in the presence of varying bathymetry. We compare some properties of these solitary waves with those of the diffractons observed in Ketcheson & Quezada de Luna (2015) and with KdV solitons. In § 4, based on the linear homogenized system, we propose a KdV-type equation, which we solve numerically and compare vs direct simulations of the shallow water equations. In § 5 we compare our KdV model with other dispersive models. Finally, we make some concluding remarks in § 6.

1.2. Solitary water waves in narrow non-rectangular channels

Dispersion of water waves due to bathymetry has been analysed and observed in a very different setting: flow in a channel like that depicted in figure 1(b) or 2. If the bottom of the channel is not flat, small-amplitude waves can be shown to obey an effective equation that includes an additional dispersion, distinct from the inherent dispersion of water waves and depending purely on the channel geometry (Peters 1966; Peregrine 1968; Shen 1969; Teng & Wu 1997; Chassagne et al. 2019). Just as in the presence of periodic bathymetry, here, the effective dispersion also increases with the amount of variation in the bathymetry. In fact, due to symmetry of the solution over periodic bathymetry, the two settings are equivalent. For example, consider the periodic bathymetry in figure 1(a) and the non-rectangular channel in figure 1(b). The solution in the periodic channel restricted to the domain extending from the middle of segment A to the middle of segment B is identical to the solution in the non-rectangular channel (with slip boundary conditions at the sides of the channel). Flow in channels with other shapes, such as that shown in figure 2, can also be viewed equivalently as flow over an infinite periodic bathymetry, as long as the channel does not have sloping sides that extend above the water surface. Thus, bathymetric solitary waves arise also in narrow non-rectangular channel flow. This equivalence also establishes a connection between our homogenized effective equations for the infinite domain and the effective equations for channels derived in Chassagne et al. (2019).

Figure 2. Channel with inclined walls. The channel is infinitely long. (a) Front view. (b) Isometric view.

Before moving on, let us mention a few additional works on how bathymetry influences wave propagation. Rosales & Papanicolau analysed the case of weakly nonlinear shallow water waves with small bathymetry changes (Rosales & Papanicolaou 1983), while Nachbin & Papanicolau studied the case of small (linear) waves with large variations in bathymetry (Nachbin & Papanicolaou 1992). Both works focused on waves in one horizontal space dimension. An extension to two dimensions was presented in Craig et al. (2005), with similar results. Based on the Green–Naghdi model (Green & Naghdi 1976), Chassagne et al. (2019) studied the dispersive effects of bores in non-rectangular channels, with an emphasis on the influence of sloping banks.

Thus far, we have dealt with the shallow water model, which neglects important physical effects such as dispersion, dissipation and surface tension. In § 5 we concentrate on waves in non-rectangular channels, with the goal of determining whether it is feasible to experimentally observe solitary waves that are created primarily by bathymetric dispersion. To do so we must go beyond the shallow water model. The shallow water equations neglect the inherent dispersive nature of water waves. If solitary waves over periodic bathymetry are observed experimentally, can one distinguish them from other solitary waves arising from other dispersive effects (such as KdV solitons)?

The model most directly relevant to the present work is perhaps that of Peregrine (Peregrine 1968) (see also Peters 1966; Shen 1969; Teng & Wu 1997), which leads to an effective KdV equation in which the dispersion coefficient is modified by the channel geometry. In §§ 4.1 and 5 we compare the dispersion in these three models (KdV, Peregrine's and ours) and show that, in certain regimes, bathymetric dispersion can be the dominant effect. This answers the question above in the affirmative. In § 4.1, we also show that our model and Peregrine's model, despite being based on different assumptions, agree well in a certain physical regime.

The code and instructions to create every plot and all the results in this work are available at https://github.com/manuel-quezada/water_wave_diffractons_RR.

1.3. Objectives and our contribution

We have two objectives in this work. The first objective, which we tackle in §§ 2 and 3, is to extend the results in Ketcheson & Quezada de Luna (2015) to obtain bathymetric solitary waves and study some of their properties. In § 4 we go beyond the work in that reference and obtain a KdV-type equation valid for weakly nonlinear regimes. The second objective, which is addressed in § 5, is to assess the feasibility of observing these waves in physical experiments. Our contribution is to

1. (i) derive an effective model for bathymetry-induced dispersion of waves in two horizontal dimensions;

2. (ii) connect two distinct areas of study: hyperbolic equations with periodic coefficients and water wave models in channels with non-flat bathymetry;

3. (iii) show that this (bathymetric) dispersion alone can lead to the formation of solitary waves;

4. (iv) derive a KdV-type equation that models these solitary waves; and

5. (v) show that bathymetric dispersion can be the dominant source of dispersion for certain classes of waves.

2. Effective dispersion due to periodic bathymetry

Water waves are inherently dispersive, and this is represented through dispersive terms or non-hydrostatic pressure terms in models such as the KdV equation (Korteweg & De Vries 1895), the Boussinesq equations (Boussinesq 1872) and the Green–Naghdi model (Green & Naghdi 1976). In this section, we demonstrate that small-amplitude shallow water waves propagating over periodic bathymetry undergo an effective dispersion. In order to clearly distinguish this source of dispersion from the dispersion present over flat bathymetry, we focus on the – dispersionless – shallow water equations over variable bathymetry

(2.1a)\begin{gather} h_t + (hu)_x + (hv)_y = 0, \end{gather}

(2.1b)\begin{gather}(hu)_t + (hu^{2})_x+gh(h+b)_x + (huv)_y = 0, \end{gather}

(2.1c)\begin{gather}(hv)_t + (huv)_x + (hv^{2})_y+gh(h+b)_y = 0. \end{gather}

Here, $h$ is the height of the column of water, $u$ and $v$ are the $x$- and $y$-velocities respectively, $g$ is the magnitude of the gravitational force and $b(x,y)$ is the periodic bathymetry. The reference point $z=0$ is chosen to coincide with the lowest point of the bathymetry; see figure 3. Unless otherwise noted, we use $g=9.8\ ({\textrm {m}}\,{\textrm {s}^{-2}})$; hereafter, we do not explicitly reference units of measure but use SI units throughout. We consider a domain that extends infinitely in both $x$ and $y$, with bathymetry periodic in $y$. We let $\varOmega$ denote the period and use it as the unit of measurement so that $\varOmega =1$.

Figure 3. Notation for shallow water equations with variable bathymetry. The reference point $z=0$ is located at the bottom of the bathymetry. The surface elevation is denoted by $\eta (x,y)$ and the undisturbed surface elevation is denoted by ${{\eta ^{*}}}$. The water depth is denoted by $h(x,y)=\eta (x,y)-b(y)$. The bathymetry $b(y)$ varies along the axis that points into the page. Figure 1(a) presents a three-dimensional view for one specific bathymetric profile. In this diagram, the bathymetry $b$ is piecewise constant; however, in general we assume it to be $y$-periodic.

The analysis and results of this section are similar to those presented in Quezada de Luna & Ketcheson (2014), which treated the acoustic wave equation in a periodic medium. These results are also connected with the work by Chassagne et al. (2019), where the authors derived an effective model for the shallow water equations that captures the dispersive effects when one-dimensional waves travel over non-flat channels. The main differences between the second reference and the results we present in this section is that the bathymetry that we consider is assumed to be changing periodically in one direction and that our effective equations are valid for propagation in two dimensions.

2.1. Linearization and homogenization

We aim to obtain a constant-coefficient homogenized system that approximates (2.1) for small-amplitude long-wavelength perturbations. To do this we follow Quezada de Luna & Ketcheson (2014) and references therein and perform a homogenization, which is valid for small-amplitude waves. We consider waves with characteristic wavelength $\lambda$ propagating in the presence of periodic bathymetry with period $\varOmega$, where $\varOmega \ll \lambda$. By letting $\delta :=\varOmega /\lambda$, we introduce a fast scale $\hat {y}:=\delta ^{-1}y$ in the $y$-direction. We assume that the solution $h$, $u$ and $v$ depend also on this fast scale; i.e. $h=h(x,y,t,\hat {y})$ and similarly for $u$ and $v$. Finally, we assume that the bathymetry depends only on the fast scale; i.e. $b=b(\hat {y})$. These are key steps in the homogenization process; see e.g. Fish & Chen (2001). Note that, by these assumptions, $(\cdot )_y\mapsto (\cdot )_y +\delta ^{-1}(\cdot )_{\hat {y}}$. Now we obtain a dimensionless version of (2.1) (after the homogenization process we go back to the variables with dimensions). This can be done by introducing the following dimensionless variables:

(2.2a–i) \begin{gather} x^{\prime}=\frac{x}{\lambda}, \quad y^{\prime}=\frac{y}{\lambda}, \quad \hat{y}^{\prime}=\frac{\hat{y}}{\lambda}, \quad t^{\prime}=\frac{c}{\lambda}t, \quad \eta^{\prime}=\frac{\eta}{\eta^{*}}, \quad h^{\prime}=\frac{h}{\eta^{*}}, \quad u^{\prime}=\frac{u}{c},\nonumber\\ v^{\prime}=\frac{u}{c}, \quad b^{\prime}=\frac{b}{\eta^{*}}, \end{gather}

where $c:=\sqrt {g{{\eta ^{*}}}}$. We remark that we scale $\hat {y}$ by $\lambda$ since the fast variation in $\hat {y}$ is introduced via its definition ($\hat {y}=\delta ^{-1}y$). After considering the system in non-conservative form and, for simplicity, dropping the tildes we get

(2.3)\begin{gather} \eta_t+[(\eta-b)u]_x+[(\eta-b)v]_y + \delta^{{-}1}[(\eta-b)v]_{\hat{y}}= 0, \end{gather}

(2.4)\begin{gather}u_t+uu_x+v(u_y+\delta^{{-}1}u_{\hat{y}})+\eta_x =0, \end{gather}

(2.5)\begin{gather}v_t+uv_x+v(v_y+\delta^{{-}1}v_{\hat{y}})+\eta_y + \delta^{{-}1}\eta_{\hat{y}}=0, \end{gather}

where $\eta =h+b$ denotes the surface elevation. We consider small-amplitude waves and perform an asymptotic expansion around $\eta ={{\eta ^{*}}}$, $u=0$ and $v=0$. The linear system is

(2.6a)\begin{gather} \eta_t + [({{\eta^{*}}}-b)u]_x + [({{\eta^{*}}}-b)v]_y + \delta^{{-}1}[({{\eta^{*}}}-b)v]_{\hat{y}} = 0, \end{gather}

(2.6b)\begin{gather}u_t + \eta_x = 0, \end{gather}

(2.6c)\begin{gather}v_t + \eta_y + \delta^{{-}1}\eta_{\hat{y}} = 0. \end{gather}

Now, let $\mu :=({{\eta ^{*}}}-b)u$ and $\nu :=({{\eta ^{*}}}-b)v$ denote the (linearized) $x$- and $y$-momentum, respectively. Doing so, we get

(2.7a)\begin{gather} \eta_t + \mu_x + \nu_y + \delta^{{-}1} \nu_{\hat{y}} = 0, \end{gather}

(2.7b)\begin{gather}\mu_t + ({{\eta^{*}}}-b(\hat{y})) \eta_x = 0, \end{gather}

(2.7c)\begin{gather}\nu_t + ({{\eta^{*}}}-b(\hat{y})) (\eta_y +\delta^{{-}1} \eta_{\hat{y}}) = 0. \end{gather}

We have explicitly noted the spatial dependence of the bathymetry $b$ in order to emphasize that (2.7) is a first-order linear hyperbolic system with spatially varying coefficients.

In the following equations and in §§ 2.1.1 and 2.1.2, to avoid confusion with subindices, we use $(\cdot )_{i,x}$ to denote differentiation of $(\cdot )_i$ with respect to $x$, and similarly for the other derivatives. Using the formal expansion $\eta (x,y,\hat {y},t)=\sum _{i=0}^{\infty }\delta ^{i}\eta _i(x,y,\hat {y},t)$ and similarly for $\mu$ and $\nu$, we get

(2.8a)\begin{gather} \sum_{i=0}^{\infty} \delta^{i}\eta_{i,t} + \sum_{i=0}^{\infty} \delta^{i}\mu_{i,x} + \sum_{i=0}^{\infty} \delta^{i}\nu_{i,y} + \delta^{{-}1}\sum_{i=0}^{\infty}\delta^{i}\nu_{i,\hat{y}} = 0, \end{gather}

(2.8b)\begin{gather}\sum_{i=0}^{\infty} \delta^{i}\mu_{i,t} + ({{\eta^{*}}}-b) \sum_{i=0}^{\infty}\delta^{i}\eta_{i,x} = 0, \end{gather}

(2.8c)\begin{gather}\sum_{i=0}^{\infty} \delta^{i}\nu_{i,t} + ({{\eta^{*}}}-b) \left(\sum_{i=0}^{\infty}\delta^{i}\eta_{i,y} + \delta^{{-}1} \sum_{i=0}^{\infty}\delta^{i}\eta_{i,\hat{y}}\right) = 0. \end{gather}

The next step is to equate terms of the same order. At each order we apply the average operator $\langle \cdot \rangle :=({1}/{|\varOmega |})\int _\varOmega (\cdot )\,{\textrm {d} y}=({1}/{\lambda })\int _0^{\lambda }(\cdot )\,\textrm {d}\hat {y}$ to obtain the homogenized leading-order system and corrections to it. In addition, we obtain expressions for the non-homogenized variables. In the next sections we present details for the derivation of the homogenized leading-order system and for the first correction. Then we present the results considering one more correction.

2.1.1. Homogenized ${{{O}}}(1)$ system

From (2.8), we equate terms of ${{{O}}}(\delta ^{-1})$ and conclude that $\nu _0=:\bar {\nu }_0(x,y,t)$ and $\eta _0=:\bar {\eta }_0(x,y,t)$ (where the bar is used to denote variables that are independent of the fast scale $\hat {y}$). From (2.8), we take terms of ${{{O}}}(1)$ to obtain

(2.9a)\begin{gather} \bar{\eta}_{0,t}+\mu_{0,x}+\bar{\nu}_{0,y}+\nu_{1,\hat{y}} = 0, \end{gather}

(2.9b)\begin{gather}\mu_{0,t}+({{\eta^{*}}}-b)\bar{\eta}_{0,x} = 0, \end{gather}

(2.9c)\begin{gather}\bar{\nu}_{0,t}+({{\eta^{*}}}-b)\left(\bar{\eta}_{0,y}+\eta_{1,\hat{y}}\right) = 0. \end{gather}

Divide (2.9c) by ${{\eta ^{*}}}-b$, apply the average operator to (2.9) and assume that $\nu _1$ and $\eta _1$ are $\hat {y}$-periodic to get

(2.10a)\begin{gather} \bar{\eta}_{0,t}+\bar{\mu}_{0,x}+\bar{\nu}_{0,y} = 0, \end{gather}

(2.10b)\begin{gather}\bar{\mu}_{0,t}+({{\eta^{*}}}-b)_m\bar{\eta}_{0,x} = 0, \end{gather}

(2.10c)\begin{gather}\bar{\nu}_{0,t}+({{\eta^{*}}}-b)_h\bar{\eta}_{0,y} = 0, \end{gather}

where

(2.11a,b)\begin{equation} ({\cdot})_m:=\langle\cdot\rangle=\frac{1}{|\varOmega|}\int_\varOmega({\cdot})\,{\textrm{d}y}\quad \textrm{and}\quad ({\cdot})_h:=\langle({\cdot})^{{-}1}\rangle^{{-}1}=\left[\frac{1}{|\varOmega|}\int_\varOmega ({\cdot})^{{-}1}\, {\textrm{d}y}\right]^{{-}1} \end{equation}

denote the mean and harmonic averages, respectively. System (2.10) is the homogenized leading-order system. It has the same form as the variable bathymetry system (2.7) but with effective constant bathymetry coefficients.

From (2.9b) and (2.10b) (and by choosing an appropriate initial condition for $\bar {\mu }_0$) we obtain

(2.12)\begin{equation} \mu_0=\frac{{{\eta^{*}}}-b}{({{\eta^{*}}}-b)_m}\bar{\mu}_0. \end{equation}

Now we obtain expressions for the non-averaged ${{{O}}}(1)$ terms in (2.9). To do this we make the following ansatz:

(2.13a)\begin{gather} \nu_1=\bar{\nu}_1+A(\hat{y})\bar{\mu}_{0,x}, \end{gather}

(2.13b)\begin{gather}\eta_1=\bar{\eta}_1+B(\hat{y})\bar{\eta}_{0,y}, \end{gather}

which is chosen to reduce (2.9) to a system of ordinary differential equations (ODEs). By substituting the ansatz (2.13), the homogenized leading-order system (2.10) and the relation for $\mu _0$ (2.12) into the ${{{O}}}(1)$ system (2.9) and by requiring that the fast variable coefficients vanish we get

(2.14a)\begin{gather} A_{,\hat{y}}+({{\eta^{*}}}-b)({{\eta^{*}}}-b)_m^{{-}1}-1 = 0, \end{gather}

(2.14b)\begin{gather}B_{,\hat{y}}-({{\eta^{*}}}-b)^{{-}1}({{\eta^{*}}}-b)_h+1 = 0. \end{gather}

We look for solutions of (2.14) with the normalization condition $\langle A\rangle =\langle B\rangle =0$. Note that $\langle A_{,\hat {y}}\rangle =\langle B_{,\hat {y}}\rangle =0$, which implies that $A$ and $B$ are $\hat {y}$-periodic.

2.1.2. Homogenized corrections

From (2.8) we collect ${{{O}}}(\delta )$ terms, plug the ansatz (2.13) and apply the average operator $\langle \cdot \rangle$ to get

(2.15a)\begin{gather} \bar{\eta}_{1,t}+\bar{\mu}_{1,x}+\bar{\nu}_{1,y} = 0, \end{gather}

(2.15b)\begin{gather}\bar{\mu}_{1,t}+({{\eta^{*}}}-b)_m\bar{\eta}_{1,x} ={-}\langle B({{\eta^{*}}}-b)\rangle\bar{\eta}_{0,xy}, \end{gather}

(2.15c)\begin{gather}\bar{\nu}_{1,t}+({{\eta^{*}}}-b)_h\bar{\eta}_{1,y} ={-}({{\eta^{*}}}-b)_h\langle A({{\eta^{*}}}-b)^{{-}1}\rangle\bar{\mu}_{0,xt}. \end{gather}

For the bathymetry that we consider in this work $\langle B({{\eta ^{*}}}-b)\rangle =\langle A({{\eta ^{*}}}-b)^{-1}\rangle =0$. Following similar (but considerably more algebraically intense) steps we obtain the homogenized second correction

(2.16a)\begin{align} \bar{\eta}_{2,t}+\bar{\mu}_{2,x}+\bar{\nu}_{2,y} &= 0, \end{align}

(2.16b)\begin{align} \bar{\mu}_{2,t}+({{\eta^{*}}}-b)_m\bar{\eta}_{2,x} &={-}\langle({{\eta^{*}}}-b)F\rangle \bar{\eta}_{0,xyy} - \langle({{\eta^{*}}}-b)E\rangle \bar{\eta}_{0,xxx}, \end{align}

(2.16c)\begin{align} \bar{\nu}_{2,t}+({{\eta^{*}}}-b)_h\bar{\eta}_{2,y} &= ({{\eta^{*}}}-b)_h^{2}\langle({{\eta^{*}}}-b)^{{-}1}D\rangle \bar{\eta}_{0,yyy}\notag\\ &\quad + ({{\eta^{*}}}-b)_m({{\eta^{*}}}-b)_h\langle({{\eta^{*}}}-b)^{{-}1}C\rangle \bar{\eta}_{0,xxy}, \end{align}

where $C,D,E$ and $F$ are fast variable functions that are given by the following ODEs:

(2.17a)\begin{gather} C_{,\hat{y}}-[1-({{\eta^{*}}}-b)_m^{{-}1}({{\eta^{*}}}-b)]B+A = 0, \end{gather}

(2.17b)\begin{gather}D_{,\hat{y}}-B = 0, \end{gather}

(2.17c)\begin{gather}E_{,\hat{y}}-({{\eta^{*}}}-b)_m({{\eta^{*}}}-b)^{{-}1}A = 0, \end{gather}

(2.17d)\begin{gather}F_{,\hat{y}}+B = 0, \end{gather}

with the normalization conditions $\langle C\rangle = \langle D\rangle = \langle E\rangle = \langle F\rangle = 0$. It can be easily shown that, for the periodic bathymetry that we consider, the coefficients in the right-hand side of (2.16) do not vanish.

Finally, given the homogenized leading-order system (2.10) and the homogenized corrections (2.15) and (2.16), we combine them into a single system by defining $\bar {\eta }:=\langle \bar {\eta }_0+\delta \bar {\eta }_1+\dots \rangle$ and similarly for $\bar {\mu }$ and $\bar {\nu }$. We obtain

(2.18a)\begin{gather} \bar{\eta}_{,t}+\bar{\mu}_{,x}+\bar{\nu}_{,y} = 0, \end{gather}

(2.18b)\begin{gather}\bar{\mu}_{,t}+({{\eta^{*}}}-b)_m\bar{\eta}_{,x} = \delta^{2}({{\eta^{*}}}-b)_m\left[\beta_1\bar{\eta}_{,xyy} + \beta_2 \bar{\eta}_{,xxx}\right], \end{gather}

(2.18c)\begin{gather}\bar{\nu}_{,t}+({{\eta^{*}}}-b)_h\bar{\eta}_{,y} = \delta^{2}({{\eta^{*}}}-b)_h\left[\gamma_1\bar{\eta}_{,yyy} + \gamma_2\bar{\eta}_{,xxy}\right], \end{gather}

where

(2.19a)\begin{gather} \beta_1 ={-}({{\eta^{*}}}-b)_m^{{-}1}\langle({{\eta^{*}}}-b)F\rangle,\quad \beta_2 ={-}({{\eta^{*}}}-b)_m^{{-}1} \langle({{\eta^{*}}}-b)E\rangle, \end{gather}

(2.19b)\begin{gather}\gamma_1 = ({{\eta^{*}}}-b)_h\langle({{\eta^{*}}}-b)^{{-}1}D\rangle,\quad \gamma_2 = ({{\eta^{*}}}-b)_m\langle({{\eta^{*}}}-b)^{{-}1}C\rangle. \end{gather}

The homogenized system (2.18) is an effective approximation of the two-dimensional linearized shallow water equations over $y$-periodic bathymetry whose period is much smaller than the characteristic wavelength.

If we assume also that the initial data do not depend on $y$, let $\bar {\eta }^{*}:=\langle {{\eta ^{*}}}-b\rangle$ denote the average undisturbed surface elevation, drop the bars in (2.18) and go back to the dimension variables, we obtain

(2.20)\begin{equation} \eta_{,tt}-g\bar{\eta}^{*}\eta_{,xx}={-}g\bar{\eta}^{*}\delta^{2}\beta_2\eta_{,xxxx}, \end{equation}

which models propagation only in the $x$-direction. The speed of small-amplitude long-wavelength perturbations is, as one might expect, $c_{{eff}}:=\sqrt {g\bar {\eta }^{*}}$. Equations (2.18)–(2.19) and (2.20) are valid for small-amplitude waves over arbitrary bathymetry that is periodic in $y$. Below, we will specialize to the piecewise-constant case.

2.2. Piecewise-constant bathymetry

Now we consider a specific case of study with piecewise-constant bathymetry

(2.21)\begin{equation} b(y) = \begin{cases} 0, & \text{if}\ n + \frac{1}{2} \le y/\varOmega < n+1, \\ {b_0}, & \text{if}\ n \le y/\varOmega < n + \frac{1}{2}, \end{cases} \end{equation}

where $n$ is any integer. This bathymetry profile is depicted in figure 1(a). The coefficients (2.19) are then

(2.22a)\begin{gather} \beta_1=\frac{{b_0}^{2}}{48(2\eta^{*}-{b_0})^{2}}\lambda^{2}, \quad \beta_2=\frac{-{b_0}^{2}}{192{{\eta^{*}}}({{\eta^{*}}}-{b_0})}\lambda^{2}, \end{gather}

(2.22b)\begin{gather}\gamma_1={-}\beta_1,\quad \gamma_2={-}\beta_2. \end{gather}

The term appearing on the right-hand side of (2.20) is dispersive; the coefficient of dispersion is in this case given by

(2.23)\begin{equation} -g\bar{\eta}^{*} \delta^{2}\beta_2 = g\bar{\eta}^{*}\left[\frac{{b_0}^{2}\varOmega^{2}}{192{{\eta^{*}}}({{\eta^{*}}}-{b_0})}\right]. \end{equation}

It is evident that this dispersion is purely an effect of the bathymetric variation; notice that it increases as ${b_0}$ grows, and vanishes as ${{\eta ^{*}}}\rightarrow \infty$ (keeping the bathymetry fixed). We remark that the dispersion in (2.20) is a consequence of changes in the bathymetry and not due to non-hydrostatic pressure effects. These dispersive effects are different from those present in dispersive water wave models, like the KdV equation, in which dispersion is present even when the bathymetry profile is flat.

In figure 4 we compare the solution of the (nonlinear) shallow water system with variable bathymetry (2.1) to that of the homogenized linear system (2.20). We take initial data

(2.24a,b)\begin{equation} \eta(x,y,t=0) ={{\eta^{*}}}+\epsilon \exp{\left(-\frac{x^{2}}{2\alpha}\right)},\quad u(x,y,0) =v(x,y,0)=0 \end{equation}

with $\epsilon =0.001, \alpha =2$ and ${{\eta ^{*}}}=0.75$. The bathymetry is given by (2.21) with

(2.25)\begin{equation} {b_0}=0.5. \end{equation}

The dispersion predicted by the linearized, homogenized model is also clearly evident in the nonlinear, variable-coefficient solution. Both models are solved to very high precision; the differences between the solutions are primarily due to the nonlinear effects that are neglected in (2.20). The shallow water equations (2.1) are solved using the finite volume code PyClaw (Ketcheson et al. 2012), with the Riemann solver developed in George (2008). The mesh resolution is $\Delta x=\Delta y=1/128$. The linear homogenized equations (2.20) are solved using a Fourier spectral collocation method in space and a fourth-order Runge–Kutta method in time; see Trefethen (2000).

Figure 4. Solution of the shallow water equations (2.1) with periodic bathymetry vs solution of the linearized homogenized approximation (2.20). We show the surface elevation $\eta$ at (from left to right) $t=20, 120$ and $t=200$. We plot two different slices of the solution at different $y$ values; however, because $\delta \ll 1$ there is almost no phase difference and the two plots are nearly exactly aligned.

3. Bathymetric solitary waves

Let us now study solutions of the nonlinear, variable-coefficient shallow water model (2.1) in a more strongly nonlinear regime. We repeat the experiment above, taking (2.24a,b) and (2.25) but with a much larger perturbation given by $\epsilon =0.05$ and $\alpha =2$ or $\alpha =10$. We again solve the equations using the finite volume solver PyClaw. The results are shown in figure 5. The mass of the initial pulse determines the number of solitary waves created. In the rest of this section we use the solitary waves shown in figure 5 and (following Ketcheson & Quezada de Luna 2015) study some properties for these solitary waves. In particular, we investigate the long-time stability and shape evolution, the scaling and speed–amplitude relations and the interaction of bathymetric solitary waves. The results from these experiments suggest that, although bathymetric solitary waves are similar to KdV solitons, they possess fundamental differences. We explore these similarities and differences in § 4. Moreover, we derive a KdV-type equation that approximates the solution of (2.1) for $x$-propagation of plane waves over periodic bathymetry like the one depicted in figure 1(a).

Figure 5. Bathymetric solitary waves at $t=340$. The initial condition is a Gaussian pulse given by (2.24a,b) with $\epsilon =0.05$, ${{\eta ^{*}}}=0.75$ and (a,b) $\alpha =2$ and (c,d) $\alpha =10$. In (a,c) we show the surface elevation plots (where the dashed line represents the location of the jump in the bathymetry) and in (b,d) we show slices along $y=0.25$ (in dashed blue) and $y=-0.25$ (in solid red).

3.1. Long-time stability and shape evolution

We first investigate the long-time behaviour and the shape evolution of these solitary waves. To do this we isolate the first solitary wave in figure 5(a,b), which corresponds to $t=340$, use it as initial condition for a new simulation and propagate it by itself until a final time of $t=1000$. Let $X(t)$ denote the location of the solitary wave's peak at time $t$; we compute

(3.1)\begin{equation} \Delta \eta =\max_{t\in[1,\dots,1000]}\left\{ \frac{\| \eta\left(x-X(0),y,t=0\right)-\eta\left(x-X(t),y,t\right) \|_{2(x,y)}}{\| \eta\left(x-X(0),y,t=0\right) \|_{2(x,y)}} \right\}, \end{equation}

which is the largest relative difference between the shape of the solution at $t\in [1,\ldots ,1000]$ and the initial condition. We perform this experiment on a grid with $\Delta x = \Delta y = 1/64$ and obtain $\Delta \eta \approx 7.66\times 10^{-4}$. Afterwards, we refine the grid so that $\Delta x = \Delta y = 1/128$ and obtain $\Delta \eta \approx 1.79\times 10^{-4}$. The results indicate that these solutions are indeed solitary waves that propagate with a fixed shape, up to numerical errors.

3.2. Scaling and speed–amplitude relations

Many solitary waves, including the diffractons found in Ketcheson & Quezada de Luna (2015), have a shape that is exactly or nearly that of a $\text {sech}^{2}$ function. Here, we investigate the shape of bathymetric solitary waves. Although these are two-dimensional waves, they vary more strongly with respect to $x$ than $y$. Our expectation (based on Ketcheson & Quezada de Luna 2015) is that the cross-section of a bathymetric solitary wave for a fixed value of $y$ should be close to a $\text {sech}^{2}$ function.

To illustrate the two-dimensional structure of these waves, in figure 6(a) we plot, for the first five solitary waves from figure 5(c,d), the peak amplitude as a function of $y$; i.e. $A(y):=\max _x \{\eta (x,y)-{{\eta ^{*}}}\}$. Observe that the variation in $y$ is stronger for larger-amplitude solitary waves; for the smallest ones, the wave is nearly uniform in $y$. This suggests that the $y$-variation is a nonlinear effect. To strengthen this argument, we plot in figure 6(b) the variation in amplitude $\Delta A:=\max _y A(y)-\min _yA(y)$ vs the mean amplitude $\bar A:=({1}/{\varOmega })\int _{-\varOmega /2}^{\varOmega /2}A(y)\,{\textrm {d}y}$. In addition, we plot a quadratic least-squares curve fitted to the data and constrained to pass through the known value $(0,0)$. The profiles plotted here are qualitatively similar to what is predicted by the model derived in Peregrine (1969), although they differ in magnitude.

Figure 6. For each solitary wave in figure 5(c,d), we plot (a) the amplitude as a function of $y$, and (b) the variation $\Delta A:=\max _y A(y)-\min _yA(y)$ as a function of mean amplitude $\bar A:=({1}/{\varOmega })\int _{-\varOmega /2}^{\varOmega /2}A(y)\,{\textrm {d}y}$. The black dashed line in (b) is a quadratic least-squares curve fitted to the data and constrained to pass through the known value $(0,0)$. In both cases, the amplitude is measured relative to the undisturbed water level ${{\eta ^{*}}}$.

We consider again the first three solitary waves in figure 5(c,d), computing the $y$-averaged surface height

(3.2)\begin{equation} \eta_m(x,t=340)\approx\frac{1}{\varOmega}\int_{-\varOmega/2}^{\varOmega/2} [\eta(x,y,t=340) - {{\eta^{*}}}]\, {\textrm{d}y}, \end{equation}

and rescaling each wave by its amplitude

(3.3)\begin{equation} \hat{\eta}_m(\hat{x}) := \frac{\eta_m\left(\hat{x}\right)}{A_m}, \end{equation}

where $A_m=\max _x \eta _m(x)$, $\hat {x}=\sqrt {A_m}(x-x_m)$, with $x_m=\text {argmax}_x \eta _m(x)$. In figure 7(a) we show the three non-scaled, $y$-averaged solitary waves (given by (3.2)) and in (b) we show the same averaged solitary waves after the scaling defined by (3.3). After this scaling, all of the bathymetric solitary waves look similar, which is a common property of many solitary waves. The black dashed line in (b) is a $\text {sech}^{2}$ function with amplitude and width fitted (in a least-squares sense) to all of the solitary waves after the scaling. In § 4 we derive a KdV-type equation that approximates the right-going part of the solution of (2.1) with variable bathymetry (2.21). The dotted cyan line in figure 7(b) is the solution of the aforementioned KdV-type equation with $A_m=1$; i.e. it is a soliton that approximates the bathymetric solitary waves for the configuration that we consider in figure 5(c,d).

Figure 7. Scaling relation for bathymetric solitary waves. In (a) we show the first three $y$-averaged solitary waves (given by (3.2)) from figure 5(c,d) centred at the origin. In (b) we show the same solitary waves after the scaling given by (3.3). The black dashed line in (b) is a $\text {sech}^{2}$ function fitted to the data and the dotted cyan line is a soliton solution of a KdV-type equation that we derive in § 4.

The KdV solitons have a linear speed–amplitude relation (Korteweg & De Vries 1895; Zabusky & Kruskal 1965). This is also true for other solitary waves. For example, stegotons, which are solitary waves created due to effective dispersion introduced by reflections in periodic media, also have a linear speed–amplitude relation (LeVeque & Yong 2003). As we discussed before, for a given bathymetric solitary wave, the amplitude is $y$-dependent. Therefore, to define the speed–amplitude relation, we must first define the amplitude of a bathymetric solitary wave. If we use the $y$-averaged wave peak amplitude, we obtain a nearly linear relation, as shown by the blue circles in figure 8. This relation also agrees well with the predicted speed–amplitude relationship for small-amplitude soliton solutions of a KdV-type model that we derive in § 4; this relation is shown with a solid purple line. If we use instead the minimum or maximum amplitude (which occur at $y=0.25$ and $y=-0.25$, respectively), we obtain nonlinear relationships, as shown also in figure 8.

Figure 8. Speed–amplitude relation for bathymetric solitary waves. We measure the amplitude based on (3.2). The solid purple line is based on a KdV-type equation that we derive in § 4 and the dashed black line is the quadratic least-squares fitted curve to the cyan circles and constrained to pass through the known value $c_{{eff}}=\sqrt {g\left \langle {{\eta ^{*}}}-b\right \rangle }$ for zero-amplitude waves. The amplitude is measured relative to the undisturbed water level ${{\eta ^{*}}}$.

3.3. Interaction of bathymetric solitary waves

Another well-known property of many solitary waves is their tendency to interact with one another only through a phase shift. In this section we study the interaction of two solitary waves that are propagating in either the same direction or opposite directions. In both situations, the solitary waves are taken from the results shown in figure 5(a,b). In all plots we show slices of the surface elevation along $y=0.25$ and $y=-0.25$.

To produce a counter-propagating collision we negate the velocity of the shorter wave. The initial condition is shown in figure 9(a). Here, the taller wave propagates to the right while the smaller one moves to the left. We show the solution at different times during and after the interaction. As a reference, we propagate the taller wave by itself and superimpose the solution using dashed lines. After the interaction, very small oscillations are visible (see figure 9d), suggesting that the interaction is not elastic. This has been reported before for diffractons (Ketcheson & Quezada de Luna 2015); however, in this case, the oscillations in the tail are much weaker. Note that there is an almost unnoticeable change in the phase for the taller solitary wave with respect to the propagation without interaction; this is due to the relatively short time of interaction. Although the phase shift is also small for diffractons (Ketcheson & Quezada de Luna 2015), the phase shift in our numerical experiments with bathymetric solitary waves is much smaller.

Figure 9. Counter-propagating collision. We show (in different colour) slices at the middle of each bathymetry section. As a reference, we plot the propagation of the taller solitary wave by itself. In (d) we zoom in on the tails of the solitary waves after the interaction to notice the oscillations and the slight change in phase.

Now we consider a collision where both waves move in the same direction. The initial condition is shown in figure 10(a). Here, both waves move to the right. Since the taller wave moves faster (see § 3.2) it eventually reaches and passes the smaller one. Again, we show the solution at different times during and after the interaction. As a reference, we propagate the taller wave by itself and superimpose the solution using dashed lines. In this case, there are no visible oscillations after the collision, suggesting an elastic collision. In contrast to the counter-propagating collision, the interaction time is larger, which leads to a noticeable shift in phase. This is a common feature for other solitary waves and for diffractons.

Figure 10. Co-propagating collision at different times. We show (in different colour) slices at the middle of each bathymetry section. As a reference, we plot the propagation of the taller solitary wave by itself. In (f) we zoom in on the tails of the solitary waves after the interaction.

4. KdV model for weakly nonlinear bathymetric solitary waves

In the previous section we explored some properties of bathymetric solitary waves. The shape of any given of these waves is not $y$-independent; i.e. these are truly two-dimensional waves that travel in one direction. However, their shape is closely connected with a one-dimensional $\text {sech}^{2}$ function. In addition, bathymetric solitary waves travel without significant change in the shape and interact similar to KdV solitons. All these properties strongly suggest that one might model small-amplitude bathymetric solitary waves via a KdV equation. We do not expect such a model to be completely accurate for large-amplitude bathymetric solitary waves since, as the amplitude increases, the waves behave less like solitons; see § 3.2. In this section, we obtain a KdV-type equation that accounts for bathymetric dispersion and is valid for weakly nonlinear waves. Before doing that, however, in the next section we compare the shallow water solutions with Peregrine's equation that accounts for additional sources of water wave dispersion.

4.1. Bathymetric solitary waves via an inherently dispersive water wave model

Another water wave model that accounts for dispersive effects due to changes in the bathymetry has been derived and analysed in Peregrine (1968) and Teng & Wu (1997). That model for flow in non-rectangular channels takes the form

(4.1)\begin{equation} \eta_t+\sqrt{g\bar{\eta}^{*}}\eta_x +\frac{3}{2}\sqrt{\frac{g}{\bar{\eta}^{*}}}(\eta-\bar{\eta}^{*})\eta_x +\frac{\kappa^{2}}{6}(\bar{\eta}^{*})^{2}\sqrt{g\bar{\eta}^{*}}\eta_{xxx}=0. \end{equation}

Equation (4.1) is a KdV-type equation with a modified dispersion coefficient. Here, $\kappa ^{2}$ is called the channel shape factor; it depends on the cross-section of the channel and is given by

(4.2)\begin{equation} \kappa^{2}=\frac{3}{(\bar{\eta}^{*})^{2}}\left[\frac{1}{|D|}\int_D\varPsi(y,z)\,\textrm{d}y\,\textrm{d}z-\frac{1}{|L|}\int_L\varPsi(y,{{\eta^{*}}})\,{\textrm{d}y}\right], \end{equation}

where $D\subset \mathbb {R}^{2}$ is the cross-section of the channel, $L\subset D$ is the top boundary of the cross-section (i.e. the undisturbed free surface) and $\varPsi$ is the solution of the following elliptic boundary value problem:

(4.3a)\begin{gather} \varPsi_{yy}+\varPsi_{zz} = 1, \end{gather}

(4.3b)\begin{gather}\varPsi_z|_{z={{\eta^{*}}}} = \bar{\eta}^{*}, \end{gather}

(4.3c)\begin{gather}\varPsi_{\boldsymbol{n}} = 0. \end{gather}

Here, $\boldsymbol {n}$ is the unit vector normal to the boundary of the cross-section of the channel. For a rectangular channel, $\kappa =1$ so (4.1) becomes the standard KdV equation.

It is natural to ask how Peregrine's model (4.1) compares with solutions of the shallow water equations in a non-rectangular channel – or equivalently, over the kind of periodic bathymetry we have considered. To answer this, we first observe that we cannot expect agreement between the models if the shallow water model leads to shock formation; this can occur if the initial data are very large or if the bathymetry variation is small (i.e. if $\kappa \approx 1$) (Ketcheson & Quezada de Luna 2020).

On the other hand, if the bathymetry variation is relatively large and the initial data are not too large, the two models can be in relatively close agreement for fairly long times. An example of this is shown in figure 11(a–c). Here, we take bathymetry given by (2.21) with ${b_0}=0.01$. The initial data are given by (2.24a,b) with $\eta ^{*}=0.015$, $\epsilon =0.001$ and $\alpha =2$. Note that, for this case, $\kappa ^{2}\approx 214.14$. We solve (4.1) using a Fourier pseudospectral collocation method. We see very close agreement between the solutions, even after the waves have travelled a distance equal to hundreds of times the channel width.

Figure 11. In green, solution of Peregrine's model (4.1). In black, $y$-averaged solution of the shallow water equations (2.1) with periodic bathymetry. The bathymetry is given by (2.21) with $b_0=0.01$ for (a–c) and $b_0=0.5$ for (d–f). The initial condition is given by (2.24a,b) with $\alpha =2$, $\eta^*=0.015$ and $\epsilon=0.001$ for (a–c) and $\alpha=2$, $\eta^*=0.75$ and $\epsilon=0.05$ for (d–f).

In figure 11(d–f) we show another example, in which bathymetric dispersion is much less dominant. The bathymetry and the amplitude of the initial data are 50 times larger, with ${b_0}=0.5$, $\eta ^{*}=0.75$ and $\epsilon =0.05$. Note that, for this case, $\kappa ^{2}\approx 1.58$. We see that the resulting solutions are completely different, even from a qualitative perspective. This is expected, since the model leading to (4.1) includes additional dispersion beyond that caused by bathymetry variation. The shallow water model neglects that additional dispersion and may therefore be less accurate in any regime (like that of the latter scenario) where it is important. In § 5, we compare the strength of these two types of dispersion independently and compare each against that predicted by Peregrine's model; see figure 15.

Remark 4.1 About the amplitude of the initial data

The shallow water equations (2.1) with the initial condition (2.24a,b) model the propagation of a left- and a right-going wave. On the other hand, the KdV-type equation (4.1) models the propagation of a right-going wave. For this reason, the amplitude $\epsilon$ of the initial condition (2.24a,b) differs by a factor of two between (2.1) and (4.1). For simplicity, we report the amplitude used for the shallow water equations, for the KdV-type equations we use half of that amplitude.

4.2. KdV-type equation with purely bathymetric dispersion

Now we derive a KdV-type equation whose dispersive effects are only a consequence of changes in bathymetry. This equation models small-amplitude bathymetric solitary waves. In the work by LeVeque & Yong (2003), the authors considered a one-dimensional nonlinear system and derived homogenized equations with a dispersive correction. Later, in Ketcheson & Quezada de Luna (2015), the authors extended the results to two dimensions. In these references, the nonlinearity was chosen to facilitate the homogenization process. Ideally, we would aim to proceed as in § 2 and obtain a nonlinear homogenized system with constant coefficients and a dispersive correction. Unfortunately, the nonlinear terms in the shallow water equations complicate the process. Based on these references and the results in § 2, it is possible to make an ansatz for what the homogenized nonlinear system should look like. We hypothesize that

(4.4a)\begin{gather} \eta_t+\left[(\eta-b_m)u\right]_x = \delta^{2}\varPhi, \end{gather}

(4.4b)\begin{gather}u_t+uu_x+g\eta_x = g\delta^{2}\beta_2\eta_{xxx} \end{gather}

is a homogenized system that models the $x$-propagation of shallow water waves over a general $y$-periodic bathymetry. Here, $b_m=\langle b\rangle ={b_0}/2$, $\beta _2$ is the coefficient of dispersion defined in (2.19) and $\delta ^{2}\varPhi =\delta ^{2}\varPhi (\eta ,\eta _x,\eta _{xx},\eta _{xxx},u,u_{x},u_{xx},u_{xxx})$ is an ${O}(\delta ^{2})$ nonlinear term that depends not only on the solution but also on its derivatives. We do not know a closed form for $\varPhi$; however, we expect it to introduce (potentially nonlinear) dispersive effects.

From (4.4) we derive a KdV-type equation. To do this we follow e.g. Garnier, Grajales & Nachbin (2007). First we consider the Riemann invariants of the left-hand side of (4.4)

(4.5a)\begin{gather} {R^{-}} =2\sqrt{g(\eta-b_m)}-u, \quad \text{on} \ {\textrm{d}\kern0.7pt x}/\textrm{d}t=u-\sqrt{g(\eta-b_m)}, \end{gather}

(4.5b)\begin{gather}{R^{+}} =2\sqrt{g(\eta-b_m)}+u, \quad \text{on}\ {\textrm{d}\kern0.7pt x}/\textrm{d}t=u+\sqrt{g(\eta-b_m)}. \end{gather}

By plugging (4.5) into (4.4) we obtain a system of partial differential equations (PDEs) for the space–time evolution of the Riemann invariants ${R^{-}}$ and ${R^{+}}$. We focus on the right-going invariant ${R^{+}}$ and choose ${R^{-}}=2\sqrt {g\bar {\eta }^{*}}$, which is set constant to match with the solution as $|x|\rightarrow \infty$. Here, $\bar {\eta }^{*}={{\eta ^{*}}}-b_m$. Doing this, we obtain

(4.6)\begin{equation} {R^{+}}_{\!\!\!\!t}-\frac{1}{4}({R^{-}}-3{R^{+}}){R^{+}}_{\!\!\!\!x}=\frac{\delta^{2}\beta_2}{8}[({R^{-}}+{R^{+}}){R^{+}}_{\!\!\!\!xxx}+3{R^{+}}_{\!\!\!\!x}{R^{+}}_{\!\!\!\!xx}] + \delta^{2} C, \end{equation}

where $C$ is an unknown function of ${R^{-}}$, ${R^{+}}$ and the derivatives of ${R^{+}}$; i.e. $C=C({R^{-}},{R^{+}},{R^{+}}_{\!\!\!\!x},{R^{+}}_{\!\!\!\!xx},{R^{+}}_{\!\!\!\!xxx})$. By using the Riemann invariants (4.5), we go back to the physical variables and obtain

(4.7)\begin{equation} \eta_t-\sqrt{g\bar{\eta}^{*}}\eta_x+3\sqrt{g(\eta-b_m)}\eta_x=\delta^{2}\frac{\beta_2}{2}\sqrt{g(\eta-b_m)}\eta_{xxx} + \delta^{2}\hat{\varPhi}, \end{equation}

where $\hat {\varPhi }$ is an unknown nonlinear function of $\eta$ and its derivatives; i.e. $\hat {\varPhi }=\hat {\varPhi }(\eta ,\eta _x,\eta _{xx},\eta _{xxx})$. Although, we do not know the exact form of $C$ and $\hat {\varPhi }$, these terms appear as a consequence of the manipulations of $\varPhi$ via the Riemann invariants. Based on § 2.1 and the references therein, we expect and assume that $\hat {\varPhi }$ is a dispersive term. Finally, we expand the nonlinear term $\sqrt {\eta -b_m}$ around ${{\eta ^{*}}}$ and drop the terms that are of size ${O}(\delta ^{2}\varepsilon )$, where $\varepsilon \sim \eta -\eta ^{*}$. We obtain

(4.8a)\begin{equation} \eta_t+\sqrt{g\bar{\eta}^{*}}\eta_x+\frac{3}{2}\sqrt{\frac{g}{\bar{\eta}^{*}}}(\eta-{{\eta^{*}}})\eta_x+\sigma(\gamma)\sqrt{g\bar{\eta}^{*}}\eta_{xxx} = 0, \end{equation}

with

(4.8b)\begin{equation} \sigma(\gamma)=\delta^{2}\frac{|\beta_2|}{2}(1+\gamma), \end{equation}

where $\gamma$ is a constant (to be determined) that accounts for the linear dispersive part of $\delta ^{2}\hat {\varPhi }$. Equation (4.8a) is a KdV-type equation with a dispersion coefficient given by (4.8b). It models the $x$-propagation of shallow water waves over a two-dimensional $y$-periodic bathymetry. Different types of bathymetry are modelled via $\beta _2$, which is given by (2.19). Just like the linear dispersive equation derived in § 2, the KdV-type equation (4.8) only accounts for dispersive effects due to changes in the bathymetry.

It is important to remark that (4.8) is a weakly nonlinear approximation of (the right-going part of) (4.4); therefore, we expect it to be a good approximation only for small-amplitude waves. Indeed, solitary wave solutions of (4.8) travel with a speed proportional to their amplitude (see § 4.2.1); however, from § 3.2, we know that the speed–amplitude relation for bathymetric solitary waves is approximately linear only for small-amplitude waves.

4.2.1. Profile and speed of weakly nonlinear bathymetric solitary waves

We now look for a travelling wave solution of (4.8) by substituting into that equation the ansatz $\eta =\bar {\eta }^{*} f(x-Vt)$, where $V$ is the speed of the travelling wave. By doing so, we obtain an ODE that we can integrate twice to get that the shape of weakly nonlinear bathymetric solitary waves is given by

(4.9a)\begin{equation} \eta-{{\eta^{*}}} = A_m{\textrm{sech}}^{2}\left[\left(\frac{A_m}{8\sigma(\gamma)\bar{\eta}^{*}}\right)^{1/2}(x-Vt)\right], \end{equation}

where $A_m$ is the amplitude of the solitary wave and

(4.9b)\begin{equation} V=\sqrt{g\bar{\eta}^{*}}\left(1+\frac{A_m}{2\bar{\eta}^{*}}\right) \end{equation}

is its speed. Let us now estimate the correction value $\gamma$ in (4.8b). To do this we use (4.9a) with $\gamma =0$ to generate multiple initial conditions for the shallow water equations (2.1) with bathymetry given by (2.21) with ${b_0}=0.5$. The initial condition for the velocity is

(4.10)\begin{equation} u(x,y,0)=2\sqrt{g(\eta(x,y,0)-b_m)}-2\sqrt{g\bar{\eta}^{*}}, \quad v(x,y,0)=0, \end{equation}

where $b_m={b_0}/2$. Let us use six initial conditions with

(4.11) \begin{align} A_m&=6.25\times 10^{{-}5}, \quad 1.25\times 10^{{-}4},\quad 2.5\times 10^{{-}4}, \quad 5\times 10^{{-}4},\notag\\ &\quad 1\times 10^{{-}3} \quad \text{and}\quad 2\times 10^{{-}3}. \end{align}

Then we propagate the waves until a final time of $t=100$. For each simulation, the initial condition quickly evolves into a solitary wave after some mass is left behind. Initial conditions that are close to being solitary waves undergo only small changes. In figure 12, we show the evolution of two of these waves (using $A_m=6.25\times 10^{-5}$ and $A_m=2\times 10^{-3}$). Finally, for each simulation, we isolate the solitary wave at $t=100$ and compute $\gamma$ such that (4.9a) is the closest (in a least-squares sense) to the corresponding isolated solitary wave. In figure 13 we plot the value of $\gamma$ as a function of amplitude. It is clear that $\gamma \approx 0$ for arbitrarily small-amplitude waves.

Figure 12. Solution of the shallow water equations (2.1) with periodic bathymetry given by (2.21) with ${b_0}=0.5$. The initial condition is given by (4.9a) with $\gamma =0$ and (4.10). We show a slice along $y=-0.25$, with initial amplitude (a–c) $A_m=6.25\times 10^{-5}$ and (d–f) $A_m=2\times 10^{-3}$.

Figure 13. Estimated $\gamma$ such that (4.9a) is the closest (in a least-squares sense) to a given bathymetric solitary wave. The bathymetric solitary waves correspond to the solution of the shallow water equations (2.1) (at $t=100$) with periodic bathymetry given by (2.21) with ${b_0}=0.5$. The initial condition is given by (4.9a) with $\sigma =\sigma (0)$, (4.10) and amplitude $A_m$ given by (4.11).

Finally, we compare the solution of (4.8) with $\gamma =0$ vs the numerical solution of the shallow water equations (2.1). We consider the same two scenarios shown in figure 11 from § 4.1. We solve (4.8) using a Fourier pseudospectral collocation method.

First, we consider the same scenario as in figure 11(a–c); the bathymetry is given by (2.21) with ${b_0}=0.01$ and the initial condition is given by (2.24a,b) with $\alpha =2$, ${{\eta ^{*}}}=0.015$ and $\epsilon =0.001$. The results are shown in figure 14(a–c). We see even closer agreement between the two models than in figure 11(a–c).

Figure 14. In red, solution of the KdV-type equation (4.8). In black, $y$-averaged solution of the shallow water equations (2.1) with periodic bathymetry. The bathymetry is given by (2.21) with ${b_0}$ as indicated below. The initial condition is given by (2.24a,b) with $\alpha =2$ and ${{\eta ^{*}}}$ and $\epsilon$ as indicated below. (a–c) Simulation of the same scenario as in figure 11(a–c), this time comparing (4.8) with the shallow water equations (2.1). We use (2.21) with ${b_0}=0.01$ and (2.24a,b) with ${{\eta ^{*}}}=0.015$ and $\epsilon =0.001$. (d–i) Simulation of the same scenario as in figure 11(d–f), this time comparing (4.8) with the shallow water equations (2.1). We use (2.21) with ${b_0}=0.5$ and (2.24a,b) with ${{\eta ^{*}}}=0.75$ and $\epsilon =0.05$. (j–o) The same as figure 14(d–i), but with an initial wave that is twice as tall ($\epsilon =0.1$).

Next, we consider the same scenario as in figure 11(d–f), which is also used in figure 5(a,b). In this case, the bathymetry is given by (2.21) with ${b_0}=0.5$ and the initial condition is given by (2.24a,b) with $\alpha =2$, ${{\eta ^{*}}}=0.75$ and $\epsilon =0.05$. The results are shown in figure 14(d–i). Recall that for this problem, the dispersion inherently present in water waves and bathymetric dispersion are both important in Peregrine's model. But the shallow water equations and the model (4.8) both account for only bathymetric dispersion, so we see much better agreement here.

Finally, in figure 14(j–o), we consider the last scenario but with an initial wave that is twice as tall. We see that the agreement between the models is worse (and the agreement with (4.1) would be even worse). Both models (4.1) and (4.8) include only a linear dispersive term, whereas in § 3.2 we observed that the speed–amplitude relationship of bathymetric solitary waves is somewhat nonlinear. Furthermore, for very large initial data, the shallow water solution will contain shocks (Ketcheson & Quezada de Luna 2020), which cannot be represented by the models (4.1) or (4.8). We conclude that for sufficiently large initial data and long times, neither of these models will remain close to the shallow water solution.

5. Comparison of dispersive effects

In this work we have used theoretical and numerical tools to analyse a new class of solitary waves. Our objective in the remainder of the paper is to determine the feasibility of observing them in experiments. In this section we compare bathymetric dispersion with that inherently present in water wave models. To do that, we compare our model (4.8) with KdV, and determine conditions under which bathymetric dispersion should be much stronger than the dispersion in KdV. We refer to the dispersion present in the KdV equation as ‘KdV dispersion’.

Consider a flat channel and a weakly nonlinear regime in which the KdV equation (Korteweg & De Vries 1895) is applicable. The right-going KdV equation is given by

(5.1)\begin{equation} \eta_t+\sqrt{g{{\eta^{*}}}}\eta_x+\frac{3}{2}\sqrt{\frac{g}{{{\eta^{*}}}}}(\eta-{{\eta^{*}}})\eta_x+\frac{1}{6}({{\eta^{*}}})^{2}\sqrt{g{{\eta^{*}}}}\eta_{xxx}=0, \end{equation}

where, as usual, $\eta (x,t)$ is the surface elevation and ${{\eta ^{*}}}$ is the undisturbed surface elevation. The KdV model has been validated experimentally, for instance in Hammack & Segur (1974), wherein water waves (over flat bathymetry) were observed to form solitary waves of the kind predicted by (5.1) after propagating over a relatively long distance. To explore the qualitative difference between bathymetric and KdV dispersion, we study the dispersion relation of (4.8) and (5.1), which are given by

(5.2a)\begin{gather} \omega_{Hom} = \sqrt{g\bar{\eta}^{*}}k\left[1 - \sigma(0) k^{2} \right], \end{gather}

(5.2b)\begin{gather}\omega_{KdV} = \sqrt{g{{\eta^{*}}}}k\left[1 -\tfrac{1}{6}({{\eta^{*}}}k)^{2} \right], \end{gather}

respectively. Note that (5.2a) can also be obtained directly from (2.20). Here, $\sigma (0)$ is given by (4.8b) (with $\gamma =0$ for small-amplitude waves) and $\bar {\eta }^{*}$ is the average undisturbed depth for the non-rectangular channel. It is natural to take the depth of the flat channel ${{\eta ^{*}}}$ equal to $\bar {\eta }^{*}$, in which case the ${O}(k)$ terms for the two models agree. The dominant dispersive effect arises from the term of ${O}(k^{3})$ in (5.2). In figure 15(a) we compare the coefficient of this term in the two models, taking ${{\eta ^{*}}}=0.75$ and a range of bathymetry profiles given by (2.21) with ${b_0} \in [0,0.75)$. In the figure, the blue and the red plots are ${(\bar {\eta }^{*})^{2}}/{6}$ and $\sigma (0)$ respectively. As one might expect, when the value of ${b_0}$ is small, the bathymetric dispersion is small compared to the KdV dispersion. On the other hand, if ${b_0}$ is close to the mean water level, bathymetric dispersion is stronger and can be of the same order or much larger than KdV dispersion. Thus, at least for small-amplitude, long-wavelength waves, the two dispersive effects can be made comparable or either one can be made dominant depending on the parameter ${b_0}$.

Figure 15. Coefficients from the ${O}(k^{3})$ terms in the dispersion relations (5.2), and that of Peregrine's model (4.1). The blue and red curves correspond to KdV equation (5.1) and our model (4.8), respectively. The dispersive coefficient of Peregrine's model is plotted in green. (a) Regime considered in figure 11(d–f). Here, ${{\eta ^{*}}}=0.75$. (b) Regime considered in figure 11(a–c). Here, ${{\eta ^{*}}}=0.015$.

As discussed in § 4.1, neither the KdV equation (5.1) nor our model (4.8) include both types of dispersive effects. Peregrine's model (4.1) includes both sources of dispersion. In figure 15(a), we also plot (in green) the coefficient of dispersion appearing in (4.1). Note that the dispersion predicted by Peregrine's model coincides with that predicted by KdV and by our model in the limits when $b_0$ is close to $0$ and $\eta ^{*}$, respectively. This behaviour is expected: when $b_0\approx 0$, the bathymetry is almost flat so the main source of dispersion is that predicted by KdV. On the other hand, when $b_0\approx \eta ^{*}=0.75$, bathymetric dispersion is dominant; see (4.8b) and (2.22). In figure 15(a), we mark (with a dashed cyan line) $b_0=0.5$, which corresponds to the situation studied in figure 11(d–f). In this case, the dispersion in Peregrine's model and bathymetric dispersion are significantly different, leading to the different solutions depicted in figure 11(d–f).

On the other hand, for the situation we considered in figure 11(a–c) with ${b_0}=0.01$ and ${{\eta ^{*}}}=0.015$, the bathymetric dispersion is much stronger. Consequently, the dispersive effects in our model (4.8) and Peregrine's model are similar, which explains the agreement in the simulations shown in figure 11(a–c). In figure 15(b), we plot the corresponding dispersive coefficients of the three models (KdV (5.1), our model (4.8) and Peregrine's model (4.1)).

From (5.2a) and (5.2b), it is possible to find a channel configuration that leads to bathymetric dispersion that has the same magnitude as the dispersive effects appearing in the classical KdV equation. We get that if

(5.3)\begin{equation} {b_0} = \frac{8(\bar{\eta}^{*})^{2}}{\varOmega^{2}} \left[{-}4\bar{\eta}^{*} + \sqrt{16(\bar{\eta}^{*})^{2}+1}\right], \end{equation}

then the shallow water equations (2.1) with bathymetry given by (2.21) are a close approximation (up to ${O}(k^{3})$) to the classical KdV equation (5.1).

Remark 5.1 About the effect of dissipation

As we concluded in this section, either source of dispersion (KdV or bathymetric dispersion) can be made dominant. Since the dissipative effects in the propagation of water waves in flat channels do not prevent the formation of solitons (see for instance Hammack & Segur 1974), it is reasonable to believe that bathymetric solitary waves can also be observed in a physical experiment. An important difference between these two scenarios is the presence of friction at the interface between the flat sections of the bathymetry. More detailed studies (or physical experiments) are needed to determine if this extra source of dissipation might prevent the formation of bathymetric solitary waves.

6. Conclusions

We have shown that bathymetric variation in an infinite periodic domain leads to an effective dispersion of water waves, and have related this to the already-studied phenomenon of dispersion of waves in non-rectangular channels. This dispersion is distinct from the dispersion accounted for in wave models like KdV, and can on its own lead to solitary wave formation, which we call bathymetric solitary waves, even when the dominant behaviour would normally be wave breaking. Weakly nonlinear plane waves in this setting approximately satisfy a KdV-type equation. This KdV-type model leads to soliton waves that approximate small-amplitude bathymetric solitary waves. However, it is important to remark that bathymetric solitary waves are truly two-dimensional waves that travel in one direction. Bathymetric solitary waves, therefore, behave similar to the solitons emerging from the derived KdV-type model (4.8) only when the amplitude is small enough. More strongly nonlinear waves exhibit more pronounced two-dimensional structure, have a nonlinear speed–amplitude relation and evidently involve nonlinear dispersion.

We have shown in §§ 4.1 and 5 that the model by Peregrine (1968) agrees with the KdV equation (5.1), which considers the inherently dispersive nature of water waves, and with our model (4.8), which considers bathymetric dispersion, in certain asymptotic regimes. In general, however, the combination of these two types of dispersion is non-trivial. In particular, as shown in figure 15, the coefficient of dispersion in Peregrine's model (Peregrine 1968) is not simply the sum of the coefficients of dispersion in KdV and in our model. A proper analysis to identify the range of validity and agreement between Peregrine's model, KdV and our model requires the solution of the elliptic PDE (4.3). Doing so could also provide insights about how these two types of dispersive effects are combined together. This investigation is an area of future work.

Although we focus on piecewise-constant bathymetry, similar phenomena appear in computational experiments with other kinds of bathymetry. To demonstrate this we consider a channel with inclined walls, like the one shown in figure 2. A similar problem was studied in Chassagne et al. (2019), where the authors considered the dispersive Green–Naghdi model (Green & Naghdi 1976) and showed that the changes in the bathymetry enhance the dispersive effects. Here, we show that solitary waves arise in this setting also as solutions of the (dispersionless) shallow water equations (2.1). In figure 16, we show the solution at $t=200$. In figure 16(a) we show the surface plot and in figure 16(b) we plot a slice along $y=0.5$.

Figure 16. Solitary waves at $t=200$ in a channel like the one shown in figure 2. The initial condition is given by (2.24a,b) with $\epsilon =0.05$, ${{\eta ^{*}}}=0.75$ and $\alpha =2$. In the (a) we show the surface plots and in (b) we show slices along $y=0.5$.

It is natural to ask whether experimental observation of these waves is feasible. We have conducted computational experiments (not shown here) using the three-dimensional Navier–Stokes equations. Preliminary results indicate that, for scenarios like those studied in this work, wave breaking is almost entirely avoided and an initial pulse breaks into multiple peaks, which then evolve into solitary waves. Further and more detailed numerical investigation of these waves is part of ongoing investigation and will be published elsewhere.