# Solitary wave-mean flow interaction in strongly nonlinear dispersive shallow water waves

Thibault Congy, Gennady El, Sergey Gavrilyuk, Mark Hoefer, Keh-Ming Shyue

► To cite this version:

Thibault Congy, Gennady El, Sergey Gavrilyuk, Mark Hoefer, Keh-Ming Shyue. Solitary wave-mean flow interaction in strongly nonlinear dispersive shallow water waves. 2025. hal-05053277

HAL Id: hal-05053277

https://hal.science/hal-05053277v1

Preprint submitted on 1 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Solitary wave-mean flow interaction in strongly nonlinear
dispersive shallow water waves

Thibault CONGY*, Gennady EL†, Sergey GAVRILYUK‡, Mark HOEFER§,
Keh-Ming SHYUE¶

May 6, 2025

Abstract

The interaction of a solitary wave and a slowly varying mean background or flow for the Serre-Green-Naghdi (SGN) equations is studied using Whitham modulation theory. The exact form of the three SGN-Whitham modulation equations—two for the mean horizontal velocity and depth decoupled from one for the solitary wave amplitude field—are obtained exactly in the solitary wave limit. Although the three equations are not diagonalizable, the restriction of the full system to simple waves for the mean equations is diagonalized in terms of Riemann invariants. The Riemann invariants are used to analytically describe the head-on and overtaking interactions of a solitary wave with a rarefaction wave and dispersive shock wave (DSW), leading to scenarios of solitary wave trapping or transmission by the mean flow. The analytical results for overtaking interactions prove that a simpler, approximate approach based on the DSW fitting method is accurate to the second order in solitary wave amplitude, beyond the first order accurate Korteweg-de Vries approximation. The analytical results also accurately predict the SGN DSW’s solitary wave edge amplitude and speed. The analytical results are favourably compared with corresponding numerical solutions of the full SGN equations. Because the SGN equations model the bi-directional propagation of strongly nonlinear, long gravity waves over a flat bottom, the analysis presented here describes large amplitude solitary wave-mean flow interactions in shallow water waves.

# 1 Introduction

A fundamental and important problem in continuum mechanics is the interaction of waves with the inhomogeneity of the medium through which they propagate. The waves could be linear or nonlinear and, in many cases, the inhomogeneity can be modelled as externally imposed through prescribed variable coefficients to a wave-type partial differential equation (PDE) [8]. When the waves and the medium are dynamically coupled, a common occurrence in geophysical fluid dynamics [7], the

<footnote>*Department of Mathematics, Physics and Electrical Engineering, Northumbria University, Newcastle upon Tyne,
NE1 8ST, UK, thibault.congy@northumbria.ac.uk</footnote>

<footnote>†Department of Mathematics, Physics and Electrical Engineering, Northumbria University, Newcastle upon Tyne,
NE1 8ST, UK, gennady.el@northumbria.ac.uk</footnote>

<footnote>‡Aix Marseille Univ, CNRS, IUSTI, UMR 7343, Marseille, France, sergey.gavrilyuk@univ-amu.fr</footnote>

<footnote>§Department of Applied Mathematics, University of Colorado, Boulder, Colorado, USA, hoefer@colorado.edu</footnote>

<footnote>¶Ocean Center, National Taiwan University, Taipei 106, Taiwan, shyue@ntu.edu.tw</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Scenarios of the solitary wave-mean flow interaction in shallow water waves. a), c): solitary wave transmission; b), d) solitary wave trapping. Reproduced, with permission, from [3].</figcaption></figure>

problem becomes more challenging to describe analytically. A natural framework to approach this class of problems is to utilize scale separation and derive equations separating the motion of waves and averaged quantities, e.g., the mean fluid density and velocity. While the focus of analytical studies was historically on linear or weakly nonlinear waves interacting with mean flows [7], a recent body of work has emerged for the case where the waves are strongly nonlinear, i.e., solitons or solitary waves; see, e.g., the review [3]. A schematic of such wave-mean flow scenarios involving solitary wave interaction with rarefaction and dispersive shock waves is shown in Fig. 1.

For the asymptotic analysis of solitary wave-mean flow interaction, the mean flow is assumed to vary on a much slower spatial scale than the solitary wave width. In this case, the mean flow evolution decouples from solitary wave motion so that the interaction is one-way: only the mean influences the motion of the solitary wave, not vice-versa. This problem was first considered theoretically and experimentally in the unidirectional case where a solitary wave and either a rarefaction wave (RW) mean flow or a dispersive shock wave (DSW) mean flow exhibited an overtaking interaction [44]. Depending on the conditions, the solitary wave was observed to be either transmitted through or trapped by the mean flow as $t \to \infty$. The key mathematical insight was the determination of the solitary wave limit of the Whitham modulation equations—a quasi-linear system of first order equations that describe the slow evolution of nonlinear wavetrains [57]—and its subsequent diagonalization in terms of Riemann invariants. A general approach for obtaining the solitary wave limit in the case of Hamiltonian partial differential equations with a local Hamiltonian can be found in [5]. While Riemann invariants can in principle be obtained for any system of two quasi-linear, first order PDEs, they are most easily obtained for completely integrable equations such as the Korteweg-de Vries (KdV) equation, either directly [56] or using the finite gap method [18, 34].
But Whitham theory can be developed for integrable and non-integrable equations alike [14]. A recent paper has obtained the diagonalization, i.e., the Riemann invariants, for the non-integrable Benjamin-Bona-Mahony (BBM) equation [4] by deriving an exact representation of the solitary wave limit of the BBM-Whitham modulation equations [24].

One of our main results in this paper is the derivation of the modulation solitary wave limit for the apparently non-integrable Serre-Green-Naghdi (SGN) equations [50, 52, 28, 27] describing long, strongly nonlinear free surface waves on a flat bottom in a fluid of constant density. The Hamiltonian structure of this system is non-local [26, 41], so the generic results of [5] cannot be used. The SGN system shares many properties with the BBM equation: they have similar dispersion relations and similar nonlocal Hamiltonian structure [48]. Compared to the BBM equation, which is a unidirectional, scalar equation, the SGN system is a Galilean invariant system containing more physics than the BBM equation. The Whitham modulation equations for the SGN equations were obtained in [12]. Their hyperbolicity was proven in [54]. Being second order in time, the SGN equations support bi-directional wave propagation and the resulting solitary wave limit quasi-linear system of the SGN-Whitham equations is third order. In this work, we prove the existence of Riemann invariants if and only if simple waves in the decoupled mean flow equations are considered.

The analysis of unidirectional wave-mean flow interaction has been carried out for solitons and RWs, DSWs for a general class of unidirectional nonlinear dispersive wave equations in [44] by analyzing the solitary wave limit of the Whitham modulation equations. In the same paper, the approach was applied to the conduit equation and compared with experiments on viscous fluid conduits. The approach was further refined and applied to the KdV equation [3] and the modified-KdV equation [55] where, in addition, the kink (monotone, heteroclinic travelling wave solution) serves as either the wave or the mean flow. The extension to oblique, two-dimensional line solitons interacting with RWs and DSWs was developed for the Kadomtsev-Petviashvili equation in [49]. A similar analysis of soliton-mean field interaction was extended to the bi-directional case of the defocusing nonlinear Schrödinger (NLS) equation where, in addition to overtaking interactions, RWs and DSWs experience head-on interactions with solitons [51].

The soliton-mean flow interaction problem for integrable equations has also been studied using the inverse scattering transform (IST) for the KdV equation [2, 3] and the focusing NLS equation [6] where some solitons were shown to leave a trailing “wake” inside the DSW after passing through. Other analytical approaches include perturbation theory [3], the Darboux transformation [47], both for KdV soliton-RW interaction, and a Hamiltonian formulation of the problem that borrows ideas from soliton perturbation theory and Whitham theory to obtain an approximate, analytical description of solitary wave-mean flow interaction for a generalized KdV equation [36] and the NLS equation [32, 35].

In this paper, we utilise the derived SGN-Whitham equations in the solitary wave limit to analytically describe the head-on and overtaking interaction of a solitary wave and a RW mean flow. The results are then extended to solitary wave-DSW mean flow interaction. One important finding in this paper is that the original approach proposed in [44] and utilised elsewhere [36] that rely upon the DSW fitting method and its conjugate wavenumber/dispersion relation [10] to obtain the soliton Riemann invariant is only approximate. By obtaining the exact representation of the solitary wave limit of the SGN-Whitham equations and corresponding Riemann invariants for simple wave mean flows, we are able to prove that the conjugate wavenumber/dispersion yields a second order in amplitude accurate prediction for solitary wave motion through a RW. It also provides a second order accurate prediction for the DSW’s solitary wave edge. Careful numerical simulations of the SGN equations agree with these findings. While this is an improvement to the
first order accurate, weakly nonlinear KdV approximation, it identifies a limitation of the DSW fitting approach.

The plan of the paper is as follows. In section 2, we give a detailed presentation of the SGN equations and their periodic traveling wave solutions are presented in section 3. In section 4, the modulation equations are given both in the mass Lagrangian coordinates and in Eulerian coordinates. The reason is that the solitary wave limit is mathematically easier to carry out in Lagrangian coordinates, while the physical interpretation is easier in Eulerian coordinates. In section 5, we study the interaction of solitary waves with RWs and DSWs. In particular, the transmission and trapping effect is studied. The closed-form analytical results are in good agreement with the numerical ones for the SGN equations. Finally, the main technical details (the different forms of the modulation equations, the passage to the solitary limit, numerical method, etc.) are given in four Appendices.

## 2 Serre–Green–Naghdi equations

The SGN equations over a flat bottom approximating the free-surface Euler equations in the long wave limit are [50, 52, 28, 27]

$$h_t + (hu)_x = 0, \qquad (1)$$

$$u_t + uu_x + gh_x = \frac{1}{h} \left( \frac{h^3}{3} (u_{xt} + uu_{xx} - u_x^2) \right)_x, \qquad (2)$$

where $h$ is the total depth, $u$ is the depth-averaged horizontal velocity, and $g$ is the acceleration due to gravity. In what follows, we scale independent and dependent variables so that $g = 1$.

The SGN equations (1), (2) are non-integrable and represent a fully nonlinear generalization of the classical Boussinesq equations [37]. The first equation (1) is the exact equation for conservation of mass, and the second equation (2) can be manipulated into the equation for conservation of horizontal momentum

$$(hu)_t + \left(hu^2 + \frac{1}{2}h^2\right)_x = \left(\frac{h^3}{3}(u_{tx} + uu_{xx} - u_x^2)\right)_x. \qquad (3)$$

An equivalent form of the momentum equation is

$$(hu)_t + \left(hu^2 + \frac{1}{2}h^2 + \frac{1}{3}h^2 \frac{d^2h}{dt^2}\right)_x = 0, \qquad (4)$$

where $\frac{d}{dt}$ is the material derivative notation, $\frac{d}{dt} = \frac{\partial}{\partial t} + u\frac{\partial}{\partial x}$, and $\frac{d^2}{dt^2} = \frac{d}{dt}\left(\frac{d}{dt}\right)$ is the second material derivative. The conservation of energy manifests in the additional conservation law

$$\left(\frac{1}{2}h\left(h + u^2 + \frac{1}{3}h^2 u_x^2\right)\right)_t + \left(hu\left(h + \frac{1}{2}u^2 + \frac{1}{2}h^2 u_x^2 - \frac{1}{3}h^2(u_{xt} + uu_{xx})\right)\right)_x = 0. \qquad (5)$$

Finally, a fourth conservation law (the so-called Bernoulli conservation law) can be derived. It is usually written for the variable $K = u - \frac{1}{3h}(h^3 u_x)_x$ representing the tangent component of the
fluid velocity at the free surface [20]

$$ \mathcal{K}_t + \left( \mathcal{K}u + h - \frac{u^2}{2} - \frac{1}{2}h_x^2 u^2 \right)_x = 0. \qquad (6) $$

A mathematical justification of the SGN model (1), (2) can be found in [37, 45]. Recent years have seen increased activity in both the study of qualitative properties of the solutions to the SGN equations and in the development of numerical discretization techniques [38, 40, 17, 21, 23].

The Lagrangian for (1), (2), where the mass conservation law (1) is considered a constraint, is given as [26]

$$ \mathcal{L} = \int_{-\infty}^{+\infty} h \left( \frac{u^2}{2} + \frac{1}{6} \left( \frac{dh}{dt} \right)^2 - \frac{h}{2} \right) dx. \qquad (7) $$

In order to obtain analytically tractable expressions, we will make use of mass Lagrangian coordinates introduced as follows. Let $h_0(x)$ be the initial position of the free surface, and $Y$ be the classical Lagrangian coordinate. The mass Lagrangian coordinate $q$ is defined as $q = \int_0^Y h_0(s)ds$. Let $x = x(t, q)$ be the trajectories of fluid particles. Since the mass conservation law is in the form $h \frac{\partial x}{\partial q} = 1$, the Lagrangian (7) can be transformed to

$$ \tilde{\mathcal{L}} = \int_{-\infty}^{+\infty} \left( \frac{u^2}{2} - \tilde{e}(\tau, \tau_t) \right) dq, \qquad (8) $$

where $\tau = 1/h$ and the potential $\tilde{e}(\tau, \tau_t)$ is

$$ \tilde{e}(\tau, \tau_t) = \frac{1}{2\tau} - \frac{1}{6} \left( \frac{\partial}{\partial t} \left( \frac{1}{\tau} \right) \right)^2 = \frac{1}{2\tau} - \frac{\tau_t^2}{6\tau^4}. \qquad (9) $$

For a general potential $\tilde{e}(\tau, \tau_t)$, the Euler-Lagrange equations for (8) are [19]:

$$ \tau_t - u_q = 0, \quad u_t + p_q = 0, \qquad (10) $$

where the pressure $p$ is defined by

$$ p = -\frac{\delta \tilde{e}}{\delta \tau} = - \left( \frac{\partial \tilde{e}}{\partial \tau} - \frac{\partial}{\partial t} \left( \frac{\partial \tilde{e}}{\partial \tau_t} \right) \right). \qquad (11) $$

System (10) is reminiscent of the $p$-system for the barotropic Euler equations. However, in our case, the pressure $p$ defined by (11) depends not only on $\tau$, but also on its first and second temporal derivatives. The variable $\tau$ is the analogue of the specific volume for the corresponding Euler equations. As a consequence of Noether's theorem, the variational formulation implies the conservation of energy and Bernoulli equation

$$ \left( \frac{u^2}{2} + \varepsilon \right)_t + (pu)_q = 0, \quad \varepsilon = \tilde{e} - \tau_t \tilde{e}_{\tau_t}, \qquad (12) $$

$$ \left( \tau u - \tau_q \frac{\partial \tilde{e}}{\partial \tau_t} \right)_t - \left( \frac{u^2}{2} - \tau p - \tilde{e} \right)_q = 0. \qquad (13) $$

that are the analogues of eqs. (5) and (6). The conservation laws (10), (12), (13) are averaged to obtain the modulation equations in Lagrangian coordinates in Appendix A.
# 3 Periodic travelling wave solutions

We now present the form of periodic travelling wave solutions in both Eulerian and Lagrangian coordinates. Several parameterisations of the four dimensional family of periodic travelling waves will be presented.

## 3.1 Eulerian coordinates

The travelling wave solutions $h(x,t) = h(\xi) = h(x - ct)$, $u(x,t) = u(\xi) = u(x - ct)$ to the SGN equations (1), (2) are specified by

$$ \left(\frac{dh}{d\xi}\right)^2 = \frac{3}{h_1 h_2 h_3} (h - h_1)(h - h_2)(h_3 - h) = \frac{3}{h_1 h_2 h_3} P(h), \quad (14) $$

$$ u = c - \sigma \frac{\sqrt{h_1 h_2 h_3}}{h}, \quad \sigma = \pm 1, \quad (15) $$

where, similar to the defocusing nonlinear Schrödinger equation, the wave is characterised by four independent parameters $h_1, h_2, h_3$ (the roots of the third order polynomial $P(h)$), the travelling wave velocity $c$, and $\sigma = +1$ ($\sigma = -1$) corresponds to the fast (slow) waves. Here, the depth variations in the wave occur in the interval $[h_2, h_3]$, for $h_1 < h_2 < h < h_3$ with the amplitude $a = h_3 - h_2$. It is assumed that there are no vacuum points for non-trivial periodic travelling wave solutions: $h_1 > 0$. The ODE in (14) can be integrated to obtain

$$ h = h_2 + (h_3 - h_2) \mathrm{cn}^2 \left( \frac{1}{2} \sqrt{\frac{3(h_3 - h_1)}{h_1 h_2 h_3}} \xi, m \right), \quad m = \frac{h_3 - h_2}{h_3 - h_1}, \quad (16) $$

where cn is the Jacobi cosine elliptic function [1]. We express the physical parameters (amplitude $a$, wavenumber $k$, and the period averages of depth $\bar{h}$ and velocity $\bar{u}$) in terms of the basic parameter set ($h_1, h_2, h_3, c$) as

$$ \begin{gathered} a = h_3 - h_2, \quad k = \sqrt{\frac{3(h_3 - h_1)}{h_1 h_2 h_3}} \frac{\pi}{2K(m)}, \quad \bar{h} = h_1 + (h_3 - h_1) \frac{E(m)}{K(m)}, \\ \overline{h(u - c)} = -\sigma \sqrt{h_1 h_2 h_3}, \quad \bar{u} = c - \sigma \sqrt{h_1 h_2 h_3} \frac{\Pi\left(1 - \frac{h_2}{h_3}, m\right)}{h_3 K(m)}, \end{gathered} \quad (17) $$

where $K, E$ and $\Pi$ are the complete elliptic integrals of the first, second, and third kinds, respectively [1]. We introduce the phase

$$ \theta = k\xi = kx - \omega t, \quad (18) $$

so that the traveling wave is $2\pi$-periodic in $\theta$: $h(\xi) = h(\theta/k) = h((\theta + 2\pi)/k)$, $u(\xi) = u(\theta/k) = u((\theta + 2\pi)/k)$. The wavelength of the traveling wave $L$ is related to the wavenumber $k$ by $L = 2\pi/k$. The solitary wave limit of (16) is achieved when $h_2 \to h_1$ so that $m \to 1$. Its explicit form is

$$ h(x,t) = \bar{h} + a \operatorname{sech}^2 \left( \frac{\sqrt{3a}}{\bar{h}\sqrt{\bar{h}} + a} (x - ct) \right), \quad u(x,t) = \bar{u} + \sigma \sqrt{\bar{h} + a} \left( 1 - \frac{\bar{h}}{h(x,t)} \right). \quad (19) $$
Fast ($\sigma = +1$) and slow ($\sigma = -1$) elevation solitary waves propagate on the background $h = \bar{h}$, $u = \bar{u}$, and are characterised by the speed-amplitude relation

$$c = c_s(a, \bar{h}, \bar{u}) \equiv \bar{u} + \sigma \sqrt{\bar{h} + a}. \qquad (20)$$

Fast solitary waves move faster than the dispersionless long-wave velocity $V_+ = \bar{u} + \sqrt{\bar{h}}$, and slow solitary waves move slower than the dispersionless long-wave velocity $V_- = \bar{u} - \sqrt{\bar{h}}$ (see Figure 2).

<figure><img src="image_3.png" /><figcaption>Figure 2: Relations (20) between the velocity of a fast (slow) solitary wave and the corresponding long-wave velocity $V_+$ ($V_-$) are represented schematically on the left (right).</figcaption></figure>

In the opposite, harmonic, limit, $h_2 \to h_3$ so that $m \to 0$, yields a small-amplitude linear wave characterised by the dispersion relation $\omega_0(k)$ (frequency-wavenumber relation) for linear waves propagating on the background $(\bar{h}, \bar{u})$

$$\omega = kc = \omega_0(k, \bar{h}, \bar{u}) \equiv k\bar{u} + \sigma k \sqrt{\frac{\bar{h}}{1 + \bar{h}^2 k^2 / 3}}. \qquad (21)$$

## 3.2 Lagrangian coordinates

Consider now SGN's travelling wave solutions in the co-moving mass Lagrangian coordinate $\zeta = q - \tilde{c}t$ instead of the Eulerian coordinate $\xi = x - ct$: $h = h(\zeta)$, $u = u(\zeta)$. Since $\frac{d\zeta}{d\xi} = h$, the equation for traveling waves is

$$\left(\frac{dh}{d\zeta}\right)^2 = \frac{3}{h_1 h_2 h_3} \frac{P(h)}{h^2}. \qquad (22)$$

The velocity $u$ is found from the relation

$$\tilde{c} \tau + u = \text{cst}, \quad \text{with} \quad \tau = \frac{1}{\bar{h}}. \qquad (23)$$

The wave velocity $\tilde{c}$ is related to the roots of the polynomial $P(h)$ by the formula (see Appendix A)

$$\tilde{c}^2 = h_1 h_2 h_3. \qquad (24)$$
## 4 Modulation equations and their solitonic reduction

The SGN equations possess all the necessary prerequisites for the application of Whitham averaging. They support a family of $2\pi$-periodic travelling wave solutions $h(\theta/k)$, $u(\theta/k)$ specified by (16), (18) and characterized by the four independent parameters $h_1, h_2, h_3, c$ while admitting four independent conservation laws (1), (2), (5) and (6).

The modulation equations for the SGN equations in Eulerian coordinates are given in [12, 54]. According to the general averaging procedure [56], the modulation system for the SGN equations can be obtained by period averaging any three conservation laws, such as mass (1), momentum (3) and energy (5), over the periodic family (16), and augmenting them by the wave conservation equation $k_t + (kc)_x = 0$. Doing so results in

$$ \begin{gathered} \bar{h}_t + (\overline{hu})_x = 0, \\ (\overline{hu})_t + \left( \frac{1}{2}\overline{h^2} + \overline{hu^2} - \frac{1}{3}\overline{h^3} \left( (u-c)u_{\xi\xi} - u_{\xi\xi}^2 \right) \right)_x = 0, \\ \left( \frac{1}{2}\overline{h} \left( h + u^2 + \frac{1}{3}h^2 u_{\xi}^2 \right) \right)_t + \left( \overline{hu} \left( h + \frac{1}{2}u^2 + \frac{1}{2}h^2 u_{\xi}^2 - \frac{1}{3}h^2(u-c)u_{\xi\xi} \right) \right)_x = 0, \\ k_t + (kc)_x = 0. \end{gathered} \tag{25} $$

The averages in (25) are evaluated using the general definition

$$ \overline{f(h)} = \frac{\int_{h_2}^{h_3} \frac{f(h)dh}{\sqrt{P(h)}}}{\int_{h_2}^{h_3} \frac{dh}{\sqrt{P(h)}}}, \tag{26} $$

where $f(h)$ is any function of $h$.

The system (25) is consistent with the averaged Bernoulli conservation law (6):

$$ \left( \overline{u} + \frac{1}{6} \overline{h^2 u \xi \xi} \right)_t + \left( \frac{1}{2} \overline{u^2} + \overline{h} - \frac{1}{2} \overline{h^2 \left( \left( \frac{2}{3} u - c \right) u_{\xi\xi} - u_{\xi}^2 \right)} \right)_x = 0, \tag{27} $$

which can be used instead of any of the averaged conservation laws (25), e.g., instead of the wave conservation equation, yielding an equivalent modulation system. This equivalence is proved in the mass Lagrangian coordinates in Appendix A.

As a result, the system (25) can be explicitly represented in canonical quasilinear form for the state vector $b = (h_1(x, t), h_2(x, t), h_3(x, t), c(x, t))^T$; see Appendix B. However, due to the lack of integrable structure for the SGN equations, the associated modulation system cannot be reduced to a diagonal form. This makes the analysis of its properties—hyperbolicity, genuine nonlinearity, simple waves, etc.—difficult. The weakly nonlinear regime of the modulation equations was studied in [12]. The hyperbolicity of the modulation equations was proven in [54]. However, their full analytical study remains a difficult task. This is why we study here a reduction of the full modulation system in the limit of waves of large wavelength (the solitary wave limit), allowing us to present some analytical results. The difficulty of obtaining the solitary wave limit of the modulation equations is that the phase equation $k_t + (kc)_x = 0$ is degenerate in this singular limit. One possibility is to
pass to the limit in the wave action conservation law, an exact conservation law of the full SGN-Whitham system. Such a method was, in particular, exploited in [29] for the KdV equation, in [24] for the BBM equation and in [5] for Hamiltonian systems of Euler-Korteweg or dispersive Eulerian type (so-called second gradient fluids). The wave action integral for the SGN-Whitham equations was already obtained in [19] in a general setting that includes the SGN-Whitham equations by using mass Lagrangian coordinates. We use this result to first find the solitary wave limit in mass Lagrangian coordinates and then re-express the limit equations in Eulerian coordinates (see Appendix A). In the mass Lagrangian coordinates, one has only to replace the polynomial $P(h)$ by $P(h)/h^2$ when computing averages (see (22)):

$$ \frac{f(h)}{\int_{h_2}^{h_3} \frac{hf(h)dh}{\sqrt{P(h)}}} = \frac{\int_{h_2}^{h_3} \frac{hdh}{\sqrt{P(h)}}}{\int_{h_2}^{h_3} \frac{hdh}{\sqrt{P(h)}}} \quad (28) $$

To double-check our computation of the solitary wave limit in mass Lagrangian coordinates, we also compute the limit directly in Eulerian coordinates (see B) using a symbolic package.

A straightforward analysis shows that in the solitary wave limit, one has $\bar{h} = h_1$, $\bar{h}^2 = h_1^2 = \bar{h}^2$, $1/\bar{h} = 1/h_1 = 1/\bar{h}$, and the average of the mass and momentum equations in mass Lagrangian coordinates (10) are thus the classical Saint-Venant equations describing hydrostatic (dispersionless) shallow water equations

$$ \left( \frac{1}{\bar{h}} \right)_t - \bar{u}_q = 0, \quad \bar{u}_t - \left( \frac{\bar{h}^2}{2} \right)_q = 0. \quad (29) $$

This decoupling of the dispersionless limit equations for the evolution of the mean flow in the solitary wave limit is a general property of dispersive hydrodynamics [10, 31, 14]. It is also physically intuitive: the effect of a single soliton on a large-scale mean flow is negligibly small.

A non-trivial amplitude equation comes from the solitary wave limit of the wave action conservation equation. For general periodic travelling waves in mass Lagrangian coordinates, the wave action equation is (see Appendix A)

$$ \left( \tilde{c} \left( \frac{\bar{h}_{\zeta}^2}{3k} + \frac{\bar{\tau}^2 - \bar{\tau}^2}{k} \right) \right)_t + \left( \tilde{c}^2 \frac{\bar{\tau}^2 - \bar{\tau}^2}{k} \right)_q = 0. \quad (30) $$

The solitary wave limit ($k \to 0$) is singular, because $\bar{h}_{\zeta}^2 \to 0$ and $\bar{\tau}^2 - \bar{\tau}^2 \to 0$. In Appendix A, we show that equation (30) limits to

$$ F(n, \bar{h})_t + G(n, \bar{h})_q = 0, \quad (31) $$

with

$$ F(n, \bar{h}) = \frac{\bar{h}^3}{\sqrt{1-n}} \left( \frac{(6-2n)\sqrt{n}}{3(1-n)} + \ln \left( \frac{1-\sqrt{n}}{1+\sqrt{n}} \right) \right), \quad (32) $$

$$ G(n, h_1) = \sigma \frac{\bar{h}^3}{1-n} \left( \frac{\sqrt{n}}{1-n} + \frac{1}{2} \ln \left( \frac{1-\sqrt{n}}{1+\sqrt{n}} \right) \right). \quad (33) $$
The quantity $n$ is related to the solitary wave amplitude $a = h_3 - h_1$ as

$$n = \frac{a}{\bar{h} + a} \qquad (34)$$

The characteristic velocity corresponding to equation (31) is

$$\lambda = \frac{G_n}{F_n} = \sigma \frac{\bar{h}^{3/2}}{\sqrt{1-n}} = \sigma \bar{h} \sqrt{\bar{h} + a} = \tilde{c}_s. \qquad (35)$$

The velocity $\tilde{c}_s$ is related to the Eulerian velocity of solitary waves $c_s$ by (20).

In Eulerian coordinates, the system (29), (31) becomes

$$\bar{h}_t + (\bar{h}\bar{u})_x = 0, \qquad (36)$$

$$\bar{u}_t + \bar{u} \bar{u}_x + \bar{h}_x = 0, \qquad (37)$$

$$z_t + \left( \bar{u} + \sigma \sqrt{\bar{h}(1+z^2)} \right) z_x + \sigma \frac{3(z^2+1)^{3/2}}{2z} \frac{z\sqrt{z^2+1} - \sinh^{-1}(z)}{2z\sqrt{z^2+1} - \sinh^{-1}(z)} \frac{\bar{h}_x}{\sqrt{\bar{h}}} - \frac{\sqrt{z^2+1}}{2z} \frac{3z+2z^3-3\sqrt{z^2+1}\sinh^{-1}(z)}{2z\sqrt{z^2+1} - \sinh^{-1}(z)} \bar{u}_x = 0. \qquad (38)$$

where $z^2 = a/\bar{h}$ and $\sinh^{-1}(z) = \ln(z + \sqrt{z^2 + 1})$. We call the system (36)–(38) the solitonic modulation system. We stress that the solitonic modulation system is an *exact reduction* of the full SGN-Whitham modulation system.

The first two equations of the solitonic modulation system are the Saint-Venant (shallow water) equations written in Eulerian coordinates. They are completely decoupled from the third equation for $z$. Other than classical fast and slow surface waves of the shallow water equations with characteristic velocities $V_{\pm} = \bar{u} \pm \sqrt{\bar{h}}$, the characteristic velocity corresponding to the amplitude equation for the dimensionless amplitude $z^2$ also represents *fast* and *slow* solitary waves propagating with the third characteristic velocity $c_s = \bar{u} + \sigma \sqrt{\bar{h}(1+z^2)}$, the solitary wave velocity (20). By direct verification, we find that the system is strictly hyperbolic $V_- < V_+ < c$, if $\sigma = +1$, $\bar{h} > 0$, and $z > 0$ or $c < V_- < V_+$, if $\sigma = -1$, $\bar{h} > 0$, and $z > 0$ and the corresponding eigenfields are genuinely nonlinear in the sense of Lax.

The solitonic modulation system (36)–(38) inherits two Riemann invariants, those of the decoupled shallow water equations (36), (37). However, the third Riemann invariant does not exist, i.e. the system cannot be fully diagonalized (see Appendix C for the proof). Nevertheless, given the solution of (36) and (37) (obtained for example by the hodograph transform or a simple wave reduction), eq. (38) can then be integrated using the method of characteristics. Moreover, if one of the shallow water Riemann invariants is constant (the simple wave reduction), there are just two equations left, e.g. (37) and (38), and so the amplitude equation (38) can also be diagonalized in that case. The resulting extra Riemann invariant $Q(x, t)$ plays the role of an adiabatic invariant for solitary wave-mean flow interaction and determines the transmission (tunnelling) and trapping conditions for the interaction of a solitary wave with a RW or DSW mean flow generated by step initial data (Riemann) problems.
# 5 Solitary wave transmission and trapping

The expression “solitary wave tunneling” comes from quantum mechanics and conveys the possibility of a quantum particle passing through a classically impenetrable potential barrier. In the context of the SGN equations, the solitary wave plays the role of the quantum particle. More precisely, for a given solitary wave amplitude, we are able to use the solitary wave limit of the SGN-Whitham equations to analytically determine if the solitary wave passes (transmits or tunnels) through a RW or its counterpart, a DSW. Both a RW and a DSW represent a continuous transition between two constant mean flows and therefore solitary wave-mean flow transmission occurs when a solitary wave, initiated on one side of a RW or DSW, emerges on the other side. Otherwise, we say that the solitary wave has been trapped.

The system (36), (37) is independent of $z$ and can be diagonalized. The corresponding Riemann invariants are:

$$r_+ = \bar{u} + 2\sqrt{\bar{h}}, \quad r_- = \bar{u} - 2\sqrt{\bar{h}}, \qquad (39)$$

yielding

$$(r_+)_t + V_+(r_+)_x = 0, \quad V_+ = \bar{u} + \sqrt{\bar{h}} = \frac{3r_+ + r_-}{4}, \qquad (40)$$

$$(r_-)_t + V_-(r_-)_x = 0, \quad V_- = \bar{u} - \sqrt{\bar{h}} = \frac{r_+ + 3r_-}{4}. \qquad (41)$$

This system has two simple wave solutions:

$$r_{-\mu} = \bar{u} - 2\mu\sqrt{\bar{h}} = \text{cst}, \quad V_{\mu} = x/t, \qquad (42)$$

where $\mu = +1$ for the fast wave and $\mu = -1$ for the slow wave.

## 5.1 Exact simple wave Riemann invariants

To study the interaction of the solitary wave with the simple wave (42), we hold one of the Riemann invariants $r_{-\mu}$ constant, yielding a relation between $\bar{h}$ and $\bar{u}$ (42) [51]. The reduced solitonic modulation system is thus obtained by substituting (42) in (36), (38) (or equivalently (37), (38)), yielding the system of two equations:

$$\bar{h}_t + (r_{-\mu} + 3\mu\sqrt{\bar{h}})\bar{h}_x = 0, \qquad (43)$$

$$z_t + (r_{-\mu} + 2\mu\sqrt{\bar{h}} + \sigma\sqrt{\bar{h}(1+z^2)})z_x + \sigma g_{\sigma\mu}(z)\frac{\bar{h}_x}{\sqrt{\bar{h}}} = 0, \qquad (44)$$

with

$$g_{\sigma\mu}(z) = \frac{3(z^2+1)^{3/2}}{2z} \frac{z\sqrt{z^2+1} - \sinh^{-1}(z)}{2z\sqrt{z^2+1} - \sinh^{-1}(z)} - \sigma\mu \frac{1+z^2}{2z} \frac{(3z+2z^3)(1+z^2)^{-1/2} - 3\sinh^{-1}(z)}{2z\sqrt{1+z^2} - \sinh^{-1}(z)} \qquad (45)$$

Since now this is a quasi-linear system of just two equations, it is diagonalizable with one Riemann invariant manifestly $\bar{h}$. By making the simple-wave ansatz $z = z(\bar{h})$ in eqs. (43) and (44), the second Riemann invariant can be obtained as the constant of integration of the ODE

$$\frac{dz}{d\bar{h}} + \frac{g_{\sigma\mu}(z)}{\bar{h}(\sqrt{1+z^2} - \sigma\mu)} = 0. \qquad (46)$$
The integral of (46) can be written in the form

$$Q_{\sigma\mu}(\bar{h}, z) = \bar{h} \exp(f_{\sigma\mu}(z)), \quad f_{\sigma\mu}(z) = \int^{z} \frac{\sqrt{1+s^2} - \sigma\mu}{g_{\sigma\mu}(s)} ds. \qquad (47)$$

One can see that the Riemann invariant $Q_{\sigma\mu}$ (47) only depends on the sign $\sigma\mu$, which determines the interaction type. The value $\sigma\mu = +1$ corresponds to the overtaking interaction between a fast solitary wave and a fast simple wave (or equivalently slow-slow waves) whereas $\sigma\mu = -1$ corresponds to the head-on interaction between a fast solitary wave and a slow simple wave (or equivalently slow-fast waves).

Note that in the small amplitude limit $z \to 0$, we obtain the following expansion for the Riemann invariant in the case of an overtaking interaction ($\sigma\mu = +1$)

$$Q_{+}(\bar{h}, z) = \bar{h} \left( 1 + \frac{z^2}{2} - \frac{z^4}{48} - \frac{61}{1440} z^6 + O(z^8) \right). \qquad (48)$$

As we will now show, this expression is asymptotically equivalent, to order $O(z^4)$, to the Riemann invariant obtained by the DSW fitting method.

## 5.2 Approximate simple wave Riemann invariants: the DSW fitting method

An efficient approach to obtain the solitary wave limit of the Whitham equations can be deduced from the dispersive shock wave (DSW) fitting method [10]. This method enables the determination, under certain assumptions, of the harmonic and solitary wave edges of a DSW directly, bypassing the derivation and asymptotic analysis of the full Whitham modulation system. The DSW fitting method has been successfully applied to many dispersive hydrodynamic systems, both integrable and non-integrable, see e.g., [15, 12, 11, 16, 42, 31, 9, 33]. In cases when the dispersive equation is integrable such as in the KdV, NLS, and Kaup-Boussinesq equations, the method's results are consistent with the available exact modulation solutions. For non-integrable systems, when the exact theory is not available, the method yields an excellent comparison with direct numerical simulations of the dispersive equation, beyond the classical weakly nonlinear KdV approximation. Here, we shall use the part of the DSW fitting construction pertaining to the determination of the DSW's solitary wave edge.

A major assumption in the DSW fitting method is that of convexity (strict hyperbolicity and genuine nonlinearity) of the full Whitham modulation system. This assumption is not always easily verifiable for the entire range of parameters involved, so one can say that DSW fitting generally provides a convex approximation of the exact DSW edge dynamics. In particular, it is known that the SGN-Whitham modulation system exhibits non-convexity for sufficiently large amplitudes [12, 31] but nevertheless, the DSW fitting results agree with direct SGN numerics remarkably well for a broad range of amplitudes, far beyond the weakly nonlinear KdV regime.

The DSW fitting method relies on the existence of exact reductions of the Whitham modulation system in two distinguished limits: the harmonic limit $a \to 0$ and the solitary wave limit $k \to 0$. The latter one is our interest here and has the general form of an equation for the solitary wave amplitude (or a related amplitude type variable), coupled to the dispersionless limit system for the large-scale background (the mean flow) [14]. In the present context of the SGN-Whitham equations, the solitary wave reduction is given by the system (36)–(38), however, the derivation of the exact
form of the amplitude equation (38) and further, of the Riemann invariant $Q_{\sigma\mu}$ in eq. (47)—the main technical hurdle in the analysis of the solitary wave limit of the modulation system—is not required for the application of the DSW fitting method. Instead, the determination of the requisite Riemann invariant involves only the linear dispersion relation $\omega_0(k, \bar{h}, \bar{u})$ for waves on the mean background $(\bar{h}, \bar{u})$ along with the simple-wave relation $\bar{u}(\bar{h})$ for the dispersionless limit equations. All of this information is readily available with no additional analysis of the SGN-Whitham modulation system required. In the context of soliton-mean flow interaction, the Riemann invariant $Q$ acquires the meaning of an adiabatic invariant [44]. This adiabatic invariant also plays a key role in the solitary wave resolution method [13, 43], where it is an analogue of the spectral parameter in the Lax pair associated with integrable equations in the semi-classical limit. See also the recent papers [36, 35] where soliton-mean field interaction has been interpreted in terms of Hamiltonian mechanics of the motion of a soliton along a large-scale background.

Note that the DSW fitting construction applies directly to the overtaking solitary wave-mean field interaction, not the head-on interaction, due to the unidirectional nature of DSW generation in Riemann problems. The DSW fitting method was applied to undular bore theory for the SGN system in [12], so we shall take advantage of the results of that work and adapt them to the current setting of the overtaking solitary wave-mean flow interaction. To be definite, we consider the fast-fast interaction, $\sigma = \mu = +1$.

For the fast simple-wave mean field we have (cf. (43))

$$ \bar{u}(\bar{h}) = r_{-} + 2\sqrt{\bar{h}}, \quad \bar{h}_{t} + V_{+}(\bar{h})\bar{h}_{x} = 0, \quad V_{+} = r_{-} + 3\sqrt{\bar{h}}, \qquad (49) $$

where $r_{-} = \bar{u} - 2\sqrt{\bar{h}} = \text{const}$ is the Riemann invariant of the dispersionless shallow-water equations. The fast branch of the dispersion relation (21) for linearised waves on the simple-wave background (49) is given by

$$ \omega_0(k, \bar{h}, \bar{u}(\bar{h})) = k \left( r_{-} + 2\sqrt{\bar{h}} \right) + k \sqrt{\frac{\bar{h}}{1 + \bar{h}^2 k^2 / 3}}. \qquad (50) $$

Solitary wave motion on the simple-wave background is defined by the conjugate dispersion relation (see [10] for details)

$$ \tilde{\omega}(\tilde{k}, \bar{h}) = -i\omega_0(\tilde{k}, \bar{h}, \bar{u}(\bar{h})) = \tilde{k}(r_{-} + 2\sqrt{\bar{h}}) + \tilde{k}\sqrt{\frac{\bar{h}}{1 - \bar{h}^2 \tilde{k}^2 / 3}}, \qquad (51) $$

where the conjugate wavenumber $\tilde{k}$ (essentially the inverse width of the solitary wave) is related to the solitary wave amplitude $a$ by

$$ \frac{\tilde{\omega}}{\tilde{k}} = c_s(a, \bar{h}, \bar{u}(\bar{h})), \qquad (52) $$

with $c_s(a, \bar{h}, \bar{u}) = \bar{u} + \sqrt{\bar{h} + a}$ the fast branch of the SGN solitary wave speed-amplitude relation (20).

Within the fitting approach to solitary wave-mean flow interaction, the amplitude modulation equation is obtained directly in diagonal form. The corresponding Riemann invariant $\tilde{Q}_+(\tilde{k}, \bar{h})$ satisfying

$$ \frac{\partial \tilde{Q}_+}{\partial t} + C(\tilde{Q}_+, \bar{h}) \frac{\partial \tilde{Q}_+}{\partial x} = 0, \quad C(\tilde{Q}_+, \bar{h}) = c_s(a(\tilde{k}, \bar{h}), \bar{h}, \bar{u}(\bar{h})) \qquad (53) $$
<figure><img src="image_4.png" /><figcaption>Figure 3: Left: comparison between the exact Riemann invariant $Q_+(47)$ (solid line) and its convex approximation $\tilde{Q}_+(57)$ (dashed line); Right: the difference $\Delta = (Q_+ - \tilde{Q}_+)/\bar{h}$.</figcaption></figure>

is found as an integral of the characteristic ODE

$$ \frac{d\tilde{k}}{d\bar{h}} = \frac{\tilde{\omega}_{\bar{h}}}{V_+(\bar{h}) - \tilde{\omega}_{\tilde{k}}} = \frac{\tilde{k}\left(\left(1 + \frac{\bar{h}^2\tilde{k}^2}{3}\right) + 2\left(1 - \frac{\bar{h}^2\tilde{k}^2}{3}\right)^{3/2} + 1\right)}{2\bar{h}\left(\left(1 - \frac{\bar{h}^2\tilde{k}^2}{3}\right)^{3/2} - 1\right)}. \quad (54) $$

Equation (54) is a convex approximation counterpart to the ODE (46) obtained by the exact evaluation of the solitary wave limit of the SGN-Whitham modulation system.

With the change of variable

$$ \alpha = \frac{1}{\sqrt{1 - \frac{\bar{h}^2 \tilde{k}^2}{3}}} = \sqrt{1 + \frac{a}{\bar{h}}} = \sqrt{1 + z^2}, \quad (55) $$

eq. (54) reduces to the separable ODE

$$ \frac{d\alpha}{d\bar{h}} = \frac{\alpha(\alpha^2 - 1)(4 - \alpha)}{2\bar{h}(1 - \alpha^3)}, \quad (56) $$

which is readily integrated with the constant of integration

$$ \tilde{Q}_+ = \frac{2^{2/5}3^{21/10}\bar{h}\alpha^{1/2}}{(\alpha + 1)^{2/5}(4 - \alpha)^{21/10}}. \quad (57) $$

The normalization coefficient in (57) is chosen by the natural requirement that in the zero-amplitude/infinite width limit $a \to 0$ of the solitary wave (19), the Riemann invariant $\tilde{Q}_+ \to \bar{h}$ so that equation (53) degenerates into the simple wave mean flow equation (49) for $\bar{h}$. Equation (57) is equivalent to eq. (49) in [12].

One can see that the expression (57) does not coincide with the rigorous asymptotic result (47) for $Q_+$ corresponding to $\mu = \sigma = 1$. However, it provides a quite accurate approximation in the physically relevant range of amplitude-depth ratios $z^2 = a/\bar{h}$. Indeed, the small amplitude expansion of the fitting approach result (57) is

$$ \tilde{Q}_+ = \bar{h} \left( 1 + \frac{z^2}{2} - \frac{z^4}{48} - \frac{37}{864} z^6 + \mathcal{O}(z^8) \right). \quad (58) $$
The first two terms in this expansion correspond to the Riemann invariant of the KdV modulation system in the soliton limit [14], i.e., is accurate to $O(a)$ ($O(z^2)$) as $a \to 0$. But the fitting approach improves upon the classical KdV prediction since (58) and the exact result (48) asymptotically agree up to $O(a^2)$ ($O(z^4)$). Even further, the difference between the two expansions at $O(a^3)$ ($O(z^6)$) is just $\sim 10^{-3}a^3$, i.e. the exact (47) and convex approximation (57) curves are in excellent agreement for physically admissible $a/\bar{h} < a_{\text{max}}/\bar{h} \approx 0.8$ [12], see Fig. 3. At the same time, the existence of such a discrepancy poses the important question: what is the cause of this $O(a^3)$ asymptotic discrepancy and what does it mean for the asymptotic validity of the DSW fitting method, particularly for non-integrable systems?

## 5.3 Riemann problem and transmission condition

We consider the Riemann problem for the system (36), (37), (38), for which the initial condition is a step function:

$$ (\bar{h}(x, 0), \bar{u}(x, 0), z(x, 0)) = \begin{cases} \left(h_{-}, u_{-}, z_{-} \equiv \sqrt{\frac{a_{-}}{h_{-}}}\right), & x < 0, \\ \left(h_{+}, u_{+}, z_{+} \equiv \sqrt{\frac{a_{+}}{h_{+}}}\right), & x > 0, \end{cases} \qquad (59) $$

with the constraint $u_{-} - 2\mu\sqrt{h_{-}} = u_{+} - 2\mu\sqrt{h_{+}}$. The solution of interest is the simple wave solution

$$ Q_{\sigma\mu}(\bar{h}, z) = Q_{\sigma\mu}(h_{-}, z_{-}) = Q_{\sigma\mu}(h_{+}, z_{+}), \quad V_{\mu} = r_{-\mu} + 3\mu\sqrt{\bar{h}} = x/t. \qquad (60) $$

Equation (60) implicitly determines $\bar{h}$ and $z^2 = a/\bar{h}$ as functions of $x/t$. It describes the fast RW with $h_{-} < h_{+}$. To simplify the discussion of the solution, we restrict our study to fast solitary waves ($\sigma = +1$). In this case, the Riemann problem for the modulation system models the interaction of an incident solitary wave with parameter $z_{-}$, initially located far to the left of the initial step $x < 0$, interacting with a RW or a DSW generated by the initial step.

It is important to note here that, although the solitary wave limit of the SGN-Whitham equations (36)–(38) do not admit DSW modulation solutions, they can be used to determine the nature of the interaction between a solitary wave and a DSW. The only region of space-time where equations (36)–(38) break down is within the DSW itself. Outside of the DSW, equations (36)–(38) are perfectly valid, so that the simple wave solution (60) applies outside of the DSW. If all we wish to determine is whether the solitary wave has been transmitted or trapped by the DSW and, if transmitted, the transmitted solitary wave amplitude, we can use the Riemann invariant relation from (60) with $h_{-} > h_{+}$ to ascertain this as described below. This concept of hydrodynamic reciprocity was first theoretically predicted and simultaneously experimentally observed in the interfacial fluid dynamics of a viscous fluid conduit [44]. The concept has since been applied to solitary wave-DSW interaction for other equations [51, 55, 49, 3]. In order to describe the interaction of the solitary wave within the DSW, one must appeal to multi-phase Whitham modulation theory as has been carried out for the KdV equation [3].

The first equation in (60) yields a relation between the left and right parameters ($h_{\pm}$, $z_{\pm}$) of the initial condition, i.e., a relation between the parameters of the incident and transmitted waves

$$ f_{\sigma\mu}(z_{+}) + \ln h_{+} = f_{\sigma\mu}(z_{-}) + \ln h_{-}, \quad z_{\pm}^{2} = \frac{a_{\pm}}{h_{\pm}}. \qquad (61) $$
If transmitted, the solitary wave amplitude $a_+ > 0$ is determined by solving for $z_+$ in eq. (61).

In the interaction with a fast simple wave, where $\sigma\mu = +1$, one can choose the constant of integration in (47) such that $f_+(0) = 0$. In that case, $f_+(z) \ge 0$ (see Fig. 4(a)) and one can find $z_+$ from Eq. (61) if

$$f_+(z_-) + \ln\left(\frac{h_-}{h_+}\right) > 0. \qquad (62)$$

This condition is called the transmission condition. It is always fulfilled in the interaction with a fast DSW where $h_- > h_+$. In the interaction with a fast RW, this condition can be written in the form:

$$z_- > z_{\min}(h_+/h_-) = f_+^{-1}\left(\ln\left(\frac{h_+}{h_-}\right)\right), \qquad (63)$$

plotted in Fig. 4(b)—solid line. If the normalized amplitude of the incident solitary wave is smaller than $z_{\min}$, then the wave becomes trapped inside the RW, similar to [44] in the case of the KdV and conduit equations, see Fig. 1(b).

If the approximate Riemann invariant $\tilde{Q}_+$ (57) is used in the transmission relation (61) then the fitting counterpart of the exact condition (63) assumes the form

$$z_- > z_{\min}(h_+/h_-) = \tilde{f}_+^{-1}\left(\frac{h_+}{h_-}\right), \qquad (64)$$

where

$$\tilde{f}_+(z) = \frac{2^{2/5}3^{21/10}\alpha^{1/2}}{(\alpha+1)^{2/5}(4-\alpha)^{21/10}}, \quad \alpha = \sqrt{1+z^2}. \qquad (65)$$

The curve (64) is plotted in Fig. 4(b) in dashed line. There is a very close agreement between the exact and convex approximation curves.

<figure><img src="image_5.png" /><figcaption>Figure 4: a) Variation of $f_+(z)$. b) The critical transmission curve $z_{\min}^2(h_+/h_-)$ for the overtaking interaction between a fast solitary wave and a fast RW ($h_- < h_+$) given by the exact formula (63) (continuous line) and the convex approximation formula (64) (dashed line). c) Relation between $z_+$ and $z_-$ for the interaction with a fast RW, ($h_-, h_+$) = (1, 1.5) (solid line) and a fast DSW, ($h_-, h_+$) = (1.5, 1) (dashed line).</figcaption></figure>

If the fast solitary wave is placed in front of the fast DSW two scenarios are possible depending on the solitary wave amplitude $z_+$ relative to the DSW leading edge amplitude $z_+$ (see Equation (66) below). If $z_+ > z_+$ then the solitary wave propagates faster than the DSW and no interaction
occurs. If $z_+ < z_+^*$ then the DSW overtakes the solitary wave and the latter gets trapped inside the DSW, see Fig. 1(d) for the illustration. The analysis of the trapping dynamics characterised by the formation of travelling breathers [3, 46] is beyond the scope of this paper.

In the head on interaction with a slow simple wave, where $\sigma\mu = -1$, $f_-(z)$ is singular at $z = 0$ (see Fig. 5(a) where $f_-(1) = 0$), and the range of $f_-(z)$ spans the whole real axis. One can always find $z_-$ from Eq. (61). The transmission condition is always satisfied in this case, i.e., the incident fast solitary wave is never trapped by the RW or DSW.

<figure><img src="image_6.png" /><figcaption>Figure 5: a) Variation of $f_-(z)$. b) Relation between $z_+^2$ and $z_-^2$ for the head on interaction of a fast solitary wave with a slow DSW, ($h_-, h_+$) = (1, 1.5) (solid line) and a slow RW, ($h_-, h_+$) = (1.5, 1) (dashed line).</figcaption></figure>

## 5.4 DSW: solitary wave edge

The transmission condition (61) can also be used to approximate the SGN DSW’s solitary wave leading edge amplitude and speed. When $h_- > h_+$, a DSW is generated by the simple wave Riemann problem (59). The DSW’s solitary wave leading edge amplitude $a_+ = h_+(z_+^*)^2$ can be estimated by evaluating the Riemann invariant at the transmission/trapping bifurcation point when $z_- = 0$. In other words, a fast solitary wave of infinitesimally small amplitude gets transmitted exactly to the DSW leading edge. As follows from the transmission relation (61) this corresponds to

$$z_+^* = f_+^{-1} \left( \ln \left( \frac{h_-}{h_+} \right) \right). \qquad (66)$$

The DSW solitary wave amplitude is $a_+ = h_+(z_+^*)^2$ and edge speed is therefore $c_s(a_+, h_+, u_+) = u_+ + \sqrt{h_+ + a_+}$.

The SGN DSW leading solitary wave amplitude prediction from eq. (66) is an approximation of the amplitude that results from integrating the appropriate integral curve (in this case, for $h_- > h_+$, a 3-wave) of the full SGN-Whitham equations. The approximation rests upon the hypothesis of hydrodynamic reciprocity—in which the solitary wave Riemann invariant $Q_{+1}(\bar{h}, z)$ (47) is the same behind and ahead of a DSW [44]—and a more subtle, convexity condition that the full Whitham modulation system remains strictly hyperbolic along the entire solitary wave trajectory for all so-called admissible states; see, Section 3.2 in reference [55] for a discussion. Although these conditions are exactly true, within the context of Whitham modulation theory, for a fast (slow) solitary wave interacting with a fast (slow) rarefaction wave mean-field, they are assumed to be true for the fast (slow) DSW mean-field. Because of this, we refer to eq. (66) as the mean-field DSW approximation.
<figure><img src="image_7.png" /><figcaption>Figure 6: DSW lead solitary wave amplitude. Solid line: mean-field DSW solution (66); dashed line: DSW fitting formula (67); circles: SGN numerical simulations.</figcaption></figure>

If instead of the exact Riemann invariant (47), its convex (DSW fitting) approximation (57) is used in the transmission condition (60), the DSW leading edge amplitude $\tilde{a}_+ = h_+(\tilde{z}_+^*)^2$ is found by solving the equation [12]

$$ \frac{h_{-}/h_{+}}{((\tilde{z}_{+}^{*})^{2} + 1)^{1/4}} - \left( \frac{3}{4 - \sqrt{(\tilde{z}_{+}^{*})^{2} + 1}} \right)^{21/10} \left( \frac{2}{1 + \sqrt{(\tilde{z}_{+}^{*})^{2} + 1}} \right)^{2/5} = 0. \quad (67) $$

Equation (67) is subject to the admissibility conditions—the convexity restrictions ensuring mono-tone behaviour of the DSW edge speeds as functions of $h_-, h_+$ [31]. This yields the admissible range of initial steps $1 < h_- / h_+ \lesssim 1.43$ for which the SGN DSW fitting results are applicable [12].

The small-jump expansions of (66) and (67) give

$$ \frac{a_{+}}{h_{+}} = 2\delta + \frac{1}{6}\delta^{2} + \frac{127}{180}\delta^{3} + \mathcal{O}(\delta^{4}), \quad (68) $$

$$ \frac{\tilde{a}_{+}}{h_{+}} = 2\delta + \frac{1}{6}\delta^{2} - \frac{71}{108}\delta^{3} + \mathcal{O}(\delta^{4}), \quad (69) $$
respectively, where $\delta = \frac{h_-}{h_+} - 1 \ll 1$. To leading order, both results agree with the famous result for the KdV equation that the lead solitary wave amplitude is twice the initial DSW jump height [30]. The mean-field DSW and DSW fitting approximations agree at second order as well, consistent with the agreement between the expansions of the exact (48) and convex approximation (58) for the solitary wave Riemann invariants. The curves $a_+(h_-/h_+)$ given by (66), (67) along with the values of $a_+$ extracted from the numerical simulations are shown in Fig. 6. One can see that both the mean-field DSW equation (66) and the DSW fitting formula (67) provide very good approximations of the SGN DSW solitary wave edge amplitude within the admissible range of initial steps. The numerical results in Fig. 6 suggest that the mean-field DSW prediction somewhat improves upon the DSW fitting result.

## 5.5 Transmission: comparison between analytical and numerical results

Figure 7 shows the interaction of a fast solitary wave with fast and slow rarefaction waves. The amplitude of the transmitted solitary wave $a_+$ obtained by the direct numerical simulation of the SGN equations (1), (2) compares well with the analytical prediction (61). The details of the numerical scheme are given in Appendix D. In the case of a fast solitary wave with a slow rarefaction wave, the amplitude of the transmitted solitary wave is always larger than that of the incident solitary wave. Physically speaking, the slow rarefaction wave is formed by retracting a piston to the right, so it will accelerate the right-going solitary wave (see Figure 7). A new feature also appears in this case: a small amplitude solitary wave forms behind the transmitted large amplitude solitary wave. This is in contrast with the typical radiation emitted when a solitary wave interacts with another wave in non-integrable systems. Despite the generation of additional solitary waves in this Riemann problem due to the non-integrable nature of the SGN equations, the accuracy of the analytical predictions from Whitham theory is remarkable.

<figure><img src="image_8.png" /><figcaption>Figure 7: Left figure: comparison between the transmission relation (61) and the numerical simulation for the fast solitary wave interaction with a fast RW (dashed red line) and a slow rarefaction wave (solid blue line). Middle figure: an example of the interaction with a fast RW. Right figure: an example of the interaction with a slow RW.</figcaption></figure>

Figure 8 shows the interaction of a solitary wave with fast and slow DSWs. The analytical formulas are in agreement with direct numerical simulations. This time, small amplitude radiation is generated after the solitary wave passes through the slow DSW.
<figure><img src="image_9.png" /><figcaption>Figure 8: Left figure: comparison between the transmission relation (61) and the numerical simulation for the fast solitary wave interaction with a fast DSW (dashed red line) and a slow DSW (solid blue line). Middle figure: an example of the interaction with a fast DSW. Right figure: an example of the interaction with a slow DSW.</figcaption></figure>

# 6 Conclusion

We have studied the interaction of a solitary wave with a slowly varying mean flow for the Serre-Green-Naghdi (SGN) equations that model fully nonlinear, bidirectional shallow water gravity waves over a flat bottom. This was achieved by the determination and analysis of the exact solitary wave limit of the Whitham modulation equations for the SGN system. The derivation of the Whitham equations was performed via the averaged conservation law approach in both Eulerian and Lagrangian coordinates. The SGN system is not integrable so the SGN-Whitham equations cannot be diagonalized, which makes their analysis a nontrivial problem.

Due to the singular nature of the solitary wave limit, the appropriate choice for the modulation parameters is crucial. We utilize the wave action density, which is shown to be a particularly well-behaved conserved quantity in the solitary wave limit. The resulting solitonic modulation system consists of the shallow water equations for the mean flow coupled to an amplitude equation for the solitary wave. Although the SGN solitonic system is not diagonalizable, its restriction to simple waves for the mean flow equations admits Riemann invariants that we use to analytically describe the head-on and overtaking interactions of a solitary wave with a rarefaction wave and dispersive shock wave (DSW). These scenarios lead to solitary wave trapping or transmission by the mean flow. The analytical results are shown to be in excellent agreement with corresponding numerical solutions of the full SGN equations. This work extends previous results on solitary wave-mean flow interaction in unidirectional, non-integrable systems (BBM and conduit equations) to the bidirectional case.

One important outcome of this paper is the comparison between the exact analytical results from SGN modulation theory for overtaking interactions with the results obtained from a simpler approach that is based on the DSW fitting method. The latter method, involving several major assumptions of the full modulation system (strict hyperbolicity, convexity), has been successfully applied to many dispersive hydrodynamic equations. It yields exact results for integrable equa-
tions, relies only upon knowledge of the linear dispersion relation and the dispersionless, mean-flow equations, and provides consistently good comparisons with direct numerical simulations for non-integrable equations. However, analytical estimates of DSW fitting theory accuracy have not been available. In this paper, we show that the DSW fitting results for the SGN equation are not exact but are accurate to the second order in the solitary wave amplitude, beyond the first order accurate Korteweg-de Vries approximation. This comparison between the exact and fitting results for the SGN-Whitham equations provides an important asymptotic benchmark for the DSW fitting method in non-integrable equations.

This work paves the way for additional studies of solitary wave interactions with mean flows in bidirectional, strongly nonlinear, physically relevant systems. One of the future directions could be the generalization of the developed theory for bidirectional solitary wave-mean flow interactions to soliton gases—large, random ensembles of interacting solitary waves [53].

**Acknowledgment** TC, GE, SG and MH would like to thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the satellite program “Emergent phenomena in nonlinear dispersive waves” at Northumbria University, Newcastle, UK where work on this paper was undertaken. This work was supported by EPSRC grant no EP/R014604/1 and NSF grant no DMS-2306319. KMS thanks the Institut Mécanique et Ingénierie of Marseille and INRIA-Sophia-Antipolis-Méditerranée for financial support during his stay at IUSTI, Marseille. The authors thank M. Rodrigues for useful discussions.

## A Lagrangian and Eulerian descriptions

If $t$ is time, and $q$ is the mass Lagrangian coordinate, the spatial Eulerian coordinate is defined through the motion of the continuum $x = x(t, q)$. The mass conservation equation in the Lagrangian coordinates can be written as $\rho(t, q)x_q = \rho_0(q)$, where $\rho(t, q)$ and $\rho_0(q)$ are the actual and the reference mass densities, so that $\rho(t, q)dx = \rho_0(q)dq$. One can see that if we choose $\rho_0(q) \equiv 1$, the Lagrangian variable $q$ will effectively coincide with the mass (one-dimensional). A general conservation law in Lagrangian coordinates $(t, q)$

$$a_t + b_q = 0,$$

is expressed in Eulerian coordinates $(t, x)$ as

$$(\rho a)_t + (\rho u a + b)_x = 0.$$

When looking for traveling wave solutions that depend on $\tilde{\xi} = q - \tilde{c}t$ in Lagrangian coordinates, or $\xi = x - ct$ in Eulerian coordinates, the velocities $\tilde{c}$ and $c$ are related by the mass conservation law

$$\rho(u - c) = -\tilde{c}.$$

We will first pass to the solitary wave limit in the modulation equations expressed in terms of mass Lagrangian coordinates and then we transform the limit equations to Eulerian form.

The governing Euler-Lagrange equations derived from the Lagrangian (8) are

$$\tau_t - u_q = 0, \quad u_t + p_q = 0, \tag{70}$$
where $p$ is given in (11). Here, $\tau = 1/h$ is the specific volume, $u$ is the velocity, and $\tilde{e}(\tau, \tau_t)$ is the potential (9). The system admits the energy equation (12) and Bernoulli equation (13). Following Whitham [57], we are looking for the solution of (70) subject to (11) in the form

$$ \begin{pmatrix} \tau \\ u \end{pmatrix} = \begin{pmatrix} \tau(T, X, \theta, \epsilon) \\ u(T, X, \theta, \epsilon) \end{pmatrix}, \quad T = \epsilon t, \quad X = \epsilon q, \quad \theta = \frac{\Theta(T, X)}{\epsilon}, $$

where $0 < \epsilon \ll 1$ is a small, positive parameter. The function $\Theta(T, X)$ is called the phase. We define the local wave number $k$ and the local frequency $\omega$ by the relations:

$$ k = \Theta_X, \quad \omega = -\Theta_T. \tag{71} $$

We will suppose that the solution is $2\pi$-periodic with respect to the rapid variable $\theta$. Below, for any function $f(T, X, \theta, \epsilon)$, $\bar{f}$ is the $2\pi$-period-average with respect to $\theta$

$$ \bar{f}(T, X, \epsilon) = \frac{1}{2\pi} \int_{0}^{2\pi} f(T, X, \theta, \epsilon) \, d\theta. \tag{72} $$

Expanding $f$ in an asymptotic series in $\epsilon$

$$ f = f_0(T, X, \theta) + \epsilon f_1(T, X, \theta) + \dots, $$

and substituting this ansatz into the governing system (70), we obtain a system of ordinary differential equations with respect to $\theta$ at leading order. One has the following first integrals (the zero subscript is omitted)

$$ \tilde{c}\tau + u = \tilde{c}\bar{\tau} + \bar{u}, \quad -\tilde{c}u + p = -\tilde{c}\bar{u} + \bar{p}, \quad \tilde{c} = \frac{\omega}{k}. \tag{73} $$

The system (73) admits the following useful consequences

$$ \overline{u^2} - (\overline{u})^2 = \tilde{c}^2(\overline{\tau^2} - (\overline{\tau})^2), \tag{74} $$

$$ \overline{u\tau} - \overline{u}\overline{\tau} = -\tilde{c}(\overline{\tau^2} - (\overline{\tau})^2), \tag{75} $$

$$ \overline{pu} - \overline{p}\overline{u} = \tilde{c}^3(\overline{\tau^2} - (\overline{\tau})^2), \tag{76} $$

$$ \overline{p\tau} - \overline{p}\overline{\tau} = -\tilde{c}^2(\overline{\tau^2} - (\overline{\tau})^2). \tag{77} $$

To leading order, $O(1)$ as $\epsilon \to 0$, one has

$$ \tilde{e}(\tau, \tau_t) \approx \tilde{e}(\tau, \eta), \quad \eta = -\omega\tau_\theta, \quad p = -(\tilde{e}_\tau + \omega(\tilde{e}_\eta)_\theta). \tag{78} $$

The system (73) also admits the first integral in the form:

$$ -\tilde{c} \left( \tilde{e} - \eta \tilde{e}_\eta + \frac{u^2}{2} \right) + pu = -\tilde{c} \left( \overline{\tilde{e} - \eta \tilde{e}_\eta + \frac{u^2}{2}} \right) + \overline{pu}. \tag{79} $$

The integral (79) allows us to obtain the nonlinear dispersion relation because the solution is $2\pi$-periodic. At $O(\epsilon)$, after averaging with respect to the rapid variable $\theta$, one obtains the following
systems of five compatible conservation laws containing only the leading order terms (the zero subscript is again omitted):

$$
\begin{align*}
\bar{\tau}_T - \bar{u}_X &= 0, \\
\bar{u}_T + \bar{p}_X &= 0, \\
\left( \frac{u^2}{2} + \varepsilon \right)_T + (\bar{p}u)_X &= 0, \quad \varepsilon = \tilde{e} - \eta \tilde{e}_\eta, \\
\left( \bar{\tau}u - k\bar{\tau}_\theta \frac{\partial \tilde{e}}{\partial \eta} \right)_T - \left( \frac{u^2}{2} - \bar{\tau}p - \tilde{e} \right)_X &= 0, \\
k_T + \omega_X &= 0.
\end{align*}
$$

The last equation is just the compatibility condition coming from the definition of the phase (71).
It can be also written in the form

$$
k_T + (\tilde{c}k)_X = 0, \quad \tilde{c} = \frac{\omega}{k}.
$$

Introducing

$$
\Delta = \overline{u\tau} - \overline{u}\tau, \quad \bar{\varepsilon} = \overline{\tilde{e}} - \eta\overline{\tilde{e}}_{\eta}, \quad E = \bar{\varepsilon} + \frac{\overline{u^2} - \overline{u}^2}{2}, \quad \Sigma = \overline{\tau_{\theta}\tilde{e}}_{\eta},
$$

one can obtain a Gibbs-type identity relating the unknowns (see [19, 25] for proof)

$$
dE + \bar{p} \, d\bar{\tau} + \tilde{c} d\Delta = \omega d\Sigma. \tag{80}
$$

It can also be written in the form:

$$
d\bar{\varepsilon} + \bar{p} \, d\bar{\tau} - \frac{\bar{c}^2}{2} d\delta = \omega d\Sigma, \quad \delta = \bar{\tau}^2 - (\bar{\tau})^2. \tag{81}
$$

The Gibbs identity (80) or (81) is equivalent to an algebraic nonlinear dispersion relation coming
from (79) when we are looking for 2π-periodic solutions [57].

Using the expressions for correlations (74)–(77), the modulation equations take the conservative
form

$$
\begin{align}
\bar{\tau}_T - \bar{u}_X &= 0, \notag \\
\bar{u}_T + \bar{p}_X &= 0, \notag \\
\left(\bar{\varepsilon} + \frac{\bar{u}^2}{2} + \frac{\bar{c}^2}{2}(\bar{\tau}^2 - \bar{\tau}^2)\right)_T + \left(\bar{p}\bar{u} + \bar{c}^3(\bar{\tau}^2 - \bar{\tau}^2)\right)_X &= 0, \tag{82} \\
\left(\bar{\tau}\bar{u} - \bar{c}(\bar{\tau}^2 - \bar{\tau}^2) - k\Sigma\right)_T - \left(\frac{\bar{u}^2}{2} + \frac{\bar{c}^2}{2}(\bar{\tau}^2 - \bar{\tau}^2) - \bar{\tau}\bar{p} + \bar{c}^2(\bar{\tau}^2 - \bar{\tau}^2) - \bar{\varepsilon} + \omega\Sigma\right)_X &= 0, \notag \\
k_T + (\bar{c}k)_X &= 0. \notag
\end{align}
$$

The five conservation laws (82) form a system of modulation equations (the SGN-Whitham system) for the four unknowns *k*, *c*, *τ*, *u*. The equations are compatible because the averaged Bernoulli law is a consequence of the mass, momentum, energy and wave conservation laws. This can be proved
by direct calculations. Using the Gibbs relation (80) or (81), one can derive the sixth conservation law for Σ (an averaged entropy conservation law)

$$ \left( \Sigma + \frac{\tilde{c}\delta}{k} \right)_T + \left( \frac{\tilde{c}^2\delta}{k} \right)_X = 0. \qquad (83) $$

For the proof, one can use the energy conservation law in the form

$$ \left( \bar{\epsilon} + \frac{\bar{u}^2}{2} + \frac{\tilde{c}^2}{2} \delta \right)_T + \left( \bar{p} \bar{u} + \tilde{c}^3 \delta \right)_X = 0. $$

Expanding it and using Gibbs relation in the form (81), one obtains

$$ -\bar{p} \bar{\tau}_T + \frac{\tilde{c}^2}{2} \delta_T + \omega \Sigma_T + \bar{u} \bar{u}_T + \frac{\tilde{c}^2}{2} \delta_T + \tilde{c} \tilde{c}_T \delta + \bar{p}_X \bar{u} + \bar{p} \bar{u}_X + \tilde{c}(\tilde{c}^2 \delta)_X + \tilde{c}^2 \tilde{c}_X \delta = 0. $$

Using the averaged mass and momentum equations, one obtains

$$ \tilde{c}^2 \delta_T + \omega \Sigma_T + \tilde{c} \tilde{c}_T \delta + \tilde{c}(\tilde{c}^2 \delta)_X + \tilde{c}^2 \delta \tilde{c}_X = 0. $$

Or, dividing by $\tilde{c}$

$$ (\tilde{c}\delta)_T + k\Sigma_T + (\tilde{c}^2\delta)_X + \tilde{c}\delta\tilde{c}_X = 0. $$

Dividing by $k$, one obtains, after some algebra, the conservation law (83).

Finally, the full, compatible SGN-Whitham system (82)-(83) consists of the averaged conservation laws of mass, momentum, energy, wave conservation, the averaged Bernoulli conservation law and that for the averaged entropy $\Sigma$. The last conservation law is also called the wave action conservation law [57].

## A.1 Solitary wave limit

The solitary wave limit is $k \to 0$, which implies $\omega \to 0$. In this case, the modulation equations become purely hydrodynamic. For example, the mass and momentum equation become

$$ \bar{\tau}_T - \bar{u}_X = 0, \quad \bar{u}_T + \bar{p}_X = 0, $$

with

$$ \bar{p} = -\tilde{e}_{\tau}(\bar{\tau}, 0). \qquad (84) $$

The conservation of waves equation identically vanishes because both $k$ and $\omega$ are zero in this limit. However, in such a limit, the equation for the wave action is non-trivial. Let us first calculate $\Sigma$ for the full (non-averaged) SGN equations. With $\tau = \frac{1}{h}$, to leading order $O(1)$, one has

$$ \tilde{e} = \frac{h}{2} - \frac{h_t^2}{6} = \frac{1}{2\tau} - \frac{\tau_t^2}{6\tau^4} \sim \frac{1}{2\tau} - \frac{\eta^2}{6\tau^4}, \quad \eta = -\omega\tau_\theta. \qquad (85) $$

Then,

$$ \tilde{e}_{\eta} \sim -\frac{\eta}{3\tau^4}, \qquad (86) $$
and

$$ \Sigma = - \frac{\overline{\tau_{\theta} \eta}}{3 \tau^4} = \frac{\omega \overline{\tau_{\theta}^2}}{3 \tau^4} = \frac{\omega}{3} \overline{h_{\theta}^2}. \qquad (87) $$

Setting $\zeta = q - \tilde{c}t$ to be the traveling wave coordinate, then

$$ \frac{d}{d\zeta} = k \frac{d}{d\theta}. \qquad (88) $$

Using this, we reduce the problem to the following one. Find the limits

$$ \lim_{k \to 0} \tilde{c} \left( \frac{\left( \frac{dh}{d\zeta} \right)^2}{3k} + \frac{\delta}{k} \right) \quad \text{and} \quad \lim_{k \to 0} \tilde{c}^2 \frac{\delta}{k}. $$

Since the limit of $\tilde{c}$ is $\tilde{c}^2 = h_1 h_2 h_3 \to h_1^2 h_3$ as $k \to 0$, we need only to find

$$ \lim_{k \to 0} \left( \frac{\left( \frac{dh}{d\zeta} \right)^2}{3k} + \frac{\delta}{k} \right) \quad \text{and} \quad \lim_{k \to 0} \frac{\delta}{k}. $$

## A.2 Computation of average values in mass Lagrangian coordinates

The wavelength in mass Lagrangian coordinates is

$$ L = \frac{2|\tilde{c}|}{\sqrt{3}} \int_{h_2}^{h_3} \frac{h dh}{\sqrt{P(h)}}. \qquad (89) $$

Also,

$$ \overline{\left(\frac{dh}{d\zeta}\right)^2} = \frac{2}{L} \sqrt{\frac{3}{\tilde{c}^2}} \int_{h_2}^{h_3} \frac{\sqrt{Pdh}}{h}. \qquad (90) $$

Then, with $k = 2\pi/L$, the limit as $m \to 1$ of the ratio

$$ \begin{aligned} \frac{\left(\frac{dh}{d\zeta}\right)^2}{3k} &= \frac{1}{2\pi} \frac{2}{3} \sqrt{\frac{3}{\tilde{c}^2}} \int_{h_2}^{h_3} \frac{\sqrt{Pdh}}{h} \to \frac{1}{2\pi} \frac{4}{\sqrt{3}} \left( \frac{(3-2n)\sqrt{n}}{3(1-n)} - \ln\left(\frac{1+\sqrt{n}}{\sqrt{1-n}}\right) \right) \\ &= \frac{1}{2\pi} \frac{4}{\sqrt{3}} \left( \frac{(3-2n)\sqrt{n}}{3(1-n)} - \frac{1}{2}\ln\left(\frac{1+\sqrt{n}}{1-\sqrt{n}}\right) \right) \\ &= \frac{1}{2\pi} \frac{4}{\sqrt{3}} \left( \frac{(3-2n)\sqrt{n}}{3(1-n)} + \frac{1}{2}\ln\left(\frac{1-\sqrt{n}}{1+\sqrt{n}}\right) \right). \end{aligned} $$

Recall,

$$ m = \frac{h_3 - h_2}{h_3 - h_1}, \quad n = \frac{h_3 - h_2}{h_3}, \quad 0 < n < m < 1. \qquad (91) $$
Indeed, making the change of variables $h = h_2 + t(h_3 - h_2)$, $t \in [0, 1]$ and then $t = \cos^2 \theta$, $\theta \in [0, \pi/2]$, one obtains

$$ 
\begin{aligned}
\int_{h_2}^{h_3} \frac{\sqrt{P} dh}{h} &= \frac{(h_3 - h_2)^2 \sqrt{h_3 - h_1}}{h_3} \int_0^1 \frac{\sqrt{t(1-t)(1-m(1-t))}}{1-n(1-t)} dt \\
&= \frac{2(h_3 - h_2)^2 \sqrt{h_3 - h_1}}{h_3} \int_0^{\pi/2} \frac{\sin^2 \theta \cos^2 \theta \sqrt{1 - m \sin^2 \theta}}{1 - n \sin^2 \theta} d\theta \\
&\rightarrow \frac{2(h_3 - h_1)^{5/2}}{h_3} \int_0^{\pi/2} \frac{\sin^2 \theta \cos^2 \theta \sqrt{1 - \sin^2 \theta}}{1 - n \sin^2 \theta} d\theta, \quad \text{as } m \rightarrow 1 \\
&= \frac{2(h_3 - h_1)^{5/2}}{h_3} \left( \frac{3 - 2n}{3n^2} - \frac{1 - n}{n^{5/2}} \ln \left( \frac{1 + \sqrt{n}}{\sqrt{1 - n}} \right) \right).
\end{aligned}
$$

The next step is the expression for $\delta/k$. One has

$$ \frac{\delta}{k} = \frac{\overline{\tau^2} - (\overline{\tau})^2}{k} = \frac{1}{2\pi} \frac{2|\tilde{c}|}{\sqrt{3}} \int_{h_2}^{h_3} \frac{hdh}{\sqrt{P}(h)} \left( \frac{\int_{h_2}^{h_3} \frac{dh}{h\sqrt{P}}}{\int_{h_2}^{h_3} \frac{hdh}{h\sqrt{P}}} - \left( \frac{\int_{h_2}^{h_3} \frac{dh}{\sqrt{P}}}{\int_{h_2}^{h_3} \frac{hdh}{\sqrt{P}}} \right)^2 \right). \quad (92) $$

The expressions of the integrals

$$ \int_{h_2}^{h_3} \frac{dh}{\sqrt{P}}, \quad \int_{h_2}^{h_3} \frac{dh}{h\sqrt{P}}, \quad \int_{h_2}^{h_3} \frac{hdh}{\sqrt{P}} \qquad (93) $$

are given in [54]. They give us

$$ \frac{\delta}{k} = \frac{1}{2\pi} \frac{4|\tilde{c}|}{\sqrt{3} h_3 \sqrt{h_3 - h_1}} \left( \Pi(n, m) - \frac{K^2(m)}{(1 - n/m)K(m) + (n/m)E(m)} \right). \quad (94) $$

As usual, $K(m)$ is the complete elliptic integral of the first kind

$$ K(m) = \int_{0}^{\pi/2} \frac{d\theta}{\sqrt{1 - m \sin^2 \theta}}, \quad (95) $$

$E(m)$ is the complete elliptic integral of the second kind

$$ E(m) = \int_{0}^{\pi/2} \sqrt{1 - m \sin^2 \theta} \, d\theta, \quad (96) $$

and $\Pi(n, m)$ is the complete elliptic integral of the third kind

$$ \Pi(n, m) = \int_{0}^{\pi/2} \frac{d\theta}{(1 - n \sin^2 \theta) \sqrt{1 - m \sin^2 \theta}}. \quad (97) $$

The limit

$$ \lim_{m \to 1} \left( \Pi(n, m) - \frac{K^2(m)}{(1 - n/m)K(m) + (n/m)E(m)} \right) \quad (98) $$
is singular. It can be obtained in the form

$$ \lim_{m \to 1} \left( \Pi(n, m) - \frac{K^2(m)}{(1-n/m)K(m) + (n/m)E(m)} \right) = \frac{n}{(1-n)^2} + \frac{\sqrt{n}}{2(1-n)} \ln \left( \frac{1-\sqrt{n}}{1+\sqrt{n}} \right). \quad (99) $$

Hence,

$$ \lim_{m \to 1} \frac{\delta}{k} = \lim_{m \to 1} \frac{\overline{\tau^2} - (\overline{\tau})^2}{k} = \frac{1}{2\pi} \frac{4}{\sqrt{3}} \left( \frac{\sqrt{n}}{1-n} + \frac{1}{2} \ln \left( \frac{1-\sqrt{n}}{1+\sqrt{n}} \right) \right). \quad (100) $$

In particular, this implies

$$ \lim_{k \to 0} \left( \frac{\left( \frac{dh}{d\zeta} \right)^2}{3k} + \frac{\delta}{k} \right) = \frac{4\sqrt{3}}{2\pi} \left( \frac{(6-2n)\sqrt{n}}{3(1-n)} + \ln \left( \frac{1-\sqrt{n}}{1+\sqrt{n}} \right) \right). \quad (101) $$

Finally, the limiting SGN-Whitham modulation equation for the solitary wave field in mass Lagrangian coordinates is eq. (31) with (32) and (33). Equation (31) is the solitary wave limit of eq. (83) and expresses the conservation of wave action. In the limit $m \to 1$, the quantity $n$ becomes eq. (34).

## B Solitary wave limit of the modulation equations in Eulerian coordinates

The solitary wave limit (31) can also be obtained from the modulation equations written in Eulerian coordinates (25), obtained previously in [12, 54]. Here, we use the formulation of [54] with the notation $c$ for the phase velocity and $g = 1$. Note that in [54], the phase velocity was denoted by $D$ and the roots of the polynomial were $h_0 < h_1 < h_2$ instead of $h_1 < h_2 < h_3$ used in this work. The SGN-Whitham modulation equations in Eulerian coordinates are

$$ \bar{h}_t + (p + \bar{h}c)_x = 0 \quad (102) $$

$$ (p + \bar{h}c)_t + \left( \bar{h}c^2 + \frac{1}{2}I_2 + 2pc \right)_x = 0 \quad (103) $$

$$ \left( \frac{1}{2} \bar{h} c^2 + \frac{1}{2} I_1 \bar{h} - \frac{1}{2} I_2 + I_3 \bar{h}^{-1} + pc \right)_t \quad (104) $$

$$ + \left( \frac{1}{2} \bar{h} c^3 + \frac{1}{2} I_1 \bar{h} c + I_3 \bar{h}^{-1} c + \frac{3}{2} pc^2 + \frac{1}{2} pI_1 \right)_x = 0, \quad (105) $$

$$ k_t + (kc)_x = 0, \quad (106) $$

with

$$ I_1 = h_1 + h_2 + h_3, \quad I_2 = h_1 h_2 + h_2 h_3 + h_1 h_3, \quad I_3 = h_1 h_2 h_3, \quad p = -\sigma \sqrt{I_3}, \quad (107) $$

$$ \bar{h} = h_1 + (h_3 - h_1) \frac{E(m)}{K(m)}, \quad \bar{h}^{-1} = \frac{\Pi(n, m)}{h_3 K(m)}, \quad \frac{2\pi}{k} = L = 4 \sqrt{\frac{h_1 h_2 h_3}{3} \frac{K(m)}{\sqrt{h_3 - h_1}}}, \quad n = 1 - \frac{h_2}{h_3}. \quad (108) $$
We first rewrite the quasi-linear system above in the new variables

$$ (h_1, h_2, h_3, c) \rightarrow \mathbf{h} = (h_1, h_2, a = h_3 - h_2, u_1 = c - \sigma\sqrt{h_1 + a}), \quad (109) $$

using a symbolic computation package, e.g., Mathematica. Note that $h_1, u_1$ are the spatial averages $\bar{h}, \bar{u}$ in the solitary wave limit $h_2 \to h_1$. The expressions are particularly long and not written here for brevity. We denote the obtained system in the abstract form

$$ \mathbf{h}_t + M(\mathbf{h})\mathbf{h}_x = 0, \quad (110) $$

where $M(\mathbf{h})$ is a $4 \times 4$ matrix.

The variables $h_1, h_2, a, u_1$ remain finite in the solitary wave limit $h_2 \to h_1$ (or $m \to 1$), and the limit of the system (110) thus yields non-trivial modulation equations for $h_1, u_1$ and $a$. The coefficients of $M(\mathbf{h})$ are rational functions of $h_1, h_2, a, u_1, E(m), K(m)$ and $\Pi(n, m)$. Using the series expansions

$$ E(m) = 1 + \mathcal{O}(1-m), \quad K(m) = 2\ln(2) - \frac{\ln(1-m)}{2} + \mathcal{O}(1-m), \quad (111) $$

$$ \Pi(n, m) = \frac{\sqrt{n} \tanh^{-1}(\sqrt{n}) - 2\ln(2) + \ln(1-m)/2}{n-1} + \mathcal{O}(1-m), \quad (112) $$

where

$$ \tanh^{-1}(\sqrt{n}) = \frac{1}{2} \ln \left( \frac{1 + \sqrt{n}}{1 - \sqrt{n}} \right) = \sinh^{-1}(z), \quad n = \frac{a}{h_1 + a}, \quad z = \sqrt{\frac{a}{h_1}}, \quad (113) $$

we obtain the solitary wave limit in Eulerian coordinates

$$ h_{1,t} + (u_1 h_1)_x = 0, \quad (114) $$

$$ u_{1,t} + u_1 u_{1,x} + h_{1,x} = 0, \quad (115) $$

$$ a_t + c a_x - \sigma \frac{(2n-3)\sqrt{n} + (3-n)(1-n)\tanh^{-1}(\sqrt{n})}{(1-n)^{3/2}(2\sqrt{n} - (1-n)\tanh^{-1}(\sqrt{n}))} \sqrt{h_1} h_{1,x} \quad (116) $$

$$ -\frac{3\sqrt{n} - (3-n)\tanh^{-1}(\sqrt{n})}{2\sqrt{n} - (1-n)\tanh^{-1}(\sqrt{n})} h_1 u_{1,x} = 0, $$

$$ c = u_1 + \sigma \sqrt{h_1 + a}. \quad (117) $$

Equation (116) for the wave amplitude $a$ is equivalent to the equation for the wave action (31) written in Eulerian coordinates as

$$ (h_1 F(n, h_1))_t + (h_1 u_1 F(n, h_1) + G(n, h_1))_x = 0, \quad (118) $$

with $F$ and $G$ given by (32) and (33), respectively. Using the notation $h_1 = \bar{h}, u_1 = \bar{u}$, and replacing the amplitude $a$ by the dimensionless amplitude $z = \sqrt{\frac{a}{\bar{h}}}$, we recover from (114)-(116) equations (36)-(38).
## C Compatibility condition

**Theorem 1**

Consider an overdetermined system of equations for the unknown $w(x_1, x_2, x_3)$

$$w_{x_1} = \alpha(x_1, x_2, x_3)w_{x_2}, \quad w_{x_1} = \beta(x_1, x_2, x_3)w_{x_3}. \qquad (119)$$

The system is compatible if and only if

$$\alpha^2\beta_{x_2} - \beta^2\alpha_{x_3} + \beta\alpha_{x_1} - \alpha\beta_{x_1} = 0. \qquad (120)$$

**Proof** The compatibility condition is

$$(\partial_{x_1} - \alpha\partial_{x_2})(\partial_{x_1} - \beta\partial_{x_3})w - (\partial_{x_1} - \beta\partial_{x_3})(\partial_{x_1} - \alpha\partial_{x_2})w = 0. \qquad (121)$$

It is equivalent to:

$$\left( \frac{\alpha\beta_{x_2} - \beta_{x_1}}{\beta} - \frac{\beta\alpha_{x_3} - \alpha_{x_1}}{\alpha} \right) u_{x_1} = 0. \qquad (122)$$

The last expression is equivalent to (120).

In our case the governing equations in the solitary wave limit are

$$\tau_t - u_q = 0, \quad u_t - g/\tau^3\tau_q = 0, \quad F(n, \tau)_t + G(n, \tau)_q = 0. \qquad (123)$$

Here, we replaced $h_1^{-1}$ by $\tau$ in the functions $F(n, h_1)$ and $G(n, h_1)$ and used the same letters $F$ and $G$ for the functions of new arguments. Let us suppose that there exists the Riemann invariant $R(n, \tau, u)$ such that

$$R_t + \lambda R_q = 0, \quad \lambda = \frac{G_n}{F_n}. \qquad (124)$$

Expanding the equation for $R$, one obtains

$$R_n(n_t + \lambda n_q) + R_\tau(\tau_t + \lambda \tau_q) + R_u(u_t + \lambda u_q) = 0. \qquad (125)$$

Also, replacing the time derivatives $n_t$, $u_t$ and $\tau_t$ from (123), we obtain two equations corresponding to the vanishing coefficients in front of the derivatives $u_q$ and $\tau_q$. They are of the form

$$R_n = \alpha(n, \tau)R_\tau, \quad R_n = \beta(n, \tau)R_u, \qquad (126)$$

with coefficients $\alpha$ and $\beta$ that are independent of $u$. We take thus $x_1 = n$, $x_2 = \tau$ and $x_3 = u$. According to Theorem 1, the compatibility condition is:

$$\beta_{\tau} = \left(\frac{\beta}{\alpha}\right)_n. \qquad (127)$$

One can check with Mathematica software that this condition is not satisfied. Thus, the Riemann invariant corresponding to the eigenvalue $\lambda$ doesn't exist.
# D Numerical method

Let $p$ be the integrated fluid pressure divided by the constant density $\rho$ and defined by

$$p = \frac{1}{2}h^2 + \frac{h^2}{3}\frac{d^2h}{dt^2}.$$

To solve the one-dimensional, homogeneous SGN equations (1), (2) numerically, as in [23] for the multidimensional SGN equations over topography, we use the $\varpi$ formulation of the equations:

$$h_t + (hu)_x = 0, \qquad (128a)$$

$$(hu)_t + \left(hu^2 + \frac{1}{2}h^2\right)_x = -\varpi_x. \qquad (128b)$$

Here $\varpi = p - \frac{1}{2}h^2$ denotes the averaged non-hydrostatic part of the pressure that is obtained by solving the linear elliptic problem

$$-\frac{h^3}{3}\left(\frac{\varpi_x}{h}\right)_x + \varpi = \frac{2}{3}h^3 u_x^2 + \frac{h^3}{3}h_{xx}. \qquad (128c)$$

In the algorithm, we employ the hyperbolic-elliptic splitting approach developed previously in [38, 22, 23]. This algorithm consists of two steps. In the first, hyperbolic step, we employ a state-of-the-art method for the numerical solution of the hyperbolic systems of equations (128a), (128b) over the time step $\Delta t$. In the second, elliptic step, using the approximate solution $h$ and $u$ computed during the hyperbolic step, we numerically invert the elliptic operator (128c) for $\varpi$ with prescribed boundary conditions.

Note that in the hyperbolic step, rather than writing the equations in the conservation form $q_t + f(q)_x = 0$ with $q = (h, hu)^T$ and $f = (hu, hu^2 + \frac{1}{2}h^2 + \varpi)^T$, which is essential in the conservative first-order setting [39] but is difficult to achieve more than first order accuracy, we write it in the form of a balance law, see [22, 23] for the details of the numerical implementation.

# References

[1] NIST Digital Library of Mathematical Functions, 2024. pages 6

[2] M. J. Ablowitz, X.-D. Luo, and J. T. Cole. Solitons, the Korteweg-de Vries equation with step boundary values, and pseudo-embedded eigenvalues. J. Math. Phys., 59(9):091406, 2018. pages 3

[3] Mark J. Ablowitz, Justin T. Cole, Gennady A. El, Mark A. Hoefer, and Xu-Dan Luo. Soliton–mean field interaction in Korteweg–de Vries dispersive hydrodynamics. Stud. Appl. Math., 151(3):795–856, 2023. pages 2, 3, 15, 17

[4] T. B. Benjamin, J. L. Bona, and J. J. Mahony. Model Equations for Long Waves in Nonlinear Dispersive Systems. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences, 272(1220):47–78, 1972. pages 3

[5] Sylvie Benzoni-Gavage, Colin Mietka, and L. Miguel Rodrigues. Modulated equations of Hamiltonian PDEs and dispersive shocks. Nonlinearity, 34(1):578, 2021. pages 2, 3, 9
[6] Gino Biondini, Sitai Li, and Dionyssios Mantzavinos. Soliton trapping, transmission, and wake in modulationally unstable media. *Phys. Rev. E*, 98(4), 2018. pages 3

[7] Oliver Bühler. *Waves and Mean Flows*. Cambridge Monographs on Mechanics. Cambridge University Press, New York, NY, 2nd edition, 2014. pages 1, 2

[8] Weng Cho Chew. *Waves and Fields in Inhomogenous Media*. John Wiley & Sons, 1999. pages 1

[9] T. Congy, G. A. El, M. A. Hoefer, and M. Shearer. Dispersive Riemann problems for the Benjamin–Bona–Mahony equation. *Stud. Appl. Math.*, 147(3):1089–1145, 2021. pages 12

[10] G. A. El. Resolution of a shock in hyperbolic systems modified by weak dispersion. *Chaos*, 15:037103, 2005. pages 3, 9, 12, 13

[11] G. A. El, A. Gammal, E. G. Khamis, R. A. Kraenkel, and A. M. Kamchatnov. Theory of optical dispersive shock waves in photorefractive media. *Phys. Rev. A*, 76(5):053813, 2007. pages 12

[12] G. A. El, R. H. J. Grimshaw, and N. F. Smyth. Unsteady undular bores in fully nonlinear shallow-water theory. *Phys. Fluids*, 18(2):027104–17, 2006. pages 3, 8, 12, 13, 14, 15, 18, 27

[13] G. A. El, R. H. J. Grimshaw, and N.F. Smyth. Asymptotic description of solitary wave trains in fully nonlinear shallow-water theory. *Phys. D*, 237(19):2423–2435, 2008. pages 13

[14] G. A. El and M. A. Hoefer. Dispersive shock waves and modulation theory. *Physica D: Nonlinear Phenomena*, 333:11–65, 2016. pages 3, 9, 12, 15

[15] G. A. El, V. V. Khodorovskii, and A. V. Tyurina. Undular bore transition in bi-directional conservative wave dynamics. *Phys. D*, 206(3-4):232–251, 2005. pages 12

[16] J. G. Esler and J. D. Pearce. Dispersive dam-break and lock-exchange flows in a two-layer fluid. *J. Fluid Mech.*, 667:555–585, 2011. pages 12

[17] N. Favrie and S. Gavrilyuk. A rapid numerical method for solving Serre–Green–Naghdi equations describing long free surface gravity waves. *Nonlinearity*, 30(7):2718, 2017. pages 5

[18] H. Flaschka, M. G. Forest, and D. W. McLaughlin. Multiphase averaging and the inverse spectral solution of the Korteweg-de Vries equation. *Comm. Pure Appl. Math.*, 33:739–784, 1980. pages 2

[19] S. L. Gavrilyuk. Large amplitude oscillations and their "thermodynamics" for continua with "memory". *Eur J Mech B:Fluids*, 13(6):753–764, 1994. pages 5, 9, 23

[20] Sergey Gavrilyuk, Henrik Kalisch, and Zahra Khorsand. A kinematic conservation law in free surface flow. *Nonlinearity*, 28(6):1805–1821, 2015. pages 5

[21] Sergey Gavrilyuk and Christian Klein. Numerical study of the Serre–Green–Naghdi equations in 2D. *Nonlinearity*, 37(4):045014. pages 5

[22] Sergey Gavrilyuk, Boniface Nkonga, Keh-Ming Shyue, and Lev Truskinovsky. Stationary shock-like transition fronts in dispersive systems. *Nonlinearity*, 33(10):5477–5509, 2020. pages 30
[23] Sergey Gavrilyuk and Keh-Ming Shyue. 2D Serre-Green-Naghdi Equations over Topography: Elliptic Operator Inversion Method. J. Hydraul. Eng., 150(1):04023054. pages 5, 30

[24] Sergey Gavrilyuk and Keh-Ming Shyue. Singular solutions of the BBM equation: Analytical and numerical study. Nonlinearity, 35(1):388–410, 2021. pages 3, 9

[25] Sergey L. Gavrilyuk and Denis Serre. A Model of a Plug-Chain System Near the Thermodynamic Critical Point: Connection with the Korteweg Theory of Capillarity and Modulation Equations. In Shigeki Morioka and Leen Van Wijngaarden, editors, IUTAM Symp. Waves Liq. Liq. Two-Phase Syst., pages 419–428. Springer Netherlands. pages 23

[26] S.L. Gavrilyuk and V.M. Teshukov. Generalized vorticity for bubbly liquid and dispersive shallow water equations. Contin. Mech. Thermodyn., 13(6):365–382, 2001. pages 3, 5

[27] A.E. Green and P.M. Naghdi. A derivation of equations for wave propagation in water of variable depth. J. Fluid Mech., 78:237–246, 1976. pages 3, 4

[28] Albert Edward Green, N. Laws, and P. M. Naghdi. On the theory of water waves. Proc. R. Soc. Lond. Math. Phys. Sci., 338(1612):43–55. pages 3, 4

[29] A. V. Gurevich, A. L. Krylov, and G. A. El. Nonlinear modulated waves in dispersive hydrodynamics. Sov. Phys. JETP, 71(5):899–910, 1990. pages 9

[30] A. V. Gurevich and L. P. Pitaevskii. Nonstationary structure of a collisionless shock wave. Sov. Phys. JETP, 38(2):291–297, 1974. pages 19

[31] M. A. Hoefer. Shock Waves in Dispersive Eulerian Fluids. J. Nonlinear Sci., 24(3):525–577, 2014. pages 9, 12, 18

[32] S. K. Ivanov and A. M. Kamchatnov. Motion of dark solitons in a non-uniform flow of Bose-Einstein condensate. Chaos, 32(11):113142, 2022. pages 3

[33] S. Jamshidi and E. R. Johnson. The long-wave potential-vorticity dynamics of coastal fronts. J. Fluid Mech., 888:A19, 2020. pages 12

[34] A. M. Kamchatnov. Nonlinear Periodic Waves and Their Modulations. World Scientific Publishing Company, 1st edition edition, 2000. pages 2

[35] A. M. Kamchatnov. Hamiltonian theory of motion of dark solitons in the theory of nonlinear Schrödinger equation. Theor Math Phys, 219(1):567–575, 2024. pages 3, 13

[36] A. M. Kamchatnov and D. V. Shaykin. Propagation of generalized Korteweg–de Vries solitons along large-scale waves. Phys. Rev. E, 108(5):054205, 2023. pages 3, 13

[37] David Lannes. The Water Waves Problem, volume no. 188. of Mathematical Surveys and Monographs. American Mathematical Society, Providence, RI, 2013. pages 4, 5

[38] O. Le Métayer, S. Gavrilyuk, and S. Hank. A numerical scheme for the Green–Naghdi model. J. Comput. Phys., 229(6):2034–2045, 2010. pages 5, 30

[39] R. J. LeVeque. Finite Volume Methods for Hyperbolic Problems. Cambridge University Press, 2002. pages 30
[40] Maojun Li, Philippe Guyenne, Fengyan Li, and Liwei Xu. High order well-balanced CDG-FE methods for shallow water waves by a Green-Naghdi model. Journal of Computational Physics, 257:169-192. pages 5

[41] Yi A Li. Hamiltonian Structure and Linear Stability of Solitary Waves of the Green-Naghdi Equations. J. Nonlinear Math. Phys., 9:99-105. pages 3

[42] Nicholas K. Lowman and M. A. Hoefer. Dispersive shock waves in viscously deformable media. J. Fluid Mech., 718:524-557, 2013. pages 12

[43] M. D. Maiden, N. A. Franco, E. G. Webb, G. A. El, and M. A. Hoefer. Solitary wave fission of a large disturbance in a viscous fluid conduit. J. Fluid Mech., 883:A10, 2020. pages 13

[44] Michelle D. Maiden, Dalton V. Anderson, Nevil A. Franco, Gennady A. El, and Mark A. Hoefer. Solitonic Dispersive Hydrodynamics: Theory and Observation. Phys. Rev. Lett., 120(14):144101, 2018. pages 2, 3, 13, 15, 16, 17

[45] N. Makarenko. A second long-wave approximation in the Cauchy-Poisson problem. Dyn. Cont. Media, 77:56-72. pages 5

[46] Yifeng Mao, Sathyanarayanan Chandramouli, Wenqian Xu, and Mark A. Hoefer. Observation of Traveling Breathers and Their Scattering in a Two-Fluid System. Phys. Rev. Lett., 131(14):147201, 2023. pages 17

[47] Ana Mucalica and Dmitry E. Pelinovsky. Solitons on the rarefaction wave background via the Darboux transformation. Proc. R. Soc. A., 478(2267):20220474, 2022. pages 3

[48] Peter J. Olver. On the Hamiltonian structure of evolution equations. Math. Proc. Camb. Phil. Soc., 88(1):71-88, 1980. pages 3

[49] S. Ryskamp, M. A. Hoefer, and G. Biondini. Oblique interactions between solitons and mean flows in the Kadomtsev-Petviashvili equation. Nonlinearity, 34(6):3583-3617, 2021. pages 3, 15

[50] F. Serre. Contribution à l'étude des écoulements permanents et variables dans les canaux. Houille Blanche, (3, 6):374-388, 830-872, June-July, December, 1953. pages 3, 4

[51] P. Sprenger, M. A. Hoefer, and G. A. El. Hydrodynamic optical soliton tunneling. Phys. Rev. E, 97(3):032218, 2018. pages 3, 11, 15

[52] C.H. Su and C. S. Gardner. Korteweg-de Vries equation and generalizations. III. Derivation of the Korteweg-de Vries and Burgers equation. J. Math. Phys., 10:536-539, 1969. pages 3, 4

[53] Pierre Suret, Stephane Randoux, Andrey Gelash, Dmitry Agafontsev, Benjamin Doyon, and Gennady El. Soliton gas: Theory, numerics, and experiments. Phys. Rev. E, 109(6):061001, 2024. pages 21

[54] Sergey Tkachenko, Sergey Gavrilyuk, and Keh-Ming Shyue. Hyperbolicity of the Modulation Equations for the Serre-Green-Naghdi Model. Water Waves, 2(2):299-326, 2020. pages 3, 8, 26, 27
[55] Kiera van der Sande, Gennady A. El, and Mark A. Hoefer. Dynamic soliton–mean flow inter-action with non-convex flux. J. Fluid Mech., 928:A21, 2021. pages 3, 15, 17

[56] G. B. Whitham. Non-linear dispersive waves. Proc. Roy. Soc. Ser. A, 283:238-261, 1965. pages 2, 8

[57] G. B. Whitham. Linear and Nonlinear Waves. Wiley, New York, 1999. pages 2, 22, 23, 24