[
Abstract
The kinematic properties of unsteady highly nonlinear 3D wave groups have been investigated using a numerical wave tank. Although carrier wave speeds based on zerocrossing analysis remain within of linear theory predictions, crests and troughs locally undertake a systematic cyclical leaning from forward to backward as the crests/troughs transition through their maximum amplitude. Consequently, both crests and troughs slow down by approximately of the linear velocity, in sharp contrast to the predictions of finite amplitude Stokes steady wavetrain theory. Velocity profiles under the crest maximum have been investigated and surface values in excess of 1.8 times the equivalent Stokes velocity can be observed. Equipartitioning between depthintegrated kinetic and potential energy holds globally on the scale of the wave group. However, equipartitioning does not occur at crests and troughs (even for low amplitude Stokes waves), where the local ratio of potential to total energy varies systemically as a function of wave steepness about a mean value of 0.67.
Local properties of highly nonlinear unsteady waves]On the local properties of highly nonlinear unsteady gravity water waves. Part 1. Slowdown, kinematics and energetics X. Barthelemy, M.L. Banner, W.L. Peirson, F. Dias, M. Allis]X. Barthelemy ^{†}^{†}thanks: Email address for correspondence: x.barthelemy@unsw.edu.au,M.L. Banner, W.L. Peirson,F. Dias and M. Allis
1 Introduction
Windgenerated ocean waves at the airsea interface propagate predominantly in group or packet structures. Group structure behaviour is not limited to ocean waves, but occurs naturally in many systems. These unsteady wave groups can exhibit a complex nonlinear life cycle, especially in focal zones where there is rapid wave energy concentration. Their structure and the associated carrier wave propagation properties of dispersion, directionality and nonlinearity, and their interplay, represents a knowledge gap beyond present analytical treatment. Understanding natural ocean wave propagation is not only relevant academically but is potentially important when studying oceanatmosphere exchanges, as well as the design of structures such as open ocean platforms. Crest speeds are a key consideration in the assessment of wave impact loading on structures, where flowinduced forces are proportional to the square of the water velocity. In other applications, measurements of initial whitecap speeds have been proposed as a method for inferring ocean wavelengths (Phillips (1985)) assuming that the Stokes dispersion relationship is applicable. Any systematic crest speed slowdown just prior to breaking onset will significantly alias wavelength determinations from the quadratic relationship between Stokes wavelength and speed.
Stokes classical deep water wave theory ((Stokes (1847)) was developed for a steady, uniform train of twodimensional (2D) nonlinear, deepwater waves of smalltointermediate mean steepness , where is the wave amplitude and is the wavelength. The wave speed c increases with as found to 5th order in wave steepness by Fenton (1985):
(1.0) 
where is the wave speed of linear (infinitesimally steep) waves.
Increasing the steepness to the Stokes limit () in eq (1) gives the limiting speed of maximallysteep steady waves of . Increased wave steepness has long been associated with higher wave speeds. Historically, Stokes theory has been the default theoretical approach for describing ocean waves (Kinsman (1965) and Wiegel (1964)). Contemporary examples are Beyá et al. (2012) who found close agreement between monochromatic wave nearsurface velocities and predictions based on 5th order Stokes theory. The observational study of Grue et al. (2003) using PIV to measure velocities under extreme and breaking crests showed that inviscid predictions are in good agreement with measured particle velocities. Johannessen (2010) showed that at the crest, nearsurface velocities agree with secondorder Stokes theory, which is capable of describing the free surface kinematics very accurately provided that the local underlying regime of free waves can be identified. A principal constraint of the Stokes wavetrain framework is the imposed uniform spatial and temporal periodicity, which provides the basis for the extensive analytic treatment in the literature of this intrinsically nonlinear free surface problem.
However, ocean waves propagate naturally in unsteady groups that evolve dynamically. Relaxing the periodicity constraints allows additional degrees of freedom to develop as shown in a suite of laboratory studies. Melville (1983) found a wave speed variation of between 17% and +32% around the linear phase velocity. Shemer (2013) demonstrated that for a wave group with a broad spectral bandwidth, the crest propagation velocity of the dominant waves differs significantly from both their phase and group velocities. Deepwater breaking wave studies (e.g. Rapp & Melville (1990), Stansell & MacFarlane (2002), Jessup & Phadnis (2005), Melville & Matusov (2002), Melville et al. (2002)) all found a significant (O(20%)) breaking crest speed slowdown relative to the expected linear phase velocity (Eq. 1). For focusing wave packets, Johannessen & Swan (2001) and Johannessen & Swan (2003) reported a slowdown in the expected crest velocity at focus of around 10%.
Understanding wave crest slowdown behaviour is central to both refining present knowledge on waterwave propagation and dynamics, and effective implementation of Phillips’ spectral framework for breaking waves (Phillips (1985), Gemmrich et al. (2008), Kleiss & Melville (2010)). The groupiness of natural ocean waves gives rise to a suite of differences from Stokes waves predictions. Baldock et al. (1996) report an ensemble of directional wave modes of different frequency components focusing in space and time to produce an unsteady wavegroup in the laboratory. Surface elevations and subsurface particle kinematics were compared with linear wave theory and the LonguetHiggins & Stewart (1960) secondorder solution of the wavewave interactions. Their study shows that nonlinear wavewave interactions produce a highly nonlinear wavegroup with crests higher and troughs shallower than expected from numerical simulations by LonguetHiggins (1987) and experiments by Miller et al. (1991). Sutherland et al. (1995) previously obtained very similar characterisations although they also concluded that the wave kinematics were Stokeslike. Song & Banner (2002), Banner & Peirson (2007) and Viotti et al. (2014) report complementary numerical and experimental studies that demonstrate that intragroup and interwave energy transfers cannot be neglected when characterising wave group kinematics and dynamics.
In this contribution we show that Stokes characterisations have significant limitations when applied to fully nonlinear unsteady wave groups. Specifically, important properties derived from Stokes (1847) theory and subsequent refinements (e.g. Fenton (1985)) change fundamentally when stationarity is relaxed. Under wave crest maxima, the velocity profiles, energetics and energy partitioning all differ systematically from those of a comparably steep Stokes wave. The unsteady leaning of waves (Tayfun (1986)) is shown to be a critical aspect, determining the actual wave crest speed. Our preliminary numerical studies of chirped wave packets have already been verified against observations of wave groups in the laboratory and in the field (Banner et al. (2014)). Fedele (2014) provided a theoretical explanation and further validation of the crest slowdown in terms of geometric phases.
2 Numerical simulation
For modulated waves, the nonlinear Schrödinger (NLS) equation and its higherorder extensions (Zakharov (1968)) have been widely applied. NLS models have been used extensively to study the evolution of quasimonochromatic waves and instabilities due to resonant interactions (e.g. the Benjamin–Feir instability (Benjamin & Feir (1967)). Extensions of the NLS equation include the Zakharov equation which has been applied to describe the longtime evolution of the spectrum of weakly nonlinear, dispersive waves.
There has been a growing interest in the development of threedimensional models which inherently incorporate nonlinearity and associated dispersion effects. The broadbandedness in both frequency and direction of real sea states poses significant challenges in numerical simulation. Bateman et al. (2001) demonstrated the importance of directionality and the consequent benefits of efficient wave modelling during comparison of numerical simulations with the laboratory observations of Johannessen & Swan (2001). Highorder spectral expansion approaches using efficient FFT solvers for application to 3D waves have been developed (e.g. Ducrozet et al. (2012), Dias et al. (2015)). Another option is to solve the full NavierStokes equations (Park et al. (2003)), but viscous flow solvers tend to be too dissipative and computationally timeconsuming.
Numerical models of potential 3D wave propagation can be divided into three main categories: (a) boundary element integral methods (BEM): e.g. Baker et al. (1982), Bateman et al. (2001), Clamond & Grue (2001), Grilli et al. (2001), Fochesato et al. (2007), Guyenne & Grilli (2006), Xue et al. (2001), Hou & Zhang (2002), Fructus et al. (2005); (b) finite element method (FEM) e.g. Ma et al. (2001); (c) spectral methods: e.g. Dommermuth & Yue (1987), West et al. (1987), Craig & Sulem (1993), Nicholls (1998), Bateman et al. (2001).
Spectral methods based on perturbation expansions are known to be very efficient. These methods reduce the water wave problem from one posed inside the entire fluid domain to one posed on the boundary alone, thus reducing the dimension of the formulation. This reduction can be accomplished by using integral equations over the boundary of the domain (socalled boundary integral methods) or by introducing boundary quantities which can be expanded as Taylor series for reference domain geometries (Xu & Guyenne (2009)). Both approaches have been summarised recently by Ma (2010). BEM techniques are efficient for representing wave propagation and overturning until the wave surface reconnects (Grilli & Subramanya (1996)).
The present study used a boundary element numerical wave tank code called WSIM, which is a 3D extension of the 2D code developed by Grilli et al. (1989) to solve the singlephase wave motions of a perfect fluid. It has been applied extensively to the solution of finite amplitude wave propagation and wave breaking problems (see chapter 3 of Ma (2010)). The ideal fluid assumption makes WSIM unable to simulate breaking impact subsequent to surface reconnection. However, its potential theory formulation enables it to simulate wave propagation in a CPUefficient way, without the diffusion issues of viscous numerical codes. The simulation of wave generation and development of the onset of breaking events can be carried out with great precision, as demonstrated by Fochesato & Dias (2006), Fochesato et al. (2007).
WSIM has been validated extensively and shows excellent energy conservation (Grilli et al. (1989), Grilli & Svendsen (1990), Grilli & Subramanya (1994), Grilli & Subramanya (1996), Grilli & Horrillo (1997), Grilli et al. (2001), Fochesato (2004), Fochesato & Dias (2006), Fochesato et al. (2007)). Its kinematical accuracy has been demonstrated by comparison with the analytical solutions for infinitesimal sine waves found in Phillips (1977).
2.1 Governing equations
WSIM uses a mixed EulerianLagrangian timeupdating scheme for irrotational motion described by the velocity potential , in a Cartesian coordinate system with constant pressure at the open water surface. is the vertical upward direction and the still water surface (Figure 1).
Fluid velocity is defined as . The continuity equation leads to the Laplace equation for the potential within the fluid domain (2.1). The symbols and are used to denote the entire domain boundary and the free surface respectively.
(2.0)  
(2.0)  
(2.0)  
(2.0) 
where is the position vector of a fluid particle on the free surface, the gravitational acceleration, the atmospheric pressure, the fluid density and the Lagrangian (or material) time derivative.
The free surface (domain ) is described by fullynonlinear kinematic (eq 2.1) and dynamic (2.1) equations. Recent developments have implemented a 3D snake wave paddle at one vertical face of the domain which is described subsequently. The far face of the numerical wave tank from the paddle has an absorbing beach which damps any incident wave energy as described in Grilli & Horrillo (1997). All remaining faces of the domain have a zeroflux boundary condition () (eq 2.1).
2.2 Discretisation
Boundary geometry and field variables are represented by 16node quadrilateral elements, providing bicubic local interpolation of the solution between the nodes. Highorder tangential derivatives required for the time updating are calculated in a local curvilinear coordinate system, using 25node fourthorder quadrilateral elements to retain third order accuracy in the spatial derivative. This results in a highorder integration scheme over the discretisation elements. The solver uses a fast PTMA scheme introduced by Fochesato & Dias (2006).
Once preliminary testing of WSIM was complete, two spatial resolutions were used during the present simulations to quantify sensitivity to the discretisation. These were 8 points per wavelength in (socalled) medium resolution and 16 points per wavelength in high resolution in the propogation () direction. The number of points in the transverse and depth directions was adjusted to keep element aspect ratios close to 1.
Even if WSIM internally works with a different reference system, the results in this paper have been chosen to be expressed in deep water units. Model time scales are nondimensionalised by the base frequency of the generating paddle signal, and its corresponding length scale is the deep water wavelength . These choices impose a gravitational acceleration of and a reference linear velocity .
Nondimensional lengths can be computed using and the time related quantities using .
The numerical scheme is a mixed EulerianLagrangian method. Consequently, the mesh is subject to Lagrangian drift during the passage of steady, steep wave groups, especially near the freesurface. For the transitory Class 3 wave groups defined in Song & Banner (2002) and used during this present investigation, the mesh did not deform significantly.
2.3 Wave group generation
The WSIM wavemaker was chosen as a bottomhinged flap paddle to generate deepwater waves.
The wavemaker motion requires specification of the velocity potential on the boundary as follows:
(2.0)  
(2.0) 
where is the distance of the centre of the rotation and is the angle of rotation of the wavemaker.
To ensure smoothness of the boundary conditions and to avoid numerical instabilities at created by the moving elements, the paddle oscillation motion is damped by a function which changes smoothly from to , as detailed in Grilli & Subramanya (1996) and Grilli & Horrillo (1997). Hence, the boundary condition can be written as follows:
(2.0)  
(2.0) 
Because of finite tanklength considerations, we chose to investigate the chirped wave packet (Class 3) and investigated cases comprising five, seven and nine carrier waves in the initial time wave packet.
2.3.1 Class 3 wave groups
Using the prescribed ’chirp’ motion from Song & Banner (2002), the following expression defines Class 3 wave packets, being the number of waves in the time series, is the paddle frequency and is the phase. is the chirp rate of the linear modulation:
The paddle motion has an oscillating part moderated by a damping function. To avoid discontinuities and assure a smooth numerical start in the code, there is a negative shift in time of 5 seconds, to ensure that any initial oscillations are of the order of the double precision zero.
2.3.2 3D converging wave group
If 3D waves are to be produced by the snake paddle, the phase becomes a function of which specify the linear convergence point coordinates, as defined by Dalrymple & Kirby (1988); Dalrymple (1989):
3 Postprocessing
3.1 Determining nearsurface quantities:
A significant challenge within this investigation was determining quantities close to the highlycurved free surface. Three complementary solutions have been implemented. Evaluating the interior fields is known to be a ’quasisingular’ problem, because it evaluates integration with kernels of the type , being the distance between the point where the field is calculated and a point on the boundary and is a positive exponent. The closer the evaluation points are to the boundary, the larger these terms become and they can possibly grow without bound.
3.2 Interior fields
The boundary element formulation of 3D potential problems implies:
where is the source point, is the potential; and is the derivative of along the unit outward normal at on the boundary . is the boundary of the region of interest and boundary conditions concerning and are specified on . when and when and is smooth at
The fundamental solutions and are defined by:
where and .
The flux at a point is given by the potential gradient:
(3.0) 
where
Equations are discretised on the boundary by boundary elements defined by the interpolation function. The integral kernels of these equations become nearly singular when the distance between and is small compared to the size of . In the following, we will denote by for brevity.
Reconstructing the inner field values relies on the ability to correctly estimate the sum (3.0) over the entire domain. Three complementary methods are used depending on the distance between and .
When the distance is large enough (relative to the element size), the integration on the individual element is carried out by a Riemann quadrature on a classical GaussLobatto point distribution.
When the distance becomes small, the classical quadrature does not retain enough precision, and the Telles (1987) method is used. The Telles method consists of binary subdivisions of the integration space to maximize precision where the singularity begins to manifest itself. The precision of this method is acceptable at moderate distance from the boundary, as described in Grilli & Subramanya (1994).
However, the Telles method becomes inefficient when the nearboundary singularities are too strong, specifically adjacent to the highly curved surface of steep waves. Consequently, a third method was developed and implemented called Projection and Angular and Radial Transformation (PART) (see Hayami (1990), Hayami (1991), Hayami & Matsumoto (1994),(Hayami (2005a), Hayami (2005b)). A change of variables is used to calculate the integration on a closed element.
When applied to the curved quadrilateral element , for the nearlysingular integral over a curved boundary element:
(3.0) 
A sequence of steps is needed in order to effect a change of variables and integrate:

1. Source projection: Find the nearest point on to using the NewtonRaphson method. are the local orthonormal coordinate system on the element .

2. Determine the relative position of the source projection in the local coordinate system.

3. Calculate the unit normal to the curved element at

4. Approximately project the curved element on the polygon in the plane tangent to at

5. Determine the geometry of the projected quadrilateral and four triangular regions defined by joining to the four corners of , and then determine the linear mapping matrices .

6. Introduce polar coordinates in , centered at . Eq. (3.0) becomes:
(3.0) 
where is the Jacobian of the mapping from Cartesian coordinates on to curvilinear coordinates on .

7. Apply radial variable transformation: in order to weaken the nearsingularity due to , which is essentially related to the radial variable only:

8. Apply the angular variable transformation , with being for each triangle the distance between and the foot of the perpendicular from and the edge of and the angle at of the triangle :
(3.0) 
in order to weaken the angular nearsingularity which arises from when is near the edge of the polygon . This uses the fact that:
(3.0) 

9. Use GaussLegendre’s formula to perform numerical integration of (3.0) in the transformed variables and :
(3.0) 
3.3 Determination of physical quantities
Local surface elevation and its corresponding depthintegrated potential energy are obtained from the free surface data, extracted from the data using the same 3rd order spline interpolation polynomials as used in the BEM simulation (section 2.2). Velocities and the associated depthintegrated kinetic energy are calculated using free surface data and the interior fields as seen above (section 3.2).
The crest and trough locations are tracked using a twostage detection algorithm:
First, extrema and the zerocrossing positions are calculated directly from the discretising spline polynomials. These positions are then linked together between time steps to form characteristics of the motion. All physical quantities are tracked along these characteristic lines.
4 Results and Discussion
The numerical simulations undertaken are summarised in Table 1.
Identifier  Class  Number of Waves  Approx. points per wave  Paddle amplitude A  /  Comments  
C3N5A0.05  3  5  8  0.05  2D  1.1014  1.2878  Small steepness case 
C3N5A0.3  3  5  8  0.3  2D  1.1171  1.22878  Medium steepness case 
C3N5A0.508  3  5  16  0.508  2D  1.2382  1.2124  
C3N5A0.511  3  5  16  0.511  2D  1.2499  1.2028  Marginal high resolution recurrence case 
C3N5A0.513  3  5  8  0.513  2D  1.2024  1.2156  
C3N5A0.514  3  5  8/16  0.514  2D  1.2416  1.1978  Marginal high resolution breaking case 
C3N5A0.516  3  5  8/16  0.516  2D  1.2380  1.2022  
C3N5A0.518  3  5  8/16  0.518  2D  1.2126  1.2023  Marginal medium resolution recurrence case 
C3N5A0.519  3  5  8/16  0.519  2D  1.2275  1.1882  Marginal medium resolution breaking case 
C3N5A0.53  3  5  8  0.53  2D  Breaking onset too close to paddle to obtain reliable values for and  
C3N5A0.56  3  5  8  0.56  2D  Breaking onset too close to paddle to obtain reliable values for and  
C3N5A0.32X10  3  5  8  0.32  10  1.3940  1.1901  
C3N5A0.33X10  3  5  8  0.33  10  1.2538  1.2290  Marginal recurrence case 
C3N5A0.34X10  3  5  8  0.34  10  1.2472  1.2250  Marginal breaking case 
C3N5A0.35X10  3  5  8  0.35  10  1.2568  1.2202  
C3N5A0.36X10  3  5  8  0.36  10  1.2735  1.2168  
C3N7A0.41  3  7  8  0.41  2D  1.2925  1.3141  
C3N7A0.42  3  7  8  0.42  2D  1.2974  1.3155  
C3N7A0.43  3  7  8  0.43  2D  1.2964  1.3233  
C3N7A0.44  3  7  8  0.44  2D  1.2981  1.3592  
C3N7A0.45  3  7  8  0.45  2D  1.3075  1.3082  
C3N7A0.46  3  7  8  0.46  2D  1.3029  1.3102  
C3N7A0.47  3  7  8  0.47  2D  1.3002  1.3505  Marginal recurrence case 
C3N7A0.48  3  7  8  0.48  2D  1.3156  1.3085  Marginal breaking case 
C3N7A0.49  3  7  8  0.49  2D  1.3430  1.3290  
C3N7A0.5  3  7  8  0.5  2D  1.2818  1.2593  
C3N9A0.42  3  9  8  0.42  2D  1.4651  1.4521  
C3N9A0.43  3  9  8  0.43  2D  1.4771  1.4942  
C3N9A0.44  3  9  8  0.44  2D  1.4836  1.5000  
C3N9A0.45  3  9  8  0.45  2D  1.4904  1.5054  
C3N9A0.46  3  9  8  0.46  2D  1.4938  1.4551  
C3N9A0.47  3  9  8  0.47  2D  1.5229  1.5265  Marginal recurrence case 
C3N9A0.48  3  9  8  0.48  2D  1.6005  1.4276  Marginal breaking case 
C3N9A0.49  3  9  8  0.49  2D  1.2659  1.4044  
method  Interpolated  Interpolated  Interpolated  
local  5.3833  0.4141  5.2393  4.1493  0.5187  4.2357  1.2974  0.2020  1.2405 
average  5.3833  0.4141  5.2238  4.1490  0.519  4.1793  1.2974  0.2020  1.2499 
dispersion plot  5.3833  0.4141  4.1493  0.5187  1.2974  0.2020  
4.1 Crest and trough slowdown
4.1.1 Local quantities
The local behaviour of the wave kinematics calls for precise definitions of wave geometry (Figure 4): amplitude , height , wavelength and crest half wavelength . The local crest steepness is defined as . An analogous definition exists for the local trough steepness .
4.1.2 Determination of the reference wave speed
For robust comparison, the reference wave speed is required and needs to be precisely and carefully determined. A spectral analysis of both time and space wave signals is undertaken to be able to estimate the true reference wave quantities rather than those based on the paddle signal.
Figure 6 shows the amplitude spectrum of a representative paddle signal. The characteristic generating frequency is systematically different from the power spectrum peak value (vertical black line) obtained by a FFT analysis. Consequently, each investigated case has to be analysed in space and time to determine the reference power spectrum peak wavenumber and peak frequency .
Several methods were used to determine the corresponding values of and and these show convergence towards identical values for the same simulation case. The linear velocity has been estimated without making any hypothesis about the (deepwater) dispersion relationship, using the peak of the frequency spectrum and the peak of the wave number spectrum , giving the linear celerity from the ratio . The deep water dispersion relationship can be used to determine , using the peak frequency .
Table 2 shows an example of three different methods to obtain using a local method, an average method and a global method for the 2D marginal breaking C3N5A0.511 case.
The local method value has been found by computing the spatial FFT at a random time and a FFT time series centered on the position of the peak of the wave group at this given time. The average method is based on determining from the ensemble mean FFT of the time series produced by 15 virtual height probes spaced equally, between two crest (or trough) maximum events. The peak wavenumber was determined from the ensemble mean FFT of 512 surface profiles between two crest (or trough) maxima. The global method is a 2D FFT of a spacetime domain of the simulation to obtain the dispersion graph (figure 7).
The values of found in Table 2 show a very small sensitivity to the extent of the time and space domain used to determine the reference velocity: remains almost constant. These frequency and wavenumber estimates are the raw values found by taking the maximum values of the spectra at a given resolution. The length of the time series and the size of the domain do not provide sufficiently high resolution in the spectral domain (indicated by and ), and a refined method was necessary. An interpolating polynomial is constructed using the neighbouring points around the power spectral maximum and the new refined maximum is found analytically using the polynomial coefficients. These new quantities give converged values of the reference velocities and that have been tabulated for all the cases in the table 1.
Some additional remarks follow on the value of the reference velocity based on the deep water dispersion relation. Differences between and have been found, which can be explained by the nature of the wave group. Each individual wave in the group has a wavelength short enough to be considered in deep water in the simulation, but the wavelength of the wavepacket envelope is long enough not to be considered in deep water. LonguetHiggins & Stewart (1964) called this a transitional wave group. These differences are secondary and do not affect the results.
In addition to the information about the spectral peak waves, the dispersion plot shows the degree of nonlinearity obtained with the multiple higher frequency “ripples”. These are the bound waves of the group, i.e. higher harmonics of the basic free waves constituting the wave group. In general, the more harmonics that are present in the dispersion graph, the more nonlinear is the whole phenomenon.
4.2 Spacetime characteristics
This study has been carried out using chirped Class 3 wave packets as defined above, but broader experimental and field studies (see Banner et al. (2014)) confirm its generality.
Figure 8 represents the characteristics plot of crests and troughs along the centre slice of the 3D numerical wave tank of a representative Class 3 wave packet. Crests, troughs and zerocrossings have been detected and followed in this coordinate system, the abscissa is the nondimensional position and the ordinate is the nondimensional time. In this plot, the theoretical linear trajectory is represented by a longdash line as a reference, crest paths are the full lines and trough paths are represented by dashed lines. Trajectories in space and time of each crest and trough are clearly seen as nonlinear, with local variations.
This plot exhibits systematic patterns for the crest and trough trajectories during the life cycle of each wave observed. The local slope of these trajectories represents the velocity of the crest or trough. A trajectory path turning anticlockwise implies slowdown and a clockwise rotation represents an acceleration. Each of the paths has a clear systematic variation in the slope and slowdown can be observed before a reacceleration, between and . These slowdowns occur close to the crest and trough maxima of each wave inside the wave packet.
4.3 Slowdown of each wave, a generic feature.
Slowdown phenomena have been observed and mentioned previously (Johannessen & Swan (2001), Johannessen & Swan (2003), Katsardi & Swan (2011)), but have not been studied as a phaseresolved wave effect. The following plots address this gap by presenting the evolution of the crest (or trough) velocity as a function of the local steepness . The curves show the trajectory of the nondimensional horizontal velocity of each carrier wave crest in the wavegroup, normalised by the theoretical linear wavespeed defined above. Velocities are calculated using the first derivative of a 3rdorder spline fitting the spacetime trajectory. Results below are presented for crests only for illustration, and troughs are not shown because they behave in exactly the same way.
Figure 9 shows the evolution of the horizontal crest velocity of every crest of a Class 3, marginal recurrent wave packet. The abscissa measures the nondimensional velocity ratio and the ordinate shows the local steepness . Each individual crest lifecycle in the group is represented by a line trajectory. A black dot indicates where each trajectory starts.
Two important finding are represented in this plot showing the instantaneous crest velocity against the local steepness. The first conclusion, which is one of the new findings of this study, is that contrary to inferences based on Stokes wavetrain theory, the crest velocity is not strongly dependent on the steepness. Minimum crest velocity in the ensemble of waves shown in figure 9 approaches the same value systematically, with a local steepness varying from to . Previous knowledge of velocity variation of crests is represented by the dashed line, which corresponds to order Stokes crest celerity.
The second new result in this plot concerns the minimum crest velocity. In this marginal recurrent case, the velocities of each crest reduce close to a minimum of of the linear reference velocity associated with this packet. The evolution of velocity is not linear with the steepness, and loops can be seen in the curves describing the history of the crest velocity, which shows the existence of hysteresis behaviour. It is noted that the minimum speed of this ensemble of crests seems to rise a few percent (from to in this case) when the local steepness rises and the strong nonlinearities enter.
These findings are corroborated by Johannessen & Swan (2001), Johannessen & Swan (2003), Katsardi & Swan (2011). These authors show the existence of a slowdown of 10% in 2D waves and of 7% in the case of 3D waves, with scatter. These slowdowns are quantified in respect to a reference linear velocity, called , which is not clearly defined. If a guess of the reference velocity estimation based on the peak of the frequency spectra is made, then the normalising velocity becomes defined above. With this hypothesis, the minimum velocity for the dominant crest in this graph goes from to , resulting in a 20% reduction in wave speed. A key question that arises from this analysis is whether this effect is local. The above papers argue that the cause of the slowdown can be associated with the appearance in the frequency wave spectrum of a peak around to of the main wavepacket peak, producing slower waves in accordance with the theoretical linear freewave dispersion relationship. However, the present study offers a complementary argument for the cause of the slowdown in section 4.4 below.
The next figure addresses whether the slowdown is only a local crest and trough effect. Figure 10 shows the velocity of the midpoint between the two zerocrossings which adjoin (or span) each crest (or trough) studied. Each zerocrossing pair has the steepness attribute of the crest they span. The abscissa is the nondimensional velocity speed with respect to the reference linear phase velocity .
The crest velocity at each crest maximum is represented by a dot, and the bars represent the extent of the variation of the crest velocity during each individual crest lifecycle. The plot shows that the zerocrossing speed changes dynamically during the wave evolution, and can undertake large relative speed variations of the order of roughly 20%. Nonetheless, the speed values when the crest reaches its maximum (dots) are close to the linear velocity, within the range . Consequently, the slowdown is a localised crest effect, and does not affect the whole packet uniformly.
4.4 Leaning, cause of the slowdown.
This section will show that the cause of the slowdown is geometric. We define the leaning as the left/right asymmetry :
(4.0) 
where represents the crest (trough) position, is the position of the downward zerocrossing and is the position of the upward zerocrossing. The leaning L quantifies the asymmetry of the geometry of the triangle formed by the two zerocrossings spanning the crest (or trough). This parameter varies from when the backwardleaning face has a vertical asymptote to +1, when the forwardleaning face is vertical.
Figure 4 defines polar coordinates associated with a local crest (or trough). is the angle at the midpoint (halfway point between the left and right zerocrossings) between the zero water level and the crest. is the radius, the distance between the crest and the middlepoint. Using polar decomposition of the velocity, the projection on the axis of the orthoradial component (the rotation part) is and the projection on the axis of the radial component is .
Figure 11 shows the sum of the normalised frame of reference velocity and the xprojection of the normalised orthoradial component of the relative velocity at the crest maximum on the abscissa against the local steepness for C3 N5, C3N7, 2D and 3D wave groups. This plot shows a very strong correlation between the rotation component and the crest velocity, and proves that the major part of the slowdown resides in the rotation (hence the leaning) of crests about the midpoint between the zerocrossings. For all the waves in the group, each individual wave tends to lean from forward to backward, reducing the crest velocity to below the expected linear crest velocity.
As an example, figure 12 summarises several key kinematic and geometric properties of the dominant crest of a C3N5 wave group. In these graphs, there is a zeroreference time (along the abscissa) when the crest reaches its maximum elevation and the local steepness is normalised to 1 at this crest maximum. The maximum local steepness for this wave is .
This plot represents the entire history of a single carrier wave, with the first panel showing the normalised steepness and the second panel showing the leaning parameter . The third panel shows the average zerocrossing speed normalised by the linear reference velocity, the fourth panel shows the crest velocity and the last panel is proportional to the leaning velocity . The graph reproduces the evolution history of the dominant crest travelling from the rear of the wave group towards the front of the group.
The wave peaks at . There is a strong correlation between the crest leaning backwards, which produces a negative time derivative of , and the crest slowdown, as seen between and . During the same time interval, the average speed of the midpoint of the adjacent zerocrossings is almost constant, with a value close to the linear reference velocity. Some peaks appear in the velocity and the time derivative of , and are associated with a strong forward jump in at . This is just an artefact of the cresttracking process and is addressed in the section (4.5).
One important finding is that always when the crest (or trough) reaches its maximum, implying that the wave crest (or trough) shape is spatially symmetric at this instant. This leaning property could potentially be used to detect crest and trough maxima in a wave field. These findings imply that the crest slowdown is a strong kinematic effect, easily understood via the composition of velocity law: the frame of reference is the average position of the two adjacent zerocrossings travelling at roughly the linear phase velocity, and the local velocity in this frame of reference is associated with the leaning forwardtobackward tendency of the crest, inducing a negative relative crest velocity. The composition of the frame velocity and the backward leaning reduce the absolute crest velocity to below the originallyexpected linear phase velocity. Other fluctuations in the absolute velocity confirming this analysis can be seen just after , where a positive velocity of the frame of reference increases the absolute velocity above the linear level.
Figure 13 presents multiple wave profiles stacked in order to form a spacetime diagram of the typical evolution along the centreline of a C3N5 maximum recurrent wave group. The axes have the nondimensional position as the abscissa and nondimensional time as the ordinate. The wave profiles are colourcoded with the value of the leaning. This figure summarises the behaviour described above: both crests and troughs undertake a systematic evolution throughout their life cycle. It can be generalised as a generic behaviour as follows: a crest (or a trough) begins its life leaning forward (red shading), peaks (the geometry of the wave becomes symmetric (green shading)) and decays leaning backwards (blue shading). Some abrupt jumps towards the front of the wave (seen as sudden changes of shading colour) are also present and are explained in the following section (4.5).
4.5 Apparent scattering and detailed geometric analysis: crestlets or riding waves?
This section explains the apparent velocity peaks found in the velocity plots (see Fig 12 around , also apparent in Fig 9 and Fig 13). Figure 14 is a spacetime plot of the same wave packet shown in Figure 13. The abscissa represents the position, the ordinate represents the time coordinate and the depth represents the nondimensional free surface elevation of a slice along the centreline of the numerical wave tank. The plot is a time sequence of wave profiles.
In this section we explain the limitation of the cresttracking algorithm and the apparent velocity spikes. The red arrows in Figure 14 identify where the phenomenon occurs. The idea of a crest or a trough being a single maximum or minimum between two zerocrossings (which is an arbitrary reference) refers to a Stokes wave. Natural unsteady waves sometimes have multiple crests and troughs between the same two zerocrossings, which evolve dynamically. This effect can be seen in Johannessen (2010) in their observational surface elevation plots. Also, Bateman et al. (2012) describes a numerical model results where crestlets and leaning in the plots can be seen in their figures 8 and 9.
During the growth of a crest, there is an interchange of these small crests or troughs, the first one evolving and then the next one appearing from the front of the packet to become temporarily the dominant crest, and then the next one appears. These crests are initiated, in the relative frame of reference of the wave group, at the front of the group and then travel rearwards. During this crest exchange, which is rapid (roughly a tenth of a time unit over a tenth of the reference length), the tracking system jumps from the decaying crest to the next growing crest, sometimes going through a ’flattened’ area of the crest. This creates a discontinuity in the detected crest position, as it is now located further towards the front of the group. This produces the apparent high velocity spike. It should be noted that between each of these peaks, the celerity of these "crestlets" decreases, back to .
The cause of these leaning effects is not yet clear, but possible insight is provided by Bettini et al. (1983) who a discuss the radiative tail of a soliton. The radiative tail is the train of trailing waves produced by a propagating wave group. In our application, these forward jumps of the crests (or troughs) align on the characteristic plot (Figure 8). The average slope gives a velocity one half of the linear velocity, this being the definition of the radiative tail in Bettini’s article. These crestlets (or troughlets) seem to propagate backwards. This effect is relative because their velocity is one half of the velocity of the wave crests (or troughs), which causes the illusion. This energy leakage creates the radiative tail.
5 Subsurface kinematics
5.1 Velocity profiles
This section provides further details on the subsurface motion by investigating the velocity profiles. Figure 15 shows the ratio of the velocity profiles computed from the NWT simulation and the 5th order Stokes velocity profiles for different cases. Shown here are results below the local crest maxima of 2D C3N5 and C3N9 breaking and nonbreaking (low steepness to marginal recurrence) , as well as some 3D cases. The Stokes waves used for comparison were fitted according to height and using the neighbouring troughs to define the wavelength. Stokes velocity profiles have been corrected for the Eulerian current by matching velocities at the bottom of the tank.
Different zones can be identified. The first zone, around , shows a relative ratio which can grow up to . While this value appears substantial, the actual velocity differences are not very large because the absolute values of the velocities are small, due to the quasiexponential attenuation of the velocity with depth.
The second zone, above , is more important as it is close to the free surface. Here the ratio varies between for maximum recurrence cases and when the wave breaks. This change is significant and highlights an important difference between the predicted Stokes orbital velocity and the actual orbital velocity computed by the simulation. The consequence can be important when estimating the impact of a wave on a structure where the drag is proportional to the square of the velocity. The Stokes fit systematically underestimates the particle crest velocity by at least 20% for very steep waves, suggesting a nearhalving of the impulse loading.
6 Energy partitioning
One of the main properties commonly assumed for waves is equipartitioning of the mean kinetic energy (KE) and potential energy (PE). With this hypothesis, meaurements can be made that only measure wave heights (hence PE), from which the kinetic energy is deduced. This equipartitioning is theoretically true for steady linear wave trains, but does it hold for a nonlinear unsteady wave group? This section investigates the energy partitioning in a global sense, first analysing the whole wave group and then studying locally what occurs at the scale of the individual crest and troughs.
6.1 Global partitioning
Having proven in the simulation validation section that energy conservation was respected, further investigation was undertaken concerning the energy partitioning for the whole wave group. Figure (16) shows the depthintegrated energy partition for the whole wave group with respect to time, for the C3N7 breaking and nonbreaking cases. Periods when the whole group was not in the integration domain have been blanked out. This figure shows that global equipartitioning holds with a small excess in the kinetic energy at the crest maximum. The ratio oscillates globally close to , showing that the total energy (TE) is almost equally divided between potential energy (PE) and kinetic energy (KE). Breaking cases have the maximum deviation from equipartitioning with of the energy residing in the PE. The energy flux necessary to produce the jet formed as the breaking event initiates has its source in the kinetic energy.
6.2 Local energy partitioning at crests and troughs
The above section showed how the depthintegrated energy partitioning still holds for the whole unsteady wave group. We now examine the local behaviour of the depthintegrated energy partitioning at both crest and trough maxima. Figure 17 shows the energy partition plotted against local steepness for 2D and 3D breaking for the maximum recurrent wave packet cases C3N5, C3N7, C3N9 and C3N10. The dashed lines represent this ratio for a 5th order Stokes wave. Triangles represent the values for the troughs and circles show the values for the crests. Remarkably, our results show that equipartitioning breaks down locally at crest and trough maxima and the ratio stays almost constant, to for the crests and 0.7 for the troughs. Locally, the Stokes wave does not conform to the equipartitioning, but behaves differently from the unsteady group waves. The ratio for the unsteady wave groups matches the Stokes wave for very small steepness and is close to , but diverges from this level at higher steepness. The notable difference between the marginal recurrent case and the breaking case is the breaking crest acquires proportionally more kinetic energy with the partition ratio falling below . It is again evident that the source of the breaking jet resides in the kinetic energy, while the wave height (hence the potential energy) remains almost constant.
7 Conclusions and recommendations
This computational study of highly nonlinear gravity water waves unsteady wave groups has revealed new insights into the behaviour of kinematic and geometric properties, as well as the associated energetics. In particular, our findings on generic crest slowdown obtained from these numerical chirped wave group simulations have been confirmed observationally for chirped as well as for bimodal wave groups in our laboratory wave basin experiments, as well as for deep water wave groups measured in an ocean tower experiment (Banner et al. (2014)).
When the theoretical constraint of steadiness is relaxed, additional asymmetric degrees of freedom are observed in the wave shape. In particular, the kinematics reveal a forwardtobackward leaning cycle for each individual crest and trough inside the wave group, resulting in a local crest or trough slowdown of around 20% of the reference linear velocity estimated from a dispersion plot. This is a geometric effect, as the velocity of the zerocrossings remains close to the linear phase velocity and the observed relative velocity variations can be explained by the geometric changes in the waveform. Velocities profiles under the crest show a large difference between the Stokesfit velocity profile and the computed velocity profile. Computed orbital velocities under the crest can be significantly higher than for the Stokes wave of corresponding height, and the ratio between them varies between 1.2 and 1.8 for the steepest cases. These values should be taken in account when designing structures close to the surface of the ocean where the Stokes values are often assumed. The associated energetics have been carefully monitored and support, within 1%, equipartitioning of the depthintegrated energy on the scale of the wave group, although at breaking onset, a slight excess of kinetic energy occurs. This depthintegrated energy equipartitioning breaks down when investigated locally at crest and trough maxima during the evolution of the wave group. These show an energy partitioning of around 69% for the potential energy at trough maxima and a partitioning range from 60% to 64% of the potential energy at crest maxima. The crest energy partitioning at breaking drops below 60% for the potential energy, indicating the initiation of the breaking jet.
Overall, these findings contribute new insights into how actual ocean wave groups propagate when the academic constraint of steadiness is relaxed. Aside from its intrinsic interest, the generic crest slowdown explains the reduced speed of initiation of breaker fronts that has been reported occasionally in the literature and is central to assimilating whitecap data accurately into spectral seastate forecast models. Also, parameterisations of airsea fluxes of momentum and energy, which depend on the square and cube of the seasurface fluid velocity, may also be modified appreciably, as may the impact loading on offshore structures.
8 Acknowledgements
Funding for this investigation was provided by the Australian Research Council under Discovery Project DP120101701. Also, FD acknowledges partial support by the European Research Council (ERC) under the research project ERC2011AdG 290562MULTIWAVE and Science Foundation Ireland under grant number SFI/12/ERC/E2227.
References
 Baker et al. (1982) Baker, Gregory R., Meiron, Daniel I. & Orszag, Steven A. 1982 Generalized vortex methods for freesurface flow problems. Journal of Fluid Mechanics 123, 477–501.
 Baldock et al. (1996) Baldock, T. E., Swan, C. & Taylor, P. H. 1996 A laboratory study of nonlinear surface waves on water. Philosophical Transactions: Mathematical, Physical and Engineering Sciences 354 (1707), pp. 649–676.
 Banner et al. (2014) Banner, M. L., Barthelemy, X., Fedele, F., Allis, M., Benetazzo, A., Dias, F. & Peirson, W. L. 2014 Linking reduced breaking crest speeds to unsteady nonlinear water wave group behavior. Phys. Rev. Lett. 112, 114502.
 Banner & Peirson (2007) Banner, Michael L. & Peirson, William L. 2007 Wave breaking onset and strength for twodimensional deepwater wave groups. Journal of Fluid Mechanics 585, 93–115.
 Bateman et al. (2012) Bateman, W.J.D., Katsardi, V. & Swan, C. 2012 Extreme ocean waves. part i. the practical application of fully nonlinear wave modelling. Applied Ocean Research 34 (0), 209 – 224.
 Bateman et al. (2001) Bateman, W.J.D., Swan, C. & Taylor, P.H. 2001 On the efficient numerical simulation of directionally spread surface water waves. Journal of Computational Physics 174 (1), 277 – 305.
 Benjamin & Feir (1967) Benjamin, T. Brooke & Feir, J. E. 1967 The disintegration of wave trains on deep water part 1. theory. Journal of Fluid Mechanics 27, 417–430.
 Bettini et al. (1983) Bettini, Alessandro, Minelli, Tullio A. & Pascoli, Donatella 1983 Solitons in undergraduate laboratory. American Journal of Physics 51 (11), 977–984.
 Beyá et al. (2012) Beyá, J.F., Peirson, W.L. & Banner, M.L. 2012 Turbulence beneath finite amplitude water waves. Experiments in Fluids 52, 1319–1330.
 Clamond & Grue (2001) Clamond, Didier & Grue, John 2001 A fast method for fully nonlinear waterwave computations. Journal of Fluid Mechanics 447, 337–355.
 Craig & Sulem (1993) Craig, W. & Sulem, C. 1993 Numerical simulation of gravity waves. Journal of Computational Physics 108 (1), 73 – 83.
 Dalrymple (1989) Dalrymple, R. A. 1989 Directional wavemaker theory with sidewall reflection. Journal of Hydraulic Research 27 (1), 23–34.
 Dalrymple & Kirby (1988) Dalrymple, Robert A. & Kirby, James T. 1988 Models for very wideangle water waves and wave diffraction. Journal of Fluid Mechanics Digital Archive 192 (1), 33–50.
 Dias et al. (2015) Dias, F., Brennan, J., Ponce de Leon, S., Clancy, C. & Dudley, J. 2015 Local analysis of wave fields produced from hindcasted rogue wave sea states. In Proceedings of the ASME 2015 34th International Conference on Ocean, Offshore and Arctic Engineering OMAE2015.
 Dommermuth & Yue (1987) Dommermuth, Douglas G. & Yue, Dick K. P. 1987 A highorder spectral method for the study of nonlinear gravity waves. Journal of Fluid Mechanics 184, 267–288.
 Ducrozet et al. (2012) Ducrozet, Guillaume, Bonnefoy, Felicien, Le Touze, David & Ferrant, Pierre 2012 A modified highorder spectral method for wavemaker modeling in a numerical wave tank. European Journal of Mechanics  B/Fluids 34, 19–34.
 Fedele (2014) Fedele, Francesco 2014 Geometric phases of water waves. EPL (Europhysics Letters) 107 (6), 69001.
 Fenton (1985) Fenton, J. 1985 A fifthorder stokes theory for steady waves. Journal of Waterway, Port, Coastal, and Ocean Engineering 111 (2), 216–234.
 Fochesato (2004) Fochesato, Christophe 2004 Modèles numériques pour les vagues et les ondes internes. PhD thesis, CMLA / Ecole Normale Supérieure de CACHAN.
 Fochesato & Dias (2006) Fochesato, Christophe & Dias, Frédéric 2006 A fast method for nonlinear threedimensional freesurface waves. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science 462 (2073), 2715–2735.
 Fochesato et al. (2007) Fochesato, Christophe, Grilli, Stéphan & Dias, Frédéric 2007 Numerical modeling of extreme rogue waves generated by directional energy focusing. Wave Motion 44 (5), 395 – 416.
 Fructus et al. (2005) Fructus, Dorian, Clamond, Didier, Grue, John & Kristiansen, Oyvind 2005 An efficient model for threedimensional surface wave simulations: Part i: Free space problems. Journal of Computational Physics 205 (2), 665 – 685.
 Gemmrich et al. (2008) Gemmrich, Johannes R., Banner, Michael L. & Garrett, Chris 2008 Spectrally resolved energy dissipation rate and momentum flux of breaking waves. J. Phys. Oceanogr. 38 (6), 1296–1312.
 Grilli & Svendsen (1990) Grilli, S.T. & Svendsen, I.A. 1990 Corner problems and global accuracy in the boundary element solution of nonlinear wave flows. Engineering Analysis with Boundary Elements 7 (4), 178 – 195.
 Grilli et al. (2001) Grilli, Stéphan T., Guyenne, Philippe & Dias, Frédéric 2001 A fully nonlinear model for threedimensional overturning waves over an arbitrary bottom. International Journal for Numerical Methods in Fluids 35 (7), 829–867.
 Grilli & Horrillo (1997) Grilli, Stéphan T. & Horrillo, Juan 1997 Numerical generation and absorption of fully nonlinear periodic waves. Journal of Engineering Mechanics 123 (10), 1060–1069.
 Grilli et al. (1989) Grilli, S. T., Skourup, J. & Svendsen, I. A. 1989 An efficient boundary element method for nonlinear water waves. Engineering Analysis with Boundary Elements 6 (2), 97 – 107.
 Grilli & Subramanya (1994) Grilli, Stéphan T. & Subramanya, Ravishankar 1994 Quasisingular integrals in the modeling of nonlinear water waves in shallow water. Engineering Analysis with Boundary Elements 13 (2), 181 – 191.
 Grilli & Subramanya (1996) Grilli, S. T. & Subramanya, R. 1996 Numerical modeling of wave breaking induced by fixed or moving boundaries. Computational Mechanics 17 (6), 374 – 391.
 Grue et al. (2003) Grue, John, Clamond, Didier, Huseby, Morten & Jensen, Atle 2003 Kinematics of extreme waves in deep water. Applied Ocean Research 25 (6), 355 – 366.
 Guyenne & Grilli (2006) Guyenne, P. & Grilli, S. T. 2006 Numerical study of threedimensional overturning waves in shallow water. Journal of Fluid Mechanics 547, 361–388.
 Hayami (1990) Hayami, K 1990 A robust numerical integration method for threedimensional boundary element analysis. In Bounday Elements XII: Applications in Stress Analysis, Potential and Diffusion (ed. Honma Tanaka, Brebbia), Bounday Elements XII, vol. 1. SpringerVerlag.
 Hayami (1991) Hayami, K 1991 A projection transformation method for nearly singular surface boundary element integrals. PhD thesis, Computational mechanics institute, Wessex institute of technology, Southampton.
 Hayami (2005a) Hayami, K. 2005a Variable transformations for nearly singular integrals in the boundary element method. Tech. Rep.. NII Technical Reports.
 Hayami (2005b) Hayami, K 2005b Variable transformations for nearly singular integrals in the boundary element method. ublications of Research Institute for Mathematical Sciences 41, 821–842.
 Hayami & Matsumoto (1994) Hayami, Ken & Matsumoto, Hideki 1994 A numerical quadrature for nearly singular boundary element integrals. Engineering Analysis with Boundary Elements 13 (2), 143 – 154.
 Hou & Zhang (2002) Hou, T.Y.a & Zhang, P.b 2002 Convergence of a boundary integral method for 3d water waves. Discrete and Continuous Dynamical Systems  Series B 2 (1), 1–34, cited By (since 1996)9.
 Jessup & Phadnis (2005) Jessup, A T & Phadnis, K R 2005 Measurement of the geometric and kinematic properties of microscale breaking waves from infrared imagery using a piv algorithm. Measurement Science and Technology 16 (10), 1961.
 Johannessen (2010) Johannessen, Thomas B. 2010 Calculations of kinematics underneath measured time histories of steep water waves. Applied Ocean Research 32 (4), 391 – 403.
 Johannessen & Swan (2001) Johannessen, T. B. & Swan, C. 2001 A laboratory study of the focusing of transient and directionally spread surface water waves. Proceedings: Mathematical, Physical and Engineering Sciences 457 (2008), pp. 971–1006.
 Johannessen & Swan (2003) Johannessen, T. B. & Swan, C. 2003 On the nonlinear dynamics of wave groups produced by the focusing of surfacewater waves. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 459 (2032), 1021–1052.
 Katsardi & Swan (2011) Katsardi, V. & Swan, C. 2011 The evolution of large nonbreaking waves in intermediate and shallow water. i. numerical calculations of unidirectional seas. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science 467 (2127), 778–805.
 Kinsman (1965) Kinsman, Blair 1965 Wind waves : their generation and propagation on the ocean surface. Englewood Cliffs, N.J.: Englewood Cliffs, N.J. : PrenticeHall.
 Kleiss & Melville (2010) Kleiss, Jessica M. & Melville, W. Kendall 2010 Observations of wave breaking kinematics in fetchlimited seas. J. Phys. Oceanogr. 40 (12), 2575–2604.
 LonguetHiggins & Stewart (1964) LonguetHiggins, M.S. & Stewart, R.w. 1964 Radiation stresses in water waves; a physical discussion, with applications. Deep Sea Research and Oceanographic Abstracts 11 (4), 529 – 562.
 LonguetHiggins (1987) LonguetHiggins, M. S. 1987 The propagation of short surface waves on longer gravity waves. Journal of Fluid Mechanics 177, 293–306.
 LonguetHiggins & Stewart (1960) LonguetHiggins, M. S. & Stewart, R. W. 1960 Changes in the form of short gravity waves on long waves and tidal currents. Journal of Fluid Mechanics 8, 565–583.
 Ma (2010) Ma, Qingwei. 2010 Advances in numerical simulation of nonlinear water waves. Hackensack, NJ: World Scientific.
 Ma et al. (2001) Ma, Q. W., Wu, G. X. & Eatock Taylor, R. 2001 Finite element simulation of fully nonlinear interaction between vertical cylinders and steep waves. part 1: methodology and numerical procedure. International Journal for Numerical Methods in Fluids 36 (3), 265–285.
 Melville (1983) Melville, W. K. 1983 Wave modulation and breakdown. Journal of Fluid Mechanics 128, 489–506.
 Melville & Matusov (2002) Melville, W. Kendall & Matusov, Peter 2002 Distribution of breaking waves at the ocean surface. Nature 417 (6884), 58–63.
 Melville et al. (2002) Melville, W. Kendall, Veron, Fabrice & White, Christopher J. 2002 The velocity field under breaking waves: coherent structures and turbulence. Journal of Fluid Mechanics 454, 203–233.
 Miller et al. (1991) Miller, Sarah J., Shemdin, Omar H. & LonguetHiggins, Michael S. 1991 Laboratory measurements of modulation of shortwave slopes by long surface waves. Journal of Fluid Mechanics 233, 389–404.
 Nicholls (1998) Nicholls, David P. 1998 Traveling water waves: Spectral continuation methods with parallel implementation. Journal of Computational Physics 143 (1), 224 – 240.
 Park et al. (2003) Park, J. C., Kim, M. H., Miyata, H. & Chun, H. H. 2003 Fully nonlinear numerical wave tank (nwt) simulations and wave runup prediction around 3d structures. Ocean Engineering 30 (15), 1969 – 1996.
 Phillips (1977) Phillips, Owen 1977 The dynamics of the upper ocean, 2nd edn. Cambridge University Press Cambridge ; New York.
 Phillips (1985) Phillips, O. M. 1985 Spectral and statistical properties of the equilibrium range in windgenerated gravity waves. Journal of Fluid Mechanics 156, 505–531.
 Rapp & Melville (1990) Rapp, R. J. & Melville, W. K. 1990 Laboratory measurements of deepwater breaking waves. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 331 (1622), 735–800.
 Shemer (2013) Shemer, L. 2013 On kinematics of very steep waves. Natural Hazards and Earth System Science 13 (8), 2101–2107.
 Song & Banner (2002) Song, JinBao & Banner, Michael L. 2002 On determining the onset and strength of breaking for deep water waves. part i: Unforced irrotational wave groups. Journal of Physical Oceanography 32 (9), 2541–2558.
 Stansell & MacFarlane (2002) Stansell, Paul & MacFarlane, Colin 2002 Experimental investigation of wave breaking criteria based on wave phase speeds. Journal of Physical Oceanography 32 (5), 1269–1283.
 Stokes (1847) Stokes, G. G. 1847 On the theory of oscillatory waves. Trans. Camb. Phil. Soc. 8, 441–455.
 Sutherland et al. (1995) Sutherland, J., Greated, C.A. & Easson, W.J. 1995 Variations in the crest kinematics of wave groups. Applied Ocean Research 17 (1), 55 – 62.
 Tayfun (1986) Tayfun, M. Aziz 1986 On narrowband representation of ocean waves: 1. theory. Journal of Geophysical Research: Oceans 91 (C6), 7743–7752.
 Telles (1987) Telles, J.C.F 1987 A selfadaptive coordinate transformation for efficient numerical evaluation of general boundary element integrals. Int. J. Numer. Meth. Engng. 24 (5), 959–973.
 Viotti et al. (2014) Viotti, Claudio, Carbone, Francesco & Dias, Frédéric 2014 Conditions for extreme wave runup on a vertical barrier by nonlinear dispersion. Journal of Fluid Mechanics 748, 768–788.
 West et al. (1987) West, Bruce J., Brueckner, Keith A., Janda, Ralph S., Milder, D. Michael & Milton, Robert L. 1987 A new numerical method for surface hydrodynamics. Journal of Geophysical Research: Oceans 92 (C11), 11803–11824.
 Wiegel (1964) Wiegel, R. L. 1964 Oceanographical Engineering. Englewood Cliffs, N.J.: Englewood Cliffs, N.J., PrenticeHall.
 Xu & Guyenne (2009) Xu, Liwei & Guyenne, Philippe 2009 Numerical simulation of threedimensional nonlinear water waves. Journal of Computational Physics 228 (22), 8446 – 8466.
 Xue et al. (2001) Xue, Ming, Xü, Hongbo, Liu, Yuming & Yue, Dick K. P. 2001 Computations of fully nonlinear threedimensional wavebody interactions. part 1. dynamics of steep threedimensional waves. Journal of Fluid Mechanics 438 (1), 11–39.
 Zakharov (1968) Zakharov, V.E. 1968 Stability of periodic waves of finite amplitude on the surface of a deep fluid. Sov. Phys. J. Appl. Mech. Tech. Phys. 4, 190–194.