Energy budget in a dispersive model for undular bores

Energy conservation properties of weak bores in free-surface flows are considered. The energy loss in the shallow-water theory for an undular bore is thought to be due to upstream oscillations that carry away the energy lost at the front of the bore. Using a higher-order dispersive model equation, this expectation is confirmed through a quantitative study, which shows that there is no energy loss if dispersion is accounted for.


INTRODUCTION
Bores and hydraulic jumps arise in open channel flow as the transition between two uniform streams with different flow depths.River bores are generally tidal phenomena, appearing if the tidal range is large enough and the geometry of the river estuary is favourable.Even though bores are commonly observed in nature, few field measurements are available, and quantitative information is most easily obtained from laboratory experiments.The study of Favre [13] focused on bores propagating into undisturbed water.Taking as a definition of the bore strength the ratio α = a 0 /h 0 of the initial difference in the surface elevation a 0 to the undisturbed water depth h 0 , he found that depending on the value of α, bores fall into one of three categories.For values of α below 0.28, the bore is purely undular.This means that there are oscillations in the downstream part, but none of these trailing waves are breaking.When the ratio α is between 0.28 and 0.75, the bore still features oscillations, but one or a few waves are breaking.For values above 0.75, the transition region is completely turbulent.
In the classical theory of bores of Rayleigh [22], conservation of mass and momentum is used in the shallow-water equations to predict the loss of energy at the bore front.It is now more or less accepted that the energy loss in strong bores is due to turbulent dissipation.Although it has not been confirmed by quantitative evidence, the energy loss in weak bores is believed to be caused by the oscillations behind the bore front.In this article, a higher-order system of dispersive evolution equations and the associated energy integral are used to study the energy balance across an undular bore.The dispersive system to be used in our study has the following form: Here η is the deflection of the free surface from its rest position, and w represents the horizontal flow velocity at a height of 2/3h 0 above the flat bottom.The coefficient g denotes the gravitational acceleration, and h 0 is the undisturbed depth of water.The system (1) was proposed by Bona et al. [6] for small-amplitude long waves with weak transverse variation as a generalization of the regularized long-wave equation [1,20].The purpose of the present work is to use equation (1) to show that if dispersive effects are included in the description of the fluid flow, then the energy loss mentioned in the previous paragraph can be accounted for.Indeed, monitoring the mechanical energy associated with (1) shows that the rate of change of the energy contained in a control volume incorporating both the bore front and the ensuing oscillations, is exactly equal to the net influx of energy to the same volume.Thus one may indeed say that the energy loss in the shallow-water theory is due to surface oscillations which cannot be described by the shallow-water system.
In a previous work, Lemoine tried to match the energy loss due to the shallow-water approximation with the wave energy contained in the trailing oscillations [19], but obtained only moderate agreement.Benjamin and Lighthill [1] attributed the failure of Lemoine's theory to accurately match experimental results to the fact that harmonic waves were used to approximate the oscillations behind the bore front.Using travellingwave solutions of the Korteweg-de Vries (KdV) equation given in terms of cnoidal functions, Benjamin and Lighthill argued that energy is disseminated by trailing oscillations, but that an additional energy loss is required.Indeed, showing that a supercritical uniform stream cannot be matched with a cnoidal periodic wave-train unless there is either a change in the momentum flux or in the energy flux, Benjamin and Lighthill concluded that viscosity should play a role in the required loss of energy.Later, Sturtevant [23] and Byatt-Smith [7] developed a theory of boundary layers, and argued that the part of the energy that is not disseminated through oscillations is absorbed by viscous action in a boundary layer near the bore front.It is clear that in any physical channel flow, there will be a boundary layer, and molecular viscosity will play an important role.On the other hand, the shallow-water theory which originally predicted the energy loss is an approximation to the Euler equations, and not the Navier-Stokes equations.Therefore, it should be possible to explain the energy loss in the shallow-water theory without resorting to the effect of viscosity.As will become clear in the following development, the inclusion of dispersive effects is indeed sufficient to explain the energy loss predicted by the shallow-water theory.In order to get an accurate description, it appears to be vital to use a system of equations such as (1) which allows for counter-propagating waves, and to use the correct form of the mechanical energy.
To set the stage for our investigations, we briefly describe the shallow-water theory and its implications for the study of bores.If the transition between the upstream and downstream flow depths in a bore is gentle, and the undisturbed water depth is not too large, a long wave approximation is justified and the shallow-water system appears as the governing set of equations for the fluid flow.In (2), h(x,t) = h 0 + η(x,t) and both h 0 and η have the same interpretation as in (1).The variable u, however, represents a uniform horizontal velocity.
The system (2) has a well-known weak solution, accommodating the transition between the uniform flows through a discontinuity.A number of conservation equations are associated with the system (2).The most relevant here are the conservation equations for mass, momentum, and energy.For the discontinuous bore solution mentioned above, one may use the integral forms of these conservation equations to obtain the following bore conditions (see [24]).
For a given function f , the notation [ f ] 2 1 connotes the difference in function values f (x 2 ) − f (x 1 ) across the discontinuity.Considering now the discontinuous bore solution of the shallow-water system, if the upstream velocity u 2 = u(x 2 ,t) and the downstream flow depth h 1 = h 0 + a 0 are given, then equations (3) provide the speed of the bore front U and the downstream velocity u 1 .
The energy contained in a control volume of unit width above the interval [x 1 , x 2 ] is given by the energy integral E(h, u) associated with the shallow-water system (2).This integral is given by The energy flow rate across the boundaries above the endpoints x i is given by Note here that the second term in F i includes the work done by the pressure force, and therefore the quantity F 1 − F 2 represents the net influx of energy into the control volume.For the discontinuous bore solution of the shallow-water system (2) one cannot insist on conservation of energy.The loss of energy across the bore may be computed exactly, and is given by In ( 6), the energy loss is represented by the negative right-hand side.
We end this section with a few remarks about the discontinuous bore solution.One may of course use a continuous transition as initial data for the shallow-water system.In this case both mass and momentum, as well as energy are exactly conserved.However, given the nonlinear hyperbolic character of (2), it is clear that all such solutions will eventually break, and discontinuities will develop.Therefore, no generality is lost by considering solutions that are discontinuous from the start.
In the next section we will derive a dispersive system of partial differential equations along with the corresponding energy integral.In Section 3, a finite difference method for solving the dispersive system is put forward.This numerical method is employed in Section 4 to conduct some numerical experiments with focus on the energy budget in a bore.

THE DISPERSIVE SYSTEM AND ITS ASSOCIATED ENERGY
Throughout this article, it is assumed that the fluid is incompressible and inviscid.The flow is assumed irrotational and essentially two-dimensional.It is also assumed that the free surface can be described as a single-valued function η(x,t) and that no wave-breaking takes place.With these provisos the free surface problem can be written in terms of the Laplace equation for a velocity potential φ (x, z,t), accompanied by a set of nonlinear boundary conditions at the bottom and at the free surface.In the nondimensional variables defined below, the problem is as follows: If l represents a typical wavelength and a 0 denotes a typical wave-amplitude, the nondimensional variables appearing above are given as with c 0 = √ gh 0 being the limiting long-wave speed.The small parameter α = a 0 h 0 represents the strength of the bore.The second small parameter is defined by β = h 2 0 l 2 , and following Whitham [24], we assume that β = O(α).Now an ansatz on the form of the velocity potential as is made.Substituting the form (7) into the Laplace equation and using the boundary condition on the bottom, gives where f = f0 is the nondimensional velocity potential on the bottom.Substituting equation ( 8) into the free surface conditions one obtains the system These are just a variant of the Boussinesq equations [6].The quantity ṽ appearing above is the first term in the expansion of the velocity φx , which is The fluid velocity w in equations ( 1) is modelled at a nondimensional depth z = 2 3 .Denoting the nondimensional velocity at this depth by w, the following relation is found when inverting (9) Replacing ṽ in the system above with the expression (10) gives, after some further manipulations of the higher-order terms (see [6]), the equations Omitting terms of quadratic or higher order in α and β , and reverting to dimensional variables, the system (1) appears.
We now turn to the energy considerations for the system derived above.To find the mechanical energy associated with the system (1) per unit width above the interval [x 1 , x 2 ], the integrals Here E k and E p are kinetic and potential energy, respectively.For potential energy, it is found that The potential energy scales as ρgh 2 0 l, which will be used when writing the energy in nondimensional form.Note that E p above includes both the potential energy of the rest state and the potential energy due to the wave motion.The expression for kinetic energy can be shown to be Using the expansion (8) for φ and the relation (10) in the integral above, and omitting terms of quartic or higher order in α and β , the nondimensional total mechanical energy associated with the dispersive system (1) is found as Note that the potential energy of the rest state has been used to normalize the total mechanical energy.This normalization was chosen for the sake of simplicity when comparing with the energy integral (4) associated with the shallow-water system.However, this choice of normalization necessitates the inclusion of terms of second and third order in α and β in the energy integral.

NUMERICAL METHOD
In this section, a strategy for solving the equations in (1) numerically with Dirichlet boundary conditions will be given.The system in question is now written as The system (11) is accompanied with the following time-independent Dirichlet conditions The equations (11) together with the boundary conditions given above are discretized using a finite difference method.A spatial grid of uniform grid size δ x is used on the interval [x 1 , x 2 ].On this grid the spatial derivatives are approximated by a second-order central-difference scheme.Let two grid functions, which approximate the solutions at time t, be denoted by V j (t) ≈ η(x 1 + j δ x,t) and W j (t) ≈ w(x 1 + j δ x,t).The boundaries are given at j = 0 and j = N + 1. Approximating the spatial derivatives in (11) with central differences yields the semi-discrete system of equations for j = 1, 2, . . ., N. The expressions for the right-hand sides are Table 1.Convergence study.For the convergence of the spatial discretization, the fixed time step is δt = 1 × 10 −2 s and the number of grid points is N = 2 k − 1.For the convergence of the temporal discretization, the fixed number of grid points used is N = 2 18  It appears now that the semi-discrete system can be written in the form where the matrix A is the discrete form of the operator 1 − x which is already encoded in the system (12).To integrate the solution in time, the well-known explicit four-stage Runge-Kutta method is used.This method is self-starting and is accurate to the fourth order in the size of the time step δt.As the Runge-Kutta time integrator is explicit, the right-hand sides above are only evaluated at times t for which the values V j (t) and W j (t) are already known.Moreover, using an explicit time integrator is seen to decouple the problem.
At each stage of the Runge-Kutta method, one takes an internal Runge-Kutta step for V j and W j separately, then updates the right-hand sides and continues the process.The operations performed on V j (t) and W j (t) in ( 12) are exactly the same, so that they can be carried out by the same subroutine.Once the time discretization is performed, it turns out that it is more advantageous to keep the matrix A on the left-hand side of the equation.Thus one ends up solving a system of the form Ax = b eight times at each time step.Since the matrix A is symmetric and tridiagonal, advancing the solution by one time step requires only O(N) operations.
Finally, a convergence study was conducted to check if the expected convergence rates of our scheme, O(δ x 2 , δt 4 ), are actually achieved in practical computations.The errors and convergence rates are obtained by comparing a numerical solution with an exact travelling wave solution of the system (1), given for instance in [9].The results of this convergence study are summarized in Table 1.

NUMERICAL STUDY
The focus in this section is on a numerical study of the undular bore, and the main goal is to gain an understanding of the energy conservation properties of the system (1) when applied to the bore.Light will also be cast on differences in the energy conservation properties of the dispersive system (1) and the shallow-water approximation (2).
Since Dirichlet boundary conditions are used at the endpoints of the interval, the flux of energy across the boundaries and into the domain, as given by (5), will be controlled and constant.The spatial domain [x 1 , x 2 ] must be chosen large enough to contain both the bore front and any disturbances propagating away from it.With such a set-up it is reasonable to assume that the values of the functions η and w at the endpoints equal the far field conditions, and are constant both in space and time.In dimensional variables, the energy corresponding to the dispersive system (1) is given by The last term can be written as h 3 0 (ww x ) x /3, which is zero when the integral is evaluated with the far field conditions at the endpoints.Therefore, the expression for the energy of the dispersive system is reduced to the formula (4) for the shallow-water energy.To distinguish these two, the energy for the dispersive system is denoted by E(η, w), and the energy corresponding to the shallow-water system is labelled E(h, u).The initial conditions are given by and the boundary conditions are given by The downstream flow velocity u 1 is given by the the bore conditions (3), and the model parameter κ is a measure of the steepness of the initial data.For a given value of κ, it is clear that a large enough choice of the domain [x 1 , x 2 ] will yield initial data that agree with the boundary conditions at the endpoints to machine precision.In Fig. 1 the initial surface profile (dashed line) is plotted along with the developed profile (solid line) at t = 6 s.While the graph only shows the main disturbance, the computational domain is actually between x 1 = −80 m and x 2 = 80 m.At each time step, the value of the energy integral E(η, w) is computed, and the rate of change of energy corresponding to the dispersive system (1) is calculated numerically.This value is then compared to the predictions of the rate of change of energy in the shallow-water equations as given by equation (6).
Such a comparison is shown in Table 2 for a range of α-values given as tabulated in column one.The second column contains the net flux of energy into the domain across the boundaries.The third and fourth columns give the rate of change of energy in the dispersive and in the shallow-water equations, respectively.The difference between the latter and the flux F 1 − F 2 is given in the last column.This difference is solely due to the approximate nature of the shallow-water system.The difference between F 1 − F 2 and the rate of change of E(η, w) is actually less than about 2 × 10 −4 % for the results given in Table 2. Table 3 shows the result for a case where there is a strong backflow u 2 = −0.8m/s.As can be gleaned from the table, Table 2.The values in Column 1 are α = a 0 /h 0 .Shown in the next three columns are the net energy flux, the rate of change of energy in the dispersive system, and the rate of change of energy for the shallow-water approximation (all in kg m/s 3 ), respectively.The last column shows the difference between the net energy flux and the rate of change of energy given by the shallow-water theory.For the dispersive system this difference is less than 2 × 10 −4 % in these experiments.The parameters h 0 = 0.1 m and u 2 = 0 were used.For the discretization, we used N = 2 16 − 1, δt = 1 × 10 −2 s, and the domain was given by x 1 = −80 m and  3.The values in Column 1 are α = a 0 /h 0 .Shown in the next three columns are the net energy flux, the rate of change of energy in the dispersive system, and the rate of change of energy for the shallow-water approximation (all in kg m/s 3 ), respectively.The last column shows the difference between the net energy flux and the rate of change of energy given by the shallow-water theory.For the dispersive system this difference was less than 7 × 10 −4 % in these experiments.The parameters h 0 = 0.1 m and u 2 = −0.8m/s were used.For the discretization, we used N = 2 16 − 1, δt = 1 × 10 −2 s, and the domain was given by x 1 = −80 m and the energy loss predicted by the shallow-water theory is an order of magnitude greater than for the bore propagating into still water.However, the dispersive system also conserves the energy in this case.

CONCLUSIONS
In this paper a study of undular bores has been conducted.The weakly nonlinear, dispersive system (1) has been used to model this phenomenon.The focus has been on weak bores, for which the common point of view is that the excess energy is disseminated by oscillations behind the bore front.The numerical study in the previous section goes far in confirming this point of view.In the numerical experiments, the energy entering the computational domain across the boundaries is carefully monitored, and at the same time, the rate of change of energy corresponding to the dispersive system (1) is computed.As shown in Tables 2  and 3, the energy which enters the domain is accounted for when the dispersive effects are included and the system (1) is used in the numerical experiments.While it is clear from the above discussion that no viscosity is needed to explain the energy loss across a weak bore, it should be mentioned that the results obtained in this article are given in an idealized setting.Some effects are missing, which for real flows and flow parameters could be as important as the dispersive effects.Among these are viscosity and wave breaking.In real flows, physical experiments show that boundary layers develop, and viscous effects become very important in the dynamics of the flow.Likewise, for increasing bore strengths, i.e. increasing values of α, one or more of the waves in the bore will eventually break.This in turn will lead to a loss of mechanical energy due to turbulent dissipation associated with the breaking event.This is not captured by the dispersive model, since this description does not include the possibility of breaking waves.
Finally, it should be mentioned that the study of bores and hydraulic jumps is pursued with vigour by several scientific communities.There are a great many works focusing on engineering and environmental aspects of river bores, such as sediment resuspension and bank erosion.See for instance [8,14] for some recent theoretical studies.For works of an experimental nature, the reader may consult [3,18], among many others.Some inquiries that are similar in spirit to the present study can be found in [10][11][12].In these works, a refined version of Whitham's average Lagrangian method is used to derive qualitative information using a number of different systems of equations, some of which are similar to the system studied here.The work presented in [4,5,15,21] is focused on the effect of incorporating a simple diffusive term in the dispersive equations, such as the model proposed in [16,17].