The Velocity Field Underneath a Breaking Rogue Wave: Laboratory Experiments Versus Numerical Simulations

Wave breaking is the most characteristic feature of the ocean surface. Physical investigations (in the field and at laboratory scale) and numerical simulations have studied the driving mechanisms that lead to wave breaking and its effects on hydrodynamic loads on marine structures. Despite computational advances, accurate numerical simulations of the complex breaking process remain challenging. Validation of numerical codes is routinely performed against experimental observations of the surface elevation. However, it is still uncertain whether simulations can accurately reproduce the velocity field under breaking waves due to the lack of ad-hoc measurements. In the present work, the velocity field recorded with a Particle Image Velocimetry method during experiments conducted in a unidirectional wave tank is directly compared to the results of a corresponding numerical simulation performed with a Navier–Stokes (NS) solver. It is found that simulations underpredict the velocity close to the wave crest compared to measurements. Higher resolutions seem necessary in order to capture the most relevant details of the flow.


Introduction
Extreme hydrodynamic forces associated with rogue waves can jeopardise the safety of marine activities [1]. In fact, numerous ship accidents recorded at sea have been related to the impact of rogue waves [2][3][4]. A detailed description of the breaking process often associated with rogue waves and the associated hydrodynamic field is required to provide improved design criteria for the maritime industry.
Despite the inherent complexity of wave breaking, the practical engineering interest inspired numerous research activities to provide a more accurate physical description of rogue waves. Model tests are routinely employed by the maritime industry to provide reliable prediction of the response of marine structures to rogue waves [5][6][7][8][9][10][11][12]. In particular, the impact of a breaking wave has been identified as the most hazardous condition for structural survivability [13][14][15].
Experimental techniques have been used extensively in the past to investigate the flow field underneath breaking waves, from which hydrodynamic forces can be calculated [7,11,12,[16][17][18]. In the last twenty years, optical techniques, such as Particle Image Velocimetry (PIV), have allowed the non-intrusive investigation of the velocity field beneath breaking waves at laboratory scale [11,12,[19][20][21][22]. Note, however, that even in the laboratory velocity measurements are sporadic compared to recording of the free surface [12] and reference therein. Data collected during physical experiments prompted engineers to formulate semi-empirical models to predict the velocity field given the measured surface elevation to predict velocity under oceanic wave where measurements are rarely available [19,23,24]. Whereas in most conditions the models provide accurate predictions of the velocity profile, they have been found to under-predict velocities at the crest of a breaking wave [12,21].
In the last decade, the numerical solution of the Navier-Stokes (NS) equation for the water wave problem has become a feasible alternative to physical model tests in flumes and basins for the maritime industry [25][26][27][28]. The major computational challenge for NS solvers is the accurate simulation of the moving interface between air and water and of the topologically complex bubbly flow, for such a large density ratio. The treatment of the free surface differentiates the solvers commonly used to simulate water waves with some solvers adopting the Volume of Fluid (VoF) approach (e.g., [26,27,29]) in which the motion of the interface is expressed by the transport equation for the averaged volume concentration [30] and others employing the Level Set (LS) approach in which the distance from the interface is advected [31][32][33].
All solvers are routinely validated against physical experiments in respect with the surface elevation (this is routinely measured during physical experiments [34]). In this regard, satisfactory agreement between laboratory experiments and numerical simulations can be achieved for grid size of the order of 400 elements per wavelength [35][36][37][38]. It remains unclear if the NS can accurately simulate the kinematic field underneath breaking rogue waves. Dommermuth et al. [39] compare Boundary Element Simulations to Acoustic Doppler Velocimetry measures for a breaking wave generated via dispersive focussing at close distance from the wave maker which prevents the full nonlinear evolution of the wave spectrum commonly observed in ocean conditions (see [12]). The agreement between measurements and simulations is good, but a quantitative metric is missing. Skyner [40] adopted a PIV technique, but also in this case the breaking is obtained for an unstable spectral shape in proximity of the wave maker (a couple of wavelengths from the wave generation). A slight mismatch is found between the breaking time and distance and in the overall breaker profile. The differences between measurements and Boundary Element Simulations are of the order of 3.9 cm/s for maximum crest velocity of 2.3 m/s, making for a remarkable agreement. More recently, [41] compared PIV measurements to NS simulations for a breaking plunging wave in surf zone, but quantitative assessment of the accuracy is lacking and agreement seems less than optimal.
Here a direct comparison of the velocity filed obtained from numerical simulations and experimental measurements is presented. We remark that it is the wave kinematic field that determines the overall wave stresses and forces. Numerical simulations of a breaking rogue wave are performed with a Level-Set Navier-Stokes (LS-NS) solver [35,36]. To optimise the computational effort, the nonlinear wave growth up to the breaking is reproduced by numerically solving the Euler equations with the Higher Order Spectral Method (HOSM), and the solution is computed shortly before the onset of the breaking is used to initialize the NS solver [42]. Corresponding experiments are conducted in the wave flume at The University of Melbourne where PIV measurements provide the velocity field.
The paper is structured as follows. First, the numerical method employed for the simulation of the wave evolution and breaking is described together with the evolution of a breaking rogue wave. Ad-hoc experiments conducted in the wave flume at The University of Melbourne are then presented. The comparison of the velocity field obtained numerically and experimentally is presented in the discussion section together with a comparison with the velocity profile predicted by standard engineering methods. Final remarks are reported in the Conclusions.

Numerical Method
The computational time associated with two-phase Navier-Stokes simulations of water waves is usually high. The treatment of the moving boundary between fluids with very different density (air and water) represents the biggest challenge. Furthermore, small scales (turbulence) and large scales (wave-wave interactions) have to be resolved at same time. These computational challenges often limit the size of the domain that can be studied numerically. To optimise the computational resources two numerical methods (the HOSM and the LS-NS) are coupled [21,35,36,42]. We use a NS in the coupling with the aim of accurately reproducing the energy dissipation due to viscous effects which occur when the large surface curvatures develop [43] and provides a general methodology to investigate post-breaking stages of wave evolution and interactions with structures that we plan to do in future. In the following subsections the numerical methods and their coupling are described in detail.

Higher Order Spectral Method (HOSM)
The HOSM is a numerical method that efficiently solves Euler equations. It accurately reconstructs the nonlinear wave evolution over large domains. The HOSM has been independently proposed by West et al. [44] and Dommermuth and Yue [45], and since then it underwent improvements and validations, e.g., [46,47]. The mathematical problem is fully described by the pair: where ψ is the potential at the surface (ψ = φ| z=η ), W the vertical velocity at the surface W = ∂ψ/∂z and ∇ the spatial derivative in the horizontal direction only, i.e., ∇ = (∂/∂x, ∂/∂y).
The system of Equations (1) and (2) is numerically solved truncating at the nonlinear order M. For realistic random sea state conditions satisfactory results can be achieved with d < λ 0 /16 where λ 0 is the characteristic wave length and the time step must satisfy the Courant-Friedrichs-Lewy (CFL) condition [48] on the phase speed c to guarantee numerical stability.

Level-Set Navier-Stokes (LS-NS)
The NS equations are a set of conservation equations, i.e., for the mass and the momentum, and can be easily derived from the Newton's equations. The momentum conservation for a incompressible flow is: where p is the pressure, T is the deviatoric stress tensor, and f the external forcing.
In the case of water wave breaking the physical domain is a two fluid domain, i.e., air and water. The two-phase (or two-fluids) flow is modelled as that of a single incompressible fluid whose physical properties varies in the domain as a function of the instantaneous position of the interface. In the Level-Set method, the interface is detected as the zero level-set of a distance function from the interface. Assuming that particles lying on a free surface remain on it, the distance function is advected with the flow as a passive scalar [31][32][33]. In the model adopted here, density and viscosity vary smoothly across the interface and surface tension forces are spread across the same region [36].
The non-dimensional NS equations for an incompressible newtonian fluid with variable physical properties are written in curvilinear coordinates and solved by using a classical fractional step approach [25,49]. In the predictor step the momentum equation is advanced in time disregarding the pressure terms, in the corrector step the pressure is reintroduced satisfying the continuity equation. This implies the solution of a Poisson equation whose coefficients involve the local density and thus, for large density ratios, it results in a hill conditioned problem that is responsible for a large portion of the total computational effort. At the initial time the distance is assumed positive in water, negative in air and d = 0 at the interface. The smoothing function is defined as: where δ is the semi-amplitude of the smoothing zone and f the generic property. The semi-amplitude of the smoothing zone is usually 3 ≤ δ ≤ 5 grid cells. In this work, a grid cell of about 2 mm is employed and thus the semi-amplitude of the smoothing zone is about 1 cm, i.e., 5 cells. The smooth density profile, together with the spreading of the surface tension effects across the same region, introduce some spurious velocity and vorticity in the layer. A careful study on this regard was conducted in Iafrati and Campana [50] where comparisons with results obtained by single-fluid, sharp interface, models available in literature were established and it was found that, within some restrictions in the grid resolution, the spurious effects remain confined in the transition layer and do not affect the fluid dynamics outside. As anticipated, the distance function d is advected with by the flow via the equation: The above equation is solved with a three-step, low-storage, Runge-Kutta scheme. Since Equation (5) is strictly valid only at the interface, the distance from the interface has to be reinitialized [32,49]. Once the distance from the interface is reinitialized, the density and viscosity distributions are computed by equation Equation (4). This numerical method has been developed and thoroughly validated [35,36].

Coupling the HOSM and the LS-NS
Starting from an initial condition the wave evolution is simulated with the HOSM up to one period prior the breaking onset (identified as ka = 0.35 , [42,51]). We note that only the breaking onset injects vorticity in the flow field [35] and potential flow approximation accurately predicts the surface elevation even in presence of strong vorticity prior the breaking onset [52]. Surface elevation and velocity potential obtained from the HOSM at this instant are used to initialise the LS-NS that only simulates the final stages of the breaking (see Figure 1). The coupling allows to significantly reduce the computational time associated with LS-NS simulations [21]. The HOSM fully simulates the nonlinear evolution occurring over tens of wave periods in less than a 1 s CPU time. On the other hand, the LS-NS runs at 1 period of wave simulation for 1 day CPU time on the same machine, i.e., 10 5 times slower than the HOSM. The procedure is summarised here [36]: 1. The HOSM surface elevation and potential at the final instant (η(x) and ψ(x)) are interpolated over a grid compatible with the resolution needed by the LS-NS simulation. This is about 400 elements per wavelength (L 0 ) [50]. A cubic spline with periodic conditions is used to resample the surface elevation and its potential. 2. Surface elevation and velocity potential are only defined at the surface. From the boundary value problem, defined on the contour, the potential at any point in the water domain is reconstructed by using boundary integral representation. The reconstruction is done numerically. The computational domain is discretised with square cells, i.e., dimension in the vertical axis is equal to the horizontal dimension, for −0.3 < z/L 0 < 0.3. Out of this layer the vertical dimension is stretched logarithmically, i.e., rectangular cells (see Figure 2). A no slip boundary condition is used to limit the computational domain in −1 < z/L 0 < 1. 3. To obtain the velocity components in the water domain the spatial derivatives of the potential are computed using a second order finite difference scheme. The value at the centre of the cell, x C , is computed using the four points at x D = x C ± ∆/5. ∆ denotes the dimension of the cell. 4. The boundary value problem is formulated also for the air domain. The boundary of the air domain are the air-water interface (at the bottom), and the arbitrarily chosen limit of the computational domain at the top (indicated in blue in Figure 2). At the free surface the continuity on the normal velocity is enforced providing the necessary Neumann condition for the solution of the boundary value problem. The velocity field in the air is then computed adopting the same methodology used for the water side (point 2. and 3.). 5. The signed distance function d is initialised. This is needed for the initialisation of the density distribution in the Level-Set technique. The distance is measured from each grid node to the closest segment defined by two successive points on the resampled free surface. The distance is positive if the grid node is in water, negative in the air side. The Point in Polygon algorithm is used to identify if a grid node belongs to the air or the water side of the domain. 6. At the interface between the two fluids (air and water), the velocity component tangential to the free surface has a jump. The discontinuity is consistent with the potential flow theory (no stress condition at the interface) but it does not apply to the two fluid approach. The jump is smoothed by gradually varying the fluid properties (density and viscosity) across the interface. Spurious component that are generated in this layer (∼5 cells, see Figure 2) because of the baroclinic variation of the velocity and of the spreading of the surface tension forces are not propagated in the rest of the domain because the density is reinitialised at each step.

Set-Up and Initial Conditions
In the ocean, modulational instability is the main nonlinear mechanism leading to wave growth and, eventually, breaking [21,51]. This mechanism can be exemplified by the interaction between a monochromatic wave (called carrier) and infinitesimal sidebands perturbations [53]. Energy is nonlinearly exchanged between the carrier and the sidebands resulting in the growth of the perturbations [54,55]. In the physical domain, the energy exchange leads to the formation of a rogue wave from an initially monochromatic wave train [53]. The time and space needed for the instability to develop can be controlled by tuning the frequency and amplitude of the sideband perturbations [56].
In particular, the nonlinear energy exchange is triggered if the Benjamin Feir Instability Index (BFI), an index of nonlinearity, is sufficiently high, i.e., BFI > 1 [35,53,56]. Under the hypothesis of deep water conditions, the index is computed as [57,58]: Here k 0 is the wavenumber of the carrier wave, ω 0 the corresponding radial frequency, ∆k the wavenumber difference between carrier and sideband, ∆ω the radial frequency difference between carrier and sideband and A 0 the amplitude sum of carrier (a) and sidebands (b ± ), i.e., The breaking rogue wave generated via modulational instability is studied numerically. In the simulation, the bandwidth is set to ∆k/k 0 = 1/5 (equivalent to ∆ω/ω 0 = 1/10). The wave group is periodic over 5 waves in the spatial domain and 10 periods in the time domain. The wave period of the carrier is set to T 0 = 0.8 s (ω 0 = 2π/T 0 ), the wave steepness k 0 A 0 = 0.1474 and the perturbation to carrier amplitude ratio b ± /a = 0.01. The resulting BFI is 2.08.
The initial surface elevation and the corresponding velocity potential at the surface given to the HOSM are defined respectively as: ψ(x, t = 0) = a g/k 0 exp(k 0 η) · sin[k 0 x] For the LS-NS density of the air and water are respectively ρ a = 1.25 kg/m 3 and ρ w = 1000 kg/m 3 . Surface tension for pure water is σ = 0.0073 N/m and the dynamical viscosity for air and water are µ a = 1.8 × 10 −5 Pa·s and µ w = 1.0 × 10 −3 Pa·s. We recall that the spatial resolution is over 400 grid points per wavelength (2048 grid points per 5 wavelengths), denser than comparable NS simulations of breaking waves used for engineering purposes [35][36][37][38].
We note that both the HOSM and the LS-NS used here adopt periodic boundary conditions (see Figure 2). This is made possible by the periodicity of the wave group itself. The equivalence between evolution in space and time has been commonly adopted (and proved) in the past [59]. We note that formally the periodic solution does not correspond to the wave evolution in a physical wave flume and, despite being a limitation of the present study, allows to make NS simulations considerably faster.

Wave Evolution
The evolution of the maximum amplitude in the wave group (A) as obtained from the HOSM is shown in Figure 3. At breaking, which occurs at about 1200 s from the HOSM initialisation and local wave steepness of 0.35, the envelope has grown by about two times compared to the initial monochromatic wave. The left sideband (b − ) has outgrown the carrier wave (a), see Figure 3. After the switch from the HOSM to the NS, the evolution of the breaking wave is shown in Figure 4. At the breaking, a jet develops at the front of the largest wave in the group (the one shown in the plot). As expected, the largest horizontal velocities are concentrated at the top of the crest where the overturning movement is taking place (see Figure 5). The figure also highlights asymmetries in the wave shape and the associated velocity field. It is worth noticing that the velocity in air is not shown for the sake of clarity. In previous laboratory experiments, the shape of the breaking wave has been found to be affected by the presence of surfactant that might modify the surface tensions [60]. For this reason a sensitivity analysis in respect of viscosity and surface tension is performed. To exaggerate the effect variation of this parameters have on the solution, viscosity and surface tension are varied singularly by ±30%. Results are summarised in Figure 6. Figure 6. The effect of viscosity and surface tension on the velocity field under a breaking wave, the difference between horizontal velocity computed with physical parameters and modified parameters is shown. Top left surface tension is decreased by 30% (s −30% ), top right surface tension is increased by 30% (s +30% ), bottom left air and water viscosity are decreased by 30% (v −30% ), bottom right air and water viscosity are increased by 30% (v +30% ). The red line shows the profile obtained using "physical" parameters.
Viscosity and surface tension can affect velocities up to ±0.1 m/s (about 10% of the maximum crest velocity). Viscosity has a much stronger effect on both the jet shape and the flow field. Note, however, that in all the cases the wave profile, with the exception of the jet, remains unaffected by variation of the physical parameters.

Set-Up and Initial Conditions
A physical experiment is conducted in the Extreme Air Sea Interaction facility (EASI) using the same set-up described in [11,12,22]. The flume is 60 m long, 2 m wide and filled up to a depth of 0.9 m. At one end waves are generated by a cylindrical wave maker and the opposite end a linear beach absorbs the incoming wave energy. The surface elevation is monitored with a number of wave probes deployed along the tank. The Particle Image Velocimetry (PIV) setup for the fluid velocity measurements is located 34 m from the wave maker in correspondence of a glass window that guarantees optical access. The field of view is approximately 170 mm × 200 mm with a resolution of ≈ 0.1 mm/px [11,12]. The experiment is performed to replicate the breaking event generated numerically. However, the notable difference is that the numerical simulation is performed with periodic boundary conditions, on the other hand experiments are conducted in a wave flume. The imposed surface at the wave-maker is: To reproduce the numerical simulation the same dominant wave period is used (T 0 = 0.8 s). Due to the limited length of the wave flume, the initial condition in the flume starts from a more advanced stage of evolution compared to the one used during the numerical simulation. In particular, amplitudes of the sidebands and initial steepness need to be tuned to obtain the wave breaking in the PIV field of view [11,12,22]. The initial sideband to carrier amplitude ratio is set b ± /a = 0.10 (cf. 0.01 in the numerical simulation) and the corresponding initial steepness is k 0 A = 0.1734.
PIV processing is performed using the software PIVLAb [61]. Images pairs are sampled every ∆ t = 2.5 ms. To obtain robust statistical results the test is repeated 20 times.

Wave Evolution
The evolution of the maximum amplitude at different distances from the wave maker is shown in Figure 7. The partitioning of the energy between carrier and sidebands is also shown. Starting from a carrier (a) more energetic than the sidebands close to the wave maker, modulational instability leads to the growth of the left side-band (b − ) that eventually outgrows the carrier consistently with the HOSM simulation. The shape of the wave train is shown in the insets, an extreme wave forms from an initially slightly modulated wave train. At the breaking (recorded at the fifth probe) a sudden drop of energy occurs. Note that energy exchanges between wave components continue after the breaking.

Discussion
The comparison between the simulations and the experiments focusses on the horizontal velocity field under the breaking event. The major difference between physical experiments and the numerical simulations is the shape of the breaking wave, see Figure 8. In the numerical simulation (red line in Figure 8) a jet departs from the tip of the crest. In the experiment (black particles in Figure 8) the jet is far less pronounced. The surface elevation of the wave front is reproduced correctly by the simulations.
It is important to note that the interface in the Level-Set method is identified by ρ = 500 kg/m 3 (this value is the average between the density of pure air and pure water). For this reason the jet identified numerically is a mixture of air and water, i.e., density in the jet are below the density of pure water. To better resolve the foamy jet a much smaller grid size would be required making the simulations unfeasible for practical engineering purposes. However, a finer grid is prone to numerical artefacts and would increase the computational time. To account for the variable density across the interface, velocities are corrected multiplying by the density ratio ρ/ρ W . Only data points where the density ratio is higher than 0.5 are considered as belonging to the water side.  Figure 8 noting that the velocity in the crest region is highly variable [21] and PIV uncertainty up to 33% [11,12]), in the simulation velocities at the air-water interface are smaller. It must be recalled that in the LS the fluid properties jump and the surface tension forces are spread across a a finite region about he interface, which affects the accuracy of the velocity in this transition zone.
The velocity obtained numerically and experimentally under the wave crest are compared in Figure 9. Overall, experiments provide larger velocities than those obtained numerically. The difference is about 10% (dashed line in Figure 9) in the wave body where both simulations and measurements are available which is not statistically significant due to measurement uncertainties [11,12]. The inset in Figure 9 highlights how the largest velocity differences are localised in the top of the crest and the front face of the wave. Close to the interface the numerical simulation fails to accurately reproduce the velocity field due to the interface treatment but also PIV measurements become challenging [11,12].
The velocity profile beneath the wave crest is extracted from simulations and experiments. As a reference, in Figure 10, the velocity profile is plotted against engineering methods developed to obtain the velocity only from the surface elevation time history. The Grue method, for which the velocity profile can be approximated by a third order monochromatic Stokes wave with the same amplitude and period [19], and Donelan method, in which the contribution of each spectral component is iteratively added to the perturbed solution [23], are compared. These two methods have been shown to provide good agreement with experimental results [12]. The velocity profiles highlight that experiments provide larger values compared to the simulations, the difference is less than 0.1 m/s. PIV reconstructs the velocity profile up to elevations very close to the top of the wave, however reflections at the interface make the image processing challenging [11,12]. On the other hand, the LS-NS accurately identifies the surface but the velocity and density close to the interface are uncertain due to the smoothing of properties across the interface itself. For this reasons the velocity is only reported up to about 5 cm (the maximum wave elevation is 6 cm). The Grue and Donelan method used to estimate the velocity profile from the measurements in this case provide an upper and lower limit for the numerical and experimental measurements. At the tip of the crest the Grue and Donelan methods predict the same maximum velocity. In this case, the Grue method accurately predicts the experimental data (cf. [12] and [21] in which the agreement is far from optimal).

Conclusions
A breaking wave induced by modulational instability has been generated both numerically and experimentally to evaluate the numerical solver results in regard of the induced fluid motion. To speed-up the simulations, Navier-Stokes (NS) are initialised from a pre-breaking solution obtained using the Higher Order Spectral Method (HOSM). The coupling of a "fast" method (HOSM) with a "slow" method (NS) allows one to reduce the computational time and, at the same time, accurately reproduce the breaking stage starting from a fully nonlinear wave condition that better reproduces ocean conditions. Whereas the HOSM only takes few seconds to simulate the entire evolution, the present NS solver is not particularly efficient and takes about a week of single processor CPU time to simulate one wave period. The general methodology presented in the manuscript can be used when the NS are coupled with the impact on a structure or a range of engineering problems that require the post-breaking dynamics (which we plan to study in future).
Comparison of the numerical simulation with the physical experiments highlights that the wave evolution is comparable and the shape of the breaking wave is accurately reproduced, despite having used a periodic support for the numerical simulations (to reduce the domain size and speed up the numerical simulation). A notable exception is the jet, that is absent in physical experiments. It is worth noting that the shape of the jet itself has been found to highly depend on the choice of the physical parameters, i.e., surface tension and viscosity.
Overall, the topology of the flow field recorded in the experiments is in reasonable agreement with the one obtained from numerical simulations. However, significant differences are observed in the velocity close to the crest (about 10%). Despite having ensured the equivalence between spatial and time evolution (in the wave tank and in the numerical simulation respectively), it cannot be ruled out that the periodic boundary condition might have been responsible for the velocity difference between experiments and simulations together with the relatively low spatial resolution (e.g., De Vita et al. [54] uses 5 times more elements per wavelength) and to the spreading of the fluid properties. In this regard, we note that in the LS the interface is correctly identified even at relatively low resolution, but the velocity at the interface is more sensitive to the grid size. Attempts to overcome limitations highlighted in this work are currently under development. In particular, it is advisable a comparison with simulations based on a numerical wave tank configuration to verify if the periodic support is a viable approximation for engineering purposes.
Application of engineering method to predict the velocity profile beneath the crest shows that, in this case, the Grue method accurately reproduces the velocity profile measured during the laboratory experiments, while the Donelan method accurately reproduces the LS velocity profile up to the mean water level.
Author Contributions: All authors contributed to conceptualization, methodology, investigation and writing of the manuscript.
Funding: This research received no external funding