2.1. General equations

For simplicity, the problem described here is considered to be two-dimensional in the vertical plane. The current field is assumed to be steady, constant in the horizontal $(x)$ direction, and to vary linearly with depth,

(2.1)\begin{equation} U(z)=U_0 + S z. \end{equation}

where $U_0$ refers to the current velocity at the free surface, and $S$ corresponds to its variation with depth. Since 3-D effects are neglected in this study, the vorticity $\varOmega$ within the flow is constant, with $\varOmega =S$. It is straightforward to show that such a current, associated with the hydrostatic pressure $P(z) = p_0 - \rho gz$, is a solution of the Euler equations when considering a problem of constant depth, $\rho$ being the density of water, and $g$ the acceleration due to gravity. Thus, wavy perturbations can be assumed in the form of a velocity field $(u(x,z,t),v(x,z,t))$, associated with the pressure field $p(x,z,t)$. The total flow fields are then given by

(2.2)\begin{equation} \left. \begin{gathered} \tilde{u}(x,z,t) = u(x,z,t) + U(z), \\ \tilde{v}(x,z,t) = v(x,z,t) \quad \mbox{and} \\ \tilde{p}(x,z,t) = p(x,z,t).\end{gathered}\right\} \end{equation}

Using this decomposition, the Euler equations reduce to

(2.3)\begin{gather} u_t + (U+u)u_x + v U_z + v u_z ={-}\frac{p_x}{\rho} \quad \mbox{and} \end{gather}

(2.4)\begin{gather}v_t + (U+u)v_x + v v_z + g ={-}\frac{p_z}{\rho}, \end{gather}

which has to be satisfied together with the continuity equation

(2.5)\begin{equation} u_x+v_z=0. \end{equation}

As demonstrated in Simmen (Reference Simmen1984), and more recently in Nwogu (Reference Nwogu2009), the wavy perturbations propagating in such current conditions are irrotational. Indeed, since the second derivative of the background current $U_{zz}$ is zero, the vorticity conservation equation involves no source term, and the vorticity field does not exchange any vorticity with the wavy perturbations. Yet, it should be emphasized that the presence of constant non-zero vorticity in the background flow might be responsible for substantial changes in the dynamical phenomena, even when considering classical scenario of symmetric travelling water waves. An interesting illustration of such changes is provided by the appearance of the so-called Kelvin cat's eye patterns, which correspond to the occurrence of flow reversal within the orbital velocity field. Discussion of these patterns can be found in Constantin, Strauss & Varvaruca (Reference Constantin, Strauss and Varvaruca2016) and Dyachenko & Hur (Reference Dyachenko and Hur2019). However, provided the absence of exchange of vorticity between the mean flow and the wavy flow, we can introduce a velocity potential $\phi (x,z,t)$ from which to derive the perturbation induced velocities ($\boldsymbol {\nabla } \phi =(u,v)$). It is emphasized that the continuity equation (2.5) is automatically satisfied if the velocity potential is a solution of Laplace's equation

(2.6)\begin{equation} {\rm \Delta} \phi=0.\end{equation}

The kinematic free surface condition can also be expressed, and if $( X,Z )$ denotes the location of a fluid particle at the free surface, this condition might be expressed

(2.7a,b)\begin{equation} \frac{\textrm{d} X}{\textrm{d} t}=u \quad \mbox{and} \quad \frac{\textrm{d} Z}{\textrm{d} t}=v - U(\eta) \frac{\partial \eta}{\partial x}, \end{equation}

where $\textrm {d}/\textrm {d} t$ refers to the material derivative $\textrm {d}/ \textrm {d} t=\partial / \partial t + u \partial /\partial x + v \partial /\partial z$ and $Z=\eta (x,t)$.

A streamfunction $\psi$ can also be introduced, so that $( \partial \psi /\partial z, - \partial \psi /\partial x) = ( u,v )$. The Euler equations (2.3) and (2.4) can now be integrated in space, giving

(2.8)\begin{equation} \frac{\partial \phi}{\partial t} + U(z) \frac{\partial \phi}{\partial x} + \frac{\boldsymbol{\nabla} \phi^2}{2} -S \psi +g z ={-}\frac{p}{\rho}. \end{equation}

When applied to the free surface, where the pressure is constant, this equation provides the classical dynamic boundary condition. Introducing the material derivative used in the kinematic condition, this condition reduces to

(2.9)\begin{equation} \frac{\textrm{d} \phi}{\textrm{d} t} + U(\eta) \frac{\partial \phi}{\partial x} - \frac{\boldsymbol{\nabla} \phi^2}{2} - S \psi +g \eta = 0. \end{equation}

At this point, the knowledge of the streamfunction $\psi$ at the free surface is still needed. We note the relationship

(2.10)\begin{equation} \frac{\partial \psi}{\partial \tau}={-} \frac{\partial \phi}{\partial n}, \end{equation}

where $( \boldsymbol {\tau },\boldsymbol {n} )$ refer respectively to the tangential and normal vectors at the free surface. Thus, the streamfunction $\psi$ can be evaluated at the free surface as soon as the normal derivative of the velocity potential is known.

2.2. Numerical solution

The numerical approach used to solve this set of equations is based on a boundary integral element method (BIEM) coupled with a mixed Euler–Lagrange procedure. At each time step, Green's second identity is discretized to solve numerically the Laplace equation (2.6). Thus, the potential and its normal derivative are known numerically, and the streamfunction $\psi$ can be deduced by integration of (2.10) along the free surface. This numerical integration is performed in the up-wave direction, starting from the down-wave end of the basin and using zero as initial value. Then, the time stepping is performed by numerical integration of (2.7a,b) and (2.9) using a fourth-order Runge–Kutta scheme. Full details of the implementation can be found in Touboul & Kharif (Reference Touboul and Kharif2010). The method has already been implemented and used successfully in the framework of focusing wave groups in the presence of uniform current (Touboul, Pelinovsky & Kharif Reference Touboul, Pelinovsky and Kharif2007; Merkoune et al. Reference Merkoune, Touboul, Abcha, Mouazé and Ezersky2013), or uniform vorticity (Touboul & Kharif Reference Touboul and Kharif2016; Kharif et al. Reference Kharif, Abid and Touboul2017).

In these previous studies, extensive testing and validation of the methodology were performed. Nonetheless, given the precision required for the present study, a convergence test was conducted, and the results are presented in figure 1. A nearly breaking, recurrent wave packet was propagated in the absence of vorticity. The maximum value of the parameter $B_x$, as defined by (2.14), was computed. Figure 1 shows the relative error for the results obtained as a function of the number of nodes used to describe the free surface numerically. Given these results, all the simulations conducted within this study involved $1500$ nodes, with a numerical time step of $dt=10^{-2}$. This corresponds to an average of $150$ nodes per wavelength. This choice reduces the relative numerical error to less than $1\,\%$ when computing key parameters of the study.

Figure 1.

Evolution of the relative error on the maximum value of the breaking parameter $B_x$ for the maximum recurrent wave as a function of the number of nodes used in the simulation.

For the final stages of the simulations, a partially automated mesh refinement/time- splitting algorithm was used, leading the local value of ${\textrm {d} x}$ to evolve from ${\textrm {d} x}=5.10^{-2}$ to ${\textrm {d} x}=5.10^{-4}$, and the time step to reduce to $\textrm {d} t=10^{-4}$. Thus, the breaking onset was confirmed for each identified breaking case without any ambiguity. The criterion adopted for confirming breaking onset was the occurrence of a multi-valued surface elevation of the deforming crest, featuring a transient vertical tangent segment.

2.3. Hydrodynamical conditions

The present study considers chirped wave packets, as defined by Song & Banner (Reference Song and Banner2002)

(2.11)\begin{align} X_p (t) &={-} 0.25 A_p \left( 1+\tanh \left( \frac{4 \omega_pt}{N {\rm \pi}} \right) \right) \left( 1 - \tanh \left( \frac{4 \left( \omega_p t - 2 N {\rm \pi}\right)}{N{\rm \pi}} \right) \right) \nonumber\\ & \quad \sin \left( \omega_p \left( t - \frac{\omega_p C_{ch} t^2}{2} \right) \right). \end{align}

In the equation above, $A_p$ refers to the wave packet amplitude, $\omega_p$ the average angular frequency of the waves, and $N$ is the number of waves present in the packet. Each wave packet is propagated numerically, considering various background flow conditions. Since the purpose of this work is to emphasize the role of vorticity on the marginal breaking criterion, the parameter $U_0$ is disregarded, and kept constant at $0$. Thus, the current variability in this study relies on a single parameter, $S$, providing a depth varying background current distribution according to (2.1). Besides, two values of the normalization length scale, $k_p={\rm \pi}$ and $k_p={\rm \pi} /2$ are considered in this study. These values provide wave trains propagating, respectively, in deep water conditions and at intermediate water depth. The values of the associated depth parameters are respectively $kh=2{\rm \pi}$ and $kh={\rm \pi}$. Two carrier wave peak frequencies $\omega _p$ are considered, taking the values $\omega _p=2{\rm \pi}$ and $\omega _p=4{\rm \pi}$ for the sake of normalization. Thus, in both cases, the peak linear carrier wave phase velocity $c_0$ is then fixed, with the value $c_0=2$.

In this study, the values of $S$ investigated correspond to moderate values of the vorticity in relation to the wave orbital velocities. For example, the current velocity associated with the background shear current corresponding to $S = 1.0$ has a velocity at the maximal crest elevation comparable with the orbital velocity of the steepest non-breaking waves.

Finally, the parameters used within the study are $N$, a measure of the bandwidth of the chirped group, and $A_p$, its initial amplitude. For a given choice of parameters ($S, N$), the wave amplitude $A_p$ is varied incrementally in order to capture as accurately as possible the transition between the steepest recurrent wave group and the least steep breaking group. An increment of $\delta A_p = 10^{-3}$ was used to quantify the upper and lower transition boundaries separating recurrent and breaking wave packets. A list of the simulations performed is presented in table 1.

Table 1.

Hydrodynamical conditions and wave group properties simulated.

2.4. Post-processing of the data

For a given set of parameters $(S, N, A_p)$, each simulation produces a set of data describing the free surface evolution. These data are post-processed to obtain accurate estimation of the prescribed quantities. Firstly, a wave detection algorithm is introduced, identifying the time evolution of each carrier wave in the packet. The waves, and more specifically their crests, are defined by their corresponding zero up-crossing and the zero down-crossing abscissa, $x_{zuc}(t)$ and $x_{zdc}(t)$ respectively. Their crest location $x_c(t)$ and its time evolution are easily obtained, together with the wave amplitude $A_c(t)$.

The instantaneous phase velocity is defined as the time derivative of the crest location, $c_x=\textrm {d} x_c/ \textrm {d} t$, for each wave crest. This point is crucial, since the parameters used in this study depend strongly on the accuracy used to compute the local phase velocity. For instance, an elegant procedure introduced by Seiffert, Ducrozet & Bonnefoy (Reference Seiffert, Ducrozet and Bonnefoy2017), based on the Hilbert transform of the free surface, provides a local wavenumber, which they link to a local phase velocity by means of the linear dispersion equation. This approach was initially considered for this study, but while computationally efficient, it turned out to be sensitive to the filtering procedure used to estimate the Hilbert transform. Phase velocities computed through this procedure were found to differ by up to $40\,\%$ in comparison with the direct computations that we adopted.

For kinematics purposes, we introduce a mean crest location, $x_{zc}(t)$, defined by the mean location of the two zero crossings, $x_{zuc}$ and $x_{zdc}$, namely $x_{zc}=( x_{zuc} + x_{zdc} )/ 2$. Its velocity is obtained from

(2.12)\begin{equation} c_{zc}=\frac{\textrm{d} x_{zc}}{\textrm{d} t}. \end{equation}

It is also straightforward to define a local crest steepness $S_c={\rm \pi} A_c/\lambda _c$, where $\lambda _c$ is a characteristic wavelength scale, defined by the distance between the two zero crossings $x_{zuc}$ and $x_{zdc}$.

Following Barthelemy et al. (Reference Barthelemy, Banner, Peirson, Dias and Allis2016), a leaning parameter $L$ was also introduced, defined as

(2.13)\begin{equation} L=2\frac{x_c-x_{zdc}}{x_{zuc}-x_{zdc}} - 1. \end{equation}

This parameter takes the value $1$ for a wave presenting a vertical front on its forward face, and the value $-1$ for a wave presenting a vertical front on its rear face. Finally, accordingly to the definition introduced by Barthelemy et al. (Reference Barthelemy, Banner, Peirson, Fedele, Allis and Dias2018), we considered the ratio between the local energy flux normalized by the local energy density and the local crest speed. On the wave surface, this quantity simplifies to the ratio of fluid speed $u_x$ to crest speed $c_x$. Thus, at the surface, it reads

(2.14)\begin{equation} B_x=\frac{u_x}{c_x}. \end{equation}

As an aside, we note that in the framework of linear theory, this value further reduces to $B_x=S_c$.