1 Introduction

The Korteweg-de Vries (KdV) equation is a model equation describing the evolution of long waves at the surface of a body of fluid. The KdV equation was derived in 1895 by Korteweg and de Vries [18], but was already featured in earlier work by Boussinesq [7]. The main assumptions on the waves to be represented by solutions of the KdV equation are that they be of small amplitude and long wavelength when compared to the undisturbed depth of the fluid, that the wave motion be predominantly one-directional, and that transverse effects be weak. In dimensional variables, the KdV equation is given by

$$\begin{aligned} \eta_t+ c_0 \eta_x + \frac{3}{2} \frac{c_0}{h_0} \eta\eta_x + \frac{c_0 h_0^2}{6} \eta_{xxx} = 0, \end{aligned}$$
(1.1)

where η(x,t) represents the excursion of the free surface, h 0 is the undisturbed water depth, g denotes the gravitational acceleration, and \(c_{0}=\sqrt{gh_{0}}\) is the limiting long-wave speed.

The equation arises in the so-called Boussinesq scaling regime where wavelength and wave amplitude are balanced in such a way as to allow the formation of traveling-wave solutions. Denoting by ℓ a typical wavelength and by a a typical amplitude of the wavefield to be described, the number α=a/h 0 represents the relative amplitude, and \(\beta= h_{0}^{2} / \ell^{2}\) measures the relative wavenumber. The waves fall into the Boussinesq regime if both α and β are small, and of similar size. In this case, the KdV-equation arises as a simplified asymptotic model describing the wavemotion. In other words, solutions of the full water-wave problem may be approximated on certain time scales by solutions of the KdV equation. This latter fact can be understood in the sense of asymptotics, but has also been established by mathematical proof by Craig [9] and Schneider and Wayne [24]. Incorporated in the arguments of these works are existence results for the water-wave problem in the context of the full Euler equations in the appropriate scaling. Such results are now also available independently (see Lannes [19] and Wu [28, 29]). In particular, Alvarez-Samaniego and Lannes have obtained long-time existence of solutions of the water-wave problem [3] which can be applied to a number of different scaling regimes. A further significant improvement of the results of [9, 24] was achieved by Bona et al. [6] who proved more refined convergence estimates of solutions of the water-wave problem to a family of long-wave systems as well as to one-directional models such as the KdV equations. Further extensions of this method and applications to other systems can be found for instance in the work of Lannes and Bonneton [21].

One of the early drivers of research relating to the KdV equation was the discovery of elastic overtaking collisions of solitons which in some sense seemed to resemble the dynamics of a linear differential equation. The discovery of this elastic solitary-wave interaction subsequently led to the discovery of an infinite number of time-invariant integrals (Miura [23]), and the development of the inverse-scattering method which can be used to provide exact closed form solutions for a broad class of initial data (Ablowitz and Segur [1], Green et al. [13]).

Apart from being a paradigm for the use of the inverse-scattering method, the KdV equation has been used in a large number of studies in the context of wave problems in fluid dynamics. Various dynamical quantities connected with the KdV equation have appeared in the literature. However, it is difficult to find definitive expressions for and derivations of the most important quantities, such as the energy flux. It is our purpose in the present work to present a framework in which mass, momentum and energy fluxes and densities can be expressed in terms of the principal unknown η of equation (1.1)

If the equation is given in the dimensional form (1.1), then the first three conserved integrals are

$$\begin{aligned} \int_{- \infty}^{\infty} \eta\, dx, \qquad \int_{- \infty}^{\infty} \eta^2 \, dx, \quad\mbox{ and } \quad\int_{- \infty}^{\infty} \biggl( \eta^3 - \frac{h_0^3}{3} \eta_x^2 \biggr) \, dx. \end{aligned}$$
(1.2)

The first integral is found to be invariant with respect to time t as soon as it is recognized that the KdV equation can be written in the form

$$\begin{aligned} \frac{\partial}{\partial t}(h_0 + \eta) + \frac {\partial}{\partial x}\biggl( c_0 \eta+ \frac{3}{4}\frac{c_0}{h_0} \eta^2 + \frac{c_0h_0^2}{6} \eta_{xx} \biggr) = 0, \end{aligned}$$
(1.3)

where the quantity appearing under the time derivative is interpreted as a mass density, and the term appearing under the spatial derivative is the mass flux through a cross section of unit width due to the passage of a surface wave. Invariance of the second and third integrals is obtained from similar identities, but the fluxes appearing in these identities do not represent any concrete physical quantities.

The mass balance (1.3) appears in the literature (see [15]), and one may ask whether it is possible to derive expressions for the momentum and energy densities and fluxes which permit the formulation of balance laws similar to (1.3). This problem has been partially solved in the case of steady solutions of the KdV equation. Indeed, Benjamin and Lighthill [4] used the spatial invariance of the mass flux per unit span Q, the momentum flux per unit span corrected for the pressure force S and the energy per unit mass R in steady Euler flow to develop a method for deriving a time-independent KdV equation which contains the quantities Q, R, and S as parameters.

In the case of the time-dependent problem it seems that most of the work dedicated to questions such as outlined above has focused on the question of conservation of the total energy. In the full water-wave problem, the total energy of the wave system is given by the Hamiltonian functional first recorded by Zakharov [30]. As it represents the total energy of a closed physical system, this Hamiltonian is invariant with respect to time. It was shown by Craig and Groves [10] that if the derivation of simplified evolution equations such as the Boussinesq system and the KdV equation is based on approximating the Hamiltonian function of the water-wave problem, then it is possible to define the total energy of the wave system in the corresponding approximation. This theory is quite satisfactory if the total energy is sought, but it yields no information about other quantities such as energy flux or momentum flux.

The rationale used in the present work is based on requiring mass, momentum and energy conservation to the same order as the evolution equation is valid. As will be reviewed in Sect. 2, the derivation of the KdV equation yields the relation

$$\begin{aligned} \tilde{\eta}_{\tilde{t}} + \tilde{\eta }_{\tilde{x}} + \frac{3}{2} \alpha\tilde{\eta}\tilde{\eta}_{\tilde{x}} + \frac{1}{6}\beta \tilde{\eta}_{\tilde{x}\tilde{x}\tilde{x}}= {\mathcal{O}}\bigl(\alpha^2, \alpha\beta, \beta^2\bigr) \end{aligned}$$
(1.4)

in non-dimensional variables. The KdV equation appears if terms of \({\mathcal{O}}(\alpha^{2},\alpha\beta, \beta^{2})\) are disregarded. As the development of the mass balance law associated to the KdV equation in Sect. 3 will bring to light, if non-dimensional mass density \(\tilde{M}(\tilde{\eta})\) and mass flux \(\tilde{q}_{M}(\tilde{\eta})\) are defined appropriately in terms of \(\tilde{\eta}\) and its derivatives, then the relation

$$\begin{aligned} \frac{\partial}{\partial\tilde{t}} \tilde {M}(\tilde{\eta}) + \frac{\partial}{\partial\tilde{x}} \tilde{q}_{M}(\tilde{\eta}) = {\mathcal{O}}\bigl(\alpha^2, \alpha \beta, \beta^2\bigr) \end{aligned}$$
(1.5)

will hold. By omitting terms of \({\mathcal{O}}(\alpha^{2}, \alpha\beta, \beta^{2})\), the approximate mass balance equation is found, and as it will turn out, this equation is the same as the KdV equation. Following the same idea, similar balance equations are developed in Sects. 4 and 5 for the approximate conservation of momentum and energy by defining \(\tilde{I}(\tilde{\eta})\), \(\tilde{q}_{I}(\tilde{\eta})\), \(\tilde{E}(\tilde{\eta})\), and \(\tilde{q}_{E}(\tilde{\eta})\) satisfying

$$\begin{aligned} \frac{\partial}{\partial\tilde{t}} \tilde {I}(\tilde{\eta}) + \frac{\partial}{\partial\tilde{x}} \tilde{q}_{I}(\tilde{\eta}) = {\mathcal{O}}\bigl(\alpha^2, \alpha \beta, \beta^2\bigr), \end{aligned}$$
(1.6)
$$\begin{aligned} \frac{\partial}{\partial\tilde{t}} \tilde {E}(\tilde{\eta}) + \frac{\partial}{\partial\tilde{x}} \tilde{q}_{E}(\tilde{\eta}) = {\mathcal{O}}\bigl(\alpha^2, \alpha \beta, \beta^2\bigr). \end{aligned}$$
(1.7)

As a by-product of developing the energy balance law, an expression for the total head in the context of the KdV approximation can also be found. Since the previous work by Benjamin and Lighthill [4] on steady solutions of the KdV equation has been a prime motivation for the present study, a comparison of the relevant quantities is presented in Sect. 6. A natural question arising from the above discussion is whether the conservation integrals (1.2) can be meaningfully interpreted in connection with conservation of total mass, momentum and energy. This question will be addressed in Sect. 7.

The derivations presented in the present work are formal, and there is no rigorous mathematical proof of the convergence of these approximations as the small parameters α and β approach zero. The main advancement of the present work is the identification of the expressions which satisfy the balance laws (1.5), (1.6) and (1.7), and the comparison with previous asymptotic results. While a proof of the validity of (1.5), (1.6) and (1.7) might proceed along the lines of the proofs of the validity of the KdV equation as a water-wave model, as shown in [6, 9, 24], such a study is beyond the scope of the present article.

Before we leave the introduction, let us mention some further related work. In the context of steady solutions one may exploit the conservation of mass, momentum and energy in the water-wave problem. Examples are the spatial Hamiltonian approach advocated by Bridges [8], and the work of Longuet-Higgins and Fenton [22] on the solitary wave. A study which is closer to the spirit of the present article is provided by Dutykh and Dias [12] who supplemented a Boussinesq system with an energy equation which yields information about an energy density similar to the quantity appearing in Sect. 5. The present work is related to a recent study of mechanical balance laws in a family of Boussinesq systems by Ali and Kalisch [2], and a much earlier attempt to record similar quantities for the single Boussinesq equation by Keulegan an Patterson [17].

We should also mention the concept of wave action conservation which yields an additional conservation law in the water-wave problem which can be used in the case of non-uniform environments, such as background currents and stratifications. This principle which was pioneered by Whitham [26] and Hayes [16] is based on a Lagrangian description of the problem, and can also be applied in the context of model equations (Grimshaw [14], Whitham [27]).

2 Velocity Field and Pressure

The main aim of this section is to establish expressions for the velocity field and pressure in the fluid in terms of the surface excursion η. These expressions are well known by-products of the derivation of the KdV equation. Nevertheless, it will be convenient to give a brief review of this derivation in order to fix ideas regarding the geometric setup and the notation.

Consider a fluid contained in a long channel of unit width. Fix a coordinate system by aligning the x-axis with the undisturbed free surface, and suppose the fluid domain extends along the entire x-axis. It is assumed that the fluid is inviscid, incompressible and of unit density, the bottom of the channel is flat, and that wave motion transverse to the x-axis can be neglected. The geometric setup is illustrated in Fig. 1. The surface water-wave problem is generally described by the Euler equations with slip conditions at the bottom, and kinematic and dynamic boundary conditions at the free surface. The unknowns are the surface excursion η(x,t), the pressure P(x,z,t), and the horizontal and vertical fluid velocities u 1(x,z,t) and u 2(x,z,t), respectively. With the setup described above, the problem may be posed on a domain \(\{(x,z) \in{\mathbb{R}}^{2} | -h_{0} < z < \eta(x,t) \}\) which extends to infinity in the positive and negative x-direction. On this domain, the two-dimensional Euler equations are

$$\begin{aligned} {\mathbf{u}}_t + (\mathbf{u} \cdot\nabla)\mathbf{u} + \nabla P = \ \mathbf{g},\ \end{aligned}$$
(2.1)
$$\begin{aligned} \nabla\cdot{\mathbf{u}} = \ 0,\ \end{aligned}$$
(2.2)

where u=(u 1,u 2) represents the velocity field, and g=(0,−g) is the body forcing. As surface tension effects are neglected, the dynamic free-surface boundary condition calls for the fluid pressure at the surface to be equal to atmospheric pressure. In addition, the kinematic condition requires the normal velocity of the free surface to be equal to the fluid velocity normal to the surface.

Fig. 1
figure 1

The free surface is described by a function η(x,t). The undisturbed water depth is h 0, and the x-axis is aligned with the free surface at rest. The fluid is shown in light gray, and the bed is shown in dark gray. The density of the is fluid is assumed to be unity

Full size image

Assuming irrotational flow and using the incompressibility, the problem can be written in terms of the Laplace equation for a velocity potential ϕ on the unknown time-dependent domain. The surface boundary conditions are then given by

$$\begin{aligned} \eta_t+\phi_x\eta_x- \phi_z&=0, \quad\mbox{on}\ z=\eta(x,t), \end{aligned}$$
(2.3)
$$\begin{aligned} \phi_t + \frac{1}{2}\bigl( \phi^2_x+ \phi^2_z \bigr) + g\eta&=0,\quad \mbox{on}\ z=\eta(x,t). \end{aligned}$$
(2.4)

In order to bring out the different sizes of the variables, the non-dimensionalization

$$ \tilde{x}=\frac{x}{\ell}, \quad\quad\tilde{z}=\frac{z+h_0}{h_0}, \quad\quad \tilde{\eta}= \frac{\eta}{a}, \quad\quad\tilde{t}= \frac{c_0t}{\ell}, \quad\quad\tilde {\phi}=\frac{c_0}{ga\ell} \phi $$

is used. In the following, the standard method of developing the potential ϕ in an asymptotic series is employed. Using the Laplace equation and the boundary condition at the flat bottom shows that the velocity potential takes the form

$$ \begin{aligned} \tilde{\phi}=\tilde{f}- \beta\frac{\tilde{z}^2}{2} \tilde{f}_{\tilde{x}\tilde{x}} + {\mathcal{O}}\bigl(\beta^2\bigr), \end{aligned} $$
(2.5)

where the function \(\tilde{f}(\tilde{x},\tilde{t})\) is chosen in such a way the \(\tilde{f}_{\tilde{x}}\) represents the horizontal velocity at the bottom. Following the method explained in Bona et al. [5] and Whitham [27], \(\tilde{\phi}\) is substituted into the free surface boundary conditions. First note that (2.4) yields the relation

$$ \tilde{\eta}+\tilde{f}_{\tilde{t}}-\frac{\beta }{2} \tilde{f}_{\tilde{x}\tilde{x}\tilde{t}} + \frac{\alpha}{2} {\tilde{f}_{\tilde{x}}}^2={\mathcal{O}}\bigl(\alpha\beta, \beta^2\bigr). $$
(2.6)

Differentiating (2.6) with respect to \(\tilde{x}\), using (2.3) as the first equation, and denoting \(\tilde{w}= \tilde{f}_{\tilde{x}}\) yields the system

$$\begin{aligned} \begin{aligned} \tilde{\eta}_{\tilde{t}} + \tilde{w}_{\tilde{x}} + \alpha( \tilde{\eta}\tilde{w})_{\tilde{x}} + \frac{1}{6} \beta\tilde {w}_{\tilde{x}\tilde{x}\tilde{x}} = {\mathcal{O}} \bigl(\alpha\beta, \beta^2\bigr), \\ \tilde{w}_{\tilde{t}} + \tilde{\eta}_{\tilde{x}} + \alpha\tilde {w}\tilde{w}_{\tilde{x}} - \frac{1}{2} \beta\tilde{w}_{\tilde{x}\tilde{x}\tilde{t}} = {\mathcal{O}}\bigl(\alpha\beta, \beta^2\bigr). \end{aligned} \end{aligned}$$
(2.7)

Now the KdV equation can be derived from (2.7) by assuming a certain relation between the horizontal velocity \(\tilde{w}\) and the surface excursion \(\tilde {\eta}\). As explained in Whitham [27], a solution of (2.7) traveling chiefly to the right (+) or to the left (−) will necessitate the relation

$$\begin{aligned} \tilde{w}= \pm\tilde{\eta}+ \alpha A + \beta B + {\mathcal{O}}\bigl( \alpha^2, \alpha\beta,\beta^2\bigr). \end{aligned}$$
(2.8)

The functions A and B can be found by substituting (2.8) into (2.7). Requiring both equations in (2.7) to yield the same equation for \(\tilde{\eta}\), and using the first-order equivalence

$$\begin{aligned} \partial_{\tilde{t}} F(\tilde{\eta}) = \mp\partial_{\tilde{x}} F( \tilde{\eta}) + {\mathcal{O}}( \alpha,\beta), \end{aligned}$$
(2.9)

where F is a polynomial in \(\tilde{\eta}\) and its derivatives, leads to

$$\begin{aligned} \begin{aligned} A = \mp\frac{1}{4} \tilde{\eta}^2, \quad\mbox{ and } \quad B = \pm\frac{1}{3} \eta_{\tilde{x}\tilde{x}}. \end{aligned} \end{aligned}$$

Thus requiring the equations in (2.7) to be consistent leads to the non-dimensional KdV equation

$$\begin{aligned} \tilde{\eta}_{\tilde{t}} \pm\biggl( \tilde {\eta}_{\tilde{x}} + \frac{3}{2} \alpha\tilde{\eta}\tilde{\eta}_{\tilde{x}} + \frac {1}{6}\beta \tilde{\eta}_{\tilde{x}\tilde{x}\tilde{x}} \biggr) = {\mathcal {O}}\bigl(\alpha^2, \alpha\beta, \beta^2\bigr), \end{aligned}$$
(2.10)

while the velocity \(\tilde{w}\) is given by

$$\begin{aligned} \tilde{w}= \pm\biggl( \tilde{\eta}- \frac{1}{4} \alpha \tilde{\eta}^2 + \frac{1}{3} \beta\tilde{\eta}_{\tilde{x}\tilde {x}} \biggr) + {\mathcal{O}} \bigl(\alpha^2, \alpha\beta, \beta^2\bigr). \end{aligned}$$
(2.11)

From (2.5), it is plain that the velocity field \((\tilde{\phi }_{\tilde{x}},\tilde{\phi}_{\tilde{z}})\) at any non-dimensional height \(\tilde{z}\) in the fluid column is given by

$$\begin{aligned} \tilde{\phi}_{\tilde{x}}(\tilde{x},\tilde {z},\tilde{t}) &= \pm\tilde{\eta}\mp \frac{1}{4} \alpha\tilde{\eta}^2 \pm\beta\biggl( \frac{1}{3}- \frac{\tilde{z}^2}{2} \biggr) \tilde{\eta}_{\tilde {x}\tilde{x}} + {\mathcal{O}}\bigl( \alpha^2, \alpha\beta, \beta^2\bigr), \\ \tilde{\phi}_{\tilde{z}}(\tilde{x},\tilde{z},\tilde{t}) &= \mp \beta\tilde{z}\tilde{\eta}_{\tilde{x}} + {\mathcal{O}} \bigl(\alpha\beta, \beta^2\bigr). \end{aligned}$$
(2.12)

Neglecting terms of second order in α and β, and reverting to dimensional variables, the KdV equation (1.1) appears in the case of waves propagating mainly to the right. A corresponding equation with different signs appears for waves propagating mainly to the left. It will be convenient later for purposes of comparison to have available the above expressions in the case of a moving reference frame. If the problem is put into a reference frame moving at a velocity U, and the non-dimensionalization \(U = c_{0} \tilde{U}\) is chosen, then the KdV equation and the expression for the horizontal velocity appear as

$$\begin{aligned} \tilde{\eta}_{\tilde{t}} - \tilde{U}\tilde {\eta}_{\tilde{x}} \pm\biggl( \tilde{\eta}_{\tilde{x}} + \frac{3}{2} \alpha\tilde{\eta}\tilde {\eta}_{\tilde{x}} + \frac{1}{6}\beta\tilde{\eta}_{\tilde{x}\tilde{x}\tilde{x}} \biggr) = {\mathcal{O}}\bigl( \alpha^2, \alpha\beta, \beta^2\bigr), \end{aligned}$$
(2.13)
$$\begin{aligned} \tilde{\phi}_{\tilde{x}}(\tilde{x},\tilde {z},\tilde{t}) = - \frac{1}{\alpha} \tilde{U}\pm \tilde{\eta}\mp\frac{1}{4} \alpha\tilde{\eta}^2 \pm\beta\biggl( \frac{1}{3}- \frac{\tilde{z}^2}{2} \biggr) \tilde{\eta}_{\tilde {x}\tilde{x}} + {\mathcal{O}}\bigl( \alpha^2, \alpha\beta, \beta^2\bigr). \end{aligned}$$
(2.14)

Note that there is no assumption on the relative magnitude of U, and the scaling of U reflects the most relevant case U=±c 0. However, the limiting long-wave speed c 0 is far greater than the horizontal velocity of any fluid particle if the amplitude of the surface waves is small. As a consequence, the velocity of the moving reference frame enters the expression (2.14) for the horizontal velocity at the order \(\frac{1}{\alpha}\).

Attention will now be turned to the computation of the fluid pressure associated with a given surface wave. In order to find the expression to the correct order, the hydrostatic pressure needs to be excluded from the calculation. Therefore, we work with the dynamic pressure defined by

$$ P' = P - P_{atm}+ g z. $$

Since the atmospheric pressure is of a magnitude much smaller than the significant terms in the equation, it will be assumed to be zero from here on. As shown by Stoker [25], P′ can be written with the help of the Bernoulli equation in the form

$$ P' = - \phi_t-\frac{1}{2}|\nabla \phi|^2. $$

Using the non-dimensionalization \(\tilde{P}' = \frac{1}{a g}P'\), and substituting the expression (2.5) yields

$$ \tilde{P}' = -\tilde{f}_{\tilde{t}} + \beta\frac{\tilde{z}^2}{2} \tilde{f}_{\tilde{x}\tilde{x}\tilde{t}} - \frac{1}{2}\alpha\tilde {f}^2_{\tilde{x}} + {\mathcal{O}}\bigl(\beta^2,\alpha\beta\bigr). $$

As shown in Ali and Kalisch [2], the relation (2.6) may be used to find the dynamic pressure in the form

$$ \tilde{P}'= \tilde{\eta}+ \frac{1}{2}\beta\bigl({\tilde {z}}^2-1\bigr) \tilde{f}_{\tilde{x}\tilde{x}\tilde{t}} + {\mathcal{O}}\bigl(\alpha\beta,\beta^2\bigr). $$

In line with the previous computation, we use the relations (2.8) and (2.9) to find the expression

$$ \tilde{P}'= \tilde{\eta}- \frac{1}{2}\beta \bigl({ \tilde{z}}^2-1\bigr)\tilde{\eta}_{\tilde{x}\tilde{x}} + {\mathcal {O}}\bigl(\alpha\beta, \beta^2\bigr). $$
(2.15)

In the remainder of this article, we use the expression (2.15), and truncate further only when dictated by the particular balance law.

3 Mass Conservation

In this section, mass conservation properties of the KdV equation are explored. Since the surface-wave problem is invariant under a change to a moving frame, the computation is done in the most general form of a reference frame traveling at a velocity U. Using the incompressibility of the fluid, mass conservation is stated in differential form by (2.2). Using this equation and the kinematic boundary condition (2.3), one can immediately derive the relation

$$\begin{aligned} \frac{\partial}{\partial t}\int^{\eta}_{-h_0} \, dz + \frac {\partial}{\partial x}\int ^{\eta}_{-h_0} \phi_{x}(x,z,t)\, dz = 0. \end{aligned}$$

In non-dimensional form this relation becomes

$$\begin{aligned} \frac{\partial}{\partial\tilde{t}}\int^{1 + \alpha\tilde{\eta }}_{0} \, d\tilde{z}+ \alpha\frac{\partial}{\partial\tilde {x}}\int ^{1 +\alpha\tilde{\eta}}_{0} \tilde{\phi}_{\tilde{x}}(\tilde {x},\tilde{z},\tilde{t})\, d\tilde{z}= 0. \end{aligned}$$

Substituting the expression for \(\tilde{\phi}_{\tilde{x}}\) in terms of \(\tilde{\eta}\) given by (2.14), and integrating with respect to \(\tilde{z}\) leads to the approximation

$$ \frac{\partial}{\partial\tilde{t}}(1 + \alpha\tilde{\eta}) + \frac{\partial}{\partial\tilde{x}}\biggl( -\tilde{U}- \tilde {U}\alpha\tilde{\eta}\pm\tilde{\eta}\pm \frac{3}{4} \alpha^2 \tilde{\eta}^2 \pm\frac{1}{6} \alpha\beta\tilde{\eta}_{\tilde{x}\tilde{x}} \biggr) = {\mathcal {O}}\bigl(\alpha^3, \alpha^2 \beta, \alpha\beta^2\bigr). $$

One may divide by α to find

$$ \frac{\partial}{\partial\tilde{t}}\tilde{\eta }+ \frac{\partial}{\partial\tilde{x}}\biggl( - \tilde{U}\tilde {\eta}\pm\tilde{\eta}\pm \frac{3}{4} \alpha\tilde{\eta}^2 \pm\frac{1}{6} \beta \tilde{\eta}_{\tilde{x}\tilde{x}} \biggr) = {\mathcal{O}}\bigl(\alpha^2, \alpha\beta, \beta^2\bigr). $$
(3.1)

Therefore, if we denote the non-dimensional mass density by

$$\tilde{M}= 1 + \alpha\tilde{\eta}, $$

and the non-dimensional mass flux by

$$\tilde{q}_M = -\tilde{U}-\tilde{U}\alpha\tilde{\eta}\pm \biggl( \alpha \tilde{\eta}+ \frac{3}{4} \alpha^2 \tilde{\eta}^2 + \frac{1}{6} \alpha\beta\tilde{\eta}_{\tilde{x}\tilde{x}} \biggr), $$

the non-dimensional mass balance (1.5) is achieved.

Unlike the KdV equation (2.10) or the formula for the horizontal velocity (2.11), the mass flux contains terms of quadratic order in α and β. However, these terms are necessary since the mass balance equation (1.5) holds to the same order as the evolution equation (2.10). Note also that the differential mass balance equation (3.1) is the same as the non-dimensional KdV equation (2.13), and if terms of order α 2, αβ and β 2 are disregarded, the KdV equation is a mass balance equation. In other words, in the approximation which leads to the KdV equation, mass is exactly conserved.

Using the scaling \(M = h_{0} \tilde{M}\) for the mass density and \(q_{M} = h_{0} c_{0} \tilde{q}_{M}\) for the mass flux reveals that the dimensional forms of these quantities are

$$ M = h_0 + \eta $$

and

$$ q_M = - U (h_0 + \eta) \pm c_0 \biggl(\eta+ \frac{3}{4h_0} \eta^2 + \frac{h_0^2}{6} \eta_{xx} \biggr). $$
(3.2)

4 Momentum Balance

This section is devoted to finding an approximate expression for momentum density and flux, which satisfy the relation (1.6). The incompressibility condition (2.2) may be used to rewrite the first component of the vector equation (2.1), in terms of the velocity potential ϕ as

$$ \phi_{xt} + \bigl(\phi_x^2 \bigr)_x + (\phi_x \phi_z)_z + P_x = 0. $$

Integrating over the fluid column and using the kinematic boundary condition (2.3) yields the relation

$$\begin{aligned} \frac{\partial}{\partial t}\int^{\eta}_{-h_0} \phi_{x} \,dz + \frac {\partial}{\partial x} \int_{-h_0}^{\eta} \bigl\{ \phi_{x}^2 + P \bigr\} dz = 0. \end{aligned}$$

Expressing the last relation in non-dimensional variables gives

$$\begin{aligned} \alpha\frac{\partial}{\partial\tilde{t}}\int_{0}^{1 + \alpha \tilde{\eta}} \tilde{\phi}_{\tilde{x}} \,d \tilde{z}+ \frac{\partial}{\partial\tilde{x}}\int_{0}^{1 + \alpha \tilde{\eta}} \bigl\{ \alpha^2 \tilde{\phi}_{\tilde{x}}^2 + \alpha \tilde{P}' - (\tilde{z}- 1) \bigr\} \,d\tilde{z}= 0. \end{aligned}$$

Substituting \(\tilde{\phi}_{\tilde{x}}\) and \(\tilde{P}'\) as found in Sect. 2 yields the balance equation

$$\begin{aligned} &\frac{\partial}{\partial\tilde{t}}\biggl(-\tilde{U}-\tilde {U}\alpha\tilde{\eta}\pm\alpha\tilde{\eta}\pm\frac{3}{4} \alpha^2 \tilde{\eta}^2 \pm\frac{1}{6} \alpha\beta \tilde{\eta}_{\tilde{x}\tilde{x}} \biggr) + \frac{\partial }{\partial\tilde{x}}\biggl( \tilde{U}^2 + \alpha \tilde{U}^2 \tilde{\eta}+ \frac{1}{2} + \alpha\tilde{\eta} \\ &\quad{} + \frac{3}{2} \alpha^2 \tilde{\eta}^2 + \frac{1}{3} \alpha\beta\tilde{\eta}_{\tilde{x}\tilde{x}} \mp2 \alpha\tilde{U}\tilde{\eta}\mp \frac{3}{2} \alpha^2 \tilde{U}\tilde{\eta}^2 \mp \frac{1}{3} \alpha\beta\tilde{U}\tilde{\eta}_{\tilde{x}\tilde {x}} \biggr) = {\mathcal{O}}\bigl( \alpha^3,\alpha^2\beta,\alpha\beta^2\bigr) . \end{aligned}$$

Thus denoting the non-dimensional momentum density by

$$\tilde{I}=- \tilde{U}- \alpha\tilde{U}\tilde{\eta} \pm \biggl( \alpha\tilde{\eta} + \frac{3}{4} \alpha^2 \tilde{\eta}^2 + \frac{1}{6} \alpha\beta\tilde{\eta}_{\tilde{x}\tilde{x}} \biggr), $$

and the non-dimensional momentum flux by

$$\tilde{q}_I = \frac{1}{2} (1 +\alpha\tilde{\eta})^2 + ( \alpha\tilde{\eta}\mp\tilde{U})^2 + \frac{1}{3} \alpha\beta\tilde{\eta}_{\tilde{x}\tilde{x}} + \alpha\tilde{U}^2 \tilde{\eta} \mp\frac{3}{2} \alpha^2 \tilde{U}\tilde{\eta}^2 \mp\frac{1}{3} \alpha\beta\tilde{U}\tilde{\eta}_{\tilde {x}\tilde{x}}, $$

the non-dimensional momentum balance (1.6) is achieved. Using the natural scalings \(I = c_{0}h_{0} \tilde{I}\) and \(q_{I} = h_{0} c_{0}^{2} \tilde{q}_{I}\), the dimensional forms of these quantities are

$$ I = - U(h_0 +\eta) \pm c_0 \biggl( \eta+ \frac{3}{4h_0} \eta^2 + \frac{h_0^2}{6} \eta_{xx} \biggr) $$

and

$$ q_I = c_0^2 \biggl( \frac{(h_0+\eta)^2}{2h_0} + \frac{h_0^2}{3}\eta_{xx} + \frac{(\eta \mp\frac{h_0}{c_0}U)^2}{h_0} + \frac{U^2}{c_0^2}\eta\mp\frac{3}{2} \frac{U}{c_0 h_0} \eta^2 \mp\frac{h_0^2}{3c_0} U \eta_{xx} \biggr). $$
(4.1)

The reader may take note that the momentum density has the same form as the mass flux (3.2), and that the expression for the momentum flux contains the contribution of the pressure force.

5 Energy Balance

Attention is now turned to the energy balance in the fluid. Using (2.1) and (2.2), an energy equation can be written in the form

$$ \frac{\partial}{\partial t}\biggl\{ \frac{1}{2} | \nabla\phi|^2 + g (z+h_0) \biggr\} + \nabla\cdot\biggl\{ \biggl( \frac{1}{2} | \nabla\phi|^2 + g (z+h_0) + P \biggr) \nabla\phi \biggr\} =0. $$

An integration over the depth of the fluid yields

$$ \frac{\partial}{\partial t}\int_{-h_0}^{\eta} \biggl\{ \frac{1}{2} |\nabla \phi|^{2} + g (z+h_0) \biggr\} dz + \frac{\partial}{\partial x}\int _{-h_0}^{\eta} \biggl\{ \frac{1}{2} |\nabla \phi|^{2} + g (z+h_0) + P \biggr\} \phi_x \, dz = 0. $$

Converting to non-dimensional variables transforms the last relation into

$$\begin{aligned} &\frac{\partial}{\partial\tilde{t}}\int_{0}^{1 + \alpha\tilde {\eta}} \biggl\{ \frac{\alpha^2}{2} \bigl( \tilde{\phi}_{\tilde {x}}^2 + \frac{1}{\beta} \tilde{\phi}_{\tilde{z}}^2 \bigr) + \tilde{z}\biggr\} d\tilde{z} \\ &\quad{} + \alpha\frac{\partial}{\partial\tilde{x}}\int_{0}^{1 + \alpha \tilde{\eta}} \biggl\{ \frac{\alpha^2}{2} \bigl( \tilde{\phi }_{\tilde{x}}^3 + \frac{1}{\beta} \tilde{\phi}_{\tilde{z}}^2 \tilde{\phi}_{\tilde{x}} \bigr) + \tilde{z} \tilde{\phi}_{\tilde{x}} \biggr\} d\tilde{z} \\ &\quad{} + \alpha^2 \frac{\partial}{\partial\tilde{x}}\int _{0}^{1 + \alpha\tilde{\eta}} \tilde{P}' \tilde{\phi}_{\tilde{x}} \,d\tilde{z}+ \alpha\frac {\partial}{\partial\tilde{x}}\int _{0}^{1 + \alpha\tilde{\eta}} (1-\tilde{z}) \tilde{\phi}_{\tilde {x}} \,d\tilde{z}= 0. \end{aligned}$$

Substituting the expressions for \(\tilde{\phi}_{\tilde{x}}\) and \(\tilde{\phi}_{\tilde{z}}\) yields

$$\begin{aligned} &\frac{\partial}{\partial\tilde{t}}\biggl( \frac{1}{2} + \alpha \tilde{\eta}+ \alpha^2 \tilde{\eta}^2 + \frac{\tilde{U}^2}{2} + \frac{1}{2} \alpha \tilde{U}^2 \tilde{\eta}\mp\tilde{U}\alpha\tilde{\eta}\mp\frac{3}{4} \alpha^2 \tilde{U}\tilde{\eta}^2 \mp\frac{1}{6} \alpha\beta \tilde{U}\tilde{\eta}_{\tilde{x}\tilde{x}} \biggr) \\ &\quad{}+ \frac{\partial}{\partial\tilde{x}}\biggl( -\frac {1}{2} \tilde{U}^3(1+\alpha\tilde{\eta}) - \frac{5}{2} \alpha^2 \tilde{U}\tilde{\eta}^2 \pm \frac{3}{2} \alpha\tilde{U}^2 \tilde{\eta}\pm\frac{9}{8} \alpha^2 \tilde{U}^2 \tilde{\eta}^2 \pm\frac{1}{4} \alpha\beta\tilde{U}^2 \tilde{\eta}_{\tilde{x}\tilde{x}} \\ &\quad{} - 2 \alpha\tilde{U}\tilde{\eta}- \frac{1}{3} \alpha \beta\tilde{U} \tilde{\eta}_{\tilde{x}\tilde{x}} - \tilde{U}\pm\alpha\tilde {\eta}\pm\frac{7}{4} \alpha^2 \tilde{\eta}^2 \pm\frac{1}{6} \alpha\beta \tilde{\eta}_{\tilde{x}\tilde{x}} \biggr) = {\mathcal{O}}\bigl(\alpha^3, \alpha^2\beta,\alpha\beta^2\bigr). \end{aligned}$$

Defining the non-dimensional energy density and flux accordingly, the non-dimensional energy balance (1.7) is achieved. Using the scalings \(E = c_{0}^{2}h_{0} \tilde{E}\) and \(q_{E} = h_{0} c_{0}^{3} \tilde{q}_{E}\), the dimensional forms of these quantities are

$$ E = c_0^2 \biggl( \frac{h_0}{2} + \eta+ \frac{1}{h_0} \eta^2 + \frac{h_0}{2} \frac{U^2}{c_0^2} + \frac{1}{2} \frac{U^2}{c_0^2} \eta\mp\frac{U}{c_0} \eta\mp \frac{3}{4} \frac{U}{h_0 c_0} \eta^2 \mp\frac{1}{6} h_0^2 \frac{U}{c_0} \eta_{xx} \biggr), $$

and

$$\begin{aligned} q_E =& c_0^3 \biggl( - \frac{1}{2} \frac{U^3}{c_0^3}(h_0 + \eta) - \frac{5}{2h_0}\frac{U}{c_0} \eta^2 \pm\frac{3}{2} \frac{U^2}{c_0^2} \eta\pm \frac{9}{8h_0} \frac{U^2}{c_0^2} \eta^2 \\ & \pm\frac{1}{4} \frac{U^2}{c_0^2} h_0^2 \eta_{xx} - 2 \frac{U}{c_0} \eta- \frac{1}{3} \frac{U}{c_0} h_0^2 \eta_{xx} - \frac{U}{c_0} h_0 \pm\eta\pm\frac{7}{4h_0} \eta^2 \pm\frac{h_0^2}{6}\eta_{xx} \biggr). \end{aligned}$$

Note also here that the energy flux incorporates the work done by pressure forces. For later reference in the comparison with the quantity R used in Benjamin and Lighthill [4], we record that non-dimensional energy per unit mass is given by

$$\frac{\alpha^2}{2} \biggl( \tilde{\phi}_{\tilde{x}}^2 + \frac{1}{\beta} \tilde{\phi}_{\tilde{z}}^2 \biggr) + \tilde{z}= \frac{\alpha^2}{2} \tilde{f}_{\tilde{x}}^2 + \tilde{z}+ {\mathcal {O}}(\alpha^2\beta). $$

The relevant terms in this expression can be identified using the analysis of the energy balance law. Thus keeping only the terms required to achieve the energy balance, and evaluating at the free surface yields the non-dimensional total head as

$$ \tilde{H} = \frac{\tilde{U}^2}{2} + \frac{\alpha^2}{2} \tilde {\eta}^2 \mp \alpha\tilde{U}\tilde{\eta}\pm\frac{\alpha^2}{4} \tilde {U}\tilde{\eta}^2 \mp\alpha\beta \tilde{U}\biggl( \frac{1}{3} - \frac{1}{2} \biggr) \tilde{\eta }_{\tilde{x}\tilde{x}} + (1+\alpha\tilde{\eta}). $$

Using the scaling \(\tilde{H}= h_{0}H\), and multiplying by g, the energy density at the surface corresponding to the quantity R used in Benjamin and Lighthill [4] appears in dimensional form as

$$ gH = \frac{U^2}{2} + \frac{g \eta^2}{2h_0} + (h_0+ \eta) \mp\frac{gU}{c_0} \eta\pm\frac{1}{4} \frac{c_0 U}{h_0^2} \eta^2 \pm\frac{1}{6}c_0 h_0^2 U \eta_{xx}. $$
(5.1)

In some cases, it might be preferable to use a different normalization for the potential energy in order to isolate mechanical energy due to the wave motion. If the potential energy of a particle is taken to be zero at the undisturbed free surface, and it is required that the potential energy of the quiescent state be zero, then the energy balance can be defined by the equation

$$ \frac{\partial}{\partial t}\biggl\{ \int_{-h_0}^{\eta} \frac {1}{2} |\nabla \phi|^{2} \,dz + \int_{0}^{\eta} g z \,dz \biggr\} + \frac{\partial}{\partial x}\int_{-h_0}^{\eta} \biggl\{ \frac{1}{2} |\nabla\phi|^{2} + g z + P \biggr\} \phi_x dz = 0. $$

In non-dimensional form this equation becomes

$$\begin{aligned} &\frac{\partial}{\partial\tilde{t}}\biggl\{ \int_{0}^{1 + \alpha\tilde{\eta}} \frac{\alpha^2}{2} \bigl( \tilde{\phi}_{\tilde{x}}^2 + \frac{1}{\beta} \tilde{\phi }_{\tilde{z}}^2 \bigr) \,d\tilde{z}+ \int_{1}^{1 + \alpha\tilde {\eta}} (\tilde{z}-1) \,d\tilde{z}\biggr\} \\ &\quad{} + \frac{\partial}{\partial\tilde{x}}\biggl\{ \int_{0}^{1 + \alpha \tilde{\eta}} \frac{\alpha^3}{2} \bigl( \tilde{\phi}_{\tilde {x}}^3 + \frac{1}{\beta} \tilde{\phi}_{\tilde{z}}^2 \tilde{\phi }_{\tilde{x}} \bigr) \,d\tilde{z}+ \int _{0}^{1 + \alpha\tilde{\eta}} \alpha(\tilde{z}-1) \tilde{\phi }_{\tilde{x}} d\tilde{z} \biggr\} \\ &\quad{} + \frac{\partial}{\partial\tilde{x}}\biggl\{ \int_{0}^{1 + \alpha \tilde{\eta}}\alpha^2 \tilde{P}' \tilde{\phi}_{\tilde{x}} \,d\tilde{z}+ \int_{0}^{1 + \alpha\tilde{\eta}} \alpha(1-\tilde{z}) \tilde{\phi}_{\tilde{x}} \,d\tilde{z}\biggr\} = 0. \end{aligned}$$
(5.2)

Noticing a cancellation in (5.2), using the expression of \(\tilde{\phi}_{\tilde{x}}\) found in (2.12), and performing an integration with respect to \(\tilde{z}\), the equation becomes

$$\begin{aligned} &\frac{\partial}{\partial\tilde{t}}\biggl( \alpha^2 \tilde{\eta }^2 + \frac{\alpha^3}{4} \tilde{\eta}^3 + \frac{\alpha^2\beta}{6} \tilde{\eta}\tilde {\eta}_{\tilde{x}\tilde{x}} + \frac{\alpha^2\beta}{6} \tilde{\eta}_{\tilde{x}} \tilde{\eta }_{\tilde{x}} \biggr) \\ &\quad{}+ \frac{\partial}{\partial\tilde{x}}\biggl( \pm\alpha ^2 \tilde{\eta}^2 \pm \frac{5}{4}\alpha^3\tilde{\eta}^3 \pm\frac{\alpha^2\beta}{2} \tilde{\eta}\tilde{\eta}_{\tilde{x}\tilde{x}} \biggr) = {\mathcal{O}}\bigl(\alpha^4, \alpha^3\beta, \alpha^2\beta^2\bigr). \end{aligned}$$

The common factor α 2 can be omitted, and the differential energy balance equation is

$$\begin{aligned} \frac{\partial}{\partial\tilde{t}}\biggl( \tilde{\eta}^2 + \frac {\alpha}{4} \tilde{\eta}^3 + \frac{\beta}{6} \tilde{\eta}\tilde{\eta}_{\tilde{x}\tilde{x}} + \frac{\beta}{6} \tilde{\eta}_{\tilde{x}} \tilde{\eta}_{\tilde{x}} \biggr) + \frac{\partial}{\partial\tilde{x}}\biggl( \pm \tilde{\eta}^2 \pm\frac{5}{4}\alpha\tilde{\eta}^3 \pm \frac{\alpha\beta}{2}\tilde{\eta}\tilde{\eta}_{\tilde{x}\tilde {x}} \biggr) ={\mathcal{O}}\bigl( \alpha^2, \alpha\beta, \beta^2\bigr). \end{aligned}$$

As a result, the energy density should be given by

$$ \tilde{E}^{*}=\alpha^2 \tilde{\eta}^2 + \frac{\alpha^3}{4} \tilde{\eta}^3 + \frac{\alpha^2\beta}{6} \tilde{\eta}\tilde {\eta}_{\tilde{x}\tilde{x}} + \frac{\alpha^2\beta}{6} \tilde{\eta}_{\tilde{x}}^2 , $$

and the energy flux is given by

$$ \tilde{q}_{E}^{*}= \pm\alpha^2 \tilde{\eta}^2 \pm\frac{5}{4}\alpha^3\tilde{\eta}^3 \pm \frac{\alpha^2\beta}{2}\tilde{\eta}\tilde{\eta}_{\tilde {x}\tilde{x}} . $$

The dimensional quantities are given by

$$ E^{*} = c_0^2 \biggl( \frac{1}{h_0} \eta^2 + \frac{1}{4 h_0^2} \eta^3 + \frac{h_0}{6} \eta\eta_{xx} + \frac{h_0}{6} \eta_{x}^2 \biggr) $$

and

$$ q_{E}^{*} = \pm c_0^3 \biggl( \frac{1}{h_0} \eta^2 + \frac{5}{4 h_0^2} \eta^3 + \frac{h_0}{2} \eta\eta_{xx} \biggr). $$

One may wonder whether there is a relation between the total energy of the surface wave system and the invariant integrals of the KdV equation (1.2). The formula for E ∗ contains terms that look as though they might combine to a such a conservation law, but the coefficients do not line up in quite the right way. Note that the last computation was done in an fixed frame of reference in order to reach tidier expressions, and to stay in line with the initial requirement that the energy of the quiescent state be zero. As it turns out, it is possible to normalize the potential energy in such a way that the total energy in a reference frame moving at the speed U=c 0 is given by a combination of the conservation laws (1.2). This issue will be addressed in Sect. 7.

6 Comparison with Q, R and S

In the following, the formulae for q M , q I and gH derived in the previous sections are compared to the corresponding quantities Q, S and R, studied by Benjamin and Lighthill [4]. We take a periodic traveling wave propagating to the left at a speed c>0 in an inertial frame. In a reference frame also moving to the left at the velocity U=−c, the wave becomes steady, and yields a positive mass flux. The surface excursion can be described by a function ζ(x), and equation (1.1) reads

$$ (c_0 - c)\zeta^{\prime} + \frac{3}{2} \frac{c_0}{h_0} \zeta\zeta^{\prime} + \frac{c_0 h_0^2}{6} \zeta^{\prime\prime\prime} = 0. $$

The standard procedure of integrating, multiplying by ζ′, and then integrating again leads to

$$ \frac{gh_0^3}{3} \biggl( \frac{d\zeta}{dx} \biggr)^2 + g \zeta^3 +2c_0 ( c_0 - c)\zeta^2 + \mathcal{A}\zeta+ \mathcal{B} = 0, $$
(6.1)

where \(\mathcal{A}\) and \(\mathcal{B}\) are constants of integration. This differential equation has the solution

$$ \zeta(x)=\zeta_2 + (\zeta_1 - \zeta_2) \mbox{cn}^2 \biggl( \sqrt{\frac{3( \zeta_1-\zeta _3)}{4h_0^3}}x ; m \biggr), $$
(6.2)

which is given in terms of the Jacobian elliptic function cn with modulus \(m=\frac{\zeta_{1}-\zeta_{2}}{\zeta_{1}-\zeta_{3}}\). The numbers ζ 1, ζ 2 and ζ 3 are the three roots of the cubic polynomial appearing in (6.1), arranged in the order ζ 3<ζ 2<ζ 1. The constants of integration in (6.1) can be written in terms of the roots as \(\mathcal{A}= g(\zeta_{1}\zeta_{2} + \zeta_{1}\zeta_{3} +\zeta_{1}\zeta _{2})\) and \(\mathcal{B}=g \zeta_{1} \zeta_{2} \zeta_{3}\). The wavespeed is given by

$$ c= c_0 \biggl( 1 + \frac{\zeta_1+\zeta_2+ \zeta_3}{2h_0} \biggr), $$
(6.3)

and the wavelength is

$$ \lambda=K(m)\sqrt{\frac{16h_0^3}{3( \zeta_1-\zeta_3) }}, $$
(6.4)

where K(m) is the complete elliptic integral of the first kind. In the current setup, ζ 1 represents the wave crest, ζ 2 is the wave trough, and ζ 3 is a parameter which has influence only on the wavelength λ and wavespeed c. In the traveling reference frame, the quantities M, q M , I, q I , E and q E can now be computed as functions of x, and these are plotted in Fig. 2 for a particular case. Note that the mass flux q M is constant since the KdV equation features exact mass conservation. Moreover, the momentum density I is also constant since it is equal to the mass flux. The momentum flux is nearly constant, but features small variations which are visible when plotted at a finer scale.

Fig. 2
figure 2

A cnoidal wave with h 0=1.0631, ζ 1=0.4369, ζ 2=−0.1631, and ζ 3=−0.1731 (all in m), and the associated functions M, q M , I, q I , E and q E

Full size image

We now turn to the comparison of q M , q I and gH defined for the evolution problem with the corresponding quantities Q, S and R defined for the steady problem. In order to facilitate the comparison, let us briefly recall the development presented by Benjamin and Lighthill [4]. Steady periodic traveling waves are considered in a moving reference frame, in which the mass flux Q is positive, and the momentum flux S and energy per unit mass R are also given. In this case, the steady KdV equation appears as

$$ \frac{1}{3}Q^2 \biggl( \frac{d\xi}{dx} \biggr)^2 + g \xi^3 - 2R\xi^2 + 2S \xi- Q^2 = 0, $$
(6.5)

where ξ is the total flow depth. Using the same method as above, the solution is found as

$$ \xi(x)=\xi_2 + (\xi_1 - \xi_2) \mbox{cn}^2 \biggl( \sqrt{\frac{3(\xi_1- \xi_3)}{4\bigl(Q^2/g\bigr)}}x; m \biggr), $$

where the constants Q,R and S are given by

$$\begin{aligned} \begin{aligned} Q &= (g\xi_1\xi_2 \xi_3)^{\frac{1}{2}}, \\ R&= \frac{g}{2}(\xi_1+\xi_2+\xi_3), \\ S & = \frac{g}{2}(\xi_1\xi_2+ \xi_1 \xi_3 +\xi_2\xi_3), \end{aligned} \end{aligned}$$
(6.6)

and ξ 3<ξ 2<ξ 1 are the roots of the cubic polynomial gξ 3−2Rξ 2+2Sξ−Q 2.

Now since ζ represents the deflection of the fluid surface from rest while ξ is the total flow depth, the solutions of (6.1) and (6.5) must be related by

$$\zeta=\xi-h_0. $$

By the same token, the solution parameters are related by

$$\begin{aligned} \zeta_1=\xi_1-h_0,\quad \zeta_2=\xi_2-h_0,\quad \mbox{and}\quad \zeta_3=\xi_3-h_0. \end{aligned}$$
(6.7)

Moreover, by comparing the coefficients of the equations, it appears that h 0=(Q 2/g)1/3, and the total head R can be expressed in terms of the wavespeed as

$$\begin{aligned} R= c_0 \biggl( \frac{c_0}{2} + c \biggr). \end{aligned}$$

One may now freely choose ξ 1,ξ 2, and ξ 3 and calculate Q, S and R from (6.6), and then use (6.7) and (6.2) to compute the corresponding values of q M (ζ), q I (ζ) and gH(ζ) as defined in the previous sections. For example, in Fig. 3, ξ 1=1.4 m , ξ 2=1 m and ξ 3=0.95 m are chosen which give h 0=1.1 m, ζ 1=0.3 m, ζ 2=−0.1 m and ζ 3=−0.15 m. The wavelength is λ=10.04 m, and the wave amplitude is a=0.2 m. Besides showing the wave profile, Fig. 3 features a comparison of the quantities q M , q I and gH as defined by (3.2), (4.1) and (5.1), respectively with the corresponding parameters Q, S and R. As can be seen, the difference between Q and q M (ζ), the difference between S and q I (ζ), and also the difference between R and gH(ζ) are all reasonably small. In order to further quantify these differences, waves with various combinations of the parameters ξ 1,ξ 2, and ξ 3 are computed, and the differences in the above quantities are plotted as functions of the two small parameters \(\alpha= \frac{a}{h_{0}}\) and \(\beta = \frac {h_{0}^{2}}{\lambda^{2}}\). For comparison of different values of α and β, it appears most convenient to keep the water depth h 0 constant. Then the mass flux Q=c 0 h 0 is also constant. The first equation in (6.6) implies that \(\xi_{3}= \frac {h_{0}^{3}}{\xi_{1} \xi_{2}}\) so the problem now depends only on the wave crest ξ 1 and the wave trough ξ 2. The condition ξ 3<ξ 2<ξ 1 leads to

$$ \frac{h_0^3}{\xi_1 \xi_2} <\xi_2 < \xi_1. $$
(6.8)

It is easy to see from (6.8) that wave crest is always larger than the water depth (ξ 1>h 0) and the wave trough is bounded below by \(\xi_{2} > h_{0}\sqrt{\frac{h_{0}}{\xi_{1}}}\). Consequently, the total head R will be restricted as \(\frac{3g}{2} \frac{h_{0}^{3}}{\xi_{1}^{2}} < R < \frac{3g}{2} \xi_{1}\), while the momentum flux will be bounded as \(\frac{3g}{2} \frac{h_{0}^{6}}{\xi_{1}^{4}} < S < \frac{3g}{2} \xi_{1}^{2}\).

Fig. 3
figure 3

The surface excursion ζ is plotted in panel (a). Panel (b) shows the mass flux q M vs. Q, panel (c) shows the momentum flux q I vs. S, and panel (d) shows the energy per unit mass gH vs. R. The values h 0=1.1, ζ 1=0.3, ζ 2=−0.1 and ζ 3=−0.15 (all in m) were chosen

Full size image

The top panels of Fig. 4 show the error between q M and Q. The center panels of Fig. 4 show the error between q I and S. The lower panels of Fig. 4 show the error between gH and R. These errors are plotted as level curves with respect to the small parameters α and β with Q held fixed. It can be seen clearly in all cases that the error diminishes with decreasing values of α and β.

Fig. 4
figure 4

The absolute errors in the approximations of Q, R and S, in terms of the model parameters α and β (left), and in terms of the wave crest ξ 1 and wave trough ξ 2 (right)

Full size image

7 Exact Conservation

In the following, we address exact conservation of mass, momentum and energy in KdV evolution. It is assumed that solutions are smooth, and that the wave motion is localized in the sense that the function η(x,t) describing the free surface decays rapidly enough as x→±∞ so that all integrals appearing here are defined. First of all, total mass can only be defined on a finite interval, so that the most preferable form to state mass conservation is (1.3). However, one may define excess mass by \(\int_{-\infty}^{\infty} \eta\, dx\), and it clearly follows from (1.3) that this quantity is constant with respect to t.

Second, in an inertial frame with U=0, and for a localized surface disturbance, the total horizontal momentum

$$\mathcal{I} = \int_{-\infty}^{\infty} I \, dx = c_0 \int_{-\infty}^{\infty} \eta \, dx + \frac{3c_0}{4h_0} \int_{-\infty}^{\infty} \eta^2 \, dx $$

can be defined, and it follows from the first two invariant integrals in (1.2) that \(\mathcal{I}\) is conserved.

Finally, let us discuss conservation of total energy. We use the expression for the horizontal velocity in a moving frame of reference provided by (2.14), and normalize the potential energy as in the second case in Sect. 5, so that the energy density is given by

$$ \int_{-h_0}^{\eta} \frac{1}{2} |\nabla \phi|^{2} \,dz + \int_{0}^{\eta} g z \,dz. $$

Focusing on waves propagating to the right and requiring (1.7), an analysis along the lines of the argument shown in Sect. 5 yields the following expression for the non-dimensional energy density in the moving reference frame:

$$\begin{aligned} \tilde{E}^{*}_U =& \frac{1}{2} \tilde{U}^2 + \frac{\alpha}{2} \tilde{U}^2 \tilde{\eta}+ \alpha^2 \tilde{\eta}^2 + \frac{1}{4} \alpha^3 \tilde{\eta}^3 - \alpha\tilde{U}\tilde{\eta}- \frac{3}{4} \alpha^2 \tilde{U} \tilde{\eta}^2 + \frac{1}{4} \alpha^3 \tilde{U} \tilde{\eta}^3 \\ &{} - \frac{1}{2} - \frac{1}{6} \alpha\beta\tilde{U} \tilde{\eta}_{\tilde{x}\tilde{x}} + \frac{1}{6} \alpha^2 \beta \tilde{U}\tilde{\eta} \tilde{\eta}_{\tilde{x}\tilde{x}} + \frac{1}{6} \alpha^2 \beta \tilde{\eta} \tilde{\eta}_{\tilde{x}\tilde{x}} + \frac{1}{6} \alpha^2 \beta \tilde{\eta}_{\tilde{x}}^2. \end{aligned}$$

Working now in a reference frame moving at the limiting long-wave speed, we have U=c 0, so that \(\tilde{U}= 1\). In this particular case, the energy density is given in dimensional form by

$$ E^{*}_{c_0} = c_0^2 \biggl( - \frac{1}{2} \eta+ \frac{1}{4 h_0} \eta^2 + \frac{1}{2 h_0^2} \eta^3 - \frac{h_0^2}{6} \eta_{xx} + \frac{h_0}{3} \eta\eta_{xx} + \frac{h_0}{6} \eta_{x}^2 \biggr). $$

The KdV equation in the reference frame moving at the speed c 0 is

$$\begin{aligned} \eta_t+ \frac{3}{2} \frac{c_0}{h_0} \eta \eta_x + \frac{c_0 h_0^2}{6} \eta_{xxx} = 0, \end{aligned}$$

and we note that this equation has the same conserved integrals (1.2).

The total mechanical energy of a localized surface wave is given in terms of the dimensional energy density \(E^{*}_{c_{0}}\) by \(\mathcal{E} = \int_{-\infty}^{\infty} E^{*}_{c_{0}} \, dx\). Substituting the expression for \(E^{*}_{c_{0}}\) yields

$$ \mathcal{E} = - \frac{1}{2} c_0^2 \int _{-\infty}^{\infty} \eta\, dx + \frac{1}{4} \frac{c_0^2}{h_0} \int_{-\infty}^{\infty} \eta^2 \, dx - \frac{1}{6} c_0^2 h_0^2 \int_{-\infty}^{\infty} \eta_{xx} \, dx + \frac{1}{2} \frac{c_0^2}{h_0^2} \int_{-\infty}^{\infty} \bigl\{ \eta^3 - \frac{h_0^3}{3} \eta_x^2 \bigr\} \, dx. $$

After observing that the third integral in the expression above vanishes, it becomes plain that the three invariant integrals in (1.2) guarantee that the total energy in the KdV approximation is exactly conserved in this case.

8 Conclusion

In this article, expressions for mass, momentum and energy densities and fluxes which are valid in the KdV approximation have been found. The quantities have been compared to the quantities Q, R and S which were previously derived by Benjamin and Lighthill in the steady case [4]. It has also been shown that exact conservation of total mass, momentum and energy holds in special cases. For the exact conservation, the mathematical formulations of the first three conservation laws (1.2) have been used. The main result of the paper is the identification of the quantities M, q M , I, q I , E and q E in the context of the KdV approximation. However, as already mentioned in the introduction, the method used in this paper is a formal one, and the results presented here should be understood as a first step towards a mathematical procedure which will give a definite proof that the balances (1.5), (1.6) and (1.7) are valid to the same order and over the same time scales as the KdV equation (1.4) itself. In order to provide such a proof one might follow the procedure pioneered by Craig [9], and recently refined by Bona et al. [6]. In this latter work, the Hamiltonian formulation of Zakharov and a careful analysis of the Dirichlet-Neumann operator, such as defined by Craig and Sulem [11] play a prominent role. The general procedure has been further extended and applied to the justification of a variety of simplified model equations and systems, and the recent monograph by Lannes [20] contains a large variety of different cases.

While the methods for a mathematical justification of the derivation of many model equations are available, it is not entirely clear how to apply them to the justification of the associated balance laws treated in the present article. Such a study will be an interesting topic for future work.