WAVE BREAKING USING A ROLLER APPROACH IN A HYBRID FINITE-VOLUME FINITE-DIFFERENCE BOUSSINESQ-TYPE MODEL

A new scheme implementing a roller approach into a hybrid finite-volume finite-difference Boussinesq-type model is presented. The relevant mathematics are outlined and a numerical solver is described. Predictions obtained from the model are validated against physical observations, demonstrating the capabilities of the scheme to replicate the complex hydrodynamic processes that occur during wave breaking. The benefits of both the hybrid scheme and the roller approach are discussed. The results illustrate the feasibility of modelling the breaking process with a rotational roller method in a finite-volume finite-difference scheme and show the obtainable accuracies. Finally, further tests and improvements to the model are proposed.


INTRODUCTION
Boussinesq-type equations (BTEs) allow numerical models to simulate wave propagation from intermediate depths to shallow water.The onset of wave breaking in shallow water requires that schemes dissipate energy to accurately replicate the hydrodynamics across this region.Several approaches have previously been adopted to represent wave breaking within numerical models.
Recently, numerous schemes have utilised one such approach where a switch from BTEs to nonlinear shallow water equations (NSWEs) is performed at smaller depths.Here, dispersive terms are removed when the ratio of the wave height to the local water depth exceeds a threshold value.This has commonly been coupled with hybrid finite-volume finite-difference (FVFD) schemes, such as the one used in the present work (Tonelli and Petti, 2009;Shi et al., 2012).Good levels of stability have been observed when using this approach, and the need for strong numerical filtering is negated.Indeed, using a filter to remove noise introduced by higher order terms has previously led to significant numerical dissipation which has masked the physical behaviour of the system.Moreover, simulating breaking by removing higher order terms provides a less physical description of the hydrodynamics and reduces the accuracy with which the free surface profile is replicated.
The roller approach is another energy dissipation model which allows BTEs to be used throughout the breaking region.When deriving BTEs with a potential flow formulation, rotation within the fluid is no longer resolved, but greater computational efficiency is gained.While this is acceptable in deeper water where the flow is more uniform along the water column, it means energy dissipation is no longer seen where the flow is more turbulent, as is the case during wave breaking.Defining a turbulent region within the flow allows the rotational component to be introduced by injecting an appropriate level of vorticity.Additionally, by determining the rotational component of the flow, horizontal velocity profiles can be obtained.This breaking approach provides a more physical description of the flow behaviour and has previously been combined with finite-difference (FD) schemes (Veeramony and Svendsen, 2000;Briganti et al., 2004;Musumeci et al., 2005).
Although the roller approach has not previously been implemented in a FVFD scheme (Brocchini, 2013), doing so yields several benefits.The FVFD scheme maintains good stability without the need for filtering, which ensures the flow dynamics reproduced by the model are physically grounded, and not dominated by numerical effects.The breaking terms introduced when using the roller approach provide energy dissipation based on physical characteristics.By retaining the Boussinesq terms, the wave shape and hydrodynamic behaviour can be better represented.

GOVERNING EQUATIONS
The present scheme uses the weakly-nonlinear formulation with breaking terms estimated using a rotational roller approach, as outlined by Veeramony and Svendsen (2000) and Musumeci et al. (2005), and is written in conservative form to suit the finite-volume (FV) solver used by the model.This requires  the equations to be written in terms of the total water depth, d, and flow rate, q, (see Fig. 1) which provides with (3) and Here, t denotes time, x the horizontal dimension, h the still water level, g the acceleration due to gravity and the free surface elevation.The breaking terms consist of which represents the convective momentum term, which represents the pressure due to vertical acceleration caused by rotational flow, and the viscous shear stress, The term u denotes the velocity, u r the rotational velocity, ūr the depth averaged rotational velocity, ν t the eddy viscosity, and z the vertical dimension.The rotational velocity is calculated from the vorticity, ω, according to By dividing the total horizontal velocity into potential and rotational components, such that the depth averaged velocity can be written as where the dispersivity is given by µ with k denoting the wavenumber.Since u p is constant across the water column for terms up to O µ 2 , it can be noted that ūp This allows the total velocity to be determined by The model determines the distribution and evolution of the vorticity using the vorticity transport equation: A perturbation approach is used to calculate a semi-analytical solution to this equation, with Here, σ represents a new coordinate system described by the relation where z = 0 defines the still water level and ζ s the roller thickness.The boundary conditions imposed on the vorticity equation are at the lower edge of the roller and ω| σ=0 = 0 (19) at the bed (see Fig. 1).This requires the definition of the roller thickness, and an injected vorticity, which are derived from hydraulic jump observations due to their similarity with breaking waves (Veeramony and Svendsen, 2000).Here, the depths at the crest and toe of the roller are given by d c and d t respectively, ω m = 6 defines the strength of the injected vorticity and the roller coordinate system is defined by where x c and x t give the locations of the roller crest and toe respectively.The initial vorticity is stated to be zero, such that ω| t=0 = 0. ( 23)

NUMERICAL SCHEME
The hybrid method uses an FV scheme to solve the NSWE terms given in Eq. ( 1).At each timestep, the conserved variable is updated according to Here, i and n indicate spatial and temporal FV cell indices respectively.W i is determined using a tridiagonal solver and is defined by where and G (n) i is calculated using an Adams-Bashforth-Moulton predictor-corrector scheme, where is used at the predictor stage and at the corrector stages.Calculation of these terms requires the FV flux term and the FD source term to be determined and then combined to give The corrector stage is repeated until successive steps produce sufficiently similar values of i obtained following the previous iteration.

Finite-volume flux solution
The FV part of the scheme requires the solution of the flux term, with the intercell fluxes obtained using the HLL approximate Riemann solver (Harten et al., 1983).The flux at each cell boundary is evaluated according to the states of adjacent cells: The wave speeds of the left and right cells are given by and respectively and the fluxes by and Here, Finally, conserved variables are reconstructed at each cell boundary according to where ∆ Hi defines the limited slope, determined from the free surface profile and flow rate, using a MUSCL-TVD scheme (Leer, 1984;Yamamoto et al., 1998;Toro, 2009).The cell slope is calculated as where the van-Leer limiter is defined by and the limited intercell jumps are calculated as Here, the intercell jump is and the Minmod limiter is defined by

Finite-difference solution
Calculating the FD component of the hybrid scheme is comparatively straightforward.The source term, S i , given in Eq. ( 31) introduces the higher order terms and breaking terms given in Eq. ( 4) into the model, with derivatives calculated to second-order accuracy.

MODEL VALIDATION
The capabilities of the model are assessed by repeating tests previously performed in a laboratory, allowing the numerical results to be compared with the experimental data collected in each case.The tests conducted here are selected in order to evaluate the breaking performance of the model in particular.As such, the experiments performed by Hansen and Svendsen (1979) are chosen to look in detail at the wave height decay in the surf-zone.The quantity of vorticity injected into the model is also assessed for the tests of Cox et al. (1995), allowing the resulting velocity profiles to be examined alongside physical results.
The simple uniform slope shown in Fig. 2 is used in all cases, with monochromatic waves generated offshore over a horizontal bathymetry.The numerical domain includes a 5 m sponge layer at the offshore boundary and introduces waves into the system using the internal generation method described by Wei et al. (1999).The moving shoreline is handled intrinsically by the Riemann solver.

Hansen and Svendsen (1979) tests
The data provided by the physical experiments of Hansen and Svendsen (1979) allows the influence of the breaking scheme on the wave height to be investigated.The parameters used for these tests are given in Table 1.The still water depth is h = 0.36 m and the distance from the start of the slope to the initial shoreline is l = 12.33 m.Several wave generation conditions are selected to provide a range of tests with both plunging and spilling breakers.In order to quantify the dissipation introduced by the breaking scheme, the model predictions are compared to detailed experimental wave height data.Fig. 3 presents the numerical results alongside the laboratory data.The results obtained when switching to NSWE to replicate breaking are also provided to demonstrate the model performance in each case.
The results show that, when used in conjunction with a hybrid FVFD scheme, the roller method is capable of providing a good prediction of the wave height decay due to breaking.The weakly-nonlinear nature of the model reduces the level of shoaling present, but consideration of this in calibration of the breaking parameters ensures breaking initiation is correctly located.The resulting dissipation of energy is close to that seen in the physical tests.The NSWE breaking approach produces some oscillations due to the discrepancy introduced when the scheme switches, but provides a similar decrease in wave height.
For spilling breakers, the model is able to reproduce the decay in wave height observed in the physical experiments.The roller approach allows a greater amount of energy dissipation to be introduced than is possible with the transition to NSWEs.Switching from BTEs to NSWEs provides a simple breaking method, but relies solely on the inherent properties of the equations, and does not allow any calibration to refine the amount of rotational flow introduced into the system.The roller approach instead adopts physical parameters, which influence the energy removed from the fluid.The plunging breaking present during Test R highlights the limitations of the roller method.The separation of flow seen in such cases means multiple free surfaces exist, thereby limiting the applicability of a depth averaged scheme.The bulk of the turbulence is assumed to exist within the roller region, where there are high levels of air entrainment, but the entire water column is treated as continuous.Despite this, the dissipation produced means the model is still able to provide a reasonable prediction of the wave height in the surf-zone.

Cox et al. (1995) tests
The experiments performed by Cox et al. (1995) provide data on the free surface and velocity under a breaking wave.The wave period is T = 2.2 s and a still water depth of h = 0.4 m is used, with the shoreline l = 14 m from the start of the slope (Cox et al., 1996).
Since the breaking terms are calculated from the vorticity predicted by the roller approach, it is important to study the distribution of this quantity across the surf-zone.The evolution of the roller and the resulting vorticity during breaking can be seen in Fig. 4. The assumption of irrotationality in the absence of breaking means no vorticity is present without a roller.As the wave propagates into shallower water, the criterion for breaking initiation is met and a roller region is defined.The injection of vorticity along the lower edge of the roller is clearly visible, with the distribution defined in Eq. ( 21) causing the majority to be focused at the toe.The energy dissipated by this rotational motion leads to the reduction in wave height discussed previously.
By acquiring a prediction of the vorticity within the fluid, it is possible to calculate horizontal velocity profiles using Eq. ( 14).Fig. 5 illustrates this velocity at several moments during the breaking process.The largest shoreward velocity is seen at the crest of the wave, while a smaller seaward velocity can be found at the trough.The peak velocities are found along the lower edge of the roller and are seen to reduce as the breaking continues, since much of the energy is dissipated during this process.
The distribution and magnitude of the velocities obtained by the numerical model can be validated by comparing the results with physical observations.Fig. 6 includes the velocities recorded by Cox et al. (1995) using a laser Doppler velocimeter (LDV) and gives the equivalent gauge data provided by the present model.Before the breaking process has initiated, the velocity is found to be uniform along the water column, with some residual variation left in the wake of a breaking wave due to the advection of vorticity.
Since no LDV data is available in the upper region of the wave, the velocities near the roller cannot be reliably assessed, but a good agreement is seen across the rest of the flow.The model is able to accurately reproduce both the wave shape and horizontal velocity profiles.This suggests that the quantity and distribution of the vorticity injected by the breaking model is accurate, introducing an appropriate level of rotational flow and creating breaking terms that correctly dissipate energy.

CONCLUSIONS AND FURTHER WORK
The present work provides a new model incorporating a roller approach into a hybrid FVFD model.This brings several benefits over alternative methods used to simulate wave breaking using an FVFD scheme or a roller approach.The combined effect of complex breaking and Boussinesq terms present in previous models has demanded significant levels of filtering which has introduced undesirable numerical dissipation into the scheme.This masks the physical behaviour of the model and limits the hydrodynamic impact of the breaking terms.The added stability provided by the FV solver removes the need to filter d and q at each timestep, allowing a more physically realistic breaking mechanism and providing a robust model that is able to reliably simulate long wave sequences.
Previous models utilising an FVFD scheme have removed higher order Boussinesq terms in order to simulate breaking, causing the wave shape to become bore-like.Instead, the present model retains these terms throughout the breaking process, thereby establishing a more physical description of the wave and allowing a more realistic surface profile to be obtained.
Velocity data provided by the model when using a roller method is particularly valuable in the surfzone.This information reveals additional details about the hydrodynamic processes occurring in this area which significantly impact morphodynamic changes in nearshore regions.Depth averaged equations, such as those used in the present model, reduce the complexity of a numerical scheme under the assumption of irrotationality.By injecting vorticity, turbulent motion can still be represented within the flow while maintaining a reasonable level of computational efficiency.
The execution time of simulations is of particular interest when comparing models of this type.While the calculation of breaking terms requires some additional computational work, this demand is minimised by updating these terms only at the end of the predictor stage.Breaking initiation and termination, and the roller geometry are all determined from the spatial properties of the wave at the current time step.The breaking criteria for breaking with NSWEs is established from historic wave profiles, resulting in higher memory requirements and greater computational needs.Consequently, the roller approach achieves comparable computational performance to NSWE breaking.
The results presented here demonstrate the model to be capable of producing results in good agreement with laboratory experiments.Several areas of further development are suggested by discrepancies between the numerical data and physical observations.Extending the model to include higher order Boussinesq terms would provide a fully-nonlinear scheme with improved shoaling capabilities.The addition of a bottom boundary layer model could also improve accuracy.As retaining accurate wave shapes during breaking ensures the model is suited to wave reforming, more complex bathymetries could also be used.Tests involving barred beaches can result in breaking termination and reinitiation, and would therefore be an effective way to evaluate this aspect of the model.

Figure 1 .
Figure 1.Injection of vorticity along the lower edge of the roller region.

Figure 2 .
Figure 2. Experimental setup for model validation tests.
Figure 3.Comparison of wave height predicted by model simulating breaking using the rotational roller approach ( ) and by switching to NSWE ( ), and experimental results of Hansen and Svendsen (1979) ( ).