# Breakdown along an ocean wave surface prior to overturning

###### Abstract

Generation of small droplets (aerosols) by waves breaking in surf zone exerts influence over a broad range of geophysical and climatological changes in proximity of coastal lines. Droplets may be formed by bubble bursting with jet and film fragmentation, wind tearing from wave crests etc., while theoretical understanding of such processes is still incomplete. Here we disclose a new multiscale mechanism in which small-scale perturbations (ripples) on water surface are controlled by the large-scale wave dynamics in a way analogous to the adiabatic transport in quantum mechanics. Unlike classical hydrodynamic instabilities, this results in a super-exponential steepness amplification for small-scale ripples. We propose that this explosive increase of ripple steepness explains the breakdown of a smooth water surface into a spray with small droplets, commonly observed in ocean waves prior to overturning even in the absence of wind. The developed quantitative theory is shown to be in excellent agreement with numerical simulations.

Key words: wave breaking, free surface, ripples, adiabatic theorem, instability

## 1 Introduction

Wave breaking in surf zone has a significant impact on the aerosol concentration (O’Dowd & de Leeuw 2007) and air-sea flux (Wallace & Wirick 1992) near the coast with important contribution to its biogeochemical cycling and ecosystems. Physics of this phenomenon goes down to small scales, where tiny bubbles and droplets are produced by fragmentation mechanisms (Deane & Stokes 2002) typically related to hydrodynamic instabilities (Villermaux 2007; Dyachenko & Newell 2016) and finite-time singularities (Longuet-Higgins 1983; Zeff et al. 2000). Energy sustaining this process comes from the large-scale dynamics, where a wave crest becomes vertical, curls over and falls onto the trough with a violent impact (Peregrine 1983). Here we call attention that even before the vertical slope is formed, clouds of small droplets are suddenly sprayed away from the water surface (Fig. 1), which is commonly observed independently of wind conditions.

In this paper, we suggest that such breakdown along a smooth water surface prior to overturning is triggered by the super-exponential amplification of small-scale fast perturbations (ripples), which are coupled directly to the large-scale slow wave motion. This mechanism is explained using the analogy with the quantum-mechanical adiabatic transport (Born & Fock 1928; Bohm 1989). The adiabatic theorem describes the evolution of Hamiltonian eigenstates under a slow change of parameters. By extending this theorem to classical linear oscillations, we show that both position and wavelength of a small-scale ripple is approximately “frozen” relative to the material (Lagrangian) fluid surface. Next, we determine the change of ripple amplitude using the adiabatic invariant, which is equal to the ratio between energy and frequency (Landau & Lifshitz 1976). Conservation of this ratio follows from the representation of a ripple as a non-autonomous linear oscillator. Combining these results yields a quantitative theory of ripple evolution (confirmed with accurate numerical simulations), which explains the onset of explosive spray dynamics at the water surface.

## 2 Evolution of breaking wave at large and small scales

Large-scale dynamics of a smooth water surface can be described by the Euler equations

(2.0) |

where the effects of viscosity and surface tension are neglected (Landau & Lifshitz 1987); comments on capillary effects will be made later. Here is the fluid velocity satisfying the incompressibility condition , is the density and is the vector of gravitational acceleration. At the free surface, normal components of fluid and surface velocities must be equal (kinematic condition) and the pressure is atmospheric (dynamic condition). The fluid velocity is tangent to the rigid bottom.

We will focus on the two-dimensional irrotational flow, when the motion is described by a complex potential , which is a holomorphic function of in the fluid domain. Then the flow velocity is expressed as and . Let with be a conformal mapping from a horizontal strip onto the fluid domain at time . Such mapping provides the free surface parametrization as for real . With the method of complex analysis (Dyachenko et al. 1996; Zakharov et al. 2002; Ribeiro et al. 2017) one can describe the flow in the whole fluid domain in terms of just two real functions at the free-surface, and . In this description, equations of motion reduce to nonlocal differential equations in one spatial dimension and time , which we simulated numerically (see Appendix A).

Figure 2 shows a familiar overturning wave profile above a flat rigid bottom at . The initial profile and the velocity potential at the surface are motivated by the linear theory. Similar solutions were reported in many earlier studies, e.g., (Baker et al. 1982; Peregrine 1983; Grilli & Svendsen 1990; Baker & Xie 2011). These results describe well the large-scale wave dynamics, but fail to display the breakdown of an initially smooth water surface into a cloud of small droplets; see Fig. 1. In order to explain this phenomenon, we now turn to the analysis of small-scale ripples evolving on top of the steepening wave profile.

One can introduce two temporal scales distinguishing the relatively slow evolution of unperturbed wave profile from the fast oscillation of small-scale surface perturbations (ripples). For short times (an order of a few ripple oscillations), we can neglect slow changes in the unperturbed profile. Since the deep-water approximation can be used for short wavelengths, the local dispersion relation between the ripple frequency and wavenumber is given by (Landau & Lifshitz 1987, §12)

(2.0) |

Here is the effective “gravity” acting on the ripple in the local reference frame, and the subscript refers to a specific point of the wave surface. The effective gravity is defined as , where is the material acceleration at the surface, for a flow under the Euler equations (2). As a result, we have

(2.0) |

where is the surface normal vector. Recall that the pressure is constant at the free boundary and, thus, its gradient is orthogonal to the surface. The value of remains positive at all times in our simulation in Fig. 2. At the final time, it varies from the minimum at the wave tip to the maximum at the foot of the wave. The group speed is small for small-scale ripples (large ) compared to the large-scale wave speed. Therefore, ripples are approximately convected along the material (Lagrangian) trajectories on the wave surface.

## 3 Adiabatic theory for ripples along water surface

In this section, we describe how ripple properties change during the process of wave steepening. Because of the large separation of temporal scales, the fast ripple oscillations can be studied under the slow (adiabatic) change of the unperturbed wave profile. Mathematically speaking, small ripples are described by the non-autonomous equations linearized near the slowly changing unperturbed wave solution. In our theory we employ two different techniques to determine separately the changes of ripple wavelength and amplitude.

Linear oscillations in a system with slowly changing coefficients are described by the adiabatic theorem. This theorem is conventionally formulated for quantum systems governed by the Schrödinger equation, where the Hamiltonian operator changes slowly (adiabatically) with time (Bohm 1989). At short times, the solution is given approximately by the eigenstate determined by the eigenvalue (energy) and the corresponding eigenvector of the instantaneous Hamiltonian. The adiabatic theorem states that the system prepared initially in the instantaneous eigenstate will remain in this state if the Hamiltonian depends on time slowly enough, and if there is a gap between the eigenvalue and the rest of the spectrum (Born & Fock 1928). The adiabatic theorem can be extended to more general (non-Hermitian) systems (Moiseyev 2011) and, in particular, to classical linear oscillations as shown in the Appendix B.

The main idea can be demonstrated with a simple example of an oscillating string of length , which has eigenmodes with wavenumbers for harmonics . By the adiabatic theorem, the string with a slowly changing length remains in the same mode, i.e., . Hence, the change of wavenumber is given by

(3.0) |

where is the stretching ratio and the index zero denotes the initial-time value. The validity of the adiabatic theorem can be extended to the case of continuous spectrum (Maamache & Saadi 2008), which corresponds to the limit of large : in this case the relation (3) describes a temporal change of the wavenumber, e.g., for a wave packet. In a similar way, the adiabatic theorem applies to the ripple evolution. This yields the same relation (3), where is the stretching ratio of the water surface computed along a ripple’s (Lagrangian) trajectory. We can write , where the surface arc-length element is divided by its initial value . This means that ripples are approximately “frozen” in Lagrangian coordinates of a wave profile in the sense that their position and wavelength evolve with the material wave surface.

It is not convenient to use the expression for the oscillation amplitude provided by the adiabatic theory; see the Appendix B. For this reason, we employ a different argument to describe the evolution of ripple amplitude. First, we note that our wave equations are Hamiltonian (Dyachenko et al. 1996; Zakharov et al. 2002) and, second, we have argued that each ripple eigenstate behaves as an independent linear oscillator with slowly changing parameters. Systems with such properties possess an approximately conserved adiabatic invariant, which is equal to the ratio between the oscillation energy and frequency (Landau & Lifshitz 1976, §49):

(3.0) |

Considering a ripple locally as a traveling wave of amplitude and wavenumber , the deep water linear theory (Landau & Lifshitz 1987, §12) yields the potential energy per wavelength , with the same value for the kinetic energy. Using the expression for the total energy and the dispersion relation (2) in (3), we obtain the explicit expression for the oscillation amplitude as

(3.0) |

where we evaluated the constant using initial-time values. 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.

Combining expressions (3) and (3), we express the ripple steepness (the ratio of height to wavelength) as

(3.0) |

This rather simple formula encompasses the full action of the changing large-scale wave profile on its small-scale ripples with several nontrivial implications. First, the ripple steepness is fully controlled by the surface stretching ratio and the local effective gravity . Due to the large exponent of the term in (3), the stretching is dominant and the ripple steepness increases when the wave surface compresses. Note that the stretching ratio strongly varies within a wave during the breaking process; see, e.g., two very close Lagrangian points at the wave tip in Fig. 2. Second, expression (3) does not depend on a ripple wavelength predicting that all ripples steepen with the same rate. In particular, this justifies the use of formula (3) for a ripple in the form of a general small-scale modulated perturbation on top of the original wave profile.

It should be mentioned, however, that the above theory is not applicable for very short ripples when capillary effects become important. We can account for capillary effects by using the modified dispersion relation , where is the surface tension coefficient (Landau & Lifshitz 1987, §62). The transition from gravity to capillary waves occurs near the wavelength , which is on the order of two centimeters for the air-water interface. For short capillary waves with , we have and the total ripple energy per unit wavelength for a traveling linear wave becomes . Analogous derivation with (3) and (3) yields , i.e., the same adiabatic mechanism controls the steepness of capillary ripples, but following a different power-law.

## 4 Numerical results

Expression (3) depends only on the unperturbed solution in Fig. 2 and, thus, can be evaluated at every point of the wave profile. Numerical results are shown by color (in logarithmic scale) in Fig. 3, where one observes an explosive growth of ripple steepness by almost two orders of magnitude near the overhanging wave tip (red color). At the foot of the wave the steepness decreases (blue color) depleting the surface roughness. These conclusions can be examined in the Supplementary video demonstrating the ocean wave breaking at low swell conditions. Note that our results are obtained within the two-dimensional model. This implies that, in three-dimensional waves, the ripple steepness changes only in the direction of wave propagation. The steepness increases super-exponentially in time as shown in Fig. 4(a) (black line), where we also presented a similar increase for capillary ripples, (thin blue line). Such explosive behavior distinguishes the adiabatic steepening mechanism from common instabilities featuring an exponential growth.

For verification of our theory, we added small ripples in the form of Gaussian wave packets with steepness and wavenumber located at different parts of the initial profile and repeated our numerical simulations; see Appendix A. Trajectories of these packets are shown in Fig. 3 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.

Fig. 4(b) demonstrates an excellent agreement between the numerical steepness measured at the center of each packet and the theoretical (adiabatic) prediction (3). In case III of Fig. 4(b), we observe both the 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 the Gaussian packet with a twice larger wavelength in the case III; see the dashed black line in Fig. 4(b). This confirmed the independence on wavelength for the ripple evolution predicted theoretically by Eq. (3). It is instructive to provide the reader with an example using dimensional variables: assuming the depth of m, the wave height reaches m and the initial wavelength of ripples corresponding to our numerical tests is cm. Such values, which are typical for ocean waves, demonstrate the minor importance of surface tension at the initial stage of ripple steepening.

Due to the rapidly increasing steepness, small ripples reach a nonlinear regime, which should happen at the yellow-to-red region in Fig. 3. For example, considering the case III, we observe in Fig. 5(a) that singularities (angles) are about to form at the ripple crests when the steepness gets close to . At these points, the curvature radius is decreasing at least exponentially with time; Fig. 5(b). In reality, surface tension becomes important when the curvature gets large, which should eventually lead to a surface fragmentation (Villermaux 2007) observed as an explosive ejection of droplets from the water surface; see Fig. 1 and the Supplementary video. According to Fig. 5, such breakdown of wave surface can be expected around the time , i.e., prior to overturning in Fig. 2.

## 5 Conclusions

We showed that the adiabatic mechanism couples large-scale wave motion to small-scale surface ripples explaining the breakdown of water surface prior to wave overturning. Unlike the common hydrodynamic instabilities, this mechanism predicts super-exponential growth of ripple steepness in time, resulting from simultaneous decrease of wavelength and increase of amplitude, with excellent quantitative agreement between the developed theory and numerical simulations. For waves breaking at the beach, this phenomenon initiates the transition from a relatively simple and well understood wave steepening to a sophisticated dynamics accompanied by droplets ejection and formation of bubbles. Our simulations anticipate the small-scale “ripple breaking” along the water surface revealing the increasing complexity of the subsequent fragmentation process. We expect that our theory describing explosive steepening of surface ripples will stimulate and provide guidelines for further theoretical and numerical studies of highly nontrivial fragmentation stages in the wave breaking phenomenon.

## Appendix A: Numerical method

A dimensionless formulation of the potential theory equations (Whitham 2011) was considered with a unit water depth, density and gravity acceleration. We studied numerically the flow with period in the horizontal direction , over a flat rigid bottom at (Fig. 1). As described in Section 2, the flow equations can be formulated in terms of real functions , and describing the profile and velocity potential at the wave surface with conformal variables . The resulting free surface equations (Dyachenko et al. 1996; Zakharov et al. 2002; Ribeiro et al. 2017) were simulated in the canonical domain using the pseudo-spectral method in space and the fourth-order Runge–Kutta method in time. We used an adaptive uniform grid for the coordinate , starting with and ending with points as the wave overturned. To suppress numerical instability at large wavenumbers, we used the 36th-order smoothing with the Fourier filter suggested by Hou (2009) and applied for high-accuracy numerical simulations of the Euler equations, see also (Agafontsev et al. 2015). Considering different types of Fourier cutoffs, we checked that the numerical results are not affected by the filter. During the simulations, we were doubling the grid size with the Fourier extrapolation of the data, ensuring that the physical part of spectrum (not dominated by round-off errors) is not filtered. We stoped the simulations, when the transition to the grid of points was required. This guaranteed a very high accuracy of our results.

Using the deep-water linear theory (Landau & Lifshitz 1987, §12), we have chosen initial conditions for the small-scale ripples as Gaussian wave packets, and . Such packets have the width , and we considered the wavenumber as providing ripples with a small wavelength . The frequency was determined from Eq. (2). The amplitude was chosen to give a small ripple steepness at . This wave packet was shifted and superimposed on top of the unperturbed initial wave profile at different locations as specified in the caption of Fig. 3. This superposition was obtained using a linear interpolation, e.g., of the unperturbed velocity potential to the perturbed profile. In some simulations, we also used wider Gaussian packets with .

## Appendix B: Adiabatic theorem for linear oscillations

A general linear system of differential equations can be written as

(5.0) |

where is the vector of dimension and is an matrix. Eigenvalues and eigenvectors are given by the equation and define the solutions . Additionally, we consider the left eigenvectors given by with the normalization condition , where denotes the Hermitian transpose. In particular, for an -level quantum system, the Schrödinger equation is written in the form (Appendix B: Adiabatic theorem for linear oscillations) by taking and with the Hamiltonian represented by the Hermitian matrix . In this case, the eigenvalues are purely imaginary with real energies , while the left and right eigenvectors are identical, . When (Appendix B: Adiabatic theorem for linear oscillations) represents small oscillations in a classical system, the matrix is real and the vector contains phase variables (generalized coordinates and speeds). In this case the left and right eigenvectors are different in general, . However, if the system is stable and conservative, the eigenvalues are purely imaginary with real frequencies ; see, e.g., (Seyranian & Mailybaev 2003).

Let us consider oscillations in a time-dependent system with

(5.0) |

where is a given matrix function with and is the total operation time. The adiabatic limit is given by , which means that the matrix elements change slowly with time. At each time, we define the instantaneous eigenvalues and the corresponding eigenvectors and of the matrix . At short times, remains approximately constant and the solution can still be approximated with instantaneous eigenstates . Behavior at long times is described by the adiabatic theorem as we show below.

Assuming that all eigenvalues are purely imaginary and different, we use the ansatz

(5.0) |

with unknown complex coefficients . Here represents the accumulated phase for mode . Following the standard derivation of adiabatic theory, see, e.g., (Bohm 1989, Ch. 20), but distinguishing the right and left eigenvectors, equations (Appendix B: Adiabatic theorem for linear oscillations) and (Appendix B: Adiabatic theorem for linear oscillations) yield

(5.0) |

where . Integrating each term of the sum with respect to in the adiabatic limit , we obtain

(5.0) |

due to the oscillating exponential term (as in the Riemann–Lebesgue lemma). Hence, in the adiabatic approximation, the solution of (Appendix B: Adiabatic theorem for linear oscillations) is given by

(5.0) |

The formulae (Appendix B: Adiabatic theorem for linear oscillations) and (Appendix B: Adiabatic theorem for linear oscillations) confirm the statement of the adiabatic theorem extended to classical linear oscillations. Assume that the system is taken initially at a pure oscillating eigenstate with for given and for . Then the system remains predominantly at this eigenstate at later times. The phase factor in (Appendix B: Adiabatic theorem for linear oscillations) can be written simply as , where is the eigenfrequency. The exponential term in (Appendix B: Adiabatic theorem for linear oscillations) defines the accumulated change of amplitude, apart from an additional phase. Computing this integral requires the knowledge of eigenvectors, which makes it not very convenient for our purposes. We used a simpler “shortcut” in Section 3 deriving the ripple amplitude using the adiabatic invariant.

In a simple example of string oscillations, eigenmodes are given by the functions with the wavenumber , where is the string length and is the integer harmonic number. The corresponding eigenvalues are with the frequencies given by the well-known expression , where is the string wave speed. The adiabatic theorem ensures that the system remains in the instantaneous eigenstate if we change the string length in time slowly enough, which implies the explicit expression for the wavenumber (3). The adiabatic theorem can be extended to oscillations in a system with continuous spectrum. In the case of a string, the continuous spectrum can be viewed as a limit of large string length, , when the wavenumbers become dense on the real axis. The adiabatic theorem remain valid in this limit and can be formulated in terms of wave packets (Maamache & Saadi 2008) leading to the same rule (3) for the change of wavenumber.

This work was supported in part by the CNPq (grant 302351/2015-9) and the Program FAPERJ Pensa Rio (grant E-26/210.874/2014). A.A.M. is grateful to Luca Biferale for useful discussions.

- Agafontsev et al. (2015) Agafontsev, D. S., Kuznetsov, E. A. & Mailybaev, A. A. 2015 Development of high vorticity structures in incompressible 3D Euler equations. Physics of Fluids 27 (8), 085102.
- Baker et al. (1982) Baker, G. R., Meiron, D. I. & Orszag, S. A. 1982 Generalized vortex methods for free-surface flow problems. J. Fluid Mech. 123, 477–501.
- Baker & Xie (2011) Baker, G. R. & Xie, C. 2011 Singularities in the complex physical plane for deep water waves. Journal of Fluid Mechanics 685, 83–116.
- Bohm (1989) Bohm, D. 1989 Quantum theory. Dover Publications.
- Born & Fock (1928) Born, M. & Fock, V. 1928 Beweis des adiabatensatzes. Zeitschrift für Physik A 51 (3), 165–180.
- Deane & Stokes (2002) Deane, G.B. & Stokes, M.D. 2002 Scale dependence of bubble creation mechanisms in breaking waves. Nature 418 (6900), 839–844.
- Dyachenko et al. (1996) Dyachenko, A. I., Zakharov, V. E. & Kuznetsov, E. A. 1996 Nonlinear dynamics of the free surface of an ideal fluid. Plasma Physics Reports 22 (10), 829–840.
- Dyachenko & Newell (2016) Dyachenko, S. & Newell, A. C. 2016 Whitecapping. Stud. Appl. Math. 137 (2), 199–213.
- Grilli & Svendsen (1990) Grilli, S. T. & Svendsen, I. A. 1990 Corner problems and global accuracy in the boundary element solution of nonlinear wave flows. Eng. Anal. Bound. Elem. 7 (4), 178–195.
- Hou (2009) Hou, T. Y. 2009 Blow-up or no blow-up? A unified computational and analytic approach to 3D incompressible Euler and Navier–Stokes equations. Acta Numerica 18, 277–346.
- Landau & Lifshitz (1976) Landau, L. D. & Lifshitz, E. M. 1976 Mechanics. Butterworth–Heinemann.
- Landau & Lifshitz (1987) Landau, L. D. & Lifshitz, E. M. 1987 Fluid mechanics. Pergamon.
- Longuet-Higgins (1983) Longuet-Higgins, M. S. 1983 Bubbles, breaking waves and hyperbolic jets at a free surface. Journal of Fluid Mechanics 127, 103–121.
- Maamache & Saadi (2008) Maamache, M. & Saadi, Y. 2008 Adiabatic theorem and generalized geometrical phase in the case of continuous spectra. Physical Review Letters 101 (15), 150407.
- Moiseyev (2011) Moiseyev, N. 2011 Non-Hermitian quantum mechanics. Cambridge University Press.
- O’Dowd & de Leeuw (2007) O’Dowd, C.D. & de Leeuw, G. 2007 Marine aerosol production: a review of the current knowledge. Phil. Trans. R. Soc. A 365 (1856), 1753–1774.
- Peregrine (1983) Peregrine, D. H. 1983 Breaking waves on beaches. Annu. Rev. Fluid Mech. 15 (1), 149–178.
- Ribeiro et al. (2017) Ribeiro, R., Milewski, P. A. & Nachbin, A. 2017 Flow structure beneath rotational water waves with stagnation points. Journal of Fluid Mechanics 812, 792–814.
- Seyranian & Mailybaev (2003) Seyranian, A. P. & Mailybaev, A. A. 2003 Multiparameter stability theory with mechanical applications. World Scientific.
- Villermaux (2007) Villermaux, E. 2007 Fragmentation. Annu. Rev. Fluid Mech. 39, 419–446.
- Wallace & Wirick (1992) Wallace, D.W.R. & Wirick, C.D. 1992 Large air-sea gas fluxes associated with breaking waves. Nature 356 (6371), 694.
- Whitham (2011) Whitham, G. B. 2011 Linear and nonlinear waves. John Wiley & Sons.
- Zakharov et al. (2002) Zakharov, V. E., Dyachenko, A. I. & Vasilyev, O. A. 2002 New method for numerical simulation of a nonstationary potential flow of incompressible fluid with a free surface. European Journal of Mechanics-B/Fluids 21 (3), 283–291.
- Zeff et al. (2000) Zeff, B.W., Kleber, B., Fineberg, J. & Lathrop, D.P. 2000 Singularity dynamics in curvature collapse and jet eruption on a fluid surface. Nature 403 (6768), 401–404.