{'fraction_non_alphanumeric': 0.038250661162355676, 'fraction_numerical': 0.014964845513771528, 'mean_word_length': 4.452259539300158, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Third-order Finite Volume/Finite Element Solution of the Fully Nonlinear Weakly Dispersive Serre Equations

The nonlinear weakly dispersive Serre equations contain higher-order dispersive terms. This includes a mixed derivative flux term which is difficult to handle numerically. The mix spatial and temporal derivative dispersive term is replaced by a combination of temporal and spatial terms. The Serre equations are re-written so that the system of equations contain homogeneous derivative terms only. The reformulated Serre equations involve the water depth and a new quantity as the conserved variables which are evolved using the finite volume method. The remaining primitive variable, the velocity is obtained by solving a second-order elliptic equation using the finite element method. To avoid the introduction of numerical dispersion that may dominate the physical dispersion, the hybrid scheme has third-order accuracy. Using analytical solutions, laboratory flume data and by simulating the dam-break problem, the proposed scheme is shown to be accurate, simple to implement and stable for a range of problems, including discontinuous flows.


Introduction
Rapidly-varying free surface flows are characterized by large surface gradients. These gradients produce vertical accelerations of fluid particles resulting in a non-hydrostatic pressure distribution. Equations that assume that the flow has a non-hydrostatic pressure distribution contain third-order dispersive terms. A system of equations that contain dispersive terms are the Serre equations. The Serre equations are fully nonlinear weakly dispersive equations. They are applicable to waves where the water depth h 0 ≪ L is much smaller than the horizontal wave length, L and up to wave breaking [1]. Bonneton et al. [2,3] consider the Serre equations as the most appropriate system for modelling highly nonlinear weakly dispersive waves at the shoreline.
A major difficulty with solving equations that contain dispersive terms, is that the dispersive terms usually contain a mix derivative term [4]. By replacing the mix derivative term in the flux term by a combination of temporal and spatial derivative terms, the Serre equations can be written in conservation law form, where the system of homogeneous equations contains a new conserved quantity and its corresponding flux term. The conserved quantities are evolved using a standard scheme for solving conservative laws. The remaining primitive variable is obtained by solving a second-order elliptic equation. A finite volume/finite element technique is proposed for the solution of the fully nonlinear and weakly dispersive Serre equations without the need for iteration or operator splitting necessary for dealing with the mixed derivative term.
It is well known that odd-order schemes introduce numerical dispersion and even-order schemes numerical dispersion. First-order and second-order schemes were developed by Zoppou and Roberts [5], using the methodology described in this paper, to solve the Serre equation. Both schemes produced dispersive waves that accurately predict the arrival of the initial wave and its amplitude. In addition, the phase of the predicted dispersive waves is very close to the recorded wave profile. However, the diffusion introduced by the first-order scheme rapidly dampens trailing dispersive waves. The second-order model slightly overestimates the amplitude of the dispersive waves. This could be a result of the use of a second-order scheme which introduced numerical dispersion. To avoid the introduction of numerical dispersion and excessive diffusion, a third-order finite volume/finite element scheme was developed.
The performance of the proposed third-order finite volume/finite element scheme for solving the conservative form of the Serre equations is evaluated using an analytical solution to the Serre equations, laboratory flume data and by simulating the dam-break problem. With the exception of the analytical solution, which is smooth, the remaining problems involve the simulation of flows with steep gradients that produce dispersive waves.
This paper is organized with the derivation of the standard Serre equation in Section 2 followed by the derivation of the alternative form of the Serre equation in terms of the new conservative variable. The properties of the linearized form of the Serre equation are also examined in Section 2. The third-order implementation of the proposed scheme is described in detail in Section 3 and the fully discretized scheme summarized in Section 4. In Section 5, the numerical scheme is validated using an analytical solution and laboratory flume data. The stability of the proposed scheme is demonstrated by simulating the dam-break problem. Finally, the performance of the numerical scheme is discussed in Section 6.

Serre Equations
For an invicid incompressible fluid with constant density, ρ the conservation of mass and momentum are given by the Euler equations and where in two planar dimensions, x = (x, z), a fluid particle at depth ξ = z − h − z b below the water surface, where h(x, t) is the water depth and z b (x) the bed elevation, see Figure 1, is subject to fluid pressure, p(x, t) and gravitational acceleration, g = (0, g) T , has a velocity u = (u(x, t), v(x, t)), where u(x, t) is the velocity in the x-coordinate and v(x, t) is the velocity in the z-coordinate and t is time. In addition to the above equations, a number of boundary conditions must be satisfied. These are; (a) the kinematic condition at the free surface (z = h + z b ), which is the atmospheric pressure at the water surface, usually taken to be p a = 0. The Serre equations assume that the point velocity in the x-direction is uniform over the water depth, so that u(x, z, t) =ū(x, t) with From (1a) it follows that the vertical velocity at any depth z − z b is given by for a horizontal bed. The vertical velocity is a linear function of the water depth.
Integrating the point quantities in (1b) over the flow depth z b to h + z b , and satisfying (2) produces the one-dimensional equations and where h andū are the primitive variables. 3 The pressure distribution in the water column is given by Multiplying (3b) by h, adding (3a) pre-multiplied byū and making use of (3a) to obtain; and which is written in terms of the conservative variables, h andūh. The continuity equation is exact because it is based on the depth-averaged velocity, which makes no assumpltion on the distribution of u(x, z, t) with depth. However, the momentum equation relies on the assumption that u(x, z, t) is uniform with depth. The terms in the square parenthesis are the dispersive terms which contain high order spatial derivative terms and a mixed spatial and temporal derivative term. Ignoring all the dispersive terms in (5b) results in the well known nonlinear shallow water wave equations, where the pressure distribution is hydrostatic, p(ξ) = ρgξ.
Equation (5) are known as the Serre equations [6][7][8], they retain full nonlinearity in the dispersive terms [9]. They have been derived by Serre [6], Su and Gardner [10] and Seabra-Santos et al. [7] and are equivalent to the depth averaged Green and Naghadi [11] equations. They are considered to be good approximations to the full Euler equations up to a wave breaking [2,3].

Alternative Conservative Form of the Sere Equations
The flux term in the momentum equation, (5b) contains a mixed derivative term which is difficult to treat numerically. It is possible to replace the mix spatial and temporal derivative term by a combination of spatial and temporal derivative terms. Consider Rearranging then Making use of the continuity equation, (5a) and the momentum equation, (5b) becomes The momentum equation can be written in terms of a new conservative form as where the new conserved quantity, G is given by the second-order elliptic equation The temporal derivative in the momentum equation has been eliminated from the flux term. In contrast to (5), the flux term now contains spatial derivatives only. For the remaining primitive variable,ū, if the data is square integrable in a rectangular domain, then from the regularity theorem of elliptic partial differential equations [12]ū ∈ H 2 in Sobolev space of square integrable second-derivatives. The primitive variable,ū will be smooth.
The alternative form of the Serre equations can be written in vector form as where the vector of state variables and

Properties of the Linearized Serre equations
Although they are evolution-type equations, the Serre equations are neither hyperbolic or parabolic and do not have any Riemann invariants. However, it is possible to establish some properties of the Serre equations by examining the behaviour of harmonic waves of the form h(x, t) = Ae i(kx−ωt) and u(x, t) = Ue i(kx−ωt) (8) where A and U are unknown coefficients, ω is the frequency, k = 1/L is the wave number and i = √ −1 in the linearized Serre equations and Substituting (8) into (9), the linearized equations become For a non-trivial solution In this case the dispersive terms have no effect on u 0 , only on the celerity of a small disturbance.
This compares with the frequency, ω 1,2 = (ū 0 ± gh 0 )k, for the shallow water wave equations. As µ → 0, the frequency for the Serre equations are identical to that of the shallow water wave equations. When µ → ∞, ω 1,2 =ū 0 . Therefore, the frequency for the Serre equation are bounded by the wave frequency of the shallow water wave equations.
For non-dispersive waves, the phase velocity, υ p = Re(ω)/k is identical to the group velocity υ g = dRe(ω)/dk. This is not the case for the Serre equations, where the phase speed is υ p = u 0 ± gh 0 3 µ 2 + 3 and the group velocity is Both are dependent on the wave number. Since the group speed is slower than the phase speed then the Serre equations describe dispersive waves.

Numerical Scheme
In a finite volume scheme the cell average values at the nodes,q j are updated by integrating (7) over the jth cell I j = [x j−1/2 , x j+1/2 ] to obtain a semi-discrete scheme where ∆x = x j+1/2 − x j−1/2 is assumed to be constant, x j = (x j+1/2 + x j−1/2 )/2 and is the average value of the state variable q(x, t) in I j at time t, which ensures that mass is conserved in each cell and the integral formulation admits shocks in the solution.
The flux , F j+1/2 is a function of the left and right extrapolated state values q + j+1/2 and q − j+1/2 , obtained from piecewise polynomials, P j (x j+1/2 ) and P j+1 (x j+1/2 ), respectively passing through consecutive values ofq j . Therefore, The reconstruction will usually produce two different values for q(x j+1/2 ). Generally, there will be a discontinuity in the state variables at x j±1/2 .
The flux of material across the interface of a cell is estimated by solving the Riemann problem, defined by the initial value problem Once the intercell flux has been established, the cell average values can be updated by evolving the solution over a single time step by solving the semi-discrete system, (12) using an ordinary differential equation solver. The overall accuracy of the numerical scheme is dependent on the accuracy of the reconstruction method and the order of accuracy of the scheme used to evolve the solution in time.

Inter-cell Flux Evaluation
The numerical approximation of the physical flux f (q(x, t)) across the boundary of a cell, F j+1/2 is given by the explicit upwind central scheme proposed by Kurganov et al. [13] as At the interface of a cell, x j±1/2 a discontinuity in the state variable will propagate with rightand left-sided local speeds, which are estimated by where λ 1 and λ 2 are the smallest and largest eigenvalues, respectively of the Jacobian ∂ f (q)/∂q which correspond to the phase speeds.
It has been established that the eigenvalues are bounded for the Serre equations and that the eigenvalues, λ 1 ≃ū + gh, λ 2 ≃ū − gh for the shallow water wave equation can be used as estimates of the upper and lower bounds of the eigenvalues for the Serre equations.
There is a restriction on the computational time-step that can be used in all explicit schemes. Stability is satisfied when the time step ∆t satisfies the Courant-Friedrichs-Lewy, (CFL) criteria [14] ∆t < ∆x 2max(|λ i |) ∀ i.

Reconstruction
To achieve third-order, O(∆x 3 ) accuracy, it is sufficient to consider quadratic piecewise polynomial reconstruction in each cell using centered at the computational node j. Over each cell, the cell average The mass conserving quadratic is given by where −1 ≤ κ ≤ 1. Only when κ = 1/3 is the reconstruction third-order. Any other value of κ results in linear reconstruction and a second-order scheme. The reconstructed cell edge values at x j±1/2 are estimated from (15). Consider the cell interface at x j+1/2 , then (15) becomes Similarly at the left boundary of the cell, at x j−1/2 The interpolated values must be limited so that no new extremum is introduced. Consider (16) and making use ofq j+1 −q j−1 =q j+1 −q j + (q j −q j−1 ), then Dividing the term in the square brackets byq j −q j−1 then 3 is a nonlinear limiter which can be used to prevent unwanted oscillations and ensures that the results are physical (bounded) and therefore stable.
Equation (17) forms the basis of the third-order symmetrical Koren limiter [15] given by where lim r j →∞ φ + (r j ) = 2 is used in the proposed model. This limiter ensures that the scheme remains TVD [16] and third-order away from extrema, where r j ≤ 0 and φ(r j ) = 0. This occurs when the gradient changes sign indicating that an extrema has been encountered within the cell.
The reconstruction reverts to a piecewise constant reconstruction. In smooth regions, r j → 1 and φ(r j ) → 1 and the reconstruction is third-order. Similarly for the reconstructed values at the left boundary of the cell, the limited values are given by where and the corresponding nonlinear limiter is φ + (r j ) = max 0, min(2r j , (2 + r j )/3, 2) .

Nodal Values
The point values for the conservative quantities, q j are estimated from the cell averages,q(x j ) using quadratic interpolation, (15) so that Written for all computational nodes results in a tridiagonal matrix, M which can be solved for the cell averages,q(x j ) given the nodal values, q j .

Time Integration
Time integration of the semi-discrete system (12) is performed using Strong Stability Preserving (SSP) Runge-Kutta schemes. SSP schemes involve a convex combination of first-order forward Euler steps that preserve the TVD properties of the Euler scheme [17,18].
A third-order three-stage SSP Runge-Kutta scheme is given by [17,19] is the discrete form of (12) and It is also subject to the time restriction, Cr = 1.
It is relatively straight forward to establish the boundary information required for the thirdorder Runge-Kutta TVD scheme by examining equation (20). In this case boundary information is requires at the time levels; t n , t n + ∆t/2, t n + ∆t, t n + 3∆t/2 and t n + 2∆t.
Combined with quadratic reconstruction, the resulting numerical scheme is theoretically O(∆x 3 , ∆t 3 ) accurate.

Fully Discrete System
The reconstructed point values,ū ± j+1/2 , h ± j+1/2 and G ± j+1/2 at each cell interface are calculated from the cell averaged values using (17) and (18). These are used to estimate the left and right intercell flux using a finite-difference discretization of (7c). The approximate Riemann solver (13) is used to obtain the intercell flux required in (21). Using an SSP Runge-Kutta scheme, the solution is advanced in time providing updated values for the cell average conservative variables, h n+1 j andḠ n+1 j . The nodal values for the remaining primitive variableū n+1 j is obtained by solving the second-order elliptic equation, (6) given point values for h n+1 j and G n+1 j , obtained using (19) and the cell average values, using finite elements..
The third-order finite difference discretizations of f (2) in (7c) and the third-order finite element solution of (6) are described in detail in the following Section.

Discretization of the Instantaneous Flux, f
The spatial derivatives in the local flux term, (7c) are evaluated using a quadratic polynomial fitted through the values, u j and u j±1/2 , which were obtained from the finite element solution of the second-order elliptic equation, (6) . The dircretized second component of the local flux can be written as and for the first component

Finite Element Solution of the Elliptic Equation
In finite elements, the second-order elliptic equation (6) for an element In the direct approximation of the finite element method, the desired function q(x) is approximated by a weighted finite series involving the unknown nodal values, q j of the desired function. The set of locally defined piecewise functions, w are also known as basis functions and N is the total number of elements in the computational domain. In the Galerkin weighted residual method υ(x) = w(x).
For a third-order scheme, quadratic basis functions are used in each element, shown in Figure  2, spanning three nodes, j − 1/2, j and j + 1/2. The finite elements coincides with the finite volume cells.
The Taylor series of a one-dimensional function, q(ξ) and its derivatives including the coordinate transform mapping [x j−1/2 , x j , x j+1/2 ] to the ξ-space [−1, 0, 1] is given in terms of the nodal values by where x = ξ∆x/2 + x j and dx = dξ∆x/2 for a uniform grid. Solving for q(ξ) and its derivative provides The coefficient matrix represents the weights for the quadratic interpolation of a quantity and its first derivative in space by the first and second rows respectively. The first row also represents the isometric mapping between x and ξ.
Which is satisfied by the Consider each term in (23), which can be written as Q eūe , where the stiffness matrix, Q e is given in the Appendix andū e = [ū j−1/2 ,ū j ,ū j+1/2 ] T . Consider the term Which can be written as P eūe where P e = ∆x 420 The final term Evaluating the integrals yields For each element [Q e + P e ]ū e − R e = 0 13 The assembled system of equations can be written in matrix form, The conservative quantities, h and G can be discontinuous across each element, see Figure  2, whilst the unknown quantity,ū is assumed to be continuous and approximated by a piecewise quadratic function. The valuesū j andū j±1/2 are use to calculate the gradient used in the second component in (7c), shown in (22) without a loose in accuracy.

Ensembles Scheme
The third-order three-stage strong stability preserving Runge-Kutta scheme solution of the Serre equations involves the following steps Step 1: Givenq j = [h j ,Ḡ j ] the remaining primitive variableū = H[q] is obtained by solving the second-order elliptic equation, (6) using the third-order finite element scheme. The operator H involves the following steps where A is the reconstruction of the point values at the cell interface using the Koren limiter and F is the third-order finite element solution of the second-order elliptic equation for the point values u e using (25).
Step 2: Perform the reconstruction using the cell averages,q = [h j ,Ḡ j ] T and solve the local Riemann problem to obtain the flux F j±1/2 of material across a cell interface. Evolve the solution using a first-order Euler time integration for the conserved quantities,q j .

Steps 3 and 4:
Repeat Steps 1 and 2 with the values obtained from Step 2 and evolve using another first-order Euler step.
Step 5: Calculate weighted intermediate values using the initial values and the results from Step 4.
Steps 6 and 7: Repeat Steps 1 and 2 using the intermediate values.
Step 8: The solution at the next time level is obtained by averaging the initial values and the values obtained from the third Euler step. This completes the third-order strong stability preserving Runge-Kutta time integration.

Numerical Simulations
Convergence rate of the proposed schemes is determined using a known analytical solution to the Serre equations. Data from two laboratory experiments are used to validate the proposed models and the simulation of the dam-break problem is used to show that the model is stable for simulating a wide range of flow problems.

Analytical Solutions
Simulating the propagation of solitons is a common test for Boussinesq-type equations. The Serre equations, (5) has the following analytical soliton solution [9] h(x, t) = a 0 + a 1 sech 2 with κ = √ 3a 1 /(2a 0 √ a 0 + a 1 ) and c = g(a 0 + a 1 ). Solitary waves propagate at constant speed without deformation. Therefore, there is a balance between nonlinear and dispersive effects. A numerical scheme must accurately model the equilibrium between amplitude and frequency dispersion in order to simulate the propagation of the wave profile at constant shape and speed. A result of a poorly balanced numerical scheme is the simulation of trailing edge dispersion waves which cause a reduction in wave height and celerity.
The results from a numerical scheme are compared to the corresponding analytical solution by using the non-dimensionless L 1 norm where, q j is the simulated values of q(x, t) at x j , and q(x j ) is the corresponding analytical solution. The L 1 norm is calculated using all the computational nodes, j = 1, . . . , m.
The prototypical example is a solitary wave predicted by (26) with, a 0 = 10m, an amplitude of a 1 = 1.0m. The soliton has a celerity, c = 10.387974m/s and κ = 0.026112/m.
The boundary conditions imposed on the models are a 0 (x) = 10m andū 0 (x) = 0m/s at the upstream and downstream boundaries. The model parameters are; ∆x = 0.1m, Cr = 0.1 and ∆t = Cr/∆xs. Using these parameters, the initial soliton profile and velocity, the analytical and the simulated water depth and velocity at t = 100s are shown in Figure 3 for the third-order scheme. The soliton amplitude is accurately predicted and the soliton speed is captured correctly.
Performing the simulation for a range of ∆x and keeping Cr = 0.1, the L 1 norm between the simulated and analytical solution was calculated for the water depth and fluid velocity. Plotting the log 10 L 1 against log 10 ∆x reveals that the proposed strategy for solving the Serre equations is third-order accurate, see Figure 4. Here, the usual convergence behaviour is observed. For coarse grids, the numerical scheme is not able to resolve the oscillations and the true convergence rate of the scheme is not realized. As the grid is refined convergence is confirmed and for very small grid spacing the error deteriorates because it is dominated by the accumulation of roundoff errors.
Clearly, for the simulation of the smooth soliton problem, the third-order schemes is capable of predicting the soliton speed and its amplitude.

Labratory Experiments
The frictionless horizontal flume experiment from Hammack and Segur [20], involving a negative amplitude rectangular wave and the more recent surge propagation experiment conducted by Chanson [21], are used to validate the hird-order schemes described in Section 3. Both produce highly dispersive waves from an abrupt change in the initial flow conditions. In these experiments the non-hydrostatic terms cannot be neglected in the momentum equation.

Undular Bore
An undular bore was created in a large tilting flume at the Civil Engineering Department, University of Queensland. The channel is 0.5m wide, 12m in length and the undular bore was created in the horizontal flume, which has a smooth PVC bed and glass walls. A radial gate located at the downstream end of the flume, x = 11.9m controls the water depth in the flume. The radial gate is used during the experiments to produce steady subcritical flow in the flume which remains constant for the duration of the experiment. Steady flow condition are established for 15 minutes prior to an experiment. Adjacent to the radial gate is a rapidly closing Tainter gate at, x = 11.15m that spans the full width of the flume. An undular bore is generated by the rapid closure of the Tainter gate, which is estimated to take less than 0.2s, when water accumulates at the Tainer gate forming an upstream progressing undular bore. The experiment ceases when the bore reaches the upstream intake structure to avoid any interference from wave reflection. Acoustic displacement meters, located at the flume centerline at; x = 8.0, 6.0, 5.0, 4.55, 4.0 and 3.0m record the progress of the bore and dispersive waves with time. Data acquisition starts 30 seconds prior to the closure of the Tainter gate.
The recorded water surface profile at the acoustic displacement meters over time are shown in Figure 5 along with the simulated water surface profile predicted by the model. At all locations in the flume, the dispersive waves are symmetrical about the predicted water level. The simulated results show that the simulated bore speed is slightly slower than the observed bore speed. This is the theoretical observation, where the group and phase speed of waves for the Serre equations are slower than for the shallow water wave equations. These results show that the third-order scheme has accurately predicted the arrival of the bore. In addition, it has accurately predicted the amplitude of the dispersive waves which have a slightly longer wavelength than the actual dispersive waves.

Rectangular Initial Wave
In the laboratory experiment conducted by Hammack and Segur [20], a wave maker consists of a rectangular piston 61cm in length at the end of a wave tank spans the full width of the tank. The tank is 31.6m in length, 61cm deep and 39.4cm wide, horizontal with vertical sides and is constructed from glass. The piston moved monotonically from its initial position, which is flush with the tank bed to its final elevation. It can be displaced vertically up or down. The upstream wall of the wave tank adjacent to the wave maker is a plane of symmetry. The length of the piston, b = 61cm represents the half-length of a hypothetical piston occupying the region The symmetrical problem is simulated using the numerical schemes. A rectangular wave propagates following a sudden downward 1cm (h 0 = 9cm) movement of the piston. The quiescent water depth, h 1 is fixed at 10cm. The water elevation is recorded at the fixed locations; x/h 1 = 0, x/h 1 = 50, x/h 1 = 100, x/h 1 = 150, and x/h 1 = 200, where x/h 1 = 0 is the downstream edge of the piston.
The upstream and downstream boundary conditions remain constant at; h 1 = 10cm and u 1 = 0m/s. In the simulations ∆x = 0.0005m, Cr = 0.2, ∆t = Cr∆x/ 0.1g and the solution is terminated at t = 50s.
There are several wave train following the dominant first wave train, see Figure 6. The amplitude of the dispersive waves in these is much smaller than the amplitude of the dispersive waves that immediately follow the shock.
Results from the third-order scheme, are also shown in Figure 6. There is excellent agreement between the simulated and observed results. The rarefaction wave, shock speed and the phase of the dispersive waves are faithfully reproduced by the numerical scheme. In addition, secondary wave trains are also reproduced by the numerical scheme, see for example Figure 6(c) at t g/h 1 − x/h 1 ≈ 110. The amplitude of the dispersive waves, though are slightly overestimated by the numerical scheme. This is not surprising since the Serre equations do not contain terms that represent the internal viscosity of the fluid.

Dam-break
The dam-break problem is a standard test for models used to solve the shallow water wave equations, which has a known analytical solution (see, for example Zoppou and Roberts [22]). It has been chosen to demonstrate the flexibility of the proposed model for simulating both subcritical and supercritical problems.
The dam-break problem is solved using the proposed third-order solution to the Serre equations. The simulated results have been plotted against the analytical solution to the shallow water wave equations for the dam break problem, which is used as reference data. The dam-break occurs in a frictionless rectangular channel, 1000m in length where the initial velocity of the water u = 0m/s and the water depth upstream of the dam, which is located at x = 500m is given by h 1 and downstream of the dam by h 0 . In all the models, ∆x = 0.1m, Cr = 0.2, ∆t = Cr∆x/ gh 1 and the solution is terminated at t = 30s. Three cases are considered; h 1 = 10m with h 0 = 1m, h 1 = 10m with h 0 = 2m and h 1 = 1.8m with h 0 = 1m. These have as their maximum Froude numbers; Fr = u/ gh = 1.18, 0.81 and 0.29 respectively, which were obtained from the analytical solution to the shallow water wave equations. The three problems involve supercritical flows, near critical flow and subcritical flows.
The simulated results using are shown in Figures 7. The arrival of the shock is accurately captured, as is the rarefaction fan and the shock height. The results for the third-order scheme, shown in Figure 7(c) are very similar to those obtained by El et al. [9], who used a second-order Lax-Wendroff scheme and the second-order finite volume model described by Zoppou et al. [23] to solve the Serre equations.
An interesting feature of the results shown in Figure 7 is that the oscillations are bounded. They are restricted to the minimum and maximum initial water depth. The simulated water velocity is also bounded. The Serre equations also conserve energy. Since energy for the Serre equations is uniformly bounded [24] then the solution is also bounded by the energy of the initial conditions.

Conclusions
By replacing the mix derivative term in the flux term by a combination of temporal and spatial derivative terms, the Serre equations can be written in conservation law form, where the system of homogeneous equations contains a new conserved quantity and its corresponding flux term. The water depth and the new conserved quantity is evolved using a finite volume scheme. The water velocity, which is the remaining primitive variable is obtained by solving a second-order elliptic equation using finite elements.
Using analytical solutions, laboratory flume data and by simulating the dam-break problem, the proposed third-order hybrid finite volume/finite element scheme is shown to be simple to implement and stable for a range of problems including rapidly varying flows. It accurately predicts the phase, arrival of the dispersive waves and their amplitude that are associated with rapidly varying flows.
The proposed scheme is currently being extended to problems with bathymetry, to twodimensional problems and to other system of dispersive equations.

Appendix A. Stiffness and Load Marix in the Finite Elemen Scheme
The coefficients in the stiffness matrix Q e are given by

23
The stiffness matrix, A, the unknowns,ū and the stiffness matrix, R for the finite element solution to the second-order elliptic equation, (6) are given by  where the superscript denotes the element number and the coefficients α and β are provided by the essential Dirichlet boundary conditions.