1 Introduction
This study reports on a theoretical criterion for particles to surf water waves. Here, a particle travelling near the phase velocity of the underlying wave that experiences an increase in its horizontal velocity due to the wave is said to surf. The theory is based on an analysis of the John equation (John Reference John1953) for zero-stress free-surface flows under the action of gravity. Employing this equation, it is shown that particles in the region of the forward face of the crest of a very steep or breaking wave, travelling at speeds near the phase velocity of the underlying wave, increase in speed. This is the main result of this paper.
As any experienced surfer knows, the ideal location to be on a wave is near the forward face of the crest of an unbroken wave approaching breaking. This is commonly referred to as the ‘curl’. Staying in the curl maximizes distance travelled, which is correlated with the ride enjoyment. A similar phenomenon occurs when white water from a breaking wave rides the underlying wave, which is an important process in both shallow and deep-water upper-ocean dynamics.
Shallow-water breaking waves transport mass, generate currents and mix the surf zone (Inman, Tait & Nordstrom Reference Inman, Tait and Nordstrom1971; Peregrine Reference Peregrine1983; Battjes Reference Battjes1988). Wave breaking in the surf zone radically alters along-shore currents, versus the traditional model of gradients in the radiation stress (Longuet-Higgins & Stewart Reference Longuet-Higgins and Stewart1964; Feddersen et al. Reference Feddersen, Olabarrieta, Guza, Winters, Raubenheimer and Elgar2016). Mixing is crucial for regulating the thermodynamics (Sinnett & Feddersen Reference Sinnett and Feddersen2014) and chemistry (Bresnahan et al. Reference Bresnahan, Wirth, Martz, Andersson, Cyronak, DAngelo, Pennise, Melville, Lenain and Statom2016) in the near-shore region, while currents transport pollutants along the coastline (Feddersen Reference Feddersen2007; Clark, Feddersen & Guza Reference Clark, Feddersen and Guza2010). A better understanding of breaking is important for characterizing near-shore processes.
The Lagrangian drift, or mass transport, of deep-water surface waves is an important feature of upper-ocean dynamics (Phillips Reference Phillips1977). Traditionally, the drift is associated with permanent progressive waves and for linear waves this is known as Stokes drift. Stokes drift is believed to be essential for generating Langmuir circulations (Craik & Leibovich Reference Craik and Leibovich1976; Leibovich Reference Leibovich1983; McWilliams, Restrepo & Lane Reference McWilliams, Restrepo and Lane2004), as well as transporting mass, which might strongly modulate global ocean circulation (McWilliams & Restrepo Reference McWilliams and Restrepo1999). However, Stokes drift is very difficult to measure in the ocean and its form for realistic sea states is still an active area of research (Smith Reference Smith2006).
Alternatively, wave breaking, which is intermittent in space and time, transports mass and generates ocean currents (Rapp & Melville Reference Rapp and Melville1990; Melville, Veron & White Reference Melville, Veron and White2002). Direct numerical simulations (Sullivan, McWilliams & Melville Reference Sullivan, McWilliams and Melville2004) and large eddy simulations (Sullivan, McWilliams & Melville Reference Sullivan, McWilliams and Melville2007) show that wave breaking, using a simple body force model for the effects of breaking on the mean flow, is an important component of the dynamics of the upper ocean. This motivated a closer examination of the integral properties of the breaking induced flow (Pizzo & Melville Reference Pizzo and Melville2013; Pizzo, Melville & Deike Reference Pizzo, Melville and Deike2016) in the context of impulsive forces and vortex dynamics.
Similarly, very recently, numerical work has been undertaken to study the Lagrangian transport due to breaking induced flow (Deike, Pizzo & Melville Reference Deike, Pizzo and Melville2017). The authors performed a direct numerical simulation of a breaking deep-water wave packet by using a dispersive focusing technique (Rapp & Melville Reference Rapp and Melville1990). Particles are tracked during the breaking process and two of these paths are shown in figure 1. Figure 1(a) shows a particle that does not surf the breaking wave and has a relatively small horizontal displacement of approximately
$0.3$
m. Figure 1(b) shows a particle, which is initially close to the particle in figure 1(a), which surfs the breaking wave and is transported around
$2.2$
m, which is more than seven times greater than the horizontal displacement in figure 1(a). These images motivated this research.

Figure 1. Particle displacements for a breaking deep-water wave packet, generated by a direct numerical simulation from Deike et al. (Reference Deike, Pizzo and Melville2017, see Deike, Popinet & Melville Reference Deike, Popinet and Melville2015; Pizzo et al.
Reference Pizzo, Melville and Deike2016 for further discussion of the numerical methods). The linear prediction of the maximum slope at breaking (Drazen, Melville & Lenain Reference Drazen, Melville and Lenain2008) is 0.4, while the central frequency of the packet is
$0.88$
Hz corresponding to a central wavelength of approximately
$2.02$
m. The particle trajectory is drawn by the black line, while the time evolution of the particle is indicated by the colour. Dispersive focusing is used to localize energy in space and time (Rapp & Melville Reference Rapp and Melville1990) leading to breaking that occurs around
$t=2~\text{s}$
. Two different particles are shown. In (a), the particle does not catch the wave and its displacement is relatively small, approximately
$0.3$
m. In (b), the particle surfs the wave and is displaced a relatively large distance of about
$2.2$
m, more than seven times further than the non-breaking case.
To this end the John equation, which describes the motion of particles on a zero-stress free surface, is considered. Recently, the structure of these equations has been closely examined (Fedele, Chandre & Farazmand Reference Fedele, Chandre and Farazmand2016), following the development of these equations in three dimensions (Sclavounos Reference Sclavounos2005). A generalization of this equation has been used to examine the geometry of the surface of breaking waves by Longuet-Higgins (Reference Longuet-Higgins1980, Reference Longuet-Higgins1982).
In this study, an asymptotic analysis of the behaviour of a perturbed particle obeying the John equation is undertaken. This is found to elucidate a mechanism that is present for particles travelling near the phase velocity of the underlying wave. The analysis yields a criterion for a particle to surf an underlying wave based on the geometry of the free surface and the kinematics of the particle. A generalization of this result is then provided and both results are examined theoretically and numerically.
The outline of this manuscript is as follows. In § 2, the equation of John and a criterion for particles to surf an underlying wave are presented. In § 3 a breaking deep-water wave is considered and the results are discussed. Finally, in § 4 implications of this work are considered.
2 Criterion for surfing particles
The equation of John (Reference John1953) is presented. Instead of following the original derivation of John, who used the language of complex variables, the work of Sclavounos (Reference Sclavounos2005, see also Fedele et al. Reference Fedele, Chandre and Farazmand2016) is reproduced for completeness (in two spatial dimensions). A criterion for particles to surf an underlying permanent progressive wave is derived. This criterion is then generalized and discussed.
2.1 The John equation
Recall, for a zero-stress surface the free surface is a line of constant pressure. Furthermore, curves for which
$\unicode[STIX]{x1D701}-\unicode[STIX]{x1D702}$
takes a constant value, where
$\unicode[STIX]{x1D702}$
represents the free surface and
$\unicode[STIX]{x1D701}$
is the vertical position of a particle, are level curves. Hence,
$\unicode[STIX]{x1D735}(\unicode[STIX]{x1D701}-\unicode[STIX]{x1D702})$
is normal to this level curve. Next, as the pressure gradients must be normal to the surface (the pressure is by definition the isotropic component of the stress), the following relationship holds:

The Euler equations in two spatial dimensions are

with
$\unicode[STIX]{x1D743}=(\unicode[STIX]{x1D709}(t),\unicode[STIX]{x1D701}(t))$
the horizontal and vertical positions of the particles, respectively,
$\unicode[STIX]{x1D70C}$
the density of water,
$p$
the pressure and
$\boldsymbol{g}=(0,-g)$
the acceleration due to gravity. Solving for the pressure and substituting this into (2.1) implies

with dots representing (total) derivatives with respect to time. Next, recall the kinematic boundary condition at the free surface is

Note,
$\dot{\unicode[STIX]{x1D702}}=\unicode[STIX]{x1D702}_{t}+\unicode[STIX]{x1D702}_{x}\dot{\unicode[STIX]{x1D709}}$
so that

Therefore, equation (2.3) becomes

which is the equation found by John (Reference John1953). This takes the form of a Ricatti-type equation, with a corresponding variational structure (Sclavounos Reference Sclavounos2005, see also Fedele et al. Reference Fedele, Chandre and Farazmand2016).
2.2 Surfing criterion for permanent progressive waves
A criterion for particles to surf an underlying permanent progressive wave is now derived. That is,
$\unicode[STIX]{x1D702}(x,t)=\unicode[STIX]{x1D702}(x-ct)$
for some constant phase velocity
$c$
. Note, permanent progressive waves exist in deep water (e.g. Stokes waves) and shallow water (e.g. soliton and cnoidal solutions to the Korteweg–de Vries equation). This assumption is relaxed in § 2.3, and examined numerically in § 3.
Under the above assumption, the governing equation may be rewritten as

with
$F,G,H$
corresponding to the relevant coefficients in (2.6). Up until this point, the equation still holds exactly for zero-stress free-surface permanent progressive waves under the action of gravity.
The horizontal velocity is decomposed, so that

where
$U_{0}$
is the initial velocity of the particle, a constant, and
$u=u(t)$
is a perturbation from this value. The particle location
$\unicode[STIX]{x1D709}$
can be found by integrating equation (2.8), i.e.

Formally, we take
$u\sim O(\unicode[STIX]{x1D716})$
for a small parameter
$\unicode[STIX]{x1D716}\equiv |u/U_{0}|$
while
$\dot{u}\sim O(1)$
, the later assumption being a restriction on the time scale over which this analysis is applicable.
Then, the governing equation to lowest order in
$\unicode[STIX]{x1D716}$
, namely
$O(\unicode[STIX]{x1D716}^{0})$
, is

where

which is a measure of how close the initial velocity of the particle is to the phase velocity of the underlying wave.
Substituting in the form of
$G$
and
$H$
from (2.6) gives

with primes representing differentiation with respect to
$x_{0}$
. Integrating equation (2.12) from
$0$
to a time
$t$
yields

where
$N$
is the anti-derivative of
$\unicode[STIX]{x1D702}^{\prime }/(1+\unicode[STIX]{x1D702}^{\prime 2})$
.
In the limit that the particle speed
$U_{0}$
approaches the phase velocity of the underlying wave, equation (2.13) becomes (using the fundamental theorem of calculus and the definition of the derivative)

Equation (2.14) provides a criterion for particles to surf the underlying wave. Particles speed up in regions where the right-hand side of (2.14) is positive, so that it is of interest to define the function

The term
$-\unicode[STIX]{x1D702}^{\prime }$
is positive on the forward face of a focusing wave, where the slope is negative, while the denominator is positive definite. This horizontal acceleration (orthogonal to gravity) then acts as a reduced gravity
${\mathcal{A}}$
, which is solely a function of the free-surface geometry and in particular takes on a maximum value when
$\unicode[STIX]{x1D702}^{\prime }=1$
and hence is bounded by
$\pm g/2$
.
2.3 Generalized surfing criterion
The criterion presented in the previous section elucidates the nature of the interaction between the particle and underlying waves. This analysis is valid for small
$\unicode[STIX]{x1D716}$
.
Professor Francesco Fedele (personal communication) has kindly provided a simple argument for the maximum acceleration experienced by particles. This result does not use perturbative analysis nor does it assume the underlying wave is of permanent progressive form.
Following Fedele (Reference Fedele2014), consider the collection of curves along the free surface with fixed slope
$\unicode[STIX]{x1D6FC}_{0}$
. That is, from Longuet-Higgins (Reference Longuet-Higgins1957, equation (2.5.18)), let
$x$
and
$x+\text{d}x$
be the location of the prescribed slope at times
$t$
and
$t+\text{d}t$
, then

so that the speed of propagation of this fixed slope, defined as
$c=\text{d}x/\text{d}t$
, is given by

where
$c=c(t;\unicode[STIX]{x1D6FC}_{0})$
is the speed of propagation of the free surface with fixed slope
$\unicode[STIX]{x1D6FC}_{0}$
. For example, taking
$\unicode[STIX]{x1D6FC}_{0}=0$
restricts us to following the path of crests or troughs of the wave.
Equation (2.6) can be rewritten as

where
$\unicode[STIX]{x1D709}_{P}$
is a point on the free surface where the slope is given by
$\unicode[STIX]{x1D6FC}_{0}$
and
$U=\dot{\unicode[STIX]{x1D709}}_{P}$
while
$A=\ddot{\unicode[STIX]{x1D709}}_{P}$
. Next, the extremum of
$A$
is found by varying over all possible slopes, or equivalently over each particle velocity
$U$
. That is,

which returns the condition found in § 2.2. The maximum acceleration,
$A(U=c)\equiv A_{m}$
, is then

which is a generalized version of the acceleration given by
${\mathcal{A}}$
.
Note, under the assumption that
$\unicode[STIX]{x1D702}=\unicode[STIX]{x1D702}(x-ct)$
,

in agreement with (2.15).
3 Numerical example: deep-water perturbed Stokes wave
In order to characterize the behaviour of the acceleration found in the previous section, a perturbed deep-water Stokes wave is examined. This is facilitated through numerical integration of the fully nonlinear potential flow equations (see, for instance, § 3.1 of Phillips Reference Phillips1977). A perturbed deep-water Stokes wave (Longuet-Higgins Reference Longuet-Higgins1978a
; Tanaka Reference Tanaka1983), with
$c=1.0922$
and
$ak=0.4296$
, for
$a$
the wave amplitude and
$k$
the wavenumber as defined in Longuet-Higgins (Reference Longuet-Higgins1975), is considered. Here,
$g=1$
and the wavelength of the unperturbed Stokes wave is taken to be
$2\unicode[STIX]{x03C0}$
. The normal mode eigenvalue, and eigenfunction scaling constant, are provided in appendix A.
The evolution of the free surface, in a reference frame moving at the phase velocity of the unperturbed Stokes wave, is shown in figure 2. The leading face of the wave steepens, eventually forming a jet on the forward face of the wave, as shown in the inset of the figure.

Figure 2. The evolution of the free surface of a perturbed deep-water Stokes wave in a frame moving at the phase velocity of the unperturbed Stokes wave,
$c$
. The forward face of the wave steepens, leading to a jet forming on the forward face of the crest, which is shown in more detail in the inset. Here
$c=1.0922$
and
$ak=0.4296$
and the eigenvalue of the normal mode has value
$0.0575$
. Refer to text for full details.

Figure 3. Surface particle kinematics for an overturning perturbed deep-water Stokes wave. The horizontal particle velocities
$\dot{\unicode[STIX]{x1D709}}$
and accelerations
$\ddot{\unicode[STIX]{x1D709}}$
are shown along the free surface at various times as the wave approaches breaking. The maximum velocity of the particles is on the forward face of the crest, and increases during the focusing and breaking event. The arrow in the upper right corner of the acceleration plots has magnitude
$g$
. There are large accelerations (nearly
$4g$
at time
$t=9.63$
) just below the crest of the wave on the forward face.
Note, I have been representing the free surface as a graph
$\unicode[STIX]{x1D702}(x,t)$
. However, there is no difficulty in extending this analysis to multivalued free surfaces, parameterized by
$(X(s),Y(s))$
a function of the parameter
$s\in (0,2\unicode[STIX]{x03C0})$
.
The particle kinematics are shown in figure 3 which shows the evolution of the velocity,
$\dot{\unicode[STIX]{x1D709}}$
, and acceleration
$\ddot{\unicode[STIX]{x1D709}}$
at various times approaching breaking. The particle velocities approach
$c$
in the region where the wave face becomes vertical, with strong horizontal velocities in the jet of the breaking wave (Vinje & Brevig Reference Vinje and Brevig1981). The initial wave has maximum horizontal fluid velocity that is
$0.79c$
. As the wave focuses the speed of the particles at the crest increase reaching values exceeding
$1.2c$
, with an asymmetry forming on the forward face. The accelerations,
$\ddot{\unicode[STIX]{x1D709}}$
, are initially weak (with an initial maximum value around
$0.43g$
) and become large (nearly
$4g$
) in the breaking region as is shown in more detail in figure 5.
The theoretical predictions of the accelerations,
${\mathcal{A}}$
(red line) and
$A_{m}$
(black line), are shown in figure 4(a,b). Figure 4(a) shows the evolution of the maximum of
${\mathcal{A}}$
and
$A_{m}$
. Both theories predict a maximum acceleration near the crest of the breaking wave. Figure 4(b) shows the functional dependence of the accelerations in space. These accelerations become localized as the wave steepens and starts to overturn. Furthermore, in regions where the wave starts to rapidly change its geometry, large accelerations are found and
${\mathcal{A}}$
is no longer a good approximation and
$A_{m}$
must be employed.

Figure 4. (a) The locations of the maximum values of the acceleration predicted by the theoretical arguments in § 2 along the free surface displacement
$\unicode[STIX]{x1D702}$
, shown in grey. Both predictions show that the maximum acceleration, that is the surfing sweet spot for the wave, occurs near the crest of the wave on the forward face. (b) A comparison of
${\mathcal{A}}$
(shown in red) and
$A_{m}$
(shown in black) as the wave focuses. Note the scales of the free surface and the accelerations are not the same. The circles show the maximum values of the accelerations. As the packet focuses, both predictions become localized at the crest of the breaking wave. However,
${\mathcal{A}}$
cannot describe the large accelerations that occur during the formation of the overturning jet. In this region
$A_{m}$
must be employed.

Figure 5. Accelerations during the breaking event according to numerically computed values (grey arrows) and theoretical predictions, shown on the right in red. The gold line represents particles travelling at speeds within 10% of the phase velocity
$c$
, the maximum acceleration is indicated by the black circle and the black arrow in the upper right corner has magnitude
$g$
. Plots (b) and (d) show the acceleration along a parameter
$s$
, with the gold line in (b) and (d) highlighting the domain of interest corresponding to the gold line in (a) and (c), respectively. There is good agreement between the theory presented in § 2.3 and the numerical values. The large accelerations become spatially localized during the overturning event.
The acceleration in the breaking region is shown in figure 5. The numerical values of the accelerations at two times during the overturning event are shown in (a,c). Particles travelling within 10 % of the phase velocity of the unperturbed Stokes wave
$c$
are highlighted by the gold curves in figures 5(a) and (c), while the corresponding region along the parameter
$s$
are shown in (b) and (d), respectively. These correspond to the regions of maximum acceleration, in accordance with the theoretical results. The theoretical prediction
$A_{m}$
is shown in red, with the numerical velocities shown again by the grey arrows. There is excellent agreement between the theoretical predictions and the numerical results for this example.
The accelerations become large in this region, reaching values of nearly
$4g$
. Furthermore, these accelerations become increasing localized to this region underneath the lip of the breaking wave.
4 Conclusion
A criterion to surf an underlying wave is derived. This was found by employing the equation of John for a zero-stress free surface under the action of gravity. The optimal location for surfing to occur, for a deep-water perturbed Stokes wave approaching breaking, was established for particles travelling near the phase velocity of the underlying wave and occurs below the crest on the forward face of the wave. Although a deep-water perturbed Stokes wave was considered, this behaviour (sometimes referred to as a crest instability (Longuet-Higgins & Dommermuth Reference Longuet-Higgins and Dommermuth1997)), is characteristic of breaking induced by a variety of methods in both shallow and deep water.
Two approaches were taken. In the first an asymptotic analysis was performed for permanent progressive waves. A generalized criterion was then presented. Both predict that maximum accelerations occur near the crest of the forward face of the wave in a localized region. The asymptotic model breaks down when the free-surface geometry is rapidly varying, and the full description,
$A_{m}$
, is needed. This predicts a localized ‘sweet spot’ to optimally surf a wave.
The phenomenon described in this paper is analogous to Landau damping in plasmas (Landau Reference Landau1946; Dawson Reference Dawson1961; Mouhot & Villani Reference Mouhot and Villani2010). Landau damping occurs when charged particles encounter an electric field. Particles travelling near the phase velocity of the electric field lose or gain energy based on the initial phase of the electric field and the charged particle. In a qualitative sense then, some charged particles surf the electric field, removing energy from it. This description partially motivated the present study.
The analysis performed in this paper is important for understanding the transfer of energy from the wave field to the water column by wave breaking (Drazen et al. Reference Drazen, Melville and Lenain2008; Romero, Melville & Kleiss Reference Romero, Melville and Kleiss2012; Pizzo et al. Reference Pizzo, Melville and Deike2016). Here it was found that the perturbed particles increase in speed due to the work done on them by the underlying wave. To adequately quantify the bulk scale behaviour of a collection of particles in a breaking wave, a statistical description of these effects would be necessary. This study might also be of use in better understanding optimal ways for a surfer to surf a wave (Dally Reference Dally2001). Furthermore, it would be of interest to consider the implications of this study for the capillary rollers found in Longuet-Higgins (Reference Longuet-Higgins1992), with potential implications for a better understanding of the onset of wave breaking (Perlin, Choi & Tian Reference Perlin, Choi and Tian2013).
Finally, detailed knowledge of the kinematics of particles on the free surface during breaking is currently an active area of research (Deike et al. Reference Deike, Pizzo and Melville2017). Laboratory and field analysis of the kinematics and geometry of shallow and deep-water breaking waves is warranted to further corroborate the ideas presented in this paper.
Acknowledgements
I thank Professor F. Fedele for kindly providing the material in § 2.3. Furthermore, I gratefully acknowledge useful conversations with K. Melville and L. Deike on this subject. I also thank L. Deike for providing figure 1. The Dold and Peregrine code was kindly provided to K. Melville by Professor Peregrine (deceased). This research was supported by grants from NSF (OCE) and ONR (Physical Oceanography) to K. Melville.
Appendix A. The Lagrangian of Balk (Reference Balk1996) and the numerical integration of the governing equations
The Lagrangian formulation of the potential flow deep-water surface wave equations, due to Balk (Reference Balk1996) is employed in § 3 to find the eigenvalues and eigenvectors of steep Stokes waves (as well as the Stokes coefficients themselves). This formulation has been developed in Pizzo (Reference Pizzo2015), and a manuscript further analysing these equations is in preparation (Pizzo & Melville Reference Pizzo and Melville2017). The full details of the model are outside of the scope of the present discussion; however, the basic formulation is presented for completeness.
The dependent variables of the system are
$\{Y_{\pm 1},Y_{\pm 2},\ldots ,Y_{\pm M}\}$
where
$M$
is the resolution of the system, and each variable
$Y_{j}$
is a function of time. These variables are related to the parametric description of the surface
$(X,Y)$
. In particular,
$X+\text{i}Y=x+\sum _{k=-M}^{M}(X_{k}+\text{i}Y_{k})\text{e}^{-\text{i}ks}$
, where
$X_{k}=\text{i}\unicode[STIX]{x1D70E}_{k}Y_{k}$
with
$\unicode[STIX]{x1D70E}_{k}=(-1,0,1)$
for
$k<0,k=0,k>0$
respectively and
$s\in (0,2\unicode[STIX]{x03C0})$
. The Lagrangian of the system (where recall
$g=1$
and the wavelength is taken to be
$2\unicode[STIX]{x03C0}$
) may be written as (Balk Reference Balk1996)

where

and mass conservation is ensured by setting
$Y_{0}=-\sum _{k=-M}^{M}|k|Y_{k}Y_{-k}$
.
The equations of motion are found by varying
$L$
with respect to
$Y_{k}$
and take the form of a large set of coupled nonlinear ordinary differential equations that are low-order polynomials in the dependent variables. Note, Balk (Reference Balk1996) pointed out that the algebraic set of equations describing Stokes waves takes the natural (quadratic) form (found in a different way by Longuet-Higgins Reference Longuet-Higgins1978b
) when
$Y_{k}=\unicode[STIX]{x1D6FC}_{k}\text{e}^{\text{i}kct}$
for
$\unicode[STIX]{x1D6FC}_{k}$
the (constant) Stokes coefficients. The initial conditions were generated by solving this quadratic set of equations for
$\unicode[STIX]{x1D6FC}_{k}$
by employing a Newton–Raphson method (Longuet-Higgins Reference Longuet-Higgins1985). The linear stability analysis was performed by letting
$Y_{k}=(\unicode[STIX]{x1D6FC}_{k}+\unicode[STIX]{x0394}A_{k}(t))\text{e}^{\text{i}kct}$
and finding the Lagrangian to
$O(\unicode[STIX]{x1D6E5}^{2})$
. This will be discussed in detail in Pizzo & Melville (Reference Pizzo and Melville2017).
The time integration of these equations is efficient unless there are regions of high curvature (Meiron, Orszag & Israeli Reference Meiron, Orszag and Israeli1981; Baker & Xie Reference Baker and Xie2011), at which point the spectral coefficients converge slowly and
$M$
becomes unreasonably large. Therefore, the time evolution of the system is solved by employing the model of Dold (Reference Dold1992), which has the advantage that points cluster in regions of high curvature, making it extremely efficient. A total of
$2048$
points were used, with the arc length between points chosen based on the initial particle speed, following Tanaka et al. (Reference Tanaka, Dold, Lewy and Peregrine1987). The scheme of Dold has been verified repeatedly in the literature (see, for example, Dold & Peregrine Reference Dold and Peregrine1986; Tanaka et al.
Reference Tanaka, Dold, Lewy and Peregrine1987; Dold Reference Dold1992) and here its accuracy was tested in several ways. First, the Stokes wave considered in this study was propagated, unperturbed, and its stability over two wave periods was observed. Furthermore, the energy was conserved to one part in
$O(10^{5})$
. Next, the equations of Balk were integrated numerically and compared with the results from the scheme of Dold, and accuracy was corroborated leading up to the breaking event. Finally, multiple resolutions were chosen to ensure the scheme converged.
For the example considered in § 3,
$\unicode[STIX]{x1D706}\approx 0.0575$
and as the eigenfunction is not unique (it is scale invariant) we take
$\unicode[STIX]{x1D6E5}=0.1$
. These values are found through the scheme outlined above, and are in agreement with existing studies (Tanaka Reference Tanaka1983; Longuet-Higgins & Dommermuth Reference Longuet-Higgins and Dommermuth1997).