Wave breaking in undular bores generated by a moving weir

A river bore is an upstream-propagating transition between two different flow depths usually caused by tidal forces. Similar flows can also be realized in controlled environments such as wave flumes, and a number of laboratory studies have been conducted in order to bring to light some of the main features of bores. One of the first in-depth investigations of undular bores in the laboratory was carried out by Favre in a dedicated wave flume,¹ and Favre’s results have been examined theoretically from a number of angles. For example, the initial formation of the free-surface oscillations was under study in Refs. 2–4, and the energy balance at the bore front has been reviewed in Refs. 5–9. Dissipative effects were considered in Refs. 10 and 11, and the breaking of the leading oscillations in the bore was studied in Ref. 12.

Favre’s experimental results highlighted in particular the transition between purely undular and partially turbulent bores in terms of the bore strength defined as the ratio of downstream to upstream flow depth, and the main purpose of the present work is to explore whether the critical bore strength can be found using standard wave models such as the classical Boussinesq system.

Favre himself used the shallow-water Saint-Venant system with frictional terms to describe the flow. One weakness of the shallow-water system is that it fails to describe the transition between undular and breaking bores because in the shallow-water theory, all waves eventually break (see Ref. 13, chap. 13.11, or, Ref. 14 chap. 11.3). On the other hand, the Boussinesq models incorporate dispersion, and small-amplitude waves do not feature the typical steepening observed in the shallow-water theory.¹⁵ Therefore, wave breaking cannot be observed directly in the Boussinesq equations, but one may use a convective criterion comparing the crest speed of the leading wave to the fluid particle speed near the crest as a diagnostic for wave breaking. In these terms, a wave breaking criterion may be defined by saying that a wave starts to break when the horizontal component of the fluid velocity near the wavecrest exceeds the wavecrest velocity.

An initial study using this wave-breaking criterion in the context of Boussinesq modeling of Favre’s experiments gave good qualitative results but was quantitatively inconclusive.¹⁶ The relatively poor agreement with the experiments may be partially due to vorticity created by the discharge used to generate the undular bore. The existence of vorticity in a similar situation was also found in Refs. 17 and 18, and a mathematical inquiry into the Favre results also suggested that vorticity might be present in such flows.¹⁹ 

In the current work, we investigate the wave breaking onset in the leading wave in a bore created by moving a weir through a wave tank. Such experiments were carried out by Wilkinson and Banner,²⁰ and as shown in Fig. 1, this generating mechanism yields bores that are markedly different from the ones created by a forced inflow. Favre found that wave breaking occurred when the bore strength was larger than 1.28. By contrast, in the experiments carried out by Wilkinson and Banner, the critical bore strength was found to be about 1.38. As will be shown in the body of the present paper, if the convective wave-breaking criterion is used in connection with a Boussinesq model in the case of a bore generated by a moving weir, then the critical bore strength may be found numerically to within an error of less than 2%. The plan of the paper is as follows. In Sec. II, the Boussinesq system used in the study of the undular bore is introduced. In Sec. III, the wave breaking criterion is explained in detail. The numerical experiments are described in Sec. IV, and the results are discussed in the Conclusion.

In the experiments by Wilkinson and Banner, a moving weir is used to generate an undular bore (see Fig. 1 for an illustration of the geometric setup). The weir is placed at the bottom of a wave flume and dragged along the tank. As the weir is moving, it creates a mound of fluid propagating ahead of the weir itself. This phenomenon is well known to experimentalists and is sometimes called upstream influence.²³ The borefront develops oscillations (which in this case are often called undulations), and at first sight, the resulting undular bore seems similar to the bore featuring in the work of Favre.¹ However, as already mentioned, there are some important differences between these two flows. In particular, vorticity seems to be an important factor in the bores generated by a constant discharge, while as far as we can tell, vorticity is a minor effect in the bore generated by the moving weir. Some authors have pointed out the importance of bottom friction on the appearance of an undular bore.^10,7 Indeed, in a river bore, if the conditions are right, a nearly steady profile of undulations can be observed. However, in the present case, the main focus is on the onset of wave breaking which happens on relatively short time scales where dissipative effects do not have a major impact.

A bore may be characterized either by the bore strength already mentioned or using the Froude number. The Froude number is defined by $Fr=U/gh0$⁠, where U is the velocity of the bore front, h₀ is the undisturbed (upstream) depth, and g is the gravitational acceleration.^24–26 While the Froude number seems to be more common in the literature, in the present case, since we are aiming for a favorable comparison with the experimental work in Ref. 20, it is expedient to use the same parameter as in that work to quantify the strength of the bore. Given the average downstream flow depth a₀ + h₀, the bore strength is defined as the ratio of flow depths by $r=a0+h0h0$⁠.

As a side note, we should mention that recent work suggests that a single parameter such as the Froude number or bore strength may not be sufficient to classify bores and hydraulic jumps. For example, the authors of Ref. 27 proposed an additional parameter related to the vorticity of the flow. However, in the present case, such an extension is unnecessary.

Boussinesq systems are approximations of the full Euler equations valid for long waves of small amplitude if the fluid is incompressible and inviscid and the flow is irrotational.^28,15 Due to their relative simplicity, Boussinesq models are frequently used to model various surface wave phenomena. Since the undular bores described in Ref. 20 have small amplitudes and long wavelengths relative to the undisturbed depth, it is natural to use a Boussinesq model in the current situation.

Here, and in the rest of this article, a Cartesian coordinate system (x, z) is chosen such that z = η(x, t) describes the free surface η, the undisturbed free surface is at z = 0, and the bottom of the wave tank is located at z = −h₀. Moreover, u(x, t) is the horizontal fluid velocity at the fixed height $z=7/9−1h0$⁠. We use the Boussinesq system

to simulate an undular bore. The Boussinesq system (1) was introduced in Ref. 21. It is valid in the Boussinesq scaling regime, which applies when the parameters α = a/h₀ and $β=h02/ℓ2$ are small and of the same order of magnitude. Here a is a typical wave amplitude, ℓ is a typical wavelength, and h₀ is the undisturbed depth, and the parameters α and β describe the relative strength of nonlinear and dispersive effects in the flow.

System (1) is a coupled system of two regularized long-wave equations. Using the same rationale as the authors of Refs. 2 and 22 who introduced the Benjamin-Bona-Mahony (BBM) equation, it can be surmised that this system is more amenable to numerical methods than a number of other Boussinesq systems that were discussed in Ref. 21. In fact, the system can be solved efficiently using a Fourier-spectral collocation method.

As shown in Ref. 21, the system has exact solutions which can be used to test the convergence of numerical codes. The system can also be extended to accommodate time-dependent bathymetry. This was first done in Ref. 29, and in the present case, a time-dependent, spatially varying bottom topography h(x, t) is included in these equations to model the weir moving along the bottom of the wave tank. The precise version of (1) to be used, including all parameter-values and some details on the numerical discretization, is given in the  Appendix.

In Sec. III, it is shown how the full fluid velocity field is reconstructed from these variables, using an asymptotic expansion of the velocity potential. It will be shown that once a solution of (1) is available and the wave crest velocity is known, then a local convective criterion can be taken as an indication of whether the wave starts breaking or not. Denoting the crest speed by U and the horizontal component of the fluid particle velocity by u(x, z, t), we say that the wave breaks if u(x, η(x, t), t) > U. This convective wave breaking criterion is standard in the study of breaking waves (see Ref. 30 and references therein) but has not been tested in the current situation.

In order to understand the wavebreaking criterion, a look at the derivation of the Boussinesq system (1) is in order. Since this derivation is well known, we just summarize a few points of importance. For a more detailed discussion, the reader is referred to Refs. 16 and 13. As discussed above, since the bottom of the wave tank is constant away from the weir, and given by z = −h₀ near the leading wave crest behind the bore front, the following discussion is based on the Boussinesq system (1). It is assumed that α and β are small and of the same order so that Eq. (1) are valid and higher order terms in α and β can be neglected in the following derivation.

With the assumptions of irrotational flow of an inviscid and incompressible fluid, a velocity potential ϕ(x, z, t) can be introduced. The surface wave problem is then formulated in terms of the velocity potential as

with the free-surface boundary conditions

As discussed in Ref. 13, since the total depth is small, an ansatz of the form $ϕ(x,z,t)=∑0∞(z+h0)nfn(x,t)$ is suggested for the velocity potential. Inserting this into the Laplace equation and using the boundary condition on the bottom gives

to second order in β. Note that f = f₀ is the velocity potential at the bottom z = −h₀ and that $ϕx(x,z,t)=fx(x,t)−12(z+h0)2fxxx(x,t)$ is the horizontal velocity at a height z in the fluid.

Using this relation, it is possible to derive systems of equations such as (1) where the horizontal velocity is modeled at a prescribed depth $z=(7/9−1)h0$⁠. Letting $u(x,t)=ϕx(x,z,t)z=(7/9−1)h0$ be the horizontal velocity at $z=(7/9−1)h0$⁠, we obtain

to order β². Combining (3) and (2) gives

up to $O(β2)$⁠.

Equation (4) can be used to approximate the horizontal fluid velocity field at any point (x, z) in the fluid domain. As discussed in Ref. 16, we assume that wave breaking commences when the horizontal fluid particle velocity at the leading wave crest exceeds the local phase velocity of the wave crest, denoted by U. Evaluating (4) at z = η(x, t) yields the following breaking criterion:

A wave solution $η(x,t),u(x,t)$ of system (1) starts breaking if

This criterion is illustrated in Fig. 2. From (5), it follows that once a solution $η(x,t),u(x,t)$ of the BBM system (A2) is available, the onset of wave breaking can be predicted provided the wave propagation speed U is known.

We consider quasi-twodimensional fluid flow in a numerical wave tank of length L. The coordinate system (x, z) is as defined in Sec. II, with the positive z direction taken vertically upwards and the horizontal x-axis along the undisturbed free surface of the fluid. The bottom of the wave tank is located at z = −h₀, and the undisturbed free surface is located at z = 0, as shown in Fig. 3. The numerical domain is x ∈ [0, L], and the fluid domain is then $0,L×−h(x,t),η(x,t)$⁠, where z = −h(x, t) describes the weir and the bottom of the tank, and z = η(x, t) is the free surface as usual.

In the numerical simulation used in the present study, system (A2) is approximated using a spectral collocation method coupled with a 4-stage Runge-Kutta time integration scheme. At each time step, the resulting linear system is solved using a preconditioned conjugate gradient method. As the bore is generated by a weir moving along the bottom of the tank, we use initial values of zero surface displacement and zero initial velocity, i.e., η(x, 0) = 0 and u(x, 0) = 0.

Wilkinson and Banner conducted numerous experiments with varying undisturbed fluid depths and showed experimentally that undular bores may persist up to bore strength r = 1.377. Since the goal of this paper is to obtain favorable results in comparison with Wilkinson and Banner,²⁰ we used the same parameters as in that work. Accordingly, a weir of height 25 mm and length 150 mm was used in our simulation. As the value of h₀ is not consistently specified in Ref. 20 for a given Froude number, we have chosen to use the fluid depth h₀ = 0.06 m in all simulations shown, except the computations shown in Fig. 4. In this case, the undisturbed depth was specified in Ref. 20 to be h₀ = 0.053 m. We chose to run the simulations for t = 15 s as larger time intervals seemed to have little impact on the eventual bore strength r. The parameters used in our simulation are summarized in Table I. The function

with the parameters A = 0.025 m, κ = 0.2, m = 2, $w$ = 1.22, σ = 100, and h₀ = 0.06 m is used to model the weir moving along the bottom of the tank. The speed $Uw$ at which the weir is moving along the bottom of the tank is the only parameter which is varied, while the rest of the parameters in (6) are held fixed.

In order to test the validity of the numerical simulations, we first run a test with a purely undular bore. We choose the case shown in Fig. 3, page 372 from,²⁰ which features a Froude number Fr = 1.24, and is far from the breaking regime. In this case only, in order to obtain a good match with the experimental data, we have used the undisturbed depth h₀ = 0.053 m. Figure 3 from Ref. 20 is shown in Fig. 4 in digitalized form, alongside the numerical approximation of η obtained in our simulation. As can be seen, there is relatively good correspondence between the experimental data and numerical simulation. In particular, the wavelengths and amplitudes of the bores agree very well.

Since the bottom is described by the constant function z = −h₀ near the leading wavecrest, Eq. (4) can be applied in approximating the particle velocity there. At each time step tⁿ = nΔt in the simulation, the position xⁿ of the leading wave crest is located. From this sequence of values, the local phase velocity Uⁿ of the wavecrest is estimated by a fourth-order, backwards finite-difference formula. The fluid particle velocity at the leading wavecrest behind the bore front is estimated from (4) with z = η(xⁿ, tⁿ) at each time step.

As $Uw$ is increased manually in each new run, the wave crest velocity and the horizontal fluid particle velocity also increase. Eventually the fluid particle velocity exceeds the phase velocity of the wave, and wave breaking commences. The process is illustrated in Fig. 5. In the left panel, the value of r is 1.38 and no intersection is taking place at times t ≤ 15 s. As shown in the right panel of Fig. 5, the particle velocity curve first intersects the phase velocity curve at t = 15 s for the value r = 1.40 which is therefore just above the critical bore strength for the current setup. The corresponding weir velocity is $Uw$ = 0.724 m/s.

Table II tabulates whether the bore is breaking or not for different values of the bore strength r obtained numerically. Each row in Table II contains values of r and Fr_U obtained from a fixed value of $Uw$⁠. We find the critical value of the bore strength in Table II to be r = 1.395, between the values 1.39 which does not feature breaking, and 1.40 which does lead to breaking. As an additional comparison, we also include the Froude numbers $FrU=U/gh0$⁠, where U is the phase speed of the leading wave behind the borefront, approximated numerically at t = 15 s by the method described above.

Comparing the results in Table II with the ones in Ref. 20, we see that the numerical critical bore strength r = 1.395 is about 1.3% higher than the experimental value r = 1.377 in Ref. 20 and that Fr_U = 1.358 is about 4.5% larger than Fr = 1.3 in Ref. 20, so the results are relatively close in both cases. The data are also represented in Fig. 6, where the left panel shows data from the experiments in Ref. 20, and the right panel shows data based on our numerical simulations.

In this work, it has been shown that the transition from laminar to partially turbulent flow in an undular bore can be captured with a standard weakly nonlinear dispersive system of Boussinesq type. Reference is made to two experimental studies: the experiments of Favre¹ and the experiments of Wilkinson and Banner.²⁰ The particular quantity of interest is the bore strength r given by the ratio of downstream to upstream flow depths.

In previous work, it was shown that using a convective wave-breaking criterion leads to qualitative agreement with the experiments of Favre in the sense that the Boussinesq system features a critical ratio indicating the commencement of wave breaking.^16,31 However, the quantitative agreement between the numerical experiments in Refs. 16 and 31 was not very good. The discrepancy may possibly be ascribed to the appearance of vorticity in the undular bore, such as suggested in Refs. 17 and 19. Indeed, it was argued in Ref. 27 that a single parameter such as the bore strength r or the Froude number Fr may not be sufficient to classify bores, even in the case where the flow is quasi-two-dimensional. In the present case of a bore generated by a moving weir, it was found that the bore strength r can be used as an effective diagnostic parameter to predict whether a bore will break or not. The critical bore strength signaling the demarcation between laminar and turbulent flow was found experimentally to be about 1.377.²⁰ Our numerical simulations indicate that the leading wave breaks at a bore strength of 1.395. Thus, experiment and simulation coincide up to an error of less than 2%.

The separate comparison using the Froude number $FrU=U/gh0$ for the phase speed of the leading wave crest behind the bore front, which was found to be Fr_U = 1.358 when the bore started breaking, is also fairly close to the critical Froude number Fr = 1.3 in Ref. 20. Using this parameter as an indicator for the onset of breaking yields agreement between experiment and numerical simulation with an error of less than 5%. Additionally, Fig. 4 shows that the shape of the free surface of the bore obtained numerically for non-breaking bores is very similar to the free surface in Ref. 20. This further validates the results obtained from the model used in this work.

One may ask whether the bore strength r or the Froude number Fr is a better indicator for whether a bore features incipient breaking. While the Froude number is probably more commonly used than the bore strength, the present study found slightly better agreement between theory and experiment if the bore strength was used as a parameter. Indeed it can also be observed in the left panel of Fig. 6 that the relation between Fr and r is slightly nonlinear in the experimental measurements.

To put our findings into context, it should be mentioned that wave breaking in an undular bore is a special case of wave breaking on a planar beach, a phenomenon which has been studied widely. It is not immediately clear whether the convective breaking criterion can be used in this case as well. Recent work paints a complicated picture. The experiments described in Ref. 32 point to the convective criterion as being a good indicator for the commencement of wave breaking. On the other hand, if the wave breaking criterion is to be used in a numerical model in order to switch between a dispersive system and a hyperbolic system (in order to exploit the dissipative nature of numerical approximations to hyperbolic systems), then the convective breaking criterion should possibly be sharpened in order to give the waves time to adjust (see Refs. 36–38 and the references therein).

Indeed, the detection of wave breaking through various wave breaking criteria has been researched by a number of groups for many years. As explained in Ref. 33, there are essentially three classes of criteria. Geometric criteria are based on the shape and in particular the steepness of the waves close to breaking, while kinematic criteria such as used here are based on violation of the kinematic free surface boundary condition. As it generally seems to be understood that no dynamic insight or advance warning of imminent breaking is provided by geometric and kinematic criteria,³⁴ dynamic criteria based on evaluation of the energy flux such as proposed in Refs. 35 and 30 seem to be favored at the moment.

Recently a new parameter based on crest speed and local energy flux and density has been put forward as a diagnostic for the initiation of wave breaking. As shown in Ref. 30, using this parameter reduces to a sharpened convective criterion when evaluated at the free surface. However, it has been tested mainly in deep and intermediate water,^30,39 and it is not clear at this point whether this diagnostic will work in shallow water such as in the current situation.

It should also be noted that there have been extensive efforts to understand various aspects of wave breaking using numerical approximations of the full Euler equations. In particular, in the case of potential flow, we mention the work on overturning breakers reported in Refs. 40 and 41. On the other hand, some authors have advocated the use of conservative Boussinesq system⁴² or indeed for the use of fully nonlinear long-wave systems, such as the Serre or Green-Naghdi equations^43,44 in order to study wave breaking on a slope (see Ref. 45 for example). However, in the present case, no such extension is needed.

This research was supported by the Research Council of Norway under Grant No. 239033/F20. The numerical code used in this work was developed jointly by A. Senthilkumar, H. Kalisch, and H. Munthe-Kaas.

In order to simulate an undular bore generated by a moving weir numerically, a BBM-BBM type system is used in this work. In Ref. 29, the following system describing unsteady fluid flow in one horizontal dimension over a time dependent, spatially varying bottom topography is derived:

Here u = u(x, t) is the horizontal fluid velocity at the height z = (θ − 1)h₀, where 0 < θ < 1, η(x, t) is the free surface, and h(x, t) describes the bottom topography. The parameter θ is an arbitrary modeling parameter which has been chosen to be equal to $7/9$ in the body of this paper. The definitions of the further parameters used in the system are as follows:

where $μ,ν∈R$⁠. Taking μ = ν = 0 in (A1) gives

where $b=12θ2−13,d=121−θ2$⁠. The BBM-BBM type system (A2) with the boundary and initial conditions stated in Sec. IV is used to simulate an undular bore numerically in this work. Like system (1), this system is valid under the Boussinesq approximation. As discussed earlier, the time dependent bottom h(x, t) is included in Eq. (A2) to model the weir moving along the bottom of the wave tank.

The numerical approximation of this system is effected by a pseudo-spectral collocation method based on a Fourier basis. This discretization entails the imposition of periodic boundary conditions. In the present one-dimensional case, this choice does not pose a serious problem since the numerical domain [0, L] can be made large enough to prevent wave interactions at the far ends on either side. Rescaling the equations to the numerical interval [0, 2π] introduces a scaling factor λ = L/(2π). Using the standard collocation points $xj=2πjN$ for j = 0, 1, …, N − 1 yields the first-order collocation derivative matrix D_N and the corresponding second-order matrix D_N,2, such as detailed in Refs. 46 and 47. Defining the N × N identity matrix by I_N, and the discrete unknowns as Σ^N and U^N, the semi-discrete system is written as

where h_N is defined by $[hN(t)]j=h(xj,t)$⁠. This system of N ordinary differential equations is integrated using a 4-stage Runge-Kutta method. The most demanding operation is the inversion of the time-dependent matrices on the left-hand side of these equations which has to be performed at each time step. The inversion is performed using a preconditioned conjugate gradient (PCG) method (see Ref. 48). The scheme has been tested for convergence in Ref. 49. Alternative methods for the numerical discretization of these equations are finite-element methods²⁹ or finite-difference methods.¹⁶ These methods would allow for straightforward incorporation of various boundary conditions, but in the present case, that is not necessary.