Wave Breaking in Undular Bores with Shear Flows

It is well known that weak hydraulic jumps and bores develop a growing number of surface oscillations behind the bore front. Defining the bore strength as the ratio of the head of the undular bore to the undisturbed depth, it was found in the classic work of Favre (Ondes de Translation. Dunod, Paris, 1935) that the regime of laminar flow is demarcated from the regime of partially turbulent flows by a sharply defined value 0.281. This critical bore strength is characterized by the eventual breaking of the leading wave of the bore front. Compared to the flow depth in the wave flume, the waves developing behind the bore front are long and of small amplitude, and it can be shown that the situation can be described approximately using the well known Kortweg–de Vries equation. In the present contribution, it is shown that if a shear flow is incorporated into the KdV equation, and a kinematic breaking criterion is used to test whether the waves are spilling, then the critical bore strength can be found theoretically within an error of less than ten percent.

Saint-Venant equations is described. Favre's experimental results have been examined theoretically from a number of angles. For example, the initial formation of the free-surface oscillations was explained from a physical perspective in [37,43], the interaction of bores was considered in [9], and the breaking of the leading oscillations in the bore was considered in [10,14,22,30]. Recent laboratory experiments revealing intricate properties of undular bores are detailed in [34,35], and a numerical study of the flow structure is given in [8].
Of particular interest in the present article is the critical bore strength dividing the regime of laminar flow from the regime of partially turbulent flow. To explain the situation, assume without loss of generality that the downstream flow depth is the undisturbed depth in the wave flume, say h 0 , and the incident depth is a 0 + h 0 . Defining the bore strength by the ratio a 0 /h 0 , it was found in [21] that there are three main bore types. If the bore strength is below 0.281, the flow is laminar, and since in this case, none of the waves are breaking, this case is termed the purely undular bore. If the ratio α = a 0 /h 0 exceeds 0.281, then the leading wave behind the transition front starts to break, and while the flow still features oscillations, there is some turbulence associated with the breaking waves. If the ratio exceeds 0.75, a fully turbulent bore appears.
The main purpose of the current paper is to demonstrate that the critical ratio found by Favre [21] can also be predicted using fairly simple nonlinear model equations such as the KdV equation in connection with a kinematic (or sometimes called convective) breaking criterion which defines the onset of breaking as the point when the horizontal component of the particle velocity exceeds the crest velocity. In effect, if we let the fluid particle velocity at the leading wavecrest of the bore be U = U (x, η(x, t), t), and the local phase velocity of the wavecrest be C, then the wave starts spilling when U /C > 1. The kinematic wave breaking criterion is one of the simplest diagnostics for predicting the onset of wave breaking (see [28,49] and the references therein), and has been shown to work well in a number of situations [25,27,29].
To study an undular bore in the context of the KdV equation, one needs to be able to pose a boundary-value problem where the incident-free surface level is imposed at, say, the left end of the domain, and the undisturbed level is imposed at the right end of the domain. 1 (see Fig. 1, left panel). Such a model has been developed for instance in [14,42]. One can then evaluate the free surface numerically, and reconstruct the velocity field in the fluid column using the traditional asymptotic expansion of the velocity potential as an asymptotic Taylor series in powers of the vertical coordinate [51]. If the wave crest velocity is also evaluated numerically, then the convective criterion can be tested as an indication of whether the wave starts breaking or not.
Previous studies using this wave-breaking criterion in connection with a Boussinesq system and the KdV equation gave good qualitative results, but were not quantitatively convincing. In particular, the critical ratio was found to be a 0 /h 0 ∼ 0.379 in [10] using a Boussinesq system, and a 0 /h 0 ∼ 0.353 in the context of the KdV equation [14].
A possible improvement on these results may be obtained from the inclusion of vorticity into the description. Indeed, it is by now well known that vorticity can have a   [21]. The value α = 0.379 was found in a previous work [10], the value α = 0.353 was found in [14], and the value α = 0.307 is found in the present work significant effect on the properties of surface waves (see for example [1,30,38,41,48] and references therein). One simple configuration is the case of a linear shear flow such as used in [48]. In particular in the case of long waves, this configuration is expected to capture many features of more general flows (see [17]), and it is our aim in the present work to determine whether the agreement with the experimental results of Favre may be improved by incorporating a constant shear flow into the governing equation. Indeed, the existence of vorticity in a similar situation was found in [26], and a mathematical inquiry into Favre's results also suggested that vorticity might be present in such flows [30]. To get an idea of the strength of vorticity in Favre's experiments, we derive a KdV equation in the presence of a linear shear flow. We then run a numerical simulation of an undular bore and try to match the wavelength of the oscillations of the numerical approximation of the KdV equation with the experimental wavelengths reported by Favre. This procedure leads to an estimate for the vorticity which is then used in simulations aimed at finding the critical bore ratio by testing the leading wave for incipient breaking. The approach outlined above yields a critical bore strength a 0 /h 0 ∼ 0.307 which is within 9% of the experimental ratio of 0.281.
The plan of the paper is as follows. In the next section, the experiments of Favre are explained in some detail. Then in Sect. 3, the KdV equation with a shear flow is explained. Section 4 contains some comments about the numerical approach, and Sect. 5 describes the results of our numerical simulations. Finally, Sect. 6 contains a brief discussion which puts our findings into context with respect to some recent studies on breaking waves.  Some experiments were performed with a horizontal tank and some were performed with a slightly inclined tank. Herein, we consider only experiments in the horizontal wave tank. The tank was supplied by a water tower of constant water level through a pipe C 1 . The water tower was supplied by fixed reservoirs located in the laboratory attic through a pipe C 2 . A third pipe C 3 was used to drain off the excess of water of the tower. The water level of the tower remained constant as long as the flow from C 2 was larger than that of C 1 . The excess water was drained off via C 3 to a reservoir located in the laboratory basement, and then returned by pumps to the reservoirs located in the attic. The water flow from C 1 was adjusted by a valve operated by a servo-motor. The pipe C 2 included a device that allowed the determination of the flow from the water tower. At the end of the tank was a sluice gate that allowed control of the water discharge and full closing of the end of the tank. Beyond the sluice gate, the water was discharged in the same reservoir in the basement as mentioned above. With this setup, it was possible to tune the system so as to guarantee a constant inflow into the wave tank and adjust the inflow to create bores of varying strength. A schematic of the setup just described is shown in Fig. 2.

The Experiments of Favre
Three measuring and recording devices were used: (1) six vertical scales located in the longitudinal axial plan of the tank and distributed along the tank every 12 m were used for measuring the water level at rest or in motion. The accuracy was on the order of one to two tenths of a millimeter. (2) To record the fluctuations of the free surface at the six locations Favre used pressure head (Pitot) tubes whose meniscus positions were recorded on sensitive paper exposed to light. An optical apparatus consisting of a prism, a lens, an electric lamp and a mirror was used to record the meniscus position. Preliminary experiments were conducted to calibrate the devices. Undular bores were generated by a sudden variation of the flow at the upstream end of the tank. Favre performed a series of 30 experiments on undular bores on still water, grouped in three series corresponding to three different initial constant water depths (0.205 m, 0.155 m and 0.1075 m). A spillway plate was placed at the downstream end of the tank whose height was fixed to the desired depth, then the tank was filled through C1 until the water overflowed the spillway plate. Water filling was stopped and 20 min later the water level was stabilized at the height of the spillway plate.
The water flow was determined as a function of the valve displacement (see Figure  38 of Favre [21]). The water flowing into the tank generated an undular bore with a front consisting of a series of undulations whose number increased with the distance of propagation. To make the free surface easily visible, Favre sprinkled aluminum sawdust on the free surface and scattered confetti soaked with black ink. These measures made the free surface clearly visible on photographs. The fronts of the undular bores were photographed when the crest of the first undulation was located at 64.78 m and 64.60 m from the upstream end of the tank for experiments no. 22 and no. 23, respectively. Recall that Favre introduced the bore strength parameter a 0 /h 0 , where a 0 was the mean height of the head of the undular bore (see Figure 41 of Favre [21]) and h 0 was the initial water depth where the first crest of the undular bore was photographed. Bore strength values corresponding to experiments no. 22 and no. 23 are given in Sect. 5. Favre found that for weak values of bore strength the bore undulations were almost sinusoidal, whereas for high values they became cnoidal (experiments nos. 22,23,24).
In experiments no. 22 and no. 23, no wave breaking was observed. In experiment no. 24, the leading wave exhibited spilling breaking after traveling a considerable distance. Experiments no. 26 and no. 29 featured wave breaking after the front of the bore had traveled a shorter distance. Breaking and non-breaking cases can be clearly identified by observing the maximum height of the leading wave. Experimental results based on Favre's data shown in Figure 49 of [21] are plotted in Fig. 3. It is apparent that the maximum height of the leading wave increases with increasing bore strength, up to a maximum of 2.06 times the incident depth a 0 . That maximum occurs in experiment no. 24 which corresponds to a 0 /h 0 = 0.281. As shown in these figures, experiments with higher bore strength feature earlier breaking and much lower wave heights for the leading wave.
As already stated above, the main purpose of the present work is to explore whether the critical bore strength can be found using fairly simple wave models such as the KdV equation. In the next section, the KdV equation will be derived in the presence of a shear flow, and then numerical simulations will be presented which suggest that the critical ratio can be found to within 9% error.

The KdV Equation in the Presence of Shear Flows
Consider a fluid contained in a long channel of unit width and depth h 0 . 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. We 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 the wave motion transverse to the x-axis can be neglected.
With the assumption of irrotational flow and using the incompressibility of the fluid, the two-dimensional Euler equations can be written in terms of the Laplace equation for a velocity potential φ, and the boundary conditions at the free surface are given in terms of φ and the surface excursion η(x, t) by where g is the gravitational acceleration.
As is well known, this problem is difficult to treat both mathematically and numerically, and in practical situations, an asymptotic approximation of the Euler equations is often used. In the case at hand, we have long waves of small to moderate amplitude on a shallow fluid, and it appears that theses waves fall approximately into the Boussinesq scaling regime. Moreover, the waves travel only in one direction down the length of the wave flume, so that an appropriate asymptotic model to be used is the KdV equation, given in dimensional form by Here, c 0 = √ gh 0 is the limiting long-wave speed, whereas already mentioned, where g is the gravitational acceleration and h 0 is the undisturbed depth.
The KdV equation features both nonlinearity and dispersion, and the balance of these two effects gives rise to both solitary-wave solutions and periodic traveling waves [31]. The equation is known to give a good description of the evolution of unidirectional surface water waves in the case when the waves are long compared to the undisturbed depth h 0 of the fluid, the average amplitude of the waves is small when compared to h 0 and transverse effects are negligible [12,16,32]. In the derivation of the KdV equation, the potential is written asymptotically as where f (x, t) represents an approximation to the velocity potential evaluated at the bottom. In addition, the assumption of wave propagation in the direction of increasing x values (to the right) leads to the relation Using (1) and (2) shows that the horizontal velocity can be expressed to second order as That means if a solution η of the KdV equation is given either in closed form or numerically, the horizontal particle velocity can always be found at any point in the fluid column. For more details, the reader may consult [2,51]. In the case at hand, the velocity is evaluated at the free surface so that the breaking criterion can be tested. Using this approach, the authors of [14] found that the critical bore strength was α = 0.353, improving earlier work where a Boussinesq model was used [10]. In the present work, we will incorporate a background shear flow to further improve the comparison.
In the presence of background vorticity, the horizontal velocity field may be defined as U (x, z, t) = u(x, z, t) − (z + h 0 ) (see Fig. 4). Model equations of KdV and Boussinesq type for surface waves in the presence of background shear can be derived in a similar fashion as in the irrotational case. The reader may consult for example [15,40,46,54]. Except for a minor modification, the equation to be used in the present study was derived in [40]. For background shear such as shown in Fig. 4, the appropriate This equation is given in dimensional form, with the non-dimensional parameter C + 0 = −˜ 2 + ˜ 2 4 + 1 quantifying the strength of the background shear, and˜ being the nondimensional vorticity scaled with˜ = g h 0 . The horizontal perturbation velocity field can be found in a similar way as for the irrotational case, and is given by

Non-Dimensionalization and Numerical Discretization
To prepare for the numerical discretization, it is convenient to scale the variables appearing in the equation above as follows: Using this scaling, the non-dimensional KdV equation can be expressed as The coefficients C + 0 and C − 0 are defined in the same way as in the last section. In the scaled variables, the horizontal velocity field appears as It is well known that the KdV equation features exact solitary-wave solutions. If the coefficients of the KdV equation are defined such as in (3), the solitary wave has the form η = η 0 sech 2 with the phase speed c given in terms of the amplitude η 0 by c = C + 0 + . This exact solution will be used later to validate the implementation of the numerical algorithm.
The numerical approximation of the solutions η(x, t) and U (x, η(x), t) is based on a finite-difference method for the spatial derivatives and a hybrid Adam-Bashforth/Crank-Nicolson time integration scheme, such as previously used in [14,42]. The local phase velocity C of the leading wavecrest is computed approximately by following the crest evolution and using a second-order central difference formula.
To define an appropriate numerical discretization, boundary conditions need to be imposed. In the far field upstream and downstream of the bore, the surface elevation η approaches the stipulated values α and 0, respectively. These boundary conditions are exact up to machine precision as long as the spatial domain is large enough. In addition, a Neumann boundary condition needs to be specified in the far field downstream (see for example [13]). Using a homogeneous Neumann boundary condition and the Dirichlet conditions as indicated above yields the initial-boundary-value problem The initial data are given by η 0 (x) = 1 2 a 0 1−tanh(kx) , where k is a model parameter denoting the steepness of the initial bore slope. In an idealized setting, one may take the limit as k approaches infinity, leading to the so-called dispersive shock problem [20,23]. However in the present case, we take the finite value k = 1.
It will be expedient to rewrite the problem to obtain homogeneous boundary conditions, and we define an auxiliary function ζ(x, t) ≡ η(x, t) − η 0 (x). Using ζ , the problem can be reformulated as an inhomogeneous equation with a forcing given in terms of ζ and η 0 , but with homogeneous boundary and initial conditions. The equation for ζ can then be written as and homogeneous boundary and initial conditions are imposed. We discretize the spatial domain [−l, l] using a finite set of points, , where x 0 = −l and x N = l, and δx = 2l/N is the distance between two neighboring grid points. The time domain is also discretized uniformly by defining t n = nδt, where t 0 is equal to zero. With this notation, the approximate function value at time t n and grid point x j is defined to be v n j ≈ ζ(x j , t n ). Regarding the spatial discretization, the first and third derivatives at a point x j are approximated by the central difference formulas and Since we impose the Dirichlet conditions v 0 = 0 and v N = 0, the equation can be solved for the grid points {x j } N −1 j=1 , and that leaves us with just two points for which the third derivative approximation is not valid. Given the Neumann condition and the central difference approximation, we have This enables us to use the third derivative approximation at the grid point x N −1 in the form As there is no Neumann condition on the left boundary, we employ a forward difference formula to approximate the third derivative at the grid point x 1 as follows: These computations were run up to final time T = 1, and with a spatial grid size of δx = 0.01 The difference formulas (5), (6), (7) and (8) applied at all grid points give rise to the discrete differentiation matrices D 1 and D 3 as follows: Applying a Crank-Nicolson method to the linear terms on the right-hand side of the equation, and an Adams-Bashforth method to the nonlinear terms, the difference equation for determining the vector v n+1 is given by F(x 2 ), ..., F(x N −1 )) T . This n × n-system of equations can easily be solved for v n+1 , and then to advance the numerical approximation by one time step, only three multiplications by sparse matrices are required.
The implementation is verified by using the well-known exact solitary-wave solution of the KdV equation which is given in the form (4). In Tables 1 and 2, the equation is solved in the case = −0.2 and with the amplitude η 0 = 0.5 on the time interval t ∈ [0, 1]. It can be clearly seen that the second-order convergence is achieved in the spatial and the time discretization.  , 610] was used, and some experiments were double checked on a larger domain to make sure that there was no detrimental effect of numerical instabilities due to the treatment of the boundary conditions. A comparison of the first five waves is made at a propagation distance of about 600 depths, similar to the analysis in [30]. A plot of the surface elevation is shown in Fig. 5. Tables 3 and  4 show the wave heights and wavelengths of the five waves behind the bore front for several values of the vorticity . The experiments were generally checked by making runs with smaller spatial and temporal grid sizes to make sure that numerical errors did not contribute to the results.
With the numerical measurements of wavelengths and wave heights of the five leading waves in hand, a comparison with the measurements from Favre's experiment can be conducted. For each value of in Tables 3 and 4, an estimate of the relative error 2   600 depths, we could observe whether the waves will break or not. In Fig. 8, the crest speed C and the horizontal velocity component at the free surface are plotted base on a computation with = −0.2213 and a 0 = 0.307. The ratio between the crest speed and the horizontal particle velocity reaches U /C ≥ 1 at about 600 depths, and this means that the leading wave is starting to spill which limits the further growth of the wave. Thus, with the inclusion of a background shear flow with constant vorticity the critical ratio is found to be 0.307.

Discussion
As mentioned in the introduction, the kinematic wave breaking criterion is one of the most commonly used criteria for the prediction of the onset of wave breaking. The criterion has been shown to pinpoint the commencement of wave breaking in various situations. In particular, the criterion was shown to perform well in deep water in [29,49], and in shallow water both on flat bathymetry [25] and on a sloping beach [27,50]. In some situations where the kinematic breaking criterion performs poorly, the problem can be ascribed to the difficulty of accurately finding the phase velocity of the waves from measurements [45], and the directionality of the waves in threedimensional situations [53]. Some authors have attempted to use the kinematic criterion in connection with numerical codes to determine whether to switch from a dispersive Boussinesq-type scheme to a shallow-water scheme with numerical dissipation to capture energy loss in breaking waves. It was suggested in [5] that this approach works best for wave shoaling if the criterion is tightened by introducing a positive constant κ < 1, and to define breaking onset as the first time that U /C > κ. However, this approach requires the determination of the constant κ. If an appropriate value for κ can be found, then this breaking criterion can give excellent results [6,11,39].
Sharpening the kinematic criterion if used as a numerical switch makes sense as the numerical dissipation needs time to have an effect on the waves. In fact, tightening the kinematic criterion was already suggested earlier based on experimental evidence [45], and more recently based on studies of wave breaking in large and intermediate depth [7,18]. In these works, a new parameter B based on crest speed and local energy flux and density was put forward as a diagnostic for the initiation of wave breaking. In particular, as shown in [7], using this parameter reduces to a sharpened convective criterion when evaluated at the free surface. In a nutshell, the sharpened criterion based on the parameter B predicts wave breaking when B ∼ 0.87, though breaking may not commence until B is actually close to 1. As already alluded to, the new criterion based on the parameter B is based on a dynamic or energetic criterion based on evaluation of energy flux and density, such as put forward for example in [44].
To complete this brief review, one may add that geometric breaking criteria are popular in some quarters. Geometric criteria are based on the shape and in particular the steepness of the waves close to breaking, such as reviewed in [4]. Among them the most used is the limiting wave steepness parameter s ≈ ak max , which can be transformed into the kinematic limit u ≈ c. For a detailed review on geometric, kinematic and dynamic criteria for breaking onset one can refer to [3].
To further improve comparison with experiments using the simple kinematic criterion used in the current work, one might use fully nonlinear Boussinesq systems, such as the Serre or Green-Naghdi equations [19,33] or higher-order KdV-type equations, such as the higher-order regularized long wave equation used in [55]. One possible obstacle with this approach is that the boundary-value problem must be shown to be stable to small perturbations.