An energetic signature for breaking inception in surface gravity waves

Abstract

A dynamical understanding of the physical process of surface gravity wave breaking remains an unresolved problem in fluid dynamics. Conceptually, breaking can be described by inception and onset, where breaking inception is the initiation of unknown irreversible processes within a wave crest that precede the visible manifestation of breaking onset. In the search for an energetic indicator of breaking inception, we use an ensemble of non-breaking and breaking crests evolving within unsteady wave packets simulated in a numerical wave tank to investigate the evolution of each term in the kinetic energy balance equation. We observe that breaking onset is preceded by around one quarter of a wave period by a rapid increase in the rate of convergence of kinetic energy that triggers an irreversible acceleration of the kinetic energy growth rate. This energetic signature, which is present only for crests that subsequently break, arises when the kinetic energy growth rate exceeds a critical threshold. At this point the additional kinetic energy convergence cannot be offset by converting excess kinetic energy to potential energy or by dissipation through friction. Our results suggest that the ratio of the leading terms of the kinetic energy balance equation at the time of this energetic signature is proportional to the strength of the breaking crest. Hence this energetic inception point both predicts the occurrence of breaking onset and indicates the strength of the breaking event.

1. Introduction

The physical process of wave breaking remains one of the classical unresolved problems of fluid dynamics, yet is of fundamental importance for understanding the interaction between the atmosphere and ocean. Wave breaking significantly influences the marine wind drag (e.g. Suzuki, Hara & Sullivan 2013) and generates enhanced turbulence and energy dissipation in the ocean, modifying the ocean boundary layer over significant depths when coupled with other processes such as Langmuir turbulence (Sullivan, McWilliams & Melville 2007). The highly nonlinear nature of the breaking process presents challenges for both observational and numerical studies, prompting a range of approaches to develop an objective diagnostic breaking parameter that is valid for any wave type or water depth. Perlin, Choi & Tian (2013) provides the most recent review of progress in this field and groups diagnostic parameters into three categories that use either the geometric, kinematic or dynamic properties of the wave crest.

More recently, Derakhti et al. (2020) introduced the concept of breaking inception, which describes the initiation of an irreversible process within the crest that leads to breaking and occurs before the instant of breaking onset, i.e. when the first surface manifestation of breaking occurs at the crest. A diagnostic parameter that is able to characterise breaking inception could therefore provide advance warning of a breaking event and potentially also quantify the breaking strength and the energy dissipated thereafter. A breaking inception parameter may also have broad application to the simulation of wave fields in models where individual breaking events cannot be resolved but in which the energetic processes and dynamic consequences of breaking are important to capture accurately.

The breaking inception indicator proposed by Derakhti et al. (2020) is based on the diagnostic parameter $B$ (Barthelemy et al. 2018), which is formally the ratio of the local energy flux to the local energy density normalised by the crest speed $\boldsymbol {c}$. At the interface, this reduces to the ratio of particle velocity to crest speed $\| \boldsymbol {u} \| / \| \boldsymbol {c} \|$. Although this resembles the kinematic breaking criterion in that the value of $B$ at visible breaking initiation is close to unity, it also reveals some remarkable complexities that are yet to be fully understood. Specifically, Barthelemy et al. (2018) found that a threshold value $B_{th} = 0.855 \pm 0.05$ exists, beyond which the crest will always evolve to break. This threshold value was subsequently verified in laboratory and computational studies (Saket et al. 2017, 2018), for a variety of wave packet types (Derakhti, Banner & Kirby 2018), water depths (Seiffert & Ducrozet 2018; Derakhti et al. 2020) and in the presence of a constant shear layer (Touboul & Banner 2021).

To avoid ambiguity, Derakhti et al. (2020) defined breaking inception as the instant at which $B$ first passes through the threshold value $B_{th}$, and characterised breaking onset as the instant when visible breaking first occurs. This threshold, which we shall refer to as the kinematic threshold for breaking inception, also provides information on the strength of the subsequent breaking event. Derakhti et al. (2018) and later Na, Chang & Lim (2020) found that the normalised rate of change of $B$ as it passes through $B_{th}$, known as $\varGamma$, accurately predicts the breaking strength parameter $b$ (Phillips 1985), which has been shown to quantify the energy dissipated through breaking (e.g. Drazen, Melville & Lenain 2008; Deike, Popinet & Melville 2015; Sutherland & Melville 2015).

The use of $B_{th}$ is a robust and useful choice as an indicator of breaking inception as it clearly separates breaking and non-breaking crests and can be determined from measurements at the sea surface. However, a dynamical explanation for this threshold value remains elusive and an explanation for why some waves break and others do not requires further investigation. Given that the energetic definition of $B$ reduces to the kinematic diagnostic $\| \boldsymbol {u} \| / \| \boldsymbol {c} \|$, it is possible that $B$ is a proxy variable that accurately distinguishes breaking from non-breaking waves but does not, in itself, track the underlying dynamical cause for breaking inception.

Dynamic breaking diagnostics have generally focused on the energy growth rate integrated over some region of the crest. Schultz, Huh & Griffin (1994) found that breaking onset could be characterised by the wave integrated potential energy exceeding $52\,\%$ of the total energy of the limiting Stokes wave, which suggests that a departure from the equipartitioning of kinetic and potential energy may be a contributor to breaking inception. Song & Banner (2002) constructed an energy growth rate based on the local depth-integrated total energy and the instantaneous wavenumber. This was shown experimentally to distinguish between breaking and non-breaking waves by Banner & Peirson (2007) but technical challenges with accurately measuring the wavenumber in complex wave packets and a requirement to track the temporal evolution over multiple wave periods have limited its application (Barthelemy et al. 2018). While these dynamic breaking diagnostics show that the energy growth rate is important to the breaking process, they also demonstrate that accurately and consistently measuring the energetics of an evolving nonlinear wave is non-trivial. Integrating the energy over some subdomain of the wave does ameliorate some of these challenges; however, the energy field is highly focused near the crest tip (e.g. Perlin, He & Bernal 1996; Alberello et al. 2018) and the integration process inherently diffuses these local energetic values.

This motivates us to investigate the evolution of the local crest energy field in the time leading up to breaking onset, with the aim of identifying an energetic process that robustly signals breaking inception. Numerical simulation allows us to pursue this in much greater detail than is possible within the constraints of laboratory experiments. We investigate an ensemble of high-resolution numerical simulations of non-breaking and breaking wave crests with a range of wave packet sizes and water depths. We track the point with the largest value of local instantaneous kinetic energy, which occurs near the crest tip, and then derive a balance equation for the kinetic energy at this location. We examine the relative contributions of the individual source and sink terms and find that the convergence of kinetic energy at the crest tip provides a reliable indicator of breaking inception a fraction of a wave period prior to breaking onset. The results also indicate a relationship between the convergence of kinetic energy, the rate of change of kinetic energy and the breaking strength parameter $\varGamma$ of Derakhti et al. (2018).

2. Experimental details

2.1. Numerical approach

We use the Gerris software package (Popinet 2003) to generate a suite of numerical simulations of non-breaking, near-breaking and breaking waves across a range of wave packet configurations and grid refinements. Gerris has been extensively validated for simulations of surface gravity waves (Wroniszewski, Verschaeve & Pedersen 2014), wave breaking kinematics (Pizzo, Deike & Melville 2016; Deike, Pizzo & Melville 2017) and energy dissipation (De Vita, Verzicco & Iafrati 2018). We configure the model to numerically solve the two-dimensional ($\boldsymbol {x} = (x,z)$), incompressible, variable density Navier–Stokes equations, including the effects of viscosity and surface tension

(2.1)$$\begin{gather} \rho \frac{{\rm D} \boldsymbol{u}}{{\rm D}t} =- \boldsymbol{\nabla} p + \rho \boldsymbol{g} + \boldsymbol{f} + \boldsymbol{n} \sigma \kappa \delta_s, \end{gather}$$

(2.2)$$\begin{gather}\frac{{\rm D} \rho}{{\rm D} t} = 0, \end{gather}$$

(2.3)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} = 0. \end{gather}$$

Here, $\rho = \rho (\boldsymbol {x},t)$ is the fluid density, $\boldsymbol {u} = (u,w)$ the fluid velocity, $p$ the pressure and $\boldsymbol {g}$ the gravitational body force. Viscous energy dissipation is characterised by $\boldsymbol {f} = \boldsymbol {\nabla }\boldsymbol{\cdot}\mu (\boldsymbol {\nabla } \boldsymbol {u} + \boldsymbol {\nabla } \boldsymbol {u}^T)$, where $\mu = \mu (\rho )$ is the dynamic viscosity. The magnitude of the surface tension force is a function of the surface tension coefficient $\sigma$ and the interface curvature $\kappa$, with the force localised at the interface by the Dirac delta $\delta _s$ and the interface normal vector $\boldsymbol {n}$. Surface tension is modelled through an improved implementation of the continuum-surface-force approach (Popinet 2009) and gravity is applied using the ‘reduced gravity’ method described by Popinet (2018) and shown to minimise spurious currents at the interface (Wroniszewski et al. 2014).

A two-phase air–water flow is simulated using the volume-of-fluid method, in which the fluid phase is tracked by the conservative tracer $\mathcal {T}$ specifying the fraction of a cell containing water. The limits of $\mathcal {T} = [0,1]$ indicate that a cell contains purely air and water, respectively, while a cell with any intermediate value contains a mixture of the two. We define the location of the air–water interface as the $\mathcal {T} = 0.5$ contour and take the value of scalars on the water side of the interface as the nearest cell to the interface contour for which $\mathcal {T} = 1$.

The model is set up as a two-dimensional ($\boldsymbol {x} = (x,z)$) numerical wave tank in which waves are generated at the left-hand boundary, propagate along the tank and are absorbed at the right-hand boundary (figure 1). Previous studies have reported no significant difference in the integrated wave energetics between two- and three-dimensional simulations (Derakhti & Kirby 2016; De Vita et al. 2018), so limiting our study to two-dimensional simulations allows us to examine a wide range of parameters over a large ensemble within computational constraints, while still accurately capturing the energetic characteristics of the waves.

Figure 1. Schematic of our numerical wave tank. Waves are generated at the paddle boundary and travel down the tank in the positive $x$ direction before being absorbed by the numerical sponge layer. A typical chirped wave (enlarged for clarity) is shown, with tank dimensions normalised by the deep-water wavelength $\lambda _p$ derived from the paddle frequency $\omega _p$.

To generate the wave packet, we simulate a bottom-mounted flexible flap paddle by deriving the exact solutions for the velocity and pressure gradient forcing from wavemaker theory (Dean & Dalrymple 1991) and apply these at the fixed boundary. This method removes the necessity of simulating a moving boundary and thereby greatly increases the computational efficiency of our simulations, while still allowing us to generate a fully nonlinear wave packet. The equivalent lateral movement of the paddle is ${<}5\,\%$ of the wavelength in most cases (table 1) so the approximation of a fixed boundary has little effect on the resultant wave packet.

Table 1. Summary of experiments included in this study. The model was configured using a range of mesh refinement levels, wave packet size $N$ and water depth $d/\lambda _p$. For each configuration the amplitude of the paddle $A_p/\lambda _p$ was varied to generate an ensemble of breaking and non-breaking crest cases.

The motion of the paddle $x_p$ with time $t$ follows the chirped packet function (Song & Banner 2002)

(2.4)\begin{align} x_p(t)&=-0.25A_p\left(1+\tanh\left[\frac{4\omega_pt}{N{\rm \pi}}\right]\right) \left(1-\tanh\left[\frac{4\left(\omega_pt-2N{\rm \pi}\right)}{N{\rm \pi}}\right]\right)\nonumber\\ &\quad \sin\left(\omega_p t\left[1-\frac{\omega_pC_{ch}t}{2}\right]\right) \end{align}

where $x_p$ is a function of the paddle forcing amplitude $A_p$, the forcing frequency $\omega _p$, the number of waves in the paddle signal $N$ and the packet linear chirp rate $C_{ch}=1.0112 \times 10^{-2}$.

The numerical wave tank is configured in non-dimensional coordinates scaled by the linear deep-water wavelength $\lambda _p = 2{\rm \pi} g / \omega _p^2$ and period $T_p = 2{\rm \pi} / \omega _p$ associated with the paddle forcing frequency $\omega _p$. The height of the tank is $1.18\lambda _p$ with a total length of $23.5\lambda _p$, the final $4.7\lambda _p$ being configured as a numerical sponge layer. These dimensions allow the wave packet to evolve over at least $18T_p$ after entering the tank, with wave breaking onset generally occurring within half of this time interval.

Energy absorption at the far end of the tank is achieved through a number of complementary approaches. The final $4.7\lambda _p$ of the tank consists of a numerical sponge layer based on that derived by Clément (1996), which effectively absorbs high-frequency waves. The reflection of low-frequency waves is minimised by gradually increasing the grid spacing within the sponge layer to enhance numerical dissipation. An outflow boundary condition is also applied to the dry portion of the lateral boundary to minimise compression of the air phase caused by the paddle motion, which further improves the performance of the model's Poisson solver.

Gerris uses a quadtree mesh structure that enables efficient adaptive mesh refinement (Popinet 2003). Each level of refinement divides the parent cell into four, resulting in a maximum resolution equivalent to a uniform mesh with $2^n \times 2^n$ grid cells, for $n$ refinement levels. As our primary interest in this study is focused on the air–water interface and the water boundary layer, we determine the maximum required resolution based on the boundary layer thickness $\delta \approx {\lambda _p}/{\sqrt {Re}}$ (Batchelor 1967, (5.7.4)), where $Re=\rho c_p\lambda _p / \mu$ is the wave Reynolds number formulated with the characteristic velocity ($c_p$) and length ($\lambda _p$) scales taken from the paddle signal. To reduce computational cost we set $Re = 4 \times 10^4$ which allows us to resolve the boundary layer with approximately four cells at a refinement level of $2^{10}$ and equates to a resolution of ${{\rm d}\kern0.06em x} = \lambda _p/870$ with the scaling used. While the wave Reynolds number for a physical deep-water gravity wave is $Re \approx 1 \times 10^6$, previous studies (Deike et al. 2017; Mostert & Deike 2020) have shown that $Re = 4 \times 10^4$ is large enough that viscous effects are not dominant and all energy within the boundary layer is adequately resolved. We also conducted a limited number of experiments with a maximum refinement level of $2^{11}$ (approximately eight cells within the boundary layer, equivalent to ${{\rm d}\kern0.06em x} = \lambda _p/1740$) to confirm that the total energy of the simulation did not change (see Appendix A). For all experiments, mesh refinement criteria are configured to ensure maximum resolution at the air–water interface and in regions of large vorticity.

2.2. The crest ensemble

We use our numerical wave tank to conduct a suite of simulations across a range of wave packet configurations, water depths and grid resolutions (table 1). Each individual crest in the wave packet is tracked in space and time (as described below) with the evolution of the crest geometry and energetics recorded. To account for the finite resolution of our numerical simulations, we characterise a crest as breaking if the interface contour exceeds the vertical by a horizontal distance ${\rm d}\eta _x \geqslant 0.5\,{{\rm d}\kern0.06em x}$ over a length ${\rm d}\eta _z \geqslant {{\rm d}\kern0.06em x}$ where ${{\rm d}\kern0.06em x}$ is the finest model grid scale. We use the qualitative term ‘near-breaking’ to describe the steepest and most energetic non-breaking crests for which the local interface contour closely approaches vertical but does not exceed it. The local crest energetics are measured at the location $\boldsymbol {x_+} = [x_+(t), z_+(t)]$ where the local kinetic energy density $E_k$ has its maximum value. A crest reference location and time are set as $[x_0, t_0] = [x_+, t]$ at the instant of breaking onset for breaking crests and at the instant of maximum local $E_k$ for non-breaking crests. The evolution of the crest in space and time is then referenced to these parameters using the non-dimensional coordinates $x^* = (x - x_0) / \lambda _p$, $z^* = z / \lambda _p$ and $t^* = (t - t_0) / T_p$.

Cubic interpolation is used to determine the location of $\boldsymbol {x_+}$ and the value of the energetic quantities at this point, but the unsteady movement of the crest (e.g. Derakhti et al. 2020; Fedele, Banner & Barthelemy 2020) does lead to a level of uncertainty in the resultant time-series. This is managed by filtering the data with a running mean of width $0.15T_p$, which we find effectively removes high-frequency noise without significantly smoothing the temporal variability in the data. Peak values at breaking onset are preserved by applying the filter independently for $t^* \leqslant 0$ and $t^* > 0$ (i.e. before and after breaking onset) with the window width gradually reducing to zero for $\vert t^* \vert < 0.15T_p$. The difference between the original and the filtered data is used to estimate the $5\,\%$, $95\,\%$ confidence interval using a bootstrap method. We use these confidence intervals to objectively discard crests for which the energetic parameters presented in this study are not correctly captured by our analysis methods, which usually occurs when a crest is impacted by droplets from another breaking crest in the same wave packet. Any crests for which the relative magnitude of the confidence intervals exceed the $95$th percentile for that parameter are discarded, which accounted for $7.5\,\%$ of the total. The final ensemble consists of $581$ non-breaking and $73$ breaking crests (table 1).

In figure 2(a), we characterise this ensemble in terms of the local $E_k$ at $\boldsymbol {x_+}$ and the crest steepness $S_c = {\rm \pi}a / \lambda _c$, which captures the unsteady and time-dependant development of the crest in terms of the amplitude $a$ and zero-crossing wavelength $\lambda _c$ (Banner et al. 2014). Figure 2(a) illustrates the distinct energetic characteristics of each wave packet type: for a given $E_k$, deep-water $N=5$ crests are typically steeper than $N=9$ crests, and $N=5$ crests in intermediate depth are even steeper. However, the local $E_k$ is not sufficient to distinguish breaking from non-breaking crests, with a mix of both cases occurring in the range $0.3 < E_k < 0.35$.

Figure 2. Summary of the individual crests included in the ensemble. The model was configured using a range of mesh refinement levels ($2^x$), wave packets $N$ and water depth $d/\lambda _p$. For each configuration the amplitude of the paddle $A_p/\lambda _p$ was varied to generate an ensemble of crests with different energy levels. Non-breaking (breaking) crests are indicated by hollow (filled) symbols. The representative crests shown in subsequent figures are labelled NB1, B1 and B2. (a) The maximum kinetic energy $E_k$ as a function of local wave steepness $S_c$. (b) The value of $B$ at $t^*=0$ (non-breaking crests with $B < 0.8$ not shown); the threshold $B_{th} = 0.855 \pm 0.05$ reported by Barthelemy et al. (2018) is shown by the hatched region.

We also examine the ensemble in terms of the breaking inception parameter $B$. We find that the location $\boldsymbol {x_+}$ of maximum particle velocity and kinetic energy is found on the forward face of the crest, corresponding with previous observational (Perlin et al. 1996) and viscous numerical (Varing et al. 2020) studies. This is in contrast to the inviscid simulations of Barthelemy et al. (2018), who chose to follow the crest tip $\boldsymbol {x_c}$. We find that the location $\boldsymbol {x_+}$ may be offset by up to $0.015\lambda _c$ from $\boldsymbol {x_c}$ and with a particle velocity $\Vert \boldsymbol {u}_+ \Vert$ up to $10\,\%$ greater than the crest tip particle velocity $\boldsymbol {u_c}$ at the time of breaking onset, with these differences generally being larger for more energetic crests. We also note that in a recent $B_{th}$ validation study by Derakhti et al. (2020), the particle velocity was taken as the maximum value within ${\approx }0.03 \lambda _c$ of the crest tip. This motivates us to construct $B$ as

(2.5)\begin{equation} B = \frac{\Vert \boldsymbol{u}_+ \Vert}{ \Vert \boldsymbol{b}_+ \Vert}, \end{equation}

where the denominator $\boldsymbol {b_+}={\rm d}\kern0.7pt \boldsymbol {x_+} / {\rm d} t$ is chosen primarily to be consistent with the mathematical formulation outlined in the following sections. We find that the threshold $B_{th} = 0.855 \pm 0.05$ is replicated by our ensemble (figure 2b) using (2.5). We also found (not shown) that the same threshold value is also observed if the crest speed $\boldsymbol {c}$ is used in the denominator of (2.5), but that $B$ is underestimated if $\Vert \boldsymbol {u_c} \Vert$ is used as the numerator. Because the observed difference between $\Vert \boldsymbol {b}_+\Vert$ and $\Vert \boldsymbol {c} \Vert$ is less than ${\pm }5\,\%$ we see that $B$ is relatively insensitive to the choice of crest velocity, although we do note that $\boldsymbol {b}$ varies more smoothly in time than $\boldsymbol {c}$ as the latter can change rapidly if the crest is impacted by surface ripples that make identification of $\boldsymbol {x_c}$ challenging. With this measure, figure 2(b) shows that our ensemble covers the threshold region between breaking and non-breaking crests with both varying wave packets and water depths.

We have selected representative near-breaking (NB1) and breaking (B1, B2) crests from the deep-water $N=5$ cases that span the $E_k$ and $B$ parameter space (figure 2) and use these in subsequent sections to illustrate the key features of the crest energetics. Some initial observations of the characteristics of these crests can be made from figure 3, in which the evolution of the local $E_k$ field for each crest is shown. The NB1 crest is the near-breaking case in our ensemble that most closely approaches breaking. As the crest grows, a distinct bulge develops on the crest tip and the local interface angle becomes near-vertical. This region is also associated with a local concentration of elevated $E_k$. The B1 case is a weakly breaking crest in which the interface only briefly exceeds the vertical before again relaxing. Conversely, the B2 case illustrates a stronger breaking example in which more extensive overturning of the interface is evident.

Figure 3. The evolution of $E_k$ for the (a) NB1 (near-breaking), (b) B1 and (c) B2 (breaking) representative crests as they progress through the growing and decaying phase. Snapshots are equally spaced at intervals of $0.09T_0$ and the vertical axis is exaggerated by a factor of 7 : 1. The $+$ indicates the location $\boldsymbol {x_+}$ where $E_k$ has its maximum value, and in (d) the temporal evolution of $E_k$ at this location is shown. Periods of active breaking are indicated in (d) by the dotted lines. The time of each snapshot (A–D) corresponds with the vertical dashed lines. Snapshot $C$ occurs at the time that $E_k$ at $\boldsymbol {x_+}$ has its maximum value for the near-breaking crest, and the time of breaking onset for breaking crests.

If only the shape of these crests were considered it could be concluded that the separation between near-breaking and breaking is simply a function of a marginal increase in steepness which eventually leads to the local interface angle exceeding the vertical. However, the values of $E_k$ at the location of the maxima $\boldsymbol {x_+}$ (figure 3d) demonstrate that the energetics of these crests follow diverging paths. Until $t^* \approx -0.15$ the values of $E_k$ are similar in all three cases, but at this point the $E_k$ in the near-breaking case reaches a plateau and begins to gently decrease, while the $E_k$ in the breaking cases undergoes rapid increase up to breaking onset ($t^*=0$) and beyond. The spatial extent of this increase in $E_k$ is seen in the snapshots of the crest evolution (figure 3b,c), with the region of intensifying $E_k$ magnitude fully encompassing the formation of the crest tip bulge. It is evident from this that the convergence of $E_k$ within the crest tip is an important factor in the breaking process. This motivates our analysis presented in the follow sections, in which we mathematically describe and quantitatively track the evolution of this process.

3. Evolution of the crest energetics

3.1. Mathematical formulation

To examine the crest energetics we first construct a balance equation for the local kinetic energy density $E_k = \tfrac {1}{2} \rho \vert \boldsymbol {u}\vert ^2$. This is derived by taking the scalar product of (2.1) with the fluid velocity $\boldsymbol {u}$ and making use of $\boldsymbol {\nabla }\boldsymbol{\cdot}\boldsymbol {u} = 0$ to obtain

(3.1)\begin{equation} \frac{{\rm D} E_k}{{\rm D}t} = \frac{\partial E_k}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} E_k =-\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u} p - \rho g w + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{f} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n} \sigma \kappa \delta_s, \end{equation}

where $w$ is the vertical particle velocity. The terms on the right-hand side account for the work against pressure, the work against gravity, viscous energy dissipation and the surface tension force, respectively.

A similar balance equation for the potential energy density $E_p = \rho g \zeta$ is derived by multiplying the mass-balance equation (2.2) by $\boldsymbol {g}$ and the particle vertical displacement $\zeta$ to give

(3.2)\begin{equation} \frac{{\rm D} E_p}{{\rm D}t} = \frac{\partial E_p}{\partial t} +\boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} E_p = \rho g w. \end{equation}

The $\rho g w$ term links (3.1) and (3.2) and quantifies the rate of conversion between $E_k$ and $E_p$. That is, an upward advection of water particles converts $E_k$ to $E_p$, with the reverse process occurring for downward advection. The net rate of change of $E_k$ will then depend on the sum of $\rho g w$ and the remaining terms on the right-hand side of (3.1). Because the $\rho g w$ term depends only on the vertical component of the particle velocity, the efficiency of this energy conversion process relative to the other terms in (3.1) depends on the direction of the velocity vector $\boldsymbol {u}$, which varies in space and time as the wave evolves.

The Lagrangian balance equations (3.1) and (3.2) provide useful insights into the energetics following a fluid parcel; however, their interpretation as a function of time is not useful in this context as the fluid particles follow an orbital motion that does not correspond with the geometric evolution of the crest. A more insightful approach results from defining a location that moves and evolves with the crest geometry. This may be a specific location such as the highest point of the crest $\boldsymbol {x_c}$, or the location where a scalar quantity has its maximum value. In our case, we choose to follow the location $\boldsymbol {x_+}$ of maximum $E_k$ as it is situated on the forward face of the crest tip, where breaking onset is first observed (figure 3). This location also varies smoothly in time, which aids the interpretation of the kinetic energy evolution. When following the location $\boldsymbol {x_+}$, the change in $E_k$ has both a local and a convective component described by the operator ${\rm D}_b E_k /{\rm D}t = \partial E_k/ \partial t + \boldsymbol {b_+} \boldsymbol{\cdot}\boldsymbol {\nabla } E_k$ (Tulin 2007, (2.2)), where $\boldsymbol {b_+}={\rm d}\kern0.7pt \boldsymbol {x_+} / {\rm d} t$. Henceforth we refer to the energetic values evaluated at this location as local values. When applied to (3.1) this leads to

(3.3)\begin{equation} \frac{{\rm D}_b E_k}{{\rm D} t}= \underbrace{- \boldsymbol{\nabla}\boldsymbol{\cdot} \left(\boldsymbol{u} p + [\boldsymbol{u-b_+}]E_k \right)}_{\text{$CON${}}} \underbrace{- \rho g w}_{\text{$K2P${}}} + \underbrace{\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{f}}_{\text{friction}} + \underbrace{\boldsymbol{u\ {\cdot}\ n} \sigma \kappa \delta_s}_{\text{sfc tension}} \end{equation}

where the relevant kinetic energy flux velocity is seen to be $\boldsymbol {u - b_+}$. This result can also be derived by considering the rate of change of kinetic energy within an arbitrarily small control volume (see Appendix B, where an application of this approach is shown explicitly). We subsequently refer to the terms on the right-hand side as the convergence term ($CON$), the kinetic to potential energy conversion term ($K2P$), friction and surface tension. The use of (3.3) allows us to track these energetic quantities and relevant source terms at our location of interest on the evolving crest. While the magnitude of kinetic energy does not discriminate between breaking and non-breaking crests (figure 2), we show in the subsequent section that by tracking the location of maximum $E_k$ we also track the location where the leading terms in (3.3) have their maximum values.

3.2. Temporal evolution

We examine ${\rm D}_b E_k / {\rm D} t$ and its components in figure 4, where each term of (3.3) is represented as a kinetic energy source/sink by applying the appropriate sign to the values, so that positive (negative) values represent an increase (decrease) in kinetic energy.

Figure 4. The rate of change of the local kinetic energy density $E_k$ at the location $\boldsymbol {x_+}$ (figure 3) and its contributing terms from (3.3) for the representative near-breaking (NB1) and breaking (B1, B2) crests. Positive (negative) values indicate a source (sink) of $E_k$. Periods of wave breaking are indicated by the solid dots.

For all cases, the kinetic energy generally increases (${\rm D}_b E_k / {\rm D} t > 0$) to reach a peak value around $t^* \approx 0$ (i.e. when $E_k$ reaches its maximum value for the non-breaking cases or when breaking onset occurs for breaking cases) before decreasing (${\rm D}_b E_k / {\rm D} t < 0$) as $t^* > 0$. The evolution of ${\rm D}_b E_k / {\rm D} t$ is dominated by the convergence term $CON$ and the kinetic to potential energy conversion term $K2P$, with surface tension of negligible magnitude and friction significant only near $t^*=0$. The $CON$ and $K2P$ terms are of nearly equal magnitude and opposite sign, such that when one is acting as a source of $E_k$ the other is a sink.

The near-cancellation of these two terms is observed throughout the evolution of the near-breaking case NB1 and is also seen throughout most of the growth phase of the breaking cases B1, B2. But a striking deviation from this balance develops as the B1 and B2 crests approach breaking onset. From $t^* \approx -0.1$, a rapid increase in $CON$ is seen which is not balanced by a corresponding increase in $K2P$. As a result, $E_k$ also increases rapidly up to breaking onset. Unlike the $CON$ term, the $K2P$ term for all three example crests are of similar magnitude and trajectory. This indicates that for a given crest there is a limit to the amount of $E_k$ that can be converted to $E_p$. The excess convergence of kinetic energy results in the development of the visible $E_k$ hotspot in figure 3 and the subsequent breaking event.

The interplay between $CON$ and $K2P$ is examined in detail for the near-breaking crest NB1 (figure 5) and breaking crests B1, B2 (figures 6–7) by exploring their spatial variation within the evolving crests. The flux of kinetic energy within the crest is driven by the $\boldsymbol {u-b}$ vectors (figures 5a, 6a and 7a), which follow the approximate shape of the wave, decelerating as it moves upward and rearward from the forward side of the wave to the crest tip and then accelerating down the rear face. This flux leads to a convergence of kinetic energy on the forward side of the crest and divergence on the rear side. In the NB1 case, the local magnitude of the $CON$ field is mostly offset by the $K2P$ term (figure 5c), so that the local rate of change of kinetic energy ${\rm D}_b E_k / {\rm D}t$ is near-zero across most of the crest (figure 5).

Figure 5. The evolution of the kinetic energy balance for the NB1 crest. The dominant terms of the local balance equation (3.3) are (a) $CON$, (b) $K2P$ and (c) ${\rm D}_b E_k / {\rm D}t$. Vectors in (a) show the relative magnitude of the $\boldsymbol {u - b}$ flux velocity at each snapshot. The $+$ indicates the location $\boldsymbol {x_+}$ where $E_k$ has its maximum instantaneous value. Snapshots are equally spaced at intervals of $0.09T_0$ and the vertical axis is exaggerated by a factor of 7 : 1, with snapshot $C$ occurring at the time that $E_k$ at $\boldsymbol {x_+}$ reaches its maximum value. The temporal evolution of each term in (a), (b) and (c) at $\boldsymbol {x_+}$ have been taken from figure 4 and shown in (d) for comparison.

Figure 6. As for figure 5 showing the B1 crest. Periods of wave breaking are indicated by the solid dots.

Figure 7. As for figure 5 showing the B2 crest. Periods of wave breaking are indicated by the solid dots.

In contrast, the $CON$ fields in the B1 (figure 6) and B2 (figure 7) cases develop a local hotspot on the forward face of the crest tip, which continues to intensify up to breaking onset ($t^* = 0$). However, this is not offset by an equivalent hotspot in the local $K2P$ field, whose characteristics are unchanged from the NB1 case. The imbalance between these two leading terms results in a corresponding intensification in the local rate of change ${\rm D}_b E_k / {\rm D}t$.

Because of the large gradients in the evolving ${\rm D}_b E_k / {\rm D}t$ field (figures 5–7) the temporal evolution of this and its contributing terms will be sensitive to the choice of the location $\boldsymbol {x_+}$ that is tracked within the crest. As described in § 2.2, we track the location of maximum $E_k$; this choice was made both for practical reasons as it can be applied in a laboratory setting and because the location varies smoothly in time. But it can also be seen to correspond closely with the location of maximum $CON$ and ${\rm D}_b E_k / {\rm D}t$ (figures 5–7) and therefore captures the maximum intensity of the signal of interest.

We quantify the degree to which ${\rm D}_b E_k / {\rm D}t$ at $\boldsymbol {x_+}$ is representative of the energetics across the crest region by comparing these values with those integrated over a larger region of interest (Appendix B) and find that the rapid increase in ${\rm D}_b E_k / {\rm D}t$ and $CON$ remains observable if the values are integrated over the top $20\,\%$ of the crest for the weakly breaking B1 case and over the top $50\,\%$ for the stronger B2 case. However, if the energetic values are integrated over the full water depth this energetic signature is no longer evident. As well as validating our use of the local energetic quantities, these results also demonstrate the subtlety of this process that is not easily observed in bulk energy values.

3.3. Energy balance and breaking inception

The key feature of the kinetic energy evolution leading up to breaking onset has been shown above to be a breakdown in the approximate equilibrium between the source of kinetic energy $CON$ and the sink $K2P$, with friction becoming significant only near $t^*=0$. This phenomenon is most clearly observed in the local values at the $E_k$ maxima, but is also evident to a lesser extent when integrating over subregions of the crest tip of various sizes (Appendix B). From this we can conclude that the evolution of the local energetic quantities at $\boldsymbol {x_+}$ accurately characterises the broader dynamics of the crest tip region.

The relationship between these terms during the final wave period leading up to $t^*=0$ is further examined in figure 8, an animated version of which is also provided as a supplementary movie available at https://doi.org/10.1017/jfm.2023.134. Here, the grey dashed line denotes an equal balance between the source ($CON$) and sink ($K2P$, $\boldsymbol {u} \boldsymbol{\cdot} \boldsymbol {f}$) terms with any departure above (below) this line indicating a resultant increase (decrease) in $E_k$ (ignoring the insignificant contribution of the surface tension term in (3.3)).

Figure 8. Relationship between the two leading terms in the $E_k$ balance equation (3.3): $K2P$ $- \boldsymbol {u} \boldsymbol{\cdot} \boldsymbol {f}$ and $CON$, for the period $t^*=-1$ to $t^*=0$. Values above (below) the dashed line indicate an imbalance between these terms which leads to an increase (decrease) in ${\rm D}_b E_k / {\rm D}t$. The energetic inception time ($\star$) and the kinematic inception time (i.e. when $B=B_{th}$, indicated by a $\times$) for breaking crests are annotated, with the superposed black lines indicating the period for which $0.83 < B < 0.88$. An animated version of this figure is provided as supplementary material.

In the NB1 near-breaking crest, the initial convergence of kinetic energy is mostly offset by the sink terms, with the magnitude of these terms eventually peaking as the trajectory of the line reverses direction. Throughout this process, the magnitude of each of these terms increases at a similar rate, which keeps the distance from the ${\rm D}_b E_k / {\rm D}t = 0$ line consistent, so that the growth rate of kinetic energy stays within reasonable bounds.

The breaking case B1 initially follows a near-identical trajectory, with the convergence of kinetic energy sufficiently offset by the sink terms to ensure a steady increase in $E_k$. However, at some time after the magnitude of $CON$ begins to decrease, it suddenly experiences a rapid increase (indicated by the $\star$ symbol) that is not balanced by a corresponding increase in the magnitude of the sink terms. As a consequence ${\rm D}_b E_k / {\rm D}t$ also rapidly grows up to breaking onset. A similar evolution is also observed for the stronger B2 breaking case. Here, the imbalance between source and sink terms is larger and so the trajectory is displaced farther from the ${\rm D}_b E_k / {\rm D}t = 0$ line, indicating a faster growth in kinetic energy. As in the B1 case, the magnitude of $CON$ begins to decrease before a striking deflection is seen and $CON$ grows rapidly, driving a subsequent rapid growth in ${\rm D}_b E_k / {\rm D}t$.

In order to test the existence of this inflection point for all cases in our ensemble, we define the critical point $p^\star$ as the final local minimum of the parametric curve $(x,y) = ([K2P - \boldsymbol {u} \boldsymbol{\cdot}\boldsymbol {f}], CON)$ that occurs before $t^*=0$. Without exception, we find that all breaking crests feature a rapid increase in $CON$ commencing at the critical point $p^\star$ (figure 9, solid symbols), which subsequently leads to a rapid increase in ${\rm D}_b E_k / {\rm D}t$ up to breaking onset. This generic feature of the crest evolution therefore represents a critical energy imbalance between the kinetic energy source and sink terms.

Figure 9. The magnitude of $CON$ as a function of (a) $K2P - \boldsymbol {u}\boldsymbol{\cdot}\boldsymbol {f}$ and (b) ${\rm D}_b E_k / {\rm D}t$, at the time that the critical point $p^{\star }$ occurs. The $5\,\%$ and $95\,\%$ confidence intervals are also shown. Breaking crests are coloured by the breaking strength indicator $\varGamma$, and non-breaking crests are coloured grey. The ${\rm D}_b E_k / {\rm D}t$ threshold separating non-breaking and breaking crests is shown by the hatched region.

The occurrence of this critical point is not in itself a sufficient condition for breaking. We observe that small inflections do occur for some non-breaking waves, although these are not followed by a rapid increase in $CON$. In these non-breaking cases (figure 9, grey symbols), the magnitude of the $CON$ term at the critical point may even exceed that of some of the breaking crests and by itself this is clearly not a distinguishing feature between non-breaking and breaking waves. However, the magnitude of ${\rm D}_b E_k / {\rm D}t$ at this critical point does distinguish between the two classes (figure 9b). For non-breaking crests, ${\rm D}_b E_k / {\rm D}t$ is mostly near-zero whereas ${\rm D}_b E_k / {\rm D}t$ increases near-linearly with increasing $CON$ for breaking crests. Moreover, the values of ${\rm D}_b E_k / {\rm D}t$ for non-breaking and breaking crests are clearly separated by a threshold region ${\rm D}_b E_k / {\rm D}t = [0.198, 0.235]$. Therefore, we see that the distinguishing energetic feature separating breaking and non-breaking crests is the occurrence of a critical point in the counterbalance between $E_k$ source and sink terms when ${\rm D}_b E_k / {\rm D}t \geqslant 0.235$. Beyond this threshold value of the kinetic energy growth rate, the sink terms cannot absorb the continuing increase in $E_k$ and the wave passes an energetic point of no return that culminates in breaking onset. This process represents an energetic indicator of breaking inception and we hereafter label this critical point in the crest evolution the energetic signature for breaking inception.

3.4. Features of the energetic signature for breaking inception

We now explore the prognostic characteristics of this new energetic signature for breaking inception and look first at the breaking strength. Up to this point we have characterised this through a visual examination of the interface evolution of our representative breaking crests B1 and B2 (e.g. figure 3), with the B2 crest exhibiting more extensive overturning. To quantify this assessment, we use the breaking strength parameter (Derakhti et al. 2018)

(3.4)\begin{equation} \varGamma = T \left. \frac{{\rm D}_b B}{{\rm D}t}\right|_{B_{th}}. \end{equation}

While other methods of defining the breaking strength exist (e.g. Drazen et al. 2008), we utilise (3.4) as it is conveniently formulated with the same local energetic quantities that we are investigating. The parameter ${\rm D}_b B / {\rm D}t$ is calculated as the average rate of change over the time that $0.83 < B < 0.88$ and the wave period $T$ is defined from the deep-water relationship using the crest zero-crossing wavelength $\lambda _c$. The value of $\varGamma$ for the B1 and B2 crests is $0.85$ and $1.1$, respectively, which aligns with our initial qualitative assessment of breaking strength. In figure 9, where all breaking crests are coloured by the magnitude of $\varGamma$, we see that the largest $\varGamma$ values are associated with the largest values of $CON$, the largest imbalance between $CON$ and $\rho g w - \boldsymbol {u}\boldsymbol{\cdot}\boldsymbol {f}$ and the largest ${\rm D}_b E_k / {\rm D}t$ values at the time that the energetic inception signature occurs. These results demonstrate that the magnitude of these terms at the instant of the energetic inception signature give an indication of the strength of the subsequent breaking event.

The timing of the energetic inception signature is also of particular interest as the identification of breaking inception provides advance warning of breaking onset. The instant at which $B=B_{th}$, which we refer to as the kinematic inception threshold, is shown for the B1 and B2 cases in figure 8 by the $\times$ symbol. The period over which $0.83 < B < 0.88$ has also been coloured black to provide an indication of the rate of change in $B$ around this time. The timing of the energetic inception signature ($\star$) is clearly separated from the kinematic inception threshold and occurs earlier in both examples. This is seen to be the case for all breaking crests in our ensemble (figure 10). While the kinematic inception threshold consistently occurs around $0.05\unicode{x2013}0.1$ wave periods prior to breaking onset regardless of wave packet size or water depth, the energetic inception signature occurs much earlier, up to $0.4$ wave periods prior to breaking onset for our deep-water crests and up to $0.7$ wave periods prior for our shallow water cases.

Figure 10. The energetic inception time ($\star$) and the kinematic inception time ($\times$) for all breaking crests in our ensemble, relative to the time of breaking onset $t^*=0$. Crests are grouped by wave packet size and water depth and the representative breaking crests B1 and B2 discussed in the text are annotated.

4. Discussion and conclusions

Using an ensemble of breaking, near-breaking and non-breaking wave crests simulated with a numerical wave tank, we have examined the evolution of the kinetic energy balance as crests transition from growth to decay. Our results provide new details on the energetic processes leading up to the onset of breaking and the key difference between non-breaking and breaking crest evolution.

Relative to the crest motion, the kinetic energy field is driven by a flux velocity that moves upward and rearward to the crest tip before descending down the rearward face of the wave. On the forward side of the crest, this flux drives a net convergence of kinetic energy, the majority of which is converted to potential energy as the fluid is lifted. In non-breaking crests, the net rate of change of kinetic energy only modestly varies between growth and decay as this convergence of kinetic energy is nearly offset by conversion to potential energy. In breaking crests the rate of kinetic energy convergence is significantly larger than in non-breaking crests, but the rate of conversion to potential energy is similar for all crests. This imbalance leads to a net increase in total kinetic energy, which continues up to breaking onset. Of the remaining terms in the kinetic energy balance equation (3.3), surface tension plays a negligible role and viscous dissipation is only significant on a local scale at the crest tip.

These energetic processes are highly localised and, while the general characteristics can be seen in bulk energetic values, the full detail is only revealed when observing the local values. The maximum kinetic energy, as well as the largest values of kinetic energy convergence and rate of change, were observed to occur on the forward face of the crest tip. We confirmed that the evolution of these terms at this location is representative of the wider crest region by comparing these with the equivalent values integrated over various subregions of the crest. We found that the local variability at this hotspot is still evident even when integrating over the top $20\,\%$ of the crest for weaker breaking cases and the top $50\,\%$ for stronger breaking cases.

These local values highlight the energetic signature that distinguishes breaking crests from their non-breaking counterparts. Throughout the evolution of a non-breaking crest, the rate of change of kinetic energy at the crest tip is bounded by the interplay between the source and sink terms, with the threshold range separating non-breaking and breaking crests in our ensemble found to be ${\rm D}_b E_k / {\rm D}t = [0.198,0.235]$. But in a breaking crest, this threshold is exceeded and any further increase in kinetic energy through convergence can no longer be offset by conversion to potential energy or dissipation through friction. The result is an irreversible and rapid increase in kinetic energy that leads to breaking onset.

Our results show that this energetic signature is a robust indicator of breaking inception. For our ensemble, this typically occurs around $0.25$ wave periods prior to breaking onset, but up to $0.7$ wave periods prior for the shallow water cases that we investigated with $d/\lambda _p = 0.2$. Of fundamental interest is that this energetic inception signature occurs significantly earlier than the kinematic inception threshold based on the transition of $B$ through the value $B_{th} = 0.855$.

A number of questions are left for future studies. We anticipate that the energetic signature of breaking inception presented here will be a consistent feature regardless of wave packet type, water depth or wind forcing, but our ensemble has so far explored only a subset of these variables and further investigation is needed before this can be confirmed. A full energetic explanation for the existence of the kinematic inception threshold $B_{th}$ also remains unresolved, particularly as the time of this threshold is clearly distinct from our energetic inception signature presented here. Finally, while this study has focused on the time period leading up to breaking onset, our results also indicate a clear correlation between the breaking inception point and the strength of the breaking event, which has implications for the amount of energy dissipated. While we show a strong relationship between the kinetic energy convergence, the rate of change of kinetic energy and the breaking strength parameter $\varGamma$, we leave a full examination of this result for future work.