2.1. Numerical approach

We use the Gerris software package (Popinet Reference Popinet2003) 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 Reference Wroniszewski, Verschaeve and Pedersen2014), wave breaking kinematics (Pizzo, Deike & Melville Reference Pizzo, Deike and Melville2016; Deike, Pizzo & Melville Reference Deike, Pizzo and Melville2017) and energy dissipation (De Vita, Verzicco & Iafrati Reference De Vita, Verzicco and Iafrati2018). 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 Reference Popinet2009) and gravity is applied using the ‘reduced gravity’ method described by Popinet (Reference Popinet2018) and shown to minimise spurious currents at the interface (Wroniszewski et al. Reference Wroniszewski, Verschaeve and Pedersen2014).

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 Reference Derakhti and Kirby2016; De Vita et al. Reference De Vita, Verzicco and Iafrati2018), 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 Reference Dean and Dalrymple1991) 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 Reference Song and Banner2002)

(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 (Reference Clément1996), 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 Reference Popinet2003). 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 Reference Batchelor1967, (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. Reference Deike, Pizzo and Melville2017; Mostert & Deike Reference Mostert and Deike2020) 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. Reference Derakhti, Kirby, Banner, Grilli and Thomson2020; Fedele, Banner & Barthelemy Reference Fedele, Banner and Barthelemy2020) 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. Reference Banner, Barthelemy, Fedele, Allis, Benetazzo, Dias and Peirson2014). 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. (Reference Barthelemy, Banner, Peirson, Fedele, Allis and Dias2018) 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. Reference Perlin, He and Bernal1996) and viscous numerical (Varing et al. Reference Varing, Filipot, Grilli, Duarte, Roeber and Yates2020) studies. This is in contrast to the inviscid simulations of Barthelemy et al. (Reference Barthelemy, Banner, Peirson, Fedele, Allis and Dias2018), 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. (Reference Derakhti, Kirby, Banner, Grilli and Thomson2020), 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.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 Reference Tulin2007, (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. Reference Derakhti, Banner and Kirby2018)

(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. Reference Drazen, Melville and Lenain2008), 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.