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 (Reference Quezada de Luna and Ketcheson2014) 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 (Reference Fish and Chen2001). 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. Reference Ketcheson, Mandli, Ahmadia, Alghamdi, Quezada de Luna, Parsani, Knepley and Emmett2012), with the Riemann solver developed in George (Reference George2008). 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 (Reference Trefethen2000).

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.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 (Reference Ketcheson and Quezada de Luna2015), 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 Reference Ketcheson and Quezada de Luna2015) 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 (Reference Peregrine1969), 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 Reference Korteweg and De Vries1895; Zabusky & Kruskal Reference Zabusky and Kruskal1965). 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 Reference LeVeque and Yong2003). 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 Reference Ketcheson and Quezada de Luna2015); 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 Reference Ketcheson and Quezada de Luna2015), 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.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 (Reference Peregrine1968) and Teng & Wu (Reference Teng and Wu1997). 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 Reference Ketcheson and Quezada de Luna2020).

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.

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 (Reference LeVeque and Yong2003), the authors considered a one-dimensional nonlinear system and derived homogenized equations with a dispersive correction. Later, in Ketcheson & Quezada de Luna (Reference Ketcheson and Quezada de Luna2015), 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 (Reference Garnier, Grajales and Nachbin2007). 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 Reference Ketcheson and Quezada de Luna2020), 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.