HARMONIC RESONANCE OF SHORT-CRESTED GRAVITY WAVES ON DEEP WATER: ON THEIR PERSISTENCY

Abstract

Three-dimensional short-crested water waves are known to host harmonic resonances (HRs). Their existence depends on their sporadicity versus their persistency. Previous studies, using a unique yet hybrid solution, suggested that HRs exhibit sporadic instability, with the domain of instability exhibiting a bubble-like structure which experiences a loss of stability followed by a re-stabilization. Through the calculation of their complete multiple solution structures and normal forms, we discuss the particular harmonic resonance (2,6). The (2,6) resonance was chosen, not only because it is of lower order, and thus more likely to be significant, but also because it is representative of a fully developed three-dimensional water wave field. Its appearance, growth rate and persistency are discussed. On our converged solutions, we show that, at an incidence angle for which HR (2,6) occurs, the associated superharmonic instability is no longer sporadic. It was also found that the multiple solution operates a subcritical pitchfork bifurcation, so regardless of the value of the control parameter, the wave steepness, a stable branch of the solution always exists. As a result, the analysis reveals two competing processes that either provoke and enhance HRs, or inhibit their appearance and development.

1. Introduction to harmonic resonance of short-crested gravity waves

Short-crested water waves (SCWs) are the genuine three-dimensional waves that exhibit a triangular symmetry. They can be generated by the reflection of a two-dimensional wave train onto a vertical seawall at incidence angle $\theta $ (see Figure 1). The reflection generates three-dimensional, doubly periodic wave patterns of permanent form, which are assumed to be monophasic along the wall (Figure 1). SCWs admit two two-dimensional limits: standing waves for $\theta =0^\circ $ and progressive Stokes waves for $\theta =90^\circ $ . Chappelear [1] then Hsu et al. [4] first calculated the second and third order of deep water SCWs, respectively. Through an asymptotic method, Roberts [10] extended the solutions of Hsu et al. [4] to any arbitrary order. Later, Ioualalen and Kharif [5, 6] studied their stability and estimated the time scales of low-to-moderate wave–wave interactions (resonances, modulations).

Figure 1 An example of a well-developed three-dimensional SCW profile $\eta (X,Y)$ with two wavelengths on each horizontal direction X and Y at angle $\theta =52.3^\circ $ . The direction of propagation of the SCW field is X and the transverse direction is Y. Parameter $\theta $ is the incidence angle of a two-dimensional (infinite-crested) wave train reflecting onto the X-aligned vertical seawall, and $\eta (X,Y)$ is the interface elevation. $Lx$ and $Ly$ are the respective horizontal wavelengths.

Roberts [10] discovered that SCWs are exposed to intrinsic resonances $(I,J)$ of harmonics I and J over horizontal directions X and Y, namely, harmonic resonances (HRs). He showed that, theoretically, there are an infinite number of resonance angles within the $[0^\circ \text {--}90^\circ ]$ $\theta $ -range. HR triggers a zero-divisor at resonance angle $\theta _{HR}$ and a small divisor for angles surrounding $\theta _{HR}$ . For any given angle $\theta $ , the immediate consequence is that the radius of convergence of the SCW power series solutions will drop from the limiting wave height maximum $h_{\max }$ to some value $h_{HR}$ , h being the steepness of SCW (half of the peak-to-trough height). Later, Ioualalen and Okamura [7] used a direct Newton method to calculate SCWs. They showed that HR yields multiple solution structures for some regimes of h, whereas Roberts [10] always obtained a single solution using an asymptotic method.

Ioualalen and Kharif [5] showed that an HR corresponds to a superharmonic instability triggered by a four-wave interaction (a quartet) of class Ia as defined by Ioualalen and Kharif [6]. The calculation of their normal forms led them to suggest that these instabilities are sporadic and should occur over a narrow h-band in the vicinity of $h_{HR}$ . They first appear past a certain value of the wave steepness h, and then the SCW solutions become stable again at a nearby h-value, that is, their so-called “bubble” of instability. As the order of the solutions increases, additional angles of resonance $\theta _{HR}$ emerge, extending to infinity. Since they would likely cover the entire $[0^\circ \text {--}90^\circ ]$ range, every angle could possibly be a resonance angle, and thus the radius of convergence of SCWs would vanish everywhere. Then the existence of SCWs is not ensured. The authors suggested that high order resonant harmonics should be physically damped rapidly through viscosity, and therefore it is physically relevant to ignore them. At low order, a sporadic bubble of instability is likely to statistically inhibit the emergence of HR, while a lack of re-stabilization would favour its emergence and possible growth. The issue is crucial and we want to evaluate this matter in the present study. It is strictly related to the choice of the SCW solutions, that is, whether all branch-solutions are obtained or not. Ioualalen and Kharif [5] used the SCW solutions that were derived from the asymptotic method of Roberts [10]. They are discussed here, in the light of complete solutions that are derived from the Newton method of Okamura [9]. The general formulation of the problem is set up in Section 2. The SCW solutions are described in Section 3 along with an introduction to the phenomenon of harmonic resonances. The stability of “resonant” SCWs is performed in Section 4, then the conclusion follows in Section 5.

2. Mathematical formulation for SCWs and their superharmonic normal forms

We consider surface gravity waves of an inviscid and incompressible fluid on deep water within a nonrotating medium. The motion is assumed to be irrotational, resulting in a velocity potential $\phi (x,y,z,t)$ , where $\eta (x,y,t)$ is the sea surface elevation with respect to the water level at rest. To formulate a generalized flow, we take $1/k$ and $(gk)^{-1/2}$ as our reference length and time, where k is the wave number of the incident wave and g is the gravitational acceleration [4, 10]. A new frame of reference ( $x^*,y^*,z^*,t^*,\phi ^*)$ moving with the wave and defined in [4, 10] is used to enhance the resolution: $x^*=x-ct$ , $y^*=y$ , $z^*=z$ , $t^*=t$ and $\phi ^*=\phi -cx^*$ , where c is the wave velocity equal to $\omega /m$ , $\omega $ is the angular frequency and $m=\sin \theta $ is the wave number in the seawall in the x-direction, $n=\cos \theta $ being the wave number in the transverse y-direction. Once the asterisks are omitted for the sake of simplicity, the flow equations are

(2.1) $$ \begin{align} \phi_{xx}+\phi_{yy}+\phi_{zz} &= 0 \quad\mbox{for } z\leqslant\eta(x,y), \end{align} $$

(2.2) $$ \begin{align} \phi_z &= 0 \quad\mbox{for } z\rightarrow-\infty, \end{align} $$

(2.3) $$ \begin{align} \eta_t+\phi_x\eta_x+\phi_y\eta_y-\phi_z&= 0\quad \mbox{on } z = \eta(x,y), \end{align} $$

(2.4) $$ \begin{align} \phi_t+\eta+\tfrac{1}{2}(\phi_x^2+\phi_y^2+\phi_z^2-c^2) &= 0\quad \mbox{on } z = \eta(x,y), \end{align} $$

where (2.1) is the linear Laplace continuity equation, (2.2) is the linear nonporous bottom condition, and (2.3) and (2.4) are the nonlinear kinematic and Bernoulli dynamical conditions at the free surface. To address both the SCW and the stability issues, we require the general solutions [5, 6], with $t=t^*$ for simplicity, to be

(2.5) $$ \begin{align} \eta(x,y,t) &= \bar{\eta}(x,y)+\eta'(x,y,t), \end{align} $$

(2.6) $$ \begin{align} \phi(x,y,z,t) &= \bar{\phi}(x,y,z)+\phi'(x,y,z,t), \end{align} $$

where $\bar {\eta }(x,y)$ and $\bar {\phi }(x,y,z)$ , the SCW solutions of permanent form, are superimposed on the infinitesimal unsteady perturbations $\eta '(x,y,t)$ and $\phi '(x,y,z,t)$ , with $\eta '\ll \bar {\eta }$ and $\phi '\ll \bar {\phi }$ . After substituting forms (2.5) and (2.6) into the system of equations (2.1)–(2.4), then retaining only linear terms when perturbations and their derivatives are present, and finally ranging all components in the first two powers of the perturbation terms (powers 0 and 1), we obtain two distinct systems of equations: at the zeroth order of the perturbation elements (unperturbed wave field), we obtain the nonlinear SCWs equations for solutions ( $\bar {\eta },\bar {\phi }$ ); and at their first order, we obtain a linear system of equations solving the linear stability problem for solutions ( $\eta ', \phi '$ ), that is, the search of the normal forms [5, 6].

At the zeroth order, the system of equations (2.1)–(2.4) is

(2.7) $$ \begin{align} \bar{\phi}_{xx}+\bar{\phi}_{yy}+\bar{\phi}_{zz} &= 0 \quad\mbox{for } z\leqslant\bar{\eta}(x,y),\nonumber\\\bar{\phi}_z &= 0 \quad\mbox{for } z\rightarrow-\infty,\nonumber\\\bar{\phi}_x\bar{\eta}_x+\bar{\phi}_y\bar{\eta}_y-\bar{\phi}_z&= 0\quad \mbox{on } z = \bar{\eta}(x,y), \end{align} $$

(2.8) $$ \begin{align} \bar{\eta}+\tfrac{1}{2}(\bar{\phi}_x^2+\bar{\phi}_y^2+\bar{\phi}_z^2-c^2) &= 0\quad \mbox{on } z = \bar{\eta}(x,y). \end{align} $$

The wave steepness is defined as the half of the nondimensional peak-to-trough height:

$$ \begin{align*} h=\frac{\bar{\eta}(0,0)-\bar{\eta}(\pi/m,0)}{2}. \end{align*} $$

Two methods have been considered for the computation of the SCW solutions of permanent forms. The first is the asymptotic method, which was proposed by Roberts [10], and was rebuilt and used by Ioualalen and Kharif [5]. The derived solutions are unique and, following Roberts [10], the radius of convergence of their power series is increased beyond singularities (poles) by using the Padé approximation through Shanks’ transform [12]. It is fair to say that the solution is no longer valid around the singularity, at the pole location. Unfortunately, this is precisely the location where Ioualalen and Kharif [5] found their sporadic instability domain. The second method is a direct numerical procedure proposed by Okamura [9] who used a Newton procedure. The method has no technical limitations, so it reveals the entire set of SCW solutions.

At first order, the system of equations (2.1)–(2.4) writes

(2.9) $$ \begin{align} \phi^{\prime}_{xx}+\phi^{\prime}_{yy}+\phi^{\prime}_{zz} &= 0 \quad\mbox{for } z\leqslant\bar{\eta}(x,y), \end{align} $$

(2.10) $$ \begin{align} \phi^{\prime}_z &= 0 \quad\mbox{for } z\rightarrow-\infty, \end{align} $$

(2.11) $$ \begin{align} \eta'(\bar{\phi}_{zz}-\bar{\eta}_x\bar{\phi}_{xz}-\bar{\eta}_y\bar{\phi}_{yz})-\bar{\eta}_x\phi_x^{\prime}-\bar{\phi}_x\eta_x^{\prime} -\bar{\eta}_y\phi_y^{\prime}-\bar{\phi}_y\eta_y^{\prime}+\phi_z^{\prime}-\eta_t^{\prime}= 0\quad \mbox{on } z = \bar{\eta}(x,y), \end{align} $$

(2.12) $$ \begin{align} \phi_t^{\prime}+\bar{\phi}_x\phi_x^{\prime}+\bar{\phi}_y\phi_y^{\prime}+\bar{\phi}_z\phi_z^{\prime} +\eta'(1+\bar{\phi}_x\bar{\phi}_{xz}+\bar{\phi}_y\bar{\phi}_{yz}+\bar{\phi}_z\bar{\phi}_{zz}) =0\quad \mbox{on } z = \bar{\eta}(x,y),\end{align} $$

for which we look for nontrivial superharmonic solutions of the form,

(2.13) $$ \begin{align} \eta'(x,y,t) &= e^{-i\sigma t}\sum_{K=-\infty}^{\infty}\sum_{L=-\infty}^{\infty}a_{KL}e^{i(Kmx+Lny)}, \end{align} $$

(2.14) $$ \begin{align} \phi'(x,y,z,t) &= e^{-i\sigma t}\sum_{K=-\infty}^{\infty}\sum_{L=-\infty}^{\infty}b_{KL}e^{i(Kmx+Lny)} e^{\gamma_{KL}z}, \end{align} $$

where $\gamma _{KL}=[(Km)^2+(Ln)^2]^{1/2}$ . Complex eigenvalues ( $\sigma $ ) and their respective eigenvector components ( $a_{KL}, b_{KL}$ ) are computed using the (spectral) Galerkin method of Ioualalen and Kharif [5, 6]. For a given frequency $\Im \textrm {m}(-\sigma _{K,L})$ , $\Re e(-\sigma _{K,L})=0$ corresponds to a neutrally stable $(K, L)$ perturbation mode, $\Re e(-\sigma _{K,L})<0$ to a stable one and $\Re e(-\sigma _{K,L})>0$ applies for an instability.

3. SCW solutions. A brief review on harmonic resonances (HRs)

The semi-numerical asymptotic method of Roberts [10] and the fully numerical Newton method of Okamura [9] both look for N-truncated doubly periodic permanent solutions of the form,

$$ \begin{align*} \bar{\eta}(x,y)&=\sum_{I=0}^N\sum_{J=2-I.{\mathrm{mod}}(2)}^N\eta_{I,J}\cos(Imx)\cos(Jny),\\\bar{\phi}(x,y,z)&=\sum_{I=0}^N\sum_{J=2-I.{\mathrm{mod}}(2)}^N\phi_{I,J}\sin(Imx)\cos(Jny)e^{\gamma_{IJ}z}, \end{align*} $$

with $\gamma _{JK}=[(Jm)^2+(Kn)^2]^{1/2}$ . Harmonical coefficients $\eta _{I,J}$ and $\phi _{I,J}$ are obtained directly with the numerical method of Okamura [9], which is also detailed by Ioualalen and Okamura [7]. Those of the asymptotic method of Roberts [10] are obtained through successive powers of the small parameter h, previously defined as the wave steepness,

(3.1) $$ \begin{align} \bar{\eta}(x,y)=\sum_{r=1}^N\eta_rh^r , \quad \bar{\phi}(x,y,z)=\sum_{r=1}^N\phi_rh^r , \quad \omega=\sum_{r=0}^N\omega_rh^r. \end{align} $$

The symmetry of the wave propagation at the vertical wall imposes the first-order linear solution to be of the form,

$$ \begin{align*}\omega_0=1, \quad \eta_1=\cos(x)\cos(y) \quad \text{and} \quad \phi_1=\sin(x)\cos(y)e^z.\end{align*} $$

The flow symmetry also forces the general solutions to be of the following form, where harmonics I and J must be of same parity [10]:

$$ \begin{align*} \eta_r&=\sum_{I,J}^r\eta_{r,I,J}\cos(Imx)\cos(Jny),\\\phi_r&=\sum_{I,J}^r\phi_{r,I,J}\sin(Imx)\cos(Jny)e^{\gamma_{IJ}z}. \end{align*} $$

The fully numerical method of Okamura [9] is based on a Newton procedure. The system of equations (2.7) and (2.8) is solved at any collocation point of a computational domain covering one wavelength in both horizontal directions X and Y. The procedure starts with a first-guess solution (FGS) and Newton’s procedure operates several iterations until convergence is reached. Then, the next wanted solution, at a different h, is obtained by taking the previous converged solution as an FGS, that is, the choice of the FGS to solve the problem prescribes the future solution. This iterative procedure allows to calculate a segment/branch of the solution. However, in some circumstances, convergence is not reached, indicating the presence of a singular point. An alternative FGS is then necessary, permitting a jump from one branch solution to another around the singularity. Then, the iterative procedure continues on the second branch. This efficient iterative method may then reveal the multiple branches of the SCW solution.

Roberts [10] was at the beginning of the SCW studies, and revealed most of their characteristics and behaviours, thanks to his semi-numerical method. Harmonics ( $I,J$ ) are ranged and their respective coefficients are computed once for all for a given incident angle $\theta $ . Then, the wave expressions (3.1) are rapidly derived for any chosen value of the expansion parameter h, which is defined as the wave steepness. To any given value of h corresponds a unique solution of the wave height $\eta $ , velocity potential $\phi $ and frequency $\omega $ . It is inherent to a perturbation method [10]. The asymptotic method used by Ioualalen and Kharif [6] is totally inspired by that of Roberts [10]. He prescribed a nonsecular condition to describe the harmonic coefficients as power series of the wave height h instead of the arbitrary expansion parameter $\epsilon $ which has no physical meaning, while Ioualalen and Kharif [6] operate an inversion of the power series of the wave height h, that is, a function of $\epsilon $ , and then replace the parameter $\epsilon $ inside the coefficient expansions, to obtain a power series of h as well. The obtained solutions are then rigorously identical.

The combination of (2.7) and (2.8) leads to the well-known equation of a natural oscillator, the free surface. Roberts [10] showed that its right-hand side may contain, through nonlinear interactions, a term proportional to $\sin (IX)\cos (IJ)e^{\gamma _{IJ}z}$ that is a solution of the homogeneous equation for some critical angles of resonance $\theta _{HR}$ , when $\cos ^2\theta _{HR}=(I^4-I^2)/(J^2-I^2)$ on deep water (Table 1, up to the order of truncation $N=11$ ).

Table 1 Resonance angles $\theta _{HR}(I,J)$ on deep water (in $^\circ $ ), up to order 11 of truncation [10]. The flow symmetry forbids mixed parity for harmonics I and J [10]. Thus, values $J=5, 7, 9, 11, \ldots $ , are not allowed for $I=2$ ; similarly for values $J=10$ and $I=3$ . In deep water, harmonic resonances (HRs) appear at $J=4$ and $J=9$ for $I=2$ and $I=3$ , respectively, at angle $0^\circ $ [10].

To accurately compute HRs, Roberts [10] analysed carefully the convergence of the power series (3.1). One must ensure that the expansion parameter h is smaller than the radius of convergence $R_c$ of the series and eventually increase $R_c$ by some techniques, like the Shanks transform [12], if $h>R_c$ . To illustrate, let us consider angles $\theta =52.3^\circ $ and $\theta =53^\circ $ , which are close to the resonance angle $\theta _{HR}=52.2388^\circ $ (Table 1). At the angle $\theta _{HR}$ , a resonance occurs for a flat surface ( $h=0$ ). At angles $52.3^\circ $ and $53^\circ $ , it occurs at a nonzero wave steepness. Let us proceed with the computation of the d’Alembert criterion [2] for the frequency $\omega (h)$ , that is, the ratio between two successive coefficients of its power expansion. Since $\omega $ retains only even orders of the expansion (3.1) for symmetry reasons [10], it is developed in powers of $H = h^2$ . Thus, for simplicity, we report the ratio $(\omega _{n+2}/\omega _{n})^{1/2}$ in Figure 2.

Figure 2 Computed ratio $(\omega _{n+2}/\omega _{n})^{1/2}$ as a function of the SCW order of truncation n up to order $n = 73$ for angles $\theta =52.3^\circ $ ( $\circ $ ) and $\theta =53^\circ $ $(\vartriangle )$ ; $\omega _{n}$ being the component of the frequency at order $r=n$ (see (3.1)). The asymptotic method of Roberts [10] is used.

At angle $\theta =52.3^\circ $ , where there is no other nearby resonance angles than $\theta _{HR}=52.2388^\circ $ up to order $n=73$ , the ratio converges into the single value $h\approx 0.054728$ (Figure 2). The process begins after order $n=6$ where HR (2,6) appears. According to the criterion, it can be considered as the “local” radius of convergence. The radius of convergence of the power series may also be approximated by using the Padé approximation [3]. The pole of the (17,17) approximant is $h=0.054726$ which is in full agreement.

At angle $53^\circ $ , the ratio drops into successive thresholds due to the presence of nearby resonance angles. The first is the same $\theta _{HR}=52.2388^\circ $ , appearing at order 6 and triggering a first convergence at $h\approx 0.191135$ . The value h is larger than the value for angle $\theta =52.3^\circ $ because angle $53^\circ $ is farther from $\theta _{HR}=52.2388^\circ $ . The ratio then trends to a second value $h\approx 0.0749$ due to HR (4,26) appearing at order 26 ( $\theta _{HR}=52.9133^\circ $ , not shown in Table 1), but operating at higher order; and finally around $h\approx 0.0231$ due to HR (5,41) ( $\theta _{HR}=52.9918^\circ $ ). The successive drops may potentially result in a zero-ratio, which is the effective radius of convergence of all power series (3.1). Furthermore, we may expect that SCW solutions will not converge because of the occurrence of HRs for every angle $\theta $ . Following Roberts [10], Ioualalen and Kharif [5] explored low-order HRs ( $I < 4$ ) by using the Shanks transform [12] to increase their radius of convergence. They ignored high-order resonances ( $I> 3$ ), conjecturing that high-order HRs would be damped through viscosity in real ocean waves.

The coefficient $\phi _{2,6}$ of the resonant component (2,6) is plotted in Figure 3 versus the wave steepness h at the nontrivial angle of resonance $\theta =52.3^\circ $ , which is close to the resonance angle 52.2388 of HR (2,6) (Table 1). The fully numerical Newton method of Okamura [9] and the semi-numerical asymptotic procedure of Roberts [10] are used. Using Newton’s method, branch (3) is obtained starting from a flat surface and increasing the wave steepness h (upper panel of Figure 3, in disks). A “severe” positive jump is then obtained around the critical singularity $h_c\approx 0.0547$ , which is precisely the radius of convergence and the pole described above using the asymptotic method. Below this critical value of h, branch (3) is similar to the solution obtained with the asymptotic method (in blue and red lines in Figure 3): both are weak and positive. The two solutions differ around the pole value. As expected, the asymptotic solution, without any processing, diverges at the singularity (in dashed blue without $\varepsilon $ -algorithm). However, when the radius of convergence is extended using the Shanks transform [12], here the $\varepsilon $ -algorithm ( $\varepsilon $ -alg.), the processed solution operates a slight jump into weak negative values (red line, with $\varepsilon $ -alg.). The passage from positive to negative values happens around the pole. The asymptotic solution cannot be reliable around this point.

Figure 3 (a) Coefficient $\phi _{2,6}$ of the harmonics (2,6) of the velocity potential versus the wave steepness h at the incidence angle $\theta =52.3^\circ $ using Newton’s method. The branches of the solution are numbered (1)–(3), branches (1) and (2) being connected through a turning point (TP). Note branch (2) is negative, although not visible at this scale. Branch (3), spanning from point ( $P_1$ ) to point ( $P_2$ ), is plotted in small disks when stable, large disks when unstable and plus signs when normal form solutions have not converged. Similarly, for branches (1) and (2), there are small and large triangles and plus signs extending from point ( $P_4$ ) to point ( $P_3$ ) through the TP. For reference, there are also plotted solutions obtained with an asymptotic method with (solid red) and without (dashed blue) use of the Shanks transform ( $\varepsilon $ -algorithm). (b) Enlargement around the TP for better visibility. In particular, the passage between stability regimes for branches (1) and (3) and around TP is clearly delineated (sizes of the disks and triangles for ease to read).

Figure 4 The coefficient $\phi _{1,1}$ of the harmonies (1,1) of (a) the velocity potential and (b) the wave frequency $\omega $ versus the wave steepness h at incidence angle $\theta =52.3^\circ $ , and their three solution branches (1), (2) and (3).

Using the Newton method past the pole, the remaining branches (2) and (3) are obtained; thanks to the choice of a judicious first-guess solution in the iterative process. Past the bifurcation at $h_c$ , the complete multiple solution is thus obtained. Branch (3) escapes from branches (2) and (1), while those latter are connected through a turning point (TP). It is interesting to note the similarity between branch (2) and the processed asymptotic solution: both are weak, negative and have the same trend. Around the TP, the complete and the processed asymptotic solutions differ significantly: both solutions operate a jump from positive to negative value, but to a larger extent for the former one.

To summarize, using the fully numerical method, we obtain the effective set of solutions, that is, a single branch prior to the TP, branch (3), and a multi-branch solution subsequent to the TP, which is composed of the continuation of branch (3) and the two connected branches (1) and (2). The asymptotic method, with or without the use of the $\varepsilon $ -alg., first provides the accurate segment of branch (3) for h smaller than the radius of convergence of the coefficient expansion. Then, around the radius of convergence, that is, at the TP, when complemented with the Shanks transform, the transform connects continuously the lower segment of branch (3) (before the singularity) to an approximate of branch (1) after the singularity, resulting in an overall single solution [3, 12]. This single solution can thus be considered as “hybrid”, because the transform operates a matching between two distinct branches. The use of the Shanks transform seems efficient around and past the singularity by increasing the radius of convergence of the power series. However, it still remains an artificial procedure because it uses one particular Padé approximant, bearing in mind that Padé approximants may differ drastically [3]. Moreover, Gilewicz [3] also proposed other alternatives such as his so-called best Padé approximant (BAP) to increase radii of convergence. Thus, the use of the “hybrid” solution is an interesting alternative, but should be avoided when an exact set of solutions is available, (like that of Okamura) [9] which we use hereafter. This is precisely what was done by Roberts [10], by also developing a fully numerical method after his asymptotic one [11].

Once the resonant component (2,6) forces a multiple solution past the TP, the entire set of coefficients of both the wave heights, the velocity potential and the wave frequency behave similarly (Figure 4). The full resonance of the SCW field can thus be described as a “contamination” of the entire SCW spectrum by the resonant harmonics (2,6) through interactions between frequency components that are present in nonlinear terms of (2.7) and (2.8). Once the SCW field is obtained, its mode of occurrence and persistency is further analysed through the computation of its normal forms.

4. Stability of resonant SCWs: their normal forms

We consider the resonance (2,6) on deep water versus the wave steepness h. As described by Ioualalen and Kharif [5], a harmonic resonance $(\pm K,L)$ corresponds to a superharmonic instability that may occur only when the two eigenvalues $\sigma _{K,L}^s=\sigma _{-K,L}^{-s}$ of opposite signatures coalesce. The collision of the two eigenmodes $(\pm K,L)$ is then interpreted as Ioualalen and Kharif’s [6] class Ia $(K,L)$ instability. It corresponds to a resonance between the two eigenmodes $(\pm K,L)$ and the $2K$ modes $(1,\pm 1)$ of the basic SCW of permanent form,

$$ \begin{align*} \Omega_1 &= -\Omega_2+K\Omega_{01}+K\Omega_{02},\\ \boldsymbol{k}_1 &= \boldsymbol{k}_2+K\boldsymbol{k}_{01}+K\boldsymbol{k}_{02}, \end{align*} $$

where $\Omega _i=|{\boldsymbol {k}}_i|^{1/2}$ , $\Omega _{0i}= |{\boldsymbol {k}}_{0i}|^{1/2}$ for $i=1,2$ , ${\boldsymbol {k}}_1=(mK,nL)$ and ${\boldsymbol {k}}_2=(-mK,nL)$ being the wave numbers of the perturbation modes, and ${\boldsymbol {k}}_{01}=(m,n)$ and ${\boldsymbol {k}}_{02}=(m,-n)$ being those of the basic wave.

The stability problem is solved using the procedure of Ioualalen and Kharif [5]. The perturbation forms (2.13) and (2.14) are truncated at order M and the SCWs of permanent form $\bar {\eta }(x,y,z)$ and $\bar {\phi }(x,y)$ are computed using the Newton method of Okamura [9]. Once both expressions are substituted into the first-order system of equations (2.9)–(2.12), we obtain a generalized eigenvalue problem of the form ${\boldsymbol {A}}{\boldsymbol {u}}=i\sigma {\boldsymbol {B}}{\boldsymbol {u}}$ that we solve with the QZ algorithm [8], where $\sigma $ are the eigenvalues and ${\boldsymbol {u}}=(a_{jk},b_{jk})^{T}$ their corresponding eigenvectors. Complex matrices ${\boldsymbol {A}}$ and ${\boldsymbol {B}}$ are functions of the basic flow. The Galerkin numerical procedure of Ioualalen and Kharif [5] is used to solve the eigenvalue problem. Equations (2.11) and (2.12) are numerically integrated over one wavelength in the two horizontal directions using fast Fourier transforms (FFTs) over $\nu \times \mu $ points of coordinates $\alpha x_u=2\pi u/\nu , u=0,\ldots ,\nu -1$ and $\beta y_v=2\pi v/\mu , v=0,\ldots ,\mu -1$ . A generalized eigenvalue problem is obtained at the free surface $z=\bar {\eta }(x,y)$ ,

$$ \begin{align*} &\sum_{K=-M}^{M} \sum_{L=-M}^{M} F_{K-r,L-s}\{E_{KL}^{(1)}\} a_{KL} + \sum_{K=-M}^{M} \sum_{L=-M}^{M} F_{K-r,L-s}\{G_{KL}^{(1)}\}b_{KL} =i\sigma a_{rs}, \\ &\sum_{K=-M}^{M} \sum_{L=-M}^{M} F_{K-r,L-s}\{E_{KL}^{(2)}\} a_{KL} + \sum_{K=-M}^{M} \sum_{L=-M}^{M} F_{K-rL-s}\{G_{KL}^{(2)}\} b_{KL} \\ &\quad=i\sigma \sum_{K=-M}^{M} \sum_{L=-M}^{M} F_{K-r,L-s}\{H_{KL}^{(2)}\} b_{KL}, \end{align*} $$

with

$$ \begin{align*} E_{KL}^{(1)} &=-\bar{\phi}_{zz}+imK\bar{\phi}_x+inL\bar{\phi}_y+ \bar{\eta}_x\bar{\phi}_{xz}+ \bar{\eta}_y\bar{\phi}_{yz}, \\ G_{KL}^{(1)} &=(imK\bar{\eta}_x +inL\bar{\eta}_y-\gamma_{KL}) e^{\gamma_{KL}\bar{\eta}}, \\ E_{KL}^{(2)} &=1+\bar{\phi}_x\bar{\phi}_{xz}+\bar{\phi}_y\bar{\phi}_{yz}+\bar{\phi}_z\bar{\phi}_{zz},\\ G_{KL}^{(2)} &=(imK\bar{\phi}_x +inL\bar{\phi}_y+\gamma_{KL}\bar{\phi}_z) e^{\gamma_{KL}\bar{\eta}}, \\ H_{KL} &=e^{\gamma_{KL}\bar{\eta}}. \end{align*} $$

The functions

$$ \begin{align*} F_{K-r,L-s}\{f_{KL}\}=\sum_{u=0}^{\nu-1}\sum_{v=0}^{\mu-1}f_{KL}e^{i\alpha(K-r)x_u}e^{i\alpha(L-s)y_v}, \end{align*} $$

where $r=-M,\ldots ,M$ and $s=-M,\ldots ,M$ , are computed using two-dimensional FFTs. Convergence of the eigenvalues $\sigma $ is obtained by increasing M. In the following, all computations are carried out up to order $M=49$ if necessary.

In Figure 5, we have plotted the frequencies [ $-\Re \textrm {e}(\sigma _{\pm 2,6})$ ] of the eigenvalues $\sigma _{\pm 2,6}$ along with their growth rates [ $-\Im \textrm {m}(\sigma _{\pm 2,6})$ ] for branches (1) to (3) of the SCW field. Let us first consider branch (3). As expected, as h increases from point (P1) towards the TP (Figures 3 and 5), the frequencies of the harmonics (−2,6) and (2,6) converge to each other while remaining neutrally stable. Then, they coalesce when they both vanish and give rise to an instability beyond the TP. As described by Ioualalen and Kharif [5], a coalescence at zero-frequency is consistent with a resonance, since both harmonic components propagate at the same phase speed as the basic SCW field. They are phase-locked, which is precisely the definition of a resonance phenomenon.

Figure 5 Frequency [ $-\Re \textrm {e}(\sigma _{\pm 2,6})$ ] (dots) and growth rate [ $-\Im \textrm {m}(\sigma _{\pm 2,6})$ ] (+ signs) versus the wave steepness h for angle $\theta =52.3^\circ $ . Positive values of the frequencies and growth rates correspond to eigenvalues $\sigma _{2,6}$ and negative values correspond to $\sigma _{-2,6}$ . (a) Solution branch (3), and (b) branches (2) and (1). Points ( $P_1$ ) to ( $P_4$ ) and TP are the same as in Figure 3.

Beyond a certain value of the wave steepness h, Ioualalen and Kharif [5], using an asymptotic method, observed a re-stabilization. They described the h-range of instability as a “bubble”-like structure with a specific maximum growth rate. In contrast, with our complete SCW solution, and in particular for our branch (3), we do not observe any re-stabilization, at least as far as the normal form problem converges, which is not the case past a certain value of the wave steepness (Table 2). We even find an inflection point from $h\approx 0.10$ enhancing the growth rate instead of reaching a maximum (Figure 5).

For the other branches, starting from point (P3) towards the TP all along branch (2), the frequencies converge to each other and also give rise to an instability around the TP at the beginning of branch (1) (Figures 3 and 5). Within branch (1) from the TP to point (P4), an instability arises with no re-stabilization of the converged solutions (Table 2, Figure 5), in contrast to that observed by Ioualalen and Kharif [5]. Similarly to the unstable segment of branch (3), an inflection point is also observed from the same wave steepness $h\approx 0.10$ , with very similar growth rates.

We propose the following reason for the re-stabilization computed by Ioualalen and Kharif [5], that is, the existence of their “bubble”-like structure. Their hybrid SCW solution operates a jump from the h-left segment of our branch (3) to branch (2). As we show here, both parts are stable. Around the jump, the solution operates a jump from the beginning of the unstable branch (3) to the beginning of the unstable branch (1). The passage automatically requires a narrow h-range of instability, that is, their “bubble”. It is thus due to the hybrid but incomplete structure of their solution.

Ultimately, in our case study, the stability analysis indicates a possible passage from partially unstable segments of branches (1) and (3) onto the always-neutrally stable branch (2) for any given wave steepness h (Figures 3 and 6). Regardless of the wave steepness h, a subcritical bifurcation may always operate from an unstable branch onto a stable one (Figures 3 and 6). There exists always a pre-existing (neutrally) stable branch solution upon which a transition from an unstable segment may occur (Figure 6). Therefore, the physical appearance of HRs is still questionable, despite the possible large growth rate of their associated instability.

Table 2 Convergence of the growth rates $-\Im \textrm {m}(\sigma _{\pm 2,6})$ versus the wave steepness h at incidence angle $\theta =52.3^\circ $ . Three values of the parameter h are taken for each branch (3) ( $h_3$ ) and (1) ( $h_1$ ), $h=0.1386$ and $h=0.1129$ being respectively the highest values where convergence is reached.

Figure 6 Schematic representation of the SCW three-solution branches (1), (2) and (3) for coefficient $\phi _{2,6}$ of the velocity potential versus the wave steepness h for angle $\theta =52.3^\circ $ (Figure 3). Dashed lines apply for unstable segments while solid lines are for neutrally stable ones. Arrows represent possible jumps form unstable to stable branch segments.

5. Conclusion

Short-crested waves are the genuine three-dimensional wave fields that are generated when two two-dimensional progressive waves intersect at a certain incidence angle. It is well known since the pioneering work of Roberts [10] that they host the HR phenomenon that questions their existence, because of the problem of small divisors. Previous studies of Ioualalen and Kharif [5] identified these resonances as a particular case of modulational instabilities that occur at zero-frequency, which means that they are phase-locked with the main wave. They computed their normal forms and suggested that they are weak and sporadic, and thus they are not likely to develop in a water wave field. The asymptotic method they used to compute the solutions has some limitations, especially around and past the singularity due to HR. Their hybrid solution, although unique, is composed of two branch solutions before and beyond the pole of the power series of the wave steepness h that are connected continuously (and artificially); thanks to the use of the Shanks transform. Around the pole, their solution cannot be considered valid, and past the pole, it cannot be reliable because it is constructed artificially. Later, Ioualalen and Okamura [7] computed SCWs with a fully numerical procedure based on the Newton method, and they showed that HRs give rise to multiple solution structures that are composed of three branches.

In the present study, the HR (2,6) is investigated, because it corresponds to a low-order resonance, and thus it is more likely to appear, according to previous studies [5]. It was also selected because it exemplifies a well-developed three-dimensional SCW field. The solutions and normal forms are computed. When converged solutions are obtained, instabilities arise for both branch (1) and a segment of branch (3) of Figure 5. In that case, it is important to note that no re-stabilization was obtained, unlike that observed by Ioualalen and Kharif [5] who used a hybrid asymptotic solution. It is found that the residual branch (2) is always stable. In fact, it is shown that the asymptotic SCW solution of Ioualalen and Kharif [5] is a composite of the stable segment of branch (3) and the stable branch (2) before and beyond the pole. It is therefore obvious that if their solution is unstable near the instability region (as they found), it must necessarily exhibit a bubble-like structure.

A further interesting result follows: the solution involves a subcritical bifurcation from an unstable to a stable mode around critical incidence angle and wave steepness h. As a result, there always exists a stable state regardless of the wave steepness, thus avoiding an HR appearance. These two competing processes certainly need to be further investigated.