Solitary water waves created by variations in bathymetry

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 KdV-type equation that approximates the behavior 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 %, as well as confirm the analysis and suggest that experimental observation of these solitary waves is possible.


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 et al., 1968]. Similar conclusions regarding general linear wave propagation were reached by using Bloch expansions in [Santosa and Symes, 1991] and by using homogenization theory in [Fish and 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 and 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 and De Vries, 1895]. This behavior 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 LeVeque, 2012, Ketcheson andQuezada de Luna, 2020]. In the multidimensional setting, effective dispersion arises not only from variation in the impedance, but also from variation in the linearized sound speed [Quezada de Luna and Ketcheson, 2014]. This latter source of effective dispersion can also lead to solitary wave formation for nonlinear waves [Ketcheson and 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.
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 1a.
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 and 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 versus 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.

Solitary water-waves in narrow non-rectangular channels
Dispersion of water waves due to bathymetry has been analyzed and observed in a very different setting: flow in a channel like that depicted in Figure 1b 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 [Peregrine, 1968, Peters, 1966, Shen, 1969, Teng and 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 1a and the non-rectangular channel in Figure 1b. 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 (a) Bathymetry that changes periodically in one direction.
(b) Narrow channel with non-flat bathymetry. 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].  Before moving on, let us mention a few additional works on how bathymetry influences wave propagation. Rosales & Papanicolau analyzed the case of weakly nonlinear shallow water waves with small bathymetry changes, [Rosales and Papanicolaou, 1983], while Nachbin & Papanicolau studied the case of small (linear) waves with large variations in bathymetry, [Nachbin and 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 and Naghdi, 1976], Chassagne et. al. [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 and 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.

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 and 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: • derive an effective model for bathymetry-induced dispersion of waves in two horizontal dimensions; • connect two distinct areas of study: hyperbolic equations with periodic coefficients and water wave models in channels with non-flat bathymetry; • show that this (bathymetric) dispersion alone can lead to the formation of solitary waves; • derive a KdV-type equation that models these solitary waves; • show that bathymetric dispersion can be the dominant source of dispersion for certain classes of waves.

Effective dispersion due to periodic bathymetry
Water waves are inherently dispersive, and this is represented through dispersive terms or nonhydrostatic pressure terms in models such as the KdV equation [Korteweg and De Vries, 1895], the Boussinesq equations [Boussinesq, 1872] and the Green-Naghdi model [Green and 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: 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. Unless otherwise noted, we use g = 9.8 m 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 Ω denote the period and use it as the unit of measurement so that Ω = 1.
The analysis and results of this section are similar to those presented in [Quezada de Luna and 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.

Linearization and homogenization
We aim to obtain a constant-coefficient homogenized system that approximates (1) for small-amplitude, long-wavelength perturbations. To do this we follow [Quezada de Luna and Ketcheson, 2014] and references therein and perform a homogenization, which is valid for small-amplitude waves. We consider waves with characteristic wavelength λ propagating in the presence of periodic bathymetry with period Ω, where Ω λ. By letting δ := Ω/λ, we introduce a fast scaleŷ := δ −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,ŷ) and similarly for u and v. Finally, we assume that the bathymetry depends only on the fast scale; i.e., b = b(ŷ). These are key steps in the homogenization process; see e.g., [Fish and Chen, 2001]. Note that, by these assumptions, (·) y → (·) y + δ −1 (·)ŷ. Now we obtain a dimensionless version of (1) (after the homogenization process we go back to the variables with dimensions). This can be done by introducing the following dimensionless variables: where c := √ gη * . We remark that we scaleŷ by λ since the fast variation inŷ is introduced via its definition (ŷ = δ −1 y). After considering the system in non-conservative form and, for simplicity, 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 η(x, y) and the undisturbed surface elevation is denoted by η * . The water depth is denoted by h(x, y) = η(x, y) − b(y). The bathymetry b(y) varies along the axis that points into the page. Figure 1a presents a 3D 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.
dropping the tildes we get where η = h + b denotes the surface elevation. We consider small-amplitude waves and perform an asymptotic expansion around η = η * , u = 0 and v = 0. The linear system is Now, let µ := (η * − b)u and ν := (η * − b)v denote the (linearized) x-and y-momentum, respectively. Doing so, we get We have explicitly noted the spatial dependence of the bathymetry b in order to emphasize that (4) is a first-order linear hyperbolic system with spatially-varying coefficients. In the following equations and in Sections 2.1.1 and 2.1.2, to avoid confusion with subindices, we use use (·) i,x to denote differentiation of (·) i with respect to x, and similarly for the other derivatives. Using the formal expansion η(x, y,ŷ, t) = ∞ i=0 δ i η i (x, y,ŷ, t) and similarly for µ and ν, we get The next step is to equate terms of the same order. At each order we apply the average operator · := 1 |Ω| Ω (·)dy = 1 λ λ 0 (·)dŷ 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.
From (6b) and (7b) (and by choosing an appropriate initial condition forμ 0 ) we obtain Now we obtain expressions for the non-averaged O(1) terms in (6). To do this we make the following ansatz: which is chosen to reduce (6) to a system of ODEs. By substituting the ansatz (9), the homogenized leading-order system (7) and the relation for µ 0 (8) into the O(1) system (6) and by requiring that the fast variable coefficients vanish we get We look for solutions of (10) with the normalization condition A = B = 0. Note that A ,ŷ = B ,ŷ = 0, which implies that A and B areŷ-periodic.

Homogenized corrections
From (5) we collect O(δ) terms, plug the ansatz (9) and apply the average operator · to get η 1,t +μ 1,x +ν 1,y = 0, (11a) For the bathymetry that we consider in this work B(η * − b) = A(η * − b) −1 = 0. Following similar (but considerably more algebraic intense) steps we obtain the homogenized second correction η 2,t +μ 2,x +ν 2,y = 0, (12a) where C, D, E and F are fast variable functions that are given by the following ODEs: with the normalization conditions C = D = E = F = 0. It can be easily shown that, for the periodic bathymetry that we consider, the coefficients in the right hand side of (12) do not vanish. Finally, given the homogenized leading-order system (7) and the homogenized corrections (11) and (12), we combine them into a single system by definingη := η 0 + δη 1 + . . . and similarly forμ and ν. We obtainη ,t +μ ,x +ν ,y = 0, (14a) where The homogenized system (14) is an effective approximation of the two dimensional linearized shallow water equations over y-periodic bathymetry whose period is much smaller than the characteristic wave length.
If we assume also that the initial data does not depend on y, letη * := η * − b denote the average undisturbed surface elevation, drop the bars in (14) and go back to the dimension variables, we obtain which models propagation only in the x-direction. The speed of small-amplitude, long wavelength perturbations is, as one might expect, c eff := √ gη * . Equations (14)- (15) and (16) are valid for smallamplitude waves over arbitrary bathymetry that is periodic in y. Below, we will specialize to the piecewise-constant case.

Piecewise-constant bathymetry
Now we consider a specific case of study with piecewise-constant bathymetry where n is any integer. This bathymetry profile is depicted in Figure 1a. The coefficients (15) are then The term appearing on the right hand side of (16) is dispersive; the coefficient of dispersion is in this case given by 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 η * → ∞ (keeping the bathymetry fixed). We remark that the dispersion in equation (16) 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 (1) to that of the homogenized linear system (16). We take initial data with = 0.001, α = 2 and η * = 0.75. The bathymetry is given by (17) 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 (16). The shallow water equations (1) are solved using the finite volume code PyClaw , with the Riemann solver developed in [George, 2008]. The mesh resolution is ∆x = ∆y = 1/128. The linear homogenized equations (16) are solved using a Fourier spectral collocation method in space and a fourth order Runge-Kutta method in time; see [Trefethen, 2000].  (1) with periodic bathymetry versus solution of the linearized homogenized approximation (16). We show the surface elevation η 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 δ 1 there is almost no phase difference and the two plots are nearly exactly aligned.

Bathymetric solitary waves
Let us now study solutions of the nonlinear, variable-coefficient shallow water model (1) in a more strongly nonlinear regime. We repeat the experiment above, taking (20) and (21) but with a much larger perturbation given by = 0.05 and α = 2 or α = 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 and 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 (1) for x-propagation of plane waves over periodic bathymetry like the one depicted in Figure 1a.

Long-time stability and shape evolution
We first investigate the long-time behavior and the shape evolution of these solitary waves. To do this we isolate the first solitary wave in Figure 5a, 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 which is the largest relative difference between the shape of the solution at t ∈ [1, . . . , 1000] and the initial condition. We perform this experiment on a grid with ∆x = ∆y = 1/64 and obtain ∆η ≈ 7.66 × 10 −4 . Afterwards, we refine the grid so that ∆x = ∆y = 1/128 and obtain ∆η ≈ 1.79 × 10 −4 .
The results indicate that these solutions are indeed solitary waves that propagate with a fixed shape, up to numerical errors.

Scaling and speed-amplitude relations
Many solitary waves, including the diffractons found in [Ketcheson and Quezada de Luna, 2015], have a shape that is exactly or nearly that of a 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 and 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 sech 2 function.
To illustrate the two-dimensional structure of these waves, in Figure 6a we plot, for the first five solitary waves from Figure 5b, the peak amplitude as a function of y; i.e. A(y) := max x {η(x, y) − η * }. 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 6b the variation in amplitude ∆A := max y A(y) − min y A(y) versus the mean amplitudeĀ := 1 Ω Ω/2 −Ω/2 A(y)dy. 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], though they differ in magnitude.
We consider again the first three solitary waves in Figure 5b, computing the y-averaged surface height and rescaling each wave by its amplitudeη where In the left panel of Figure  7 we show the three non-scaled, y-averaged solitary waves (given by (23)) and in the right panel we show the same averaged solitary waves after the scaling defined by (24). 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 the right panel is a sech 2 function with amplitude and width fitted (in a least  −Ω/2 A(y)dy. 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 η * .
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 (1) with variable bathymetry (17). The dotted cyan line in the right panel of Figure 7 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 5b.
KdV solitons have a linear speed-amplitude relation [Korteweg andDe Vries, 1895, Zabusky andKruskal, 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 and 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.

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  We measure the amplitude based on (23). 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 = g η * − b for zero-amplitude waves. The amplitude is measured relative to the undisturbed water level η * .
propagating in either the same direction or opposite directions. In both situations, the solitary waves are taken from the results shown in Figure 5a. 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 the top-left panel of Figure 9. 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 the bottom-right panel), suggesting that the interaction is not elastic. This has been reported before for diffractons [Ketcheson and 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 and Quezada de Luna, 2015], the phase shift in our numerical experiments with bathymetric solitary waves is much smaller. Now we consider a collision where both waves move in the same direction. The initial condition is shown in the top-left panel of Figure 10. 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.

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 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.

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 analyzed in [Peregrine, 1968, Teng andWu, 1997]. That model for flow in non-rectangular channels takes the form Equation (25) is a KdV-type equation with a modified dispersion coefficient. Here κ 2 is called the channel shape factor; it depends on the cross section of the channel and is given by where D ⊂ R 2 is the cross section of the channel, L ⊂ D is the top boundary of the cross section (i.e., the undisturbed free surface) and Ψ is the solution of the following elliptic boundary value problem: Here n is the unit vector normal to the boundary of the cross section of the channel. For a rectangular channel, κ = 1 so (25) becomes the standard KdV equation.
It is natural to ask how Peregrine's model (25) 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 is very large or if the bathymetry variation is small (i.e. if κ ≈ 1) [Ketcheson and Quezada de Luna, 2020].
On the other hand, if the bathymetry variation is relatively large and the initial data is not too large, the two models can be in relatively close agreement for fairly long times. An example of this is shown in Figure 11a. Here we take bathymetry given by (17) with b 0 = 0.01. The initial data is given by (20) with η * = 0.015, = 0.001, and α = 2. Note that for this case κ 2 ≈ 214.14. We solve (25) using a Fourier pseudospectral collocation method. We see very close agreement between the solutions, even after the waves have traveled a distance equal to hundreds of times the channel width.
In Figure 11b 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, η * = 0.75, and = 0.05. Note that for this case κ 2 ≈ 1.58. We see that the resulting solutions are completely different, even from a qualitative perspective. This is expected, since the model leading to (25) 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 Section 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.1 (About the amplitude of the initial data). The shallow water equations (1) with the initial condition (20) model the propagation of a left-and a right-going wave. On the other hand, the KdV-type equation (25) models the propagation of a right-going wave. For this reason, the amplitude of the initial condition (20) differs by a factor of two between (1) and (25). For simplicity, we report the amplitude used for the shallow water equations, for the KdV-type equations we use half of that amplitude.

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 bathymetrtic solitary waves. In the work by [LeVeque and Yong, 2003], the authors considered a one-dimensional nonlinear system and derived homogenized equations with a dispersive correction. Later in [Ketcheson and 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 is a homogenized system that models the x-propagation of shallow water waves over general y-periodic bathymetry. Here b m = b = 1 2 b 0 , β 2 is the coefficient of dispersion defined in (15) and δ 2 Φ = δ 2 Φ(η, η x , η xx , η xxx , u, u x , u xx , u xxx ) is an O(δ 2 ) nonlinear term that depends not only on the solution but also on its derivatives. We do not know a closed form for Φ; however, we expect it to introduce (potentially nonlinear) dispersive effects.
From (28) we derive a KdV-type equation. To do this we follow e.g. [Garnier et al., 2007]. First we consider the Riemann invariants of the left hand side of (28): By plugging (29) into (28) we obtain a system of 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 √ gη * , which is set constant to match with the solution as |x| → ∞. Hereη * = η * − b m . Doing this, we obtain 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 (29), we go back to the physical variables and obtain whereΦ is an unknown nonlinear function of η and its derivatives; i.e.,Φ =Φ(η, η x , η xx , η xxx ). Although, we do not know the exact form of C andΦ, these terms appear as a consequence of the manipulations of Φ via the Riemann invariants. Based on Section 2.1 and the references therein, we expect and assume thatΦ is a dispersive term. Finally, we expand the nonlinear term √ η − b m around η * and drop the terms that are of size O(δ 2 ε), where ε ∼ η − η * . We obtain with where γ is a constant (to be determined) that accounts for the linear dispersive part of δ 2Φ . Equation (31a) is a KdV-type equation with a dispersion coefficient given by (31b). It models the x-propagation of shallow water waves over two-dimensional y-periodic bathymetry. Different types of bathymetry are modeled via β 2 , which is given by (15). Just like the linear dispersive equation derived in §2, the KdV-type equation (31) only accounts for dispersive effects due to changes in the bathymetry. It is important to remark that (31) is a weakly nonlinear approximation of (the right-going part of) (28); therefore, we expect it to be a good approximation only for small-amplitude waves. Indeed, solitary wave solutions of (31) 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.

Profile and speed of weakly nonlinear bathymetric solitary waves
We now look for a traveling wave solution of (31) by substituting into that equation the ansatz η =η * f (x − V t), where V is the speed of the traveling 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 where A m is the amplitude of the solitary wave and is its speed. Let us now estimate the correction value γ in (31b). To do this we use (32a) with γ = 0 to generate multiple initial conditions for the shallow water equations (1) with bathymetry given by (17) with b 0 = 0.5. The initial condition for the velocity is where b m = 1 2 b 0 . Let us use six initial conditions with A m = 6.25 × 10 −5 , 1.25 × 10 −4 , 2.5 × 10 −4 , 5 × 10 −4 , 1 × 10 −3 and 2 × 10 −3 .
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 × 10 −5 and A m = 2 × 10 −3 ). Finally, for each simulation, we isolate the solitary wave at t = 100 and compute γ such that (32a) is the closest (in a least squares sense) to the corresponding isolated solitary wave. In Figure 13 we plot the value of γ as a function of amplitude. It is clear that γ ≈ 0 for arbitrarily small-amplitude waves. Finally, we compare the solution of (31) with γ = 0 versus the numerical solution of the shallow water equations (1). We consider the same two scenarios shown in Figure 11 from §4.1. We solve (31) using a Fourier pseudospectral collocation method.
First we consider the same scenario as in Figure 11a; the bathymetry is given by (17) with b 0 = 0.01 and the initial condition is given by (20) with α = 2, η * = 0.015 and = 0.001. The results are shown in Figure 14a. We see even closer agreement between the two models than in Figure 11a.
Next we consider the same scenario as in Figure 11b, which is also used in Figure 5a. In this case, the bathymetry is given by (17) with b 0 = 0.5 and the initial condition is given by (20) with α = 2, η * = 0.75 and = 0.05. The results are shown in Figure 14b. Recall that for Peregrine's model (refer back to Figure 11b) in this case the dispersion inherently present in water waves and bathymetric dispersion are both important. But the shallow water equations and the model (31) both account for only bathymetric dispersion, so we see much better agreement here.
Finally, in Figure 14c, 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 (25) would be even worse). Both models (25) and (31) 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 and Quezada de Luna, 2020], which cannot be represented by the models (25) or (31). We conclude that for sufficiently large initial data and long times, neither of these models will remain close to the shallow water solution.
(a) Simulation of the same scenario as in Figure 11a, this time comparing (31) with the shallow water equations (1). We use (17)  (c) The same as Figure 14b, but with an initial wave that is twice as tall ( = 0.1). Figure 14: In red, solution of the KdV-type equation (31). In black, y-averaged solution of the shallow water equations (1) with periodic bathymetry. The bathymetry is given by (17) with b 0 as indicated in the subfigures. The initial condition is given by (20) with α = 2 and η * and as indicated in the subfigures.

Comparison of dispersive effects
In this work we have used theoretical and numerical tools to analyze 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 (31) 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 and De Vries, 1895] is applicable. The right-going KdV equation is given by where as usual η(x, t) is the surface elevation and η * is the undisturbed surface elevation. The KdV model has been validated experimentally, for instance in [Hammack and Segur, 1974], wherein water waves (over flat bathymetry) were observed to form solitary waves of the kind predicted by (35) after propagating over a relatively long distance. To explore the qualitative difference between bathymetric and KdV dispersion, we study the dispersion relation of (31) and (35), which are given by respectively. Note that (36a) can also be obtained directly from (16). Here σ(0) is given by (31b) (with γ = 0 for small-amplitude waves) andη * is the average undisturbed depth for the non-rectangular channel. It is natural to take the depth of the flat channel η * equal toη * , 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 (36).
In Figure 15a we compare the coefficient of this term in the two models, taking η * = 0.75 and a range of bathymetry profiles given by (17) with b 0 ∈ [0, 0.75). In the figure, the blue and the red plots are (η * ) 2 6 and σ(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 .
As discussed in Section 4.1, neither the KdV equation (35) nor our model (31) include both types of dispersive effects. Peregrine's model (25) includes both sources of dispersion. In Figure 15a, we also plot (in green) the coefficient of dispersion appearing in (25). 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 η * , respectively. This behavior is expected: when b 0 ≈ 0, the bathymetry is almost flat so the main source of dispersion is that predicted by KdV. On the other hand, when b 0 ≈ η * = 0.75, bathymetric dispersion is dominant; see (31b) and (18). In Figure 15a, we mark (with a dashed cyan line) b 0 = 0.5, which corresponds to the situation studied in Figure 11b. In this case, the dispersion in Peregrine's model and bathymetric dispersion are significantly different, leading to the different solutions depicted in Figure 11b.
On the other hand, for the situation we considered in Figure 11a with b 0 = 0.01 and η * = 0.015, the bathymetric dispersion is much stronger. Consequently, the dispersive effects in our model (31) and Peregrine's model are similar, which explains the agreement in the simulations shown in Figure  11a. In Figure 15b, we plot the corresponding dispersive coefficients of the three models (KdV (35), our model (31) and Peregrine's model (25)).
From equations (36a) and (36b), 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 then the shallow water equations (1) with bathymetry given by (17) are a close approximation (up to O(k 3 )) to the classical KdV equation (35).  Remark 5.0.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 and 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.

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 behavior 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 KdVtype model (31) 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 (35), which considers the inherently dispersive nature of water waves, and with our model (31), 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 (27). 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 and 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 (1). In Figure 16, we show the solution at t = 200. In the left panel we show the surface plot and in the right panel we plot a slice 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 3D 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.