1 Introduction

The Korteweg–de Vries (KdV) equation is a model equation describing the evolution of long-crested waves at the surface of a body of fluid. The equation is derived as a physical model equation under the assumption that there is an approximate balance between nonlinear steepening effects and dispersive spreading. In mathematical terms, this balance is expressed by introducing two small parameters \(\alpha \) and \(\beta \) measuring the wave amplitude and the wave length, respectively. Suppose that the undisturbed depth of the fluid is given by \(h_0\). If represents a typical wave amplitude, and \(\ell \) represents a typical wavelength, the two parameters and \(\beta = h_0^2 / \ell ^2\) should be small and of the same order. This is the Boussinesq scaling. If in addition the wave motion is predominantly in a single direction, then the KdV equation is an approximate model describing the surface wave motion [9, 30]. If the undisturbed depth \(h_0\) is taken as a unit of distance and \(\sqrt{h_0/g}\) is taken as a unit of time, then the KdV equation appears in the form

$$\begin{aligned} \eta _t + \eta _x + { \textstyle \frac{3}{2}} \eta \eta _x + { \textstyle \frac{1}{6}} \eta _{xxx} = 0. \end{aligned}$$
(1.1)

In the present article, the focus is on whether—or rather how—the KdV equation is able to describe incipient wave breaking. While the KdV equation does not admit the distinctive steepening and development of infinite gradients known from nonlinear hyperbolic equations, it will be shown here that the KdV equation features a different kind of wave breaking which is more closely related to spilling at the wave crest.

Let us recall that the KdV equation can be thought of as a combination of the simple nonlinear balance law

$$\begin{aligned} \eta _t + \eta _x + { \textstyle \frac{3}{2}} \eta \eta _x = 0, \end{aligned}$$
(1.2)

and the linear dispersive equation

$$\begin{aligned} \eta _t + \eta _x + { \textstyle \frac{1}{6}}\eta _{xxx} = 0. \end{aligned}$$
(1.3)

As indicated in Fig. 1, all waves featuring negative slope break at some point in the model (1.2), but no waves break in the model (1.3) (however, (1.3) features dispersive blow-up for some data [10, 27, 34]). The KdV Eq. (1.1) allows a balance of nonlinear steepening effects and dispersive spreading which arrests the typical hyperbolic wave breaking exhibited by (1.2), and leads to the formation of steady traveling waves, such as solitary and cnoidal waves [1, 29] which propagate without a change in the wave profile.

In the context of the full Euler equations, the water-wave problem also admits various steady traveling-wave solutions. In this case, wave breaking can be induced by a number of instabilities, such as the superharmonic instability [43, 44] and transverse instabilities [22, 25]. In the KdV equation, both the solitary and cnoidal wave solutions are known to be stable [6, 14], and therefore wave breaking triggered by various instabilities is not reproduced by the KdV model.

Solitary-wave solutions of (1.1) are of the form

$$\begin{aligned} \eta (x,t) =H \text{ sech }^2\big ( { \textstyle \frac{\sqrt{3H}}{2}}(x-x_0-ct)\big ), \end{aligned}$$
(1.4)

where the phase velocity is given by \(c=1+\frac{H}{2}.\) While these formulas define solutions of the KdV equation for all waveheights H, it turns out that solutions of large waveheight are inconsistent with the model in which the KdV equation is valid. In particular, as will be shown in Sect. 2, for waveheights exceeding the critical waveheight \(H_{\text {max solitary}} = 0.6879\), the phase velocity of the wave is smaller than the particle velocity at the crest, a fact which may be interpreted as incipient wave breaking.

Fig. 1
figure 1

Left panel: Solution of Burgers equation \(\eta _t + { \textstyle \frac{3}{2}}\eta \eta _x = 0\) which features hyperbolic wave breaking. Right panel: Solution of the Airy equation \(\eta _t + { \textstyle \frac{1}{6}}\eta _{xxx}\) which features dispersive spreading. Initial data are shown as dashed curves. If the dashed curve were used as initial data for the KdV Eq. (1.1), then this wave would propagate without change in shape

Full size image

To elaborate on this argument, note that it is shown in [46] that the derivation of the equation as a surface water-wave model enables the reconstruction of an approximation of the fluid velocity field underneath the surface. In particular, it is possible to derive relations expressing the horizontal and vertical velocity components in terms of the principal unknown variable \(\eta \) which describes the shape of the free surface. The horizontal velocity component in the fluid is given by

$$\begin{aligned} u(x,y,t) = \eta - { \textstyle \frac{1}{4}} \eta ^2 + \big ({ \textstyle \frac{1}{3}} - { \textstyle \frac{y^2}{2}} \big ) \eta _{xx}, \end{aligned}$$
(1.5)

and this relation holds to the same order in the asymptotic parameters \(\alpha \) and \(\beta \) as to which the KdV equation is valid.

Since the horizontal component of the particle velocity can be found approximately, it may be compared to the local phase velocity of the wave. This leads to one of the most fundamental breaking criteria used in the literature, the convective breaking criterion which predicts wave breaking if the particle velocity at the wavecrest exceeds the phase velocity of the wave (see for example [32, 45]). More specifically, this breaking criterion (often termed kinematic criterion) is based on the local Froude number, and predicts breaking if

$$\begin{aligned} \frac{U}{c} \ge 1, \end{aligned}$$
(1.6)

where \(U = u(x,\eta (x,t),t)\) is the horizontal component of the velocity field evaluated at the free surface, c is the local phase velocity of the wave. To illustrate this point of view, Fig. 2 shows two waves with corresponding fluid particles near the crest of the wave. In the left panel, a wave with sufficiently small amplitude is depicted. A fluid particle located near the wavecrest remains in the free surface and recedes from the crest in unison with the wave motion. However in the right panel, a wave of higher amplitude is shown, and particles located near the crest have higher horizontal velocity than the wave itself.

The utility of the convective criterion has been the subject of some discussions in the literature. For wave breaking in shallow water, such as in focus in the present article, some studies such as [47] conclude that the convective criterion does a fair job of predicting wave breaking even for three-dimensional waves. On the other hand, studies focusing on breaking of deep-water waves such as [40] observe that the kinematic criterion has been useful in some practical situations, while there is also evidence that it may not be a reliable indicator to pinpoint the onset of breaking in the case of deep water. In some of the cases where the convective breaking criterion performs badly, the difficulty of obtaining accurate readings for the phase velocity of the wave in experimental situations may be responsible [39].

Fig. 2
figure 2

Left panel: Surface profile of a regular wave, and trajectory of a particle in the free surface. Right panel: Surface profile of a breaking wave. A particle contained in the free surface is leaving the fluid domain, indicating incipient wave breaking

Full size image

Note that the criterion (1.6) is stated in sufficient generality to be applicable to both experimental and numerical work. In the case of numerical modeling, and in particular in the case of phase-resolving models, there is no difficulty in evaluating the local wave and particle velocities. The local phase velocity can usually be approximated using a variety of methods. For example, Fourier techniques have been used in [36]. What is more, in the case of solitary and cnoidal waves, the phase and particle velocities are known in closed form, and it is straightforward to test the breaking criterion (1.6). These computations will be carried out in Sects. 2 and 3.

Fig. 3
figure 3

Left panel: Breaking wave in Boussinesq model with sloping bottom. Right panel: Wave breaking due to a large discharge at left boundary

Full size image

While it is relatively easy to evaluate the breaking criterion in several situations for the Eq. (1.1), it should also be noted that neither the hyperbolic steepening nor the convective wave breaking is likely to happen for free surface waves unless some focusing effect is present, such as an underlying current [40, 48], strong three-dimensional effects [47], or specially tailored initial data [10]. On the other hand wave breaking is likely to occur in the presence of forcing in the form of an uneven bottom topography or a discharge. These two possibilities are sketched in Fig. 3. In particular the left panel of Fig. 3 indicates one of the most widely known type of wave breaking, the development of breakers on a sloping beach. In this case, depending on a number of parameters such as waveheight, wavelength and bottom slope, a variety of breaking phenomena takes place, ranging from spilling at the crest, to overturning and surging breaking [19].

In the right panel of Fig. 3, the case of a forced inflow is indicated. In this case, a so-called bore generally emerges. Sufficiently small discharges lead to the appearance of an undular bore which is a characterized by a moderately steep front followed by a leading wave and trailing smaller waves. For larger inflows, the leading wave behind the bore may break, and for large enough inflows, the bore can be entirely turbulent [15, 28]. In the final section of the present paper, the case of a forced inflow and undular bore is discussed in some detail, and the question of the transition from purely undular to undular with spilling breaking is investigated. The numerical results are compared to experimental findings of Favre [21], who determined a critical inflow condition dividing between laminar flows and flows featuring breaking waves. Our numerical simulations show incipient breaking near the critical case of Favre, but some discrepancy remains. Further details on this comparison are given in Sect. 4.

Before we enter the main development of this paper, some remarks are in order. First, it should be noted that convective wave breaking of the type discussed here is unlikely to happen in hyperbolic equations or systems. While the simple conservation law \(\eta _t + \frac{3}{2}\eta \eta _x = 0\) does not contain expressions for the local particle velocity, the nonlinear shallow-water equations

$$\begin{aligned} h_t + (uh)_x = 0, \\ u_t + h_x + uu_x = 0, \end{aligned}$$

include both the surface deflection \(\eta (x,t) = h(x,t)-1\), and an average horizontal particle velocity u(x, t) in the description. In this case, a simple wave propagating to the right is given by the Riemann invariant \(u+2 \sqrt{h}\), and has phase velocity \(u+\sqrt{h}\). As shown in [46], for a wave profile \(h=H(x)\) with undisturbed depth 1, the fluid velocity is \(u = 2 \sqrt{H}-2\), and the phase velocity of the wave is \(C = 3 \sqrt{H}-2\). Thus we see that the phase velocity is always larger than the fluid velocity in this case. On the other hand, the typical hyperbolic steepening and eventual breaking will happen for appropriate initial conditions.

Secondly, it is important to note that wave breaking can be induced by instability even in cases without any forcing. This mechanism is well known in the case of waves on deeper water. Traditionally the wave steepness has been used as a predictor for the occurrence of crest instabilities and subsequent breaking of deep water waves [40], but recent work has indicated that criteria based on the energy flux in wave groups may be successful [38]. Breaking induced by crest instabilities may also occur in shallow water, such as in the situation at hand. First of all, solitary waves exist in the fully nonlinear Euler equations up to a waveheight of \(H=0.8332\). However, in [43], it was shown that the solitary wave is unstable in the context of the full Euler potential theory for waveheights exceeding 0.78, and the authors of [44] showed how this instability can be related to wave breaking. In [25], a three-dimensional situation was considered, and transverse instabilities were found at a nondimensional waveheight \(H = 0.713\). For periodic wave in shallow water, similar results were provided in [33] and [22].

Finally, it should be mentioned that the relation (1.5) and similar expressions for the vertical velocity may also be used advantageously for a number of purposes, such as the definition of a family of Boussinesq models [8], describing particle paths beneath the surface [12], and the study of mechanical balance laws associated to the evolution equation [2, 3]. On the other hand, there is a variety of different strategies to construct the velocity field in the fluid, such as the analytic method used in [23] to find the velocities associated to a periodic cnoidal solution of the KdV equation. Particle paths under periodic traveling waves can also be found to rather high accuracy using the Lagrangian approach [18]. Concerning the notion of mechanical balance laws, we mention in particular that an analysis along these lines casts doubt on the conservation of mechanical energy in the KdV approximation [3]. The conservation of energy in both the KdV and higher-order KdV equations has also been questioned recently in a study which utilizes a Lagrangian framework [24].

The plan of the paper is as follows. In Sect. 2, the convective breaking criterion is applied to the solitary wave, and it is found that the limiting waveheight is \(H_{\text {max solitary}} = 0.6879\). In Sect. 3, the breaking criterion is applied to the cnoidal wave solutions. It is found that cnoidal wave solutions may feature breaking at a much smaller waveheight than the solitary wave. Finally, Sect. 4 is devoted to the study of wave breaking in a dynamic situation with a forced inflow. The resulting undular bore features a leading wave which may break, and the focus is on finding the critical case dividing the field of purely undular bores from partially turbulent breaking bores.

2 Convective wave breaking for solitary waves

The surface water-wave problem for an inviscid fluid is generally described by the Euler equations with an impenetrable bottom, and kinematic and dynamic boundary conditions at the free surface. Assuming weak transverse effects, the unknowns are the surface elevation \(\eta (x,t)\), the horizontal and vertical fluid velocities u(x, y, t) and v(x, y, t), respectively, and the pressure P(x, y, t). If the assumption of irrotational flow is made, then a velocity potential \(\phi (x,y,t)\) can be used. In the non-dimensional variables mentioned in the introduction, the problem may be posed on the domain \(\left\{ (x,y) \in \mathbb {R}^2 \, | \, 0< y < 1+\eta (x,t) \right\} \). The derivation of the KdV equation is based on approximating the velocity potential by an asymptotic series. In suitably scaled variables, the horizontal velocity component takes the form

$$\begin{aligned} u = \phi _x = w - \beta { \textstyle \frac{y^2}{2}} w_{xx} + O(\beta ^2), \end{aligned}$$
(2.1)

where w(x, t) represents the horizontal velocity at the bottom. On the other hand, the one-way assumption inherent in the KdV equation necessitates that the velocity at the bed is given by

$$\begin{aligned} w = \eta - { \textstyle \frac{1}{4}} \alpha \eta ^2 + { \textstyle \frac{1}{3}} \beta \eta _{xx} + O(\alpha ^2, \alpha \beta , \beta ^2). \end{aligned}$$

By combining these two relations and neglecting terms of \(O(\alpha ^2, \alpha \beta , \beta ^2)\), we obtain an expression for the horizontal component of the velocity field in terms of \(\eta \) and y:

$$\begin{aligned} u = \phi _x = \eta - { \textstyle \frac{1}{4}} \alpha \eta ^2 + \beta \big ({ \textstyle \frac{1}{3}} - { \textstyle \frac{y^2}{2}}\big )\eta _{xx}. \end{aligned}$$

In the original variables, the horizontal velocity is then expressed as

$$\begin{aligned} u(x,y,t) = \eta - { \textstyle \frac{1}{4}} \eta ^2 + \big ({ \textstyle \frac{1}{3}} - { \textstyle \frac{y^2}{2}}\big ) \eta _{xx}, \end{aligned}$$
(2.2)

and the breaking criterion (1.6) can be written as

$$\begin{aligned} \eta - { \textstyle \frac{1}{4}} \eta ^2 + \big ({ \textstyle \frac{1}{3}} - { \textstyle \frac{(1 + \eta )^2}{2}}\big ) \eta _{xx} \ge c. \end{aligned}$$
(2.3)

We now apply this breaking criterion to the solitary wave given by (1.4), where H is the waveheight, c is the wavespeed, and \(x_0\) is the initial location of the wave crest. Differentiating twice with respect to x yields

$$\begin{aligned} \eta _{xx}= & {} 3 H^2 {{\mathrm{sech}}}^2 \big ({ \textstyle \frac{\sqrt{3 H}}{2}} (x - ct - x_0) \big ) \tanh ^2 \big ( { \textstyle \frac{\sqrt{3 H}}{2}} (x - ct - x_0) \big ) \\&- { \textstyle \frac{3}{2}} H^2 {{\mathrm{sech}}}^4 \big ({ \textstyle \frac{\sqrt{3 H}}{2}} (x - ct - x_0) \big ). \end{aligned}$$

Since the solitary wave retains its shape for all time we may evaluate at \(x-ct-x_0 = 0\), in which case we obtain \(\eta = H\) and \(\eta _{xx} = -{ \textstyle \frac{3}{2}} H^2\). Substituting these values into Eq. (2.2) and evaluating at the surface yields the breaking criterion (1.6)

$$\begin{aligned} H - { \textstyle \frac{1}{4}} H^2 - { \textstyle \frac{3}{2}} \big ( { \textstyle \frac{1}{3}} - { \textstyle \frac{(1 + H)^2}{2}} \big ) H^2 \ge 1 + { \textstyle \frac{1}{2}}H. \end{aligned}$$
(2.4)

Setting the left hand side equal to the right hand side and rearranging terms yields

$$\begin{aligned} \mathcal {P}(H) = { \textstyle \frac{3}{4}} H^4 + { \textstyle \frac{3}{2}} H^3 + { \textstyle \frac{1}{2}} H - 1 = 0 \end{aligned}$$
(2.5)

for the critical value of the waveheight. We see that \(\mathcal {P}\) is a fourth order polynomial in H, and since \( \mathcal {P}'(H) = 3 H^3 + \frac{9}{2} H^2 + \frac{1}{2} > 0\) for \(H \ge 0\) and \(\mathcal {P}(0) < 0\) while \( \mathcal {P}(1) > 0\), it can have only one positive root which lies in [0, 1]. This root can be found numerically to obtain the following value of the maximum allowable waveheight for the solitary wave:

$$\begin{aligned} H_{\text {max solitary}} = 0.6879. \end{aligned}$$
(2.6)

Even though the KdV equation admits solutions with any waveheight, waves with a waveheight larger than \(H_{\text {max solitary}}\) do not describe actual surface waves since these waves already feature incipient wave breaking. Of course the waveshape of a solitary-wave solution of the KdV equation may not be close to a surface water wave already for smaller waveheights than \(H_{\text {max solitary}}\) (cf. [13]).

3 Maximum waveheight for periodic waves

The cnoidal waves are defined with the help of the incomplete elliptic integral of the first kind. For a given elliptic parameter \(0<m<1\), this integral is

$$\begin{aligned} u = \int _0^{\phi } { \textstyle \frac{d\theta }{\sqrt{1 - m\sin ^2\theta }}}. \end{aligned}$$

The elliptic functions \({{\mathrm{sn}}}(\cdot | m)\), \({{\mathrm{cn}}}(\cdot | m)\) and \({{\mathrm{dn}}}(\cdot |m)\) are defined as

$$\begin{aligned} \begin{aligned} {{\mathrm{sn}}}(u | m)&= \sin \phi , \\ {{\mathrm{cn}}}(u | m)&= \cos \phi , \\ {{\mathrm{dn}}}(u |m)&= \sqrt{1 - m\sin ^2\phi }. \end{aligned} \end{aligned}$$

On the other hand, the complete elliptic integrals of the first and second kind, respectively, are defined by

$$\begin{aligned} K(m) = \int _0^{\pi /2} { \textstyle \frac{d\theta }{\sqrt{1 - m\sin ^2\theta }}} \ \ \ \text { and } \ \ E(m) = \int _0^{\pi /2}\sqrt{1 - m\sin ^2\theta } \ d\theta . \end{aligned}$$

As shown in [20], the cnoidal travelling wave solutions of (1.1) are given in terms of three parameters \(f_1\), \(f_2\) and \(f_3\) by

$$\begin{aligned} \begin{aligned} \eta (x,t)&= f_2 + (f_1 - f_2) {{\mathrm{cn}}}^2 \Big ( { \textstyle \frac{\sqrt{3(f_1 - f_3)}}{2}} (x - ct - x_0) \big | m \Big ), \end{aligned} \end{aligned}$$
(3.1)

where \(H = f_1 - f_2\) is the wave height, \(c = 1 + \frac{1}{2}(f_1 + f_2 + f_3)\) is the wavespeed, and \(m=(f_1-f_2)/(f_1-f_3)\) is the elliptic parameter. Any cnoidal wave is completely determined as long as these three constants are fixed. If \(\sigma \) is defined by \(\sigma ^2 = \frac{4}{3(f_1 - f_3)}\), then the wavelength can be written as \(\lambda = 2 \sigma K(m)\).

Since the wave oscillates about the equilibrium level of the surface, the mean value of the surface displacement over one wavelength should be zero. This requirement can be expressed by integrating over one wavelength, viz.

$$\begin{aligned} \int _0^{\lambda } \eta \, dx = 0, \end{aligned}$$

and it can be shown that this integral is

$$\begin{aligned} \int _0^{\lambda } \eta \, dx = 2 \sigma \Big ( (f_1 - f_3) E(m) + f_3 K(m) \Big ) = 2 \sigma \Big ( { \textstyle \frac{f_1-f_2}{m}}E(m) + f_3 K(m) \Big ). \end{aligned}$$
(3.2)

Since \(H = f_1 - f_2\), we therefore have \(f_3 = -\frac{H}{m} \frac{E(m)}{K(m)}\). The parameters \(f_1\) and \(f_2\), can then also be expressed in terms of wave height H and the elliptic parameter m as follows:

$$\begin{aligned} \begin{aligned} f_1&= { \textstyle \frac{H}{m}} \left( 1 - { \textstyle \frac{E(m)}{K(m)}}\right) , \\ f_2&= { \textstyle \frac{H}{m}}\left( 1 - m - { \textstyle \frac{E(m)}{K(m)}}\right) . \end{aligned} \end{aligned}$$

Thus, fixing H and m is sufficient to specify any cnoidal wave, as long as the undisturbed water level is set to zero.

To determine whether a given cnoidal wave is a reasonable description of a water wave in the KdV approximation, we test the wave breaking criterion (2.3). By differentiating (3.1) twice with respect to x we obtain

$$\begin{aligned} \begin{aligned} \eta _{xx} =&{ \textstyle \frac{3}{2}} (f_1 - f_2)(f_1 - f_3) \Big \{ {{\mathrm{sn}}}^2 \big ({ \textstyle \frac{\zeta }{\sigma }}|m \big ) {{\mathrm{dn}}}^2 \big ( { \textstyle \frac{\zeta }{\sigma }} |m \big ) \\&+m^2 {{\mathrm{sn}}}^2\left( { \textstyle \frac{\zeta }{\sigma }} |m\right) {{\mathrm{cn}}}^2\left( { \textstyle \frac{\zeta }{\sigma }} |m\right) - {{\mathrm{cn}}}^2\big ({ \textstyle \frac{\zeta }{\sigma }} |m\big ) {{\mathrm{dn}}}^2\big ({ \textstyle \frac{\zeta }{\sigma }}|m \big ) \Big \}, \end{aligned} \end{aligned}$$

where the argument is given by \(\zeta = x-ct-x_0\). Evaluating at \(\zeta = 0\) yields

$$\begin{aligned} \begin{aligned} \eta&= f_1, \\ \eta _{xx}&= -{ \textstyle \frac{3}{2}}(f_1 - f_2)(f_1 - f_3). \end{aligned} \end{aligned}$$

Substituting these expressions into Eq. (2.3) transforms the breaking criterion (2.3) into

$$\begin{aligned} f_1 - { \textstyle \frac{1}{4}}f_1^2 - { \textstyle \frac{3}{2}}\left( { \textstyle \frac{1}{3}} - { \textstyle \frac{(1 + f_1)^2}{2}}\right) (f_1 - f_2)(f_1 - f_3) \ge 1 + { \textstyle \frac{1}{2}}(f_1 + f_2 + f_3). \end{aligned}$$
(3.3)

Defining the constants

$$\begin{aligned} a = { \textstyle \frac{1}{m}}\left( 1 - { \textstyle \frac{E(m)}{K(m)}}\right) , \ b = { \textstyle \frac{1}{m}} \left( 1 - m - { \textstyle \frac{E(m)}{K(m)}}\right) \ \text { and } \ c = { \textstyle \frac{1}{m}} { \textstyle \frac{E(m)}{K(m)}}, \end{aligned}$$

we can write the constants \(f_1\), \(f_2\) and \(f_3\) as

$$\begin{aligned} f_1 = Ha, \ f_2 = Hb \ \text { and } \ f_3 = -Hc. \end{aligned}$$

Substituting these definitions into Eq. (3.3) and setting the left and right hand side equal yields

$$\begin{aligned} \begin{aligned} Q_m(H) = \frac{3}{4}&H^4(a^4 + a^3c - a^3b - a^2bc) \\ + \ \frac{3}{2}&H^3(a^2c + a^3 - a^2b - abc) \\ + \ \frac{1}{4}&H^2(ac - bc - ab) \\ + \ \frac{1}{2}&H(a - b + c) - 1 = 0. \end{aligned} \end{aligned}$$
(3.4)

This is a fourth-order polynomial in H, and by fixing a value for m it can be solved numerically to obtain the maximum allowable wave height for the cnoidal wave, \(H_{\text {max cnoidal}}(m)\). Since the periodic wave reduces to the solitary wave in the nonlinear limit (\(m \rightarrow 1\)), and \(m \rightarrow 0^+\) is the linear limit, we search for real roots of (3.4) in the interval [0, 1].

Different values of m and corresponding values of \(H_{\text {max cnoidal}}(m)\) are listed in Table 1. The table also shows corresponding values of wave length, \(\lambda \) and the dimensionless parameters \(\alpha \) and \(\beta \). In addition to Stokes’ number \(\mathcal {S}\), defined by the ratio \(\alpha /\beta \) is tabulated.

Figure 4 shows a plot of \(H_{\text {max cnoidal}}(m)\) obtained using a finer resolution in m than displayed in the table above. It can be seen that the maximum wave height seems to approach zero as m approaches zero, which fact is related to a going ’outside’ of the domain of validity of the Boussinesq approximation, where \(\alpha \) and \(\beta \) are assumed to be of the same order of magnitude. It is also worth noting that the limiting values of a, b and c as \(m \rightarrow 1\) are:

$$\begin{aligned} a = 1, \ b = 0, \ \text { and } \ c = 0 \end{aligned}$$

Applying this to Eq. (3.4) gives the breaking criterion for the solitary wave, Eq. (2.5). This is also easy to see if we plot the maximum wave height of the cnoidal wave as a function of wavelength \(\lambda \).

Table 1 Critical waveheight for the cnoidal solution of the KdV-equation, calculated for various values of the elliptic parameter m. Corresponding values of wavelength and the dimensionless parameters \(\alpha \) and \(\beta \), in addition to the Stokes number \(\mathcal {S}\) are also listed
Full size table
Fig. 4
figure 4

Left: Maximum allowable wave height for the cnoidal solution as a function of m. Right: Maximum allowable wave height for the cnoidal solution as a function of wavelength

Full size image

4 Numerical study of wave breaking in undular bores

A bore or traveling hydraulic jump is the transition between two uniform flows with different depths. Assuming that one of these flow depths is the undisturbed water level \(h_0\), and the incident water level is given by \(h_0 + a_0\), we denote by \(\alpha =a_0 / h_0\) the strength of the bore. With the nondimensionalization explained in the introduction, the undisturbed depth is 1, and the incident water level is \(\alpha \). In this section we aim to find the threshold for which the bore transitions from being purely undular to one which exhibits breaking in the leading wave. This threshold is defined in terms of the ratio \(\alpha \), and was found experimentally to be \(\alpha = 0.28\) [21].

As already mentioned, if the kinematic breaking criterion is used in connection with phase-resolving numerical models, there is no difficulty in evaluating the local wave and particle velocities. When coupled with the Boussinesq ansatz for the horizontal velocity, the convective criterion allowed for a somewhat accurate description of the breaking of the leading wave behind an undular bore [7]. In this study, numerical integration of a Boussinesq system for undular bores indicated that breaking should occur if the bore strength is greater than 0.379. While this result is qualitatively in the right direction, the critical bore strength is overestimated by about \(30\%\). In the following, it will be shown how this result can be improved.

In order to find the critical bore strength in the numerical KdV model, we will therefore use the value of \(\alpha =0.25\) as a starting point, and then proceed with small increments until the breaking criterion is reached, i.e. when the horizontal particle speed of the fluid exceeds the phase speed. The initial position of the bore front will be set to be the origin, and the bore will then travel downstream in the positive x-direction. Initial data are given by

$$\begin{aligned} \eta _0 (x) = { \textstyle \frac{1}{2}} a_0 \big [ 1 - \tanh (kx) \big ], \end{aligned}$$
(4.1)

where k is a model parameter denoting the steepness of the initial bore slope. The numerical approach to obtaining approximate solutions of (1.1) with the appropriate boundary conditions is detailed in the appendix.

In order to test the convective breaking criterion (1.6), at each time step the horizontal particle velocity and the phase speed of the bore front have to be calculated numerically. First, given the numerical solution \(\eta (x_j,t_n) = \eta _j^n\), we approximate the second spatial derivative of the surface, by the second-order central-difference formula

$$\begin{aligned} \eta _{xx}(x_j,t) \approx \frac{\eta ^n_{j-1} - 2\eta ^n_j + \eta ^n_{j+1}}{\delta x^2}, \end{aligned}$$

where \(\delta x\) is the spatial grid size. We then define the approximate horizontal particle velocity at a grid point \(x_j\) and time \(t_n\), and evaluated at the free surface as \(u^n_j = u(x_j,\eta _j^n,t_n)\), such that Eq. (2.2) yields

$$\begin{aligned} u^n_j = \eta _j^n - { \textstyle \frac{1}{4}} (\eta _j^n)^2 + \Big ( { \textstyle \frac{1}{3}} + { \textstyle \frac{(1+\eta _j^n)^2}{2}} \Big ) { \textstyle \frac{\eta _{j-1}^n - 2\eta _j^n + \eta _{j+1}^n}{\delta x^2}}. \end{aligned}$$
(4.2)

The phase speed of the first wave behind the bore front can be approximated numerically. If the location of the maximum of the wave, say \(X_n\) is found numerically and recorded at each time \(t_n\), the velocity \(c_n\) of the leading wave can be computed approximated using a finite-difference approximation of the first derivative: \(c_n \sim (X_n - X_{n-1})/\delta t\). The left panel of Fig. 5 shows the initial profile with \(a_0 = 0.25\) and \(k = 1\) as a dashed curve, and the resulting bore profile at \(t=6\) as a solid curve. The figure also shows the horizontal component of the velocity field underneath the first wavecrest, and it is apparent that this bore does not break as the particle velocity remains smaller than the propagation velocity of the leading wave.

After increasing \(\alpha \) by increments of 0.001 breaking is obtained for \(\alpha _{crit} = 0.353\). A computation for \(\alpha = 0.5\) (which is already past breaking) is shown in the right panel of Fig. 5. It can be seen that the horizontal particle velocity is greater than the front velocity after the solution has developed for some time. As expected, the critical bore strength is larger than the one found experimentally, but is fairly close to the one found in [7].

Fig. 5
figure 5

The surface profiles and the horizontal velocities below the wave-crests at nondimensional time \(t=6\) for two sets of initial data indicated with dashed-dotted curves. Left panel: Phase velocity \(c=1.1\) is greater than u at the surface; no breaking. Right panel: \(c=1.35\) is smaller than u at the surface: breaking

Full size image

The maximum waveheight \(H_{max}\) of the solution may be extracted at the particular time step where breaking first occurs. If this waveheight is compared to the critical waveheight of the solitary wave from Sect. 2, it is found that the leading wave behind the bore breaks at a waveheight slightly larger than the solitary wave. However, the values obtained for maximum wave height of the bore seem to approach the maximum height of the solitary wave as the bore strength decreases. Since smaller bore strengths will have larger lead times before the largest wave reaches the point of breaking, in these cases, the leading wave has already separated as an individual solitary wave, and the waveheight at breaking matches closely the value obtained for \(H_{\text {max solitary}}\).

The figure below shows a plot of the calculated values of \(H_{\text {max}}\) as a function of time. A curve of the form \(y_{\text {fit}} = At^B + C\) has also been fitted to the data.

Fig. 6
figure 6

Left panel: Maximum wave height as a function of time. Right panel: Bore strength as a function of corresponding breaking times

Full size image

The parameters A, B and C in the fitted curve are given by

$$\begin{aligned} A = 0.1764, \ B = -0.8582 \ \text { and } \ C = 0.6863, \end{aligned}$$

and the residual sum of squares is \(= 2.54\cdot 10^{-5}\). The curve fits the data fairly well, and it is evident from the form of the equation that it approaches the value \(C = 0.6863\) for large t. This value is very close to the maximum possible height of the solitary wave \(H_{\text {max solitary}}=0.6879\). We can therefore confirm that the bore eventually disintegrates into several solitary waves, but as the bore strength, \(\alpha \), increases, the bore will break before this happens and thus cannot be described in terms of solitary waves. This observation also points to the possibility of obtaining a lower value of the critical bore strength \(\alpha _{\text {crit}}\) if a similar curve fit is used for this parameter. To this end, we plot the bore strength as a function of corresponding breaking times. If a curve of the same form as above is used, the parameters A, B and C in the fitted curve are given by

$$\begin{aligned} A = 0.5751, \ B = -0.6908 \ \text { and } \ C = 0.3283. \end{aligned}$$

The residual is \(= 7.74\cdot 10^{-5}\), and visual inspection of Fig. 6 indicates that the fitted curve fits the data rather well also in this case. The numerical value of C indicates an asymptotic critical bore strength of \(= 0.3283\) for large values of t, and this value is somewhat lower than the value of \(\alpha _{\text {crit}}\) obtained previously.

In order to assess whether the values found are reasonable, the results should be compared to experimental findings. One set of controlled experiments was run by H. Favre, and is described in his monograph [21]. In these experiments, it was found that the critical bore strength was close to \(\alpha = a_0 / h_0 = 0.28\). Thus given a certain undisturbed depth \(h_0\), if the bore strength was chosen in such a way that the ratio \(a_0\) to \(h_0\), we less than 0.28, the result was a purely undular bore. If the bore strength was chosen larger than that, the leading wave behind the bore front would break, leading to a much lower waveheight behind the bore front. It was also found that much larger bore strengths led to turbulent bores, but since the KdV equation is only valid for laminar flows, this regime is beyond the scope of the present investigation.

5 Conclusion

It has been shown that the derivation of the KdV equation as a model for surface water waves admits the definition of a local particle velocity field. The horizontal component of this velocity field may be evaluated near the crest of a wave. For solitary and cnoidal waves, evaluating the local Froude number leads to a convective wave breaking criterion which shows that there is a limiting waveheight beyond which the waves feature incipient wave breaking. Thus even though the solitary wave and the cnoidal waves are solutions of the KdV equation for any waveheight, they do not represent valid approximate surface waves for waveheights beyond the critical waveheight.

The critical waveheight for solitary waves is \(H_{\text {max solitary}} = 0.6879\) is quite large, and it is not surprising that these waves do not represent reasonable approximations to real surface waves. For the cnoidal waves, the situation is more complex. For waves near Stokes number 1, the limiting waveheight according to the convective breaking criterion is in the range 0.15 to 0.3 (see Table 1), and these waves are in the range where it is generally thought that solutions of the KdV equation give a faithful approximation of real surface waves. As the elliptic parameter m approaches 1, the critical waveheight of the cnoidal waves approaches \(H_{\text {max solitary}} = 0.6879\). On the other hand, the limit in the case when \(m \rightarrow 0\) is the linear case in which solutions have negligible amplitude, and the critical waveheight in this case is indeed 0.

Fig. 7
figure 7

The critical amplitude \(a_0\) dividing the purely undular case from the case where the leading wave features wave breaking graphed against the undisturbed depth \(h_0\). The slope of the curves is the non-dimensional critical bore strength. The value \(\alpha = 0.28\) was found in Favre’s experimental work [21]. The value \(\alpha = 0.379\) was found in previous work [7]. The value \(\alpha = 0.3283\) is found in the current work

Full size image

The convective breaking criterion has also been used to determine the critical bore strength \(\alpha _{\text {crit}}\) for which an undular bore first exhibits breaking. As can be seen in Fig. 7, the result \(\alpha = 0.323\) is still somewhat higher than the experimentally determined value \(\alpha = 0.28\), which was obtained using wave tank experiments [21], but lower than the ratio obtained using a similar numerical model in [7]. The discrepancy with the experimental findings could be ascribed to the assumptions made on the physical system, i.e. the absence of transverse motion of the fluid and an irrotational and inviscid flow. While recent studies have highlighted the importance of these effects [15, 35], it is not clear that the inclusion of any of these will improve the quantitative agreement between the current analysis and the experimental findings of [21]. For example frictional damping due to molecular viscosity and boundary effects are more likely to delay the onset of breaking than to facilitate it. It therefore seems likely that the limiting procedure used in obtaining the model equations is responsible for some of the discrepancy, and a possible avenue for improving the quantitative description of wave breaking in an undular bore would be to include some of the higher-order terms which were neglected in the KdV model. Such higher-order models including both nonlinear and dispersive effects are well known (see [24] and the references therein).

Another possible improvement is the relaxation of the breaking criterion recently advocated in [4, 5]. In these works, several breaking criteria were used in connection with numerical Boussinesq models, and with the aim of quantifying the onset of breaking in shoaling waves. It was found that the convective criterion (1.6) did not perform as well as other criteria based on local waveheight [26, 42], and a relaxation constant was introduced in order to improve the analysis. However, since this constant has to be determined empirically from a given set of experiments, we have chosen not to pursue such a development here.