2 Lagrangian dynamical features at the onset of wave breaking

Our goal is to study dynamical features for a small ripple riding on the front face of a nonlinear breaking wave. First we will describe, through numerical simulations, Lagrangian properties (of material points) on the surface of a breaking wave. Then we consider the added ripple.

The phenomenon of gravity wave breaking is described by the Euler equations

(2.1) $$\begin{eqnarray}\boldsymbol{v}_{t}+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}=-\unicode[STIX]{x1D735}p/\unicode[STIX]{x1D70C}+\boldsymbol{g},\quad \text{div}\,\boldsymbol{v}=0,\end{eqnarray}$$

where the effects of viscosity and surface tension are neglected; comments on capillary effects will be made later. In the two-dimensional formulation, $(x,y)$ are the horizontal and vertical coordinates, $\boldsymbol{v}=(u,v)$ is the fluid velocity satisfying the incompressibility condition, $\unicode[STIX]{x1D70C}$ is the constant density and $\boldsymbol{g}$ is the vector of gravitational acceleration. We will consider the potential flow over a flat rigid bottom, which is bounded from above by a free surface $y=F(x,t)$ . A dimensionless formulation of the potential theory equations (Whitham Reference Whitham1974; Landau & Lifshitz Reference Landau and Lifshitz1987) will be considered with unit water depth, density and gravitational acceleration. For numerical convenience, we assume that the flow is periodic in the horizontal direction with period $2\unicode[STIX]{x03C0}$ , and set the rigid bottom coordinate at $y=-1$ .

We write the (kinematic and dynamic) boundary conditions at the free surface as

(2.2a-c ) $$\begin{eqnarray}y=F(x,t):\quad v=F_{t}+F_{x}u,\quad p=P_{atm},\end{eqnarray}$$

where $P_{atm}$ is a constant atmospheric pressure, and the subscripts $t$ , $x$ and $y$ are used in this section to denote partial derivatives. At the bottom, we have

(2.3a,b ) $$\begin{eqnarray}y=-1:\quad v=0.\end{eqnarray}$$

For potential incompressible flow, we can introduce the complex potential $\unicode[STIX]{x1D6F7}=\unicode[STIX]{x1D711}+\text{i}\unicode[STIX]{x1D713}$ , which is a holomorphic function of $z=x+\text{i}y$ in the fluid domain. The real potential function $\unicode[STIX]{x1D711}(x,y,t)$ and the streamfunction $\unicode[STIX]{x1D713}(x,y,t)$ are related to the fluid velocities as

(2.4a,b ) $$\begin{eqnarray}u=\unicode[STIX]{x1D711}_{x}=\unicode[STIX]{x1D713}_{y},\quad v=\unicode[STIX]{x1D711}_{y}=-\unicode[STIX]{x1D713}_{x}.\end{eqnarray}$$

Expressing velocities from (2.4) and the pressure from the Bernoulli equation, the boundary conditions (2.2)–(2.3) can be written as (Whitham Reference Whitham1974; Landau & Lifshitz Reference Landau and Lifshitz1987)

(2.5a-c ) $$\begin{eqnarray}\displaystyle y=F(x,t):\quad \unicode[STIX]{x1D711}_{y}=F_{t}+F_{x}\unicode[STIX]{x1D711}_{x},\quad \unicode[STIX]{x1D711}_{t}+{\textstyle \frac{1}{2}}(\unicode[STIX]{x1D711}_{x}^{2}+\unicode[STIX]{x1D711}_{y}^{2})+gy=0,\end{eqnarray}$$

(2.6a,b ) $$\begin{eqnarray}\displaystyle y=-1:\quad \unicode[STIX]{x1D711}_{y}=0.\end{eqnarray}$$

Let $z(\unicode[STIX]{x1D701},t)$ with $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D709}+\text{i}\unicode[STIX]{x1D702}$ be a conformal mapping from a horizontal strip $-K\leqslant \unicode[STIX]{x1D702}\leqslant 0$ onto the fluid domain at time $t$ . Such a mapping provides the free surface parametrization as $x+\text{i}y=z(\unicode[STIX]{x1D709},t)$ , with a real coordinate $\unicode[STIX]{x1D709}$ . Note that this mapping does not require that the free surface equations can be resolved with respect to the vertical coordinate, $y=F(x,t)$ , i.e. it can be used when the free boundary has overhanging sections. With the method of complex analysis (Dyachenko, Zakharov & Kuznetsov Reference Dyachenko, Zakharov and Kuznetsov1996b ; Zakharov, Dyachenko & Vasilyev Reference Zakharov, Dyachenko and Vasilyev2002; Ribeiro, Milewski & Nachbin Reference Ribeiro, Milewski and Nachbin2017) one can describe the flow in the whole fluid domain in terms of real functions defined at the free surface; see § A.1. In this description, the equations of motion reduce to non-local differential equations in one spatial dimension $\unicode[STIX]{x1D709}$ and time $t$ . This setting is very convenient for simulating numerically the potential theory equations taking advantage of the properties of harmonic functions in a strip.

The following numerical results illustrate the wave overturning; details of the numerical method are presented in § A.2. Figure 1 shows a familiar overturning wave profile above a flat rigid bottom at $y=-1$ . Of particular interest at this stage, we demonstrate the surface compression near the tip of a breaker displayed by numerical markers. We have chosen the initial profile $y=0.35\cos x$ and the velocity potential at the surface $\unicode[STIX]{x1D711}=(0.35/\sqrt{\tanh 1})\sin x$ , motivated by linear theory. We will use this specific initial profile for all numerical simulations throughout the paper. We also performed simulations for different initial conditions (not reported here) and observed qualitatively the same results for all aspects discussed in this work.

Figure 1. Profile of a breaking wave over a flat bottom $y=-1$ at three different times: $t=0$ , $2.6$ and $3.35$ . Dot markers correspond to material points, which are distributed at equal distances at the initial time; for a better visualization we display only a few markers. The free surface is strongly compressed at the overhanging tip, as seen by the increasing marker density, while it gets stretched at the front slope on the right.

Figure 2. Temporal dependence of the minimum curvature radius $R$ along the free surface; logarithmic vertical scale.

The curvature of the profile displayed in figure 1 increases rapidly at the overhanging tip. The plot of its minimal curvature radius as a function of time is presented in a logarithmic vertical scale in figure 2, demonstrating the nearly exponential decrease at later times. Similar solutions were reported in many earlier studies, such as, for example, in Baker, Meiron & Orszag (Reference Baker, Meiron and Orszag1982), Peregrine (Reference Peregrine1983), Grilli & Svendsen (Reference Grilli and Svendsen1990), Baker & Xie (Reference Baker and Xie2011). There exist initial conditions that can be rigorously tracked until a splash singularity (e.g. intersection of the wave tip with the bottom) occurs in finite time (Castro et al. Reference Castro, Córdoba, Fefferman, Gancedo and Gómez-Serrano2012). But this strong overturning is not the main goal of our work. Considering the incipient breaking wave as the large-scale underlying flow, our focus here will be on the study of much shorter small-amplitude ripples evolving on the surface of such steepening wave profiles (see figure 3).

Figure 3. Small ripple (solid blue line) travelling on top of the unperturbed wave profile (dotted line). Motion of the ripple is approximately described by the sum of the unperturbed flow speed $U$ and the relative group/phase speeds.

Small-amplitude ripples are described, as a first approximation, with equations linearized about the time-dependent unperturbed wave solution. Let us introduce a local (arc length) coordinate $s$ along the surface, which defines the surface spatial coordinates as $x(s,t)$ and $y(s,t)$ . Then, it is convenient to consider ripple perturbations in the form of a slowly modulated wavetrain (Bretherton & Garrett Reference Bretherton and Garrett1968; Peregrine Reference Peregrine1976) with frequency $\unicode[STIX]{x1D714}(s,t)$ , a carrier wavenumber $k(s,t)$ and amplitude $a(s,t)$ . The regime of interest is such that these ripple parameters may vary with position and time, where appreciable changes are observed after many periods ( $2\unicode[STIX]{x03C0}/\unicode[STIX]{x1D714}$ ) and many wavelengths $(2\unicode[STIX]{x03C0}/k)$ . Then, in the first approximation, the underlying flow due to the wave steepening is locally constant with respect to the ripple, but accelerating in time. The frequency and wavenumber of the ripple can be derived from the phase function $\unicode[STIX]{x1D703}(s,t)$ by

(2.7a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D714}=-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}t},\quad k=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}s}.\end{eqnarray}$$

If $U(s,t)$ is the velocity at the fluid surface with respect to coordinate $s$ , then, as mentioned in the introduction, one has (Bretherton & Garrett Reference Bretherton and Garrett1968)

(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D714}=Uk+\unicode[STIX]{x1D6FA},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FA}(s,t)$ is the local intrinsic frequency of the modulated Fourier mode in the Lagrangian reference frame.

Since the deep-water approximation can be used for short wavelengths, the ripple’s dispersion relation in a first approximation is given by (Landau & Lifshitz Reference Landau and Lifshitz1987, §12)

(2.9) $$\begin{eqnarray}\unicode[STIX]{x1D6FA}=\sqrt{g_{\ast }k}.\end{eqnarray}$$

Here $g_{\ast }(s,t)$ is the effective (intrinsic) gravity acting on the ripple in the local reference frame, as its background flow is moving towards overturning. The effective gravity is defined as $\boldsymbol{g}_{\ast }=\boldsymbol{g}-\boldsymbol{a}$ , where $\boldsymbol{a}=(\text{D}\boldsymbol{v}/\text{D}t)=-(\unicode[STIX]{x1D735}p)/\unicode[STIX]{x1D70C}+\boldsymbol{g}$ is the material acceleration at the surface, under the Euler equations (2.1). As a result, we have

(2.10) $$\begin{eqnarray}\boldsymbol{g}_{\ast }=\frac{\unicode[STIX]{x1D735}p}{\unicode[STIX]{x1D70C}}=-g_{\ast }\boldsymbol{n},\end{eqnarray}$$

where $\boldsymbol{n}$ is the surface normal vector. Recall that the pressure is constant at the free boundary and, thus, its gradient is orthogonal to the surface. For the numerical results in figure 1, the value of $g_{\ast }$ along the wave profile at different times is shown in figure 4(a) with the time dependence of its minimum and maximum presented in figure 4(b). At the final time, $g_{\ast }$ varies from the minimum $g_{\ast }\approx 0.2g$ at the wave tip to the maximum $g_{\ast }\approx 3g$ at the foot of the wave. The value of $g_{\ast }$ remains positive at all times in our simulation. As proved by Wu (Reference Wu1997), the positivity of $g_{\ast }$ holds as long as the interface is non-self-intersecting. Therefore we are in the stable Rayleigh–Taylor regime. The local phase speed $c_{p}$ and group speed $c_{g}$ are defined as

(2.11a,b ) $$\begin{eqnarray}c_{p}=\frac{\unicode[STIX]{x1D6FA}}{k}=\sqrt{\frac{g_{\ast }}{k}},\quad c_{g}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}{\unicode[STIX]{x2202}k}=\frac{1}{2}\sqrt{\frac{g_{\ast }}{k}},\end{eqnarray}$$

in the Lagrangian reference frame.

Figure 4. (a) Profiles of the effective gravity $g_{\ast }$ (in units of $g$ ) along the water surface at different times; times and selected markers are the same as in figure 1. (b) Minimum and maximum of the effective gravity at the water surface as functions of time.

3 Wave action and the adiabatic approximation for ripple evolution

The consistency conditions for the second derivatives of the phase in (2.7) yield the relation (Bretherton & Garrett Reference Bretherton and Garrett1968)

(3.1) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}k}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}{\unicode[STIX]{x2202}s}=0,\end{eqnarray}$$

which can be written using (2.8) and (2.11) as the conservation law

(3.2) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}k}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}[(U+c_{p})k]=0.\end{eqnarray}$$

In the adiabatic approximation, i.e. when the temporal and spatial scales of the ripple are much smaller than the scales of the underlying flow (here the unperturbed wave), one also has the conservation law (Bretherton & Garrett Reference Bretherton and Garrett1968)

(3.3) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\frac{E}{\unicode[STIX]{x1D6FA}}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\left[(U+c_{g})\frac{E}{\unicode[STIX]{x1D6FA}}\right]=0,\end{eqnarray}$$

for the wave action density $E/\unicode[STIX]{x1D6FA}$ . Here $E$ is the ripple energy, which can be obtained by considering a linear wave of amplitude $a$ with the effective gravitational acceleration $g_{\ast }$ . Considering the time-averaged values over an oscillation period and equipartition of the kinetic and potential energies, the local energy density of the ripple is written as (Landau & Lifshitz Reference Landau and Lifshitz1987)

(3.4) $$\begin{eqnarray}E={\textstyle \frac{1}{2}}\unicode[STIX]{x1D70C}g_{\ast }a^{2}.\end{eqnarray}$$

Note that the energy of the entire system is conserved, which means that the adiabatic changes just described reflect the energy exchange between the large-scale motion of the overturning wave and small-scale ripples on the water surface.

The two conservation laws (3.2) and (3.3) with the expressions (2.9), (2.11) and (3.4) define the evolution of the local wavenumber $k$ and amplitude $a$ of the ripple. We will now show that these equations can be solved approximately for the small ripples, when the ripple wavelength $\ell =2\unicode[STIX]{x03C0}/k$ is considered small compared to the scale of the unperturbed nonlinear wave. As it follows from (2.11), such an assumption implies that the ripple phase and group velocities are small compared to the local flow speed,

(3.5a,b ) $$\begin{eqnarray}c_{p}\ll U,\quad c_{g}\ll U.\end{eqnarray}$$

Hence, equations (3.2) and (3.3) in a first approximation reduce to the form

(3.6a,b ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}k}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}(Uk)=0,\quad \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}\left(\frac{E}{\unicode[STIX]{x1D6FA}}\right)+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}\left(U\,\frac{E}{\unicode[STIX]{x1D6FA}}\right)=0.\end{eqnarray}$$

Recall that $U(s,t)$ in these equations is the local flow speed on the surface of unperturbed steepening wave; for the linearized formulation, it is not affected by a small-amplitude ripple motion.

Both equations in (3.6) have the form of the continuity equation

(3.7) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}(U\unicode[STIX]{x1D70E})=0.\end{eqnarray}$$

Consider

(3.8a,b ) $$\begin{eqnarray}t=0:\quad \unicode[STIX]{x1D70E}(x)\equiv 1\end{eqnarray}$$

to represent an initial uniform marker (material tracer) distribution along the free surface; see figure 1. In this case, the (marker density) function $\unicode[STIX]{x1D70E}(s,t)$ describes the stretching (for $\unicode[STIX]{x1D70E}<1$ ) and compression (for $\unicode[STIX]{x1D70E}>1$ ) of these material markers along the free surface in time. By using (3.6) and (3.7), one can check that

(3.9a,b ) $$\begin{eqnarray}\frac{\text{D}}{\text{D}t}\left(\frac{k}{\unicode[STIX]{x1D70E}}\right)=0,\quad \frac{\text{D}}{\text{D}t}\left(\frac{E}{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FA}}\right)=0,\end{eqnarray}$$

where the $\text{D}/\text{D}t=\unicode[STIX]{x2202}/\unicode[STIX]{x2202}t+U\unicode[STIX]{x2202}/\unicode[STIX]{x2202}s$ is the material derivative. In other words, the following quantities

(3.10a,b ) $$\begin{eqnarray}\frac{k}{\unicode[STIX]{x1D70E}}=\text{const},\quad \frac{E}{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FA}}=\text{const},\end{eqnarray}$$

are invariant along Lagrangian trajectories at the fluid surface. These two quantities, which refer to the local ripple properties, represent approximate (adiabatic) Lagrangian invariants on the free surface. Nonlinear effects of compression are now built-in to expressions (3.10).

The first relation in (3.10) implies that the ripple wavelength $\ell =2\unicode[STIX]{x03C0}/k$ changes proportionally to $1/\unicode[STIX]{x1D70E}$ , i.e.

(3.11) $$\begin{eqnarray}\frac{\ell }{\ell _{0}}=\frac{1}{\unicode[STIX]{x1D70E}},\end{eqnarray}$$

where the zero subscript denotes the initial length value at $t=0$ ; recall that $\unicode[STIX]{x1D70E}_{0}=1$ due to (3.8a,b ). This formula captures the physical feature that the ripple travels along the Lagrangian trajectory, while stretching or compressing according to the material marker’s dynamics. To interpret the second relation in (3.10), recall that $E/\unicode[STIX]{x1D6FA}$ was defined as the wave action density per unit surface length. Therefore, $E/(\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FA})$ is the conserved Lagrangian wave action density, corresponding to the unit surface length at the initial time.

The conserved wave action in (3.10) will capture the explosive instability. It can be written using the dispersion relation (2.9) for the intrinsic frequency and the energy density expression (3.4) as

(3.12) $$\begin{eqnarray}\frac{E}{\unicode[STIX]{x1D70E}\unicode[STIX]{x1D6FA}}=\frac{\unicode[STIX]{x1D70C}g_{\ast }a^{2}}{2\unicode[STIX]{x1D70E}\sqrt{g_{\ast }k}}=\frac{\unicode[STIX]{x1D70C}\ell _{0}^{1/2}}{2^{3/2}\unicode[STIX]{x03C0}^{1/2}}\frac{g_{\ast }^{1/2}a^{2}}{\unicode[STIX]{x1D70E}^{3/2}},\end{eqnarray}$$

where in the last equality we used (3.11) to express $k=2\unicode[STIX]{x03C0}/\ell =2\unicode[STIX]{x03C0}\unicode[STIX]{x1D70E}/\ell _{0}$ . Since the first factor in the last expression of (3.12) is constant, the conservation property (3.10) yields

(3.13) $$\begin{eqnarray}\frac{g_{\ast }^{1/2}a^{2}}{\unicode[STIX]{x1D70E}^{3/2}}=\text{const}.\end{eqnarray}$$

Evaluating the constant from the initial time, one obtains

(3.14) $$\begin{eqnarray}\frac{a}{a_{0}}=\left(\frac{\unicode[STIX]{x1D70E}^{3}g_{\ast 0}}{g_{\ast }}\right)^{1/4},\end{eqnarray}$$

where $a_{0}$ and $g_{\ast 0}$ are, respectively, the values of $a$ and $g_{\ast }$ at $t=0$ . Combining expressions (3.11) and (3.14), we express the ripple steepness (the ratio of height to wavelength) as

(3.15) $$\begin{eqnarray}S=\frac{2a}{\ell }=\frac{2\unicode[STIX]{x1D70E}a}{\ell _{0}}\end{eqnarray}$$

which yields our final formula

(3.16) $$\begin{eqnarray}\frac{S}{S_{0}}=\left(\frac{\unicode[STIX]{x1D70E}^{7}g_{\ast 0}}{g_{\ast }}\right)^{1/4},\end{eqnarray}$$

where $S_{0}=2a_{0}/\ell _{0}$ .

Recall that all relations (3.10)–(3.16) are deduced along Lagrangian trajectories at the water surface. As will be shown numerically, it is remarkable that these rather simple formulas encompass the full action of the changing large-scale wave profile on the passive small-scale ripple with several non-trivial implications. First, the ripple steepness is fully controlled by the surface compression ratio $\unicode[STIX]{x1D70E}$ and the local effective gravity $g_{\ast }$ . Due to the rather large exponent of the term $\unicode[STIX]{x1D70E}^{7/4}$ in (3.16), the marker density function has a strong effect. Going back to § 2, note that the compression ratio strongly varies within a wave during the breaking process; see, e.g. two Lagrangian points very close to the wave tip in figure 1. Second, expression (3.16) does not depend on the ripple wavelength, predicting that all ripples steepen at the same rate. In particular, this justifies the use of formula (3.16) for a ripple in the form of a general short-wavelength modulated perturbation on top of the original wave profile.

Note that the presented approach can be used for other dispersion relations instead of (2.9). For example, $\unicode[STIX]{x1D714}=\sqrt{g_{\ast }k+\unicode[STIX]{x1D6FE}k^{3}/\unicode[STIX]{x1D70C}}$ takes into account the surface tension coefficient $\unicode[STIX]{x1D6FE}$ (Landau & Lifshitz Reference Landau and Lifshitz1987, § 62) and may be useful in the case of wind generated ripples. In particular, for very short capillary ripples, a similar derivation yields $S/S_{0}=\unicode[STIX]{x1D70E}^{5/4}$ , i.e. the adiabatic amplification mechanism holds with a different power law.

4 Numerical results

In the simulations, we consider the ripples in the form of short Gaussian wave packets. According to relations (3.5a,b ), such wave packets follow approximately the Lagrangian fluid trajectories at the surface. Expression (3.16) for the ripple steepness depends only on the unperturbed solution in figure 1 and, thus, can be evaluated at every point of the wave profile. The profiles of the effective gravity $g_{\ast }$ were already shown in figure 4. Numerically computed profiles of the surface marker density $\unicode[STIX]{x1D70E}$ are presented in figure 5, demonstrating the strong compression (large $\unicode[STIX]{x1D70E}$ ) at the wave tip, and some stretching ( $\unicode[STIX]{x1D70E}<1$ ) at the wave foot. Together, these results provide the change in ripple steepness along the wave profile, which is depicted in a logarithmic colour scale in figure 6.

Figure 5. Profile of the surface marker density $\unicode[STIX]{x1D70E}$ versus the corresponding horizontal coordinate $x$ at different times. For $t=3.35$ the graph is multi-valued due to the overhanging wave profile; see figure 1.

Here one observes an explosive growth of ripple steepness by a factor of almost 50 near the overhanging wave tip (red colour). At the foot of the wave the steepness decreases (blue colour) depleting the surface roughness. Note that our results are obtained within the two-dimensional model. The third dimension is not crucial for the dynamics of short-wave ripples, since their speeds in the direction transverse to the wave can be neglected due to relations (3.5a,b ). The steepness increases super-exponentially at later times, as shown in figure 7(a). Such explosive behaviour distinguishes the present adiabatic steepening mechanism from common instabilities featuring an exponential growth.

Figure 6. Colour (in log scale) shows the change of ripple steepness, $S(t)/S_{0}$ , at different points on the wave profile. The steepness decreases at the wave foot (blue) and then increases by almost two orders of magnitude at the wave tip (red). Black curves on a surface show the wave profiles at times $t=0,~0.5,~1,\ldots ,3$ and the final time $3.35$ . Dashed red curves indicate the trajectories of small ripples (centres of Gaussian wave packets) located initially at (I) $x=2$ , (II) $x=2.9$ and (III) $x=3.8$ ; the simulation in case III was stopped at $t=2.575$ . The insets present the shape and steepness of these ripples at the initial time (circle at the bottom) and at different later times (circles at the top), shown with the subtracted background profile and magnified vertical scale 100:1.

Figure 7. (a) Increase of ripple steepness (maximum value within a wave), demonstrating a super-exponential growth at later times; vertical log scale. The ripple steepness increases approximately 50 times. (b) Change of steepness with time for three different Gaussian ripples, see also figure 6. Adiabatic theoretical prediction is compared with numerical simulations. The dashed line corresponds to the simulation of a ripple with a twice larger wavelength.

For verification of our theory, we added small ripples in the form of Gaussian wave packets with steepness $S\approx 0.02$ and wavenumber $k=128$ located at different parts of the initial nonlinear wave profile and repeated our numerical simulations; see § A.3 for numerical aspects of creating such perturbed initial conditions. Trajectories of these packets, computed via the centres of their envelopes, are shown in figure 6 by red dotted lines. Their shapes at different times are displayed in the insets, where the background profile was subtracted and the vertical scale was magnified with the ratio 100:1. One can see that the ripple steepness increases as a combination of two factors: the decrease of ripple wavelength and the increase of its amplitude. Recall that such ripples are small perturbations, which do not affect considerably the overturning wave profile.

Figure 7(b) demonstrates an excellent agreement between the numerical steepness measured at the centre of each packet and the theoretical (adiabatic) prediction (3.16). In case III of figure 7(b), we observe both a decrease of steepness (depleting the surface roughness) at early times followed by its sudden increase by more than one order of magnitude. We also performed the simulation for a Gaussian packet with a twice larger wavelength in case III; see the dashed black line in figure 7(b). This confirms the theoretical prediction given by (3.16), where the steepness ratio does not depend on the ripples’ wavelength.

It is instructive to provide the reader with an example of the corresponding dimensional variables for our simulations: one can take the water depth to be $5$ m, the wave height to be $3$ m and the initial ripples to have a wavelength of $0.25$ m. Thus, the ripples are short compared to the underlying wave, but not necessarily in the regime where surface tension plays a role. Nevertheless, due to the strong compression this will change later in time. Nonlinear effects also become important at later times when the ripples become steep, breaking down our adiabatic model. This should typically happen in the yellow-to-red region around the wave tip in figure 6. Considering as an example the case III, we observe in figure 8(a) that singularities (sharp angles) are about to form at the ripple crests when the steepness gets close to $S\approx 0.18$ . At these points, the curvature radius is decreasing at least exponentially with time, as shown in figure 8(b), suggesting that surface tension becomes inevitably important in the vicinity of the ripple crests. In our simulations, according to figure 8, this can be expected around the time $t\approx 2.6$ , i.e. prior to overturning in figure 1.

Figure 8. (a) Onset of angle formation at ripple crests. Shown is a segment of the wave profile at the late time $t=2.575$ in case III; see figure 6. (b) Time dependence of curvature radius at ripple crests indicated in (a); vertical log scale.

In figure 9 we present numerical results for our model, now taking surface tension into account. In this case, the dynamical (stress balance) condition for the pressure at the free surface yields

(4.1) $$\begin{eqnarray}p=P_{atm}+\unicode[STIX]{x1D6FE}/R,\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}$ is the surface tension coefficient and $R$ is the curvature radius of the free surface. The simulation producing figure 9 was performed with the dimensionless value $\unicode[STIX]{x1D6FE}=6\times 10^{-6}$ , which corresponds to a realistic value of surface tension for breaking waves of moderate height; see § A.1 for details of the numerical method. The combined effect of surface tension and nonlinearity can be seen in the magnified profile presented in figure 9(b). It reveals the secondary ‘ripple breaking’ when the curvature at the ripple crests become large, followed by the generation of the so-called ‘parasitic’ capillary waves (Ceniceros & Hou Reference Ceniceros and Hou1999). Such capillary waves are known to form near the crests of nonlinear gravity waves in a resonant manner (Longuet-Higgins Reference Longuet-Higgins1995), and could be a mechanism for whitecapping (Dyachenko & Newell Reference Dyachenko and Newell2016).

Figure 9. Ripple shape at the late time $t=2.85$ for the simulation with dimensionless surface tension coefficient $\unicode[STIX]{x1D6FE}=6\times 10^{-6}$ . The initial wave packet corresponds to the case III, where we increased by four times the ripple wavelength to improve numerical resolution. (b) Represents the magnified region of the ripple shown with the small red square in (a). The arrow indicates the direction of the (small) group speed of the ripple relative to the fluid. One can see the formation of parasitic capillary waves developing in front of the ripple crests.