# Vorticity effects on steady nonlinear periodic gravity-capillary waves in finite depth

S. Halder, M Francius, A K Dhar, S Mukherjee, H C Hsu, C Kharif

► To cite this version:

S. Halder, M Francius, A K Dhar, S Mukherjee, H C Hsu, et al. Vorticity effects on steady nonlinear periodic gravity-capillary waves in finite depth. Journal of Fluid Mechanics, 2025. hal-05308477

HAL Id: hal-05308477

https://hal.science/hal-05308477v1

Submitted on 10 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
Banner appropriate to article type will appear here in typeset article

1

2

# Vorticity effects on steady nonlinear periodic gravity-capillary waves in finite depth

3

S. Halder¹, M. Francius²†, A. K. Dhar¹, S. Mukherjee¹, H.C. Hsu³ and C. Kharif⁴

4

¹Department of Mathematics, Indian Institute of Engineering Science and Technology, Shibpur, Howrah,

5

711103, West Bengal, India

6

²Université de Toulon, Aix Marseille Université, CNRS, IRD, MIO, Toulon, France

7

³Department of Marine Environment and Engineering, National Sun Yat-Sen University, Kaohsiung 804,

8

Taiwan

9

⁴Aix-Marseille Université, CNRS, Centrale Marseille, IRPHE, UMR 7342, 13384, Marseille, France.

10

(Received xx; revised xx; accepted xx)

11 Periodic gravity-capillary waves on a fluid of finite depth with constant vorticity are studied theoretically and numerically. The classical Stokes expansion method is applied to obtain the wave profile and the interior flow up to the fourth-order of approximation, which thereby extends the works of Barakat & Houston (1968) and Hsu et al. (2016). The classical perturbation scheme possesses singularities for certain wavenumbers, whose variations with depth are shown to be affected by the vorticity. This analysis also reveals that for any given value of the physical depth, there exist a threshold value of the vorticity above which there are no singularities in the theoretical solution. The validity of the third and fourth-order solutions is examined by comparison with exact numerical results, which are obtained with a method based on conformal mapping and Fourier series expansions of the wave surface. The outcomes of this comparison are surprising as they report important differences in the internal flow structure, when compared with the third-order predictions, even though both approximations predict almost perfectly the phase velocity and the surface profiles. Usually, this occur when the wavenumber is far enough a critical value and the steepness is not too large. In these non resonant cases, it is found that the fourth-order theory is more consistent with the exact numerical results. With negative vorticity the improvement is noticeable both beneath the crest and the trough, whereas with positive vorticity the fourth-order theory does a better job either beneath the crest or beneath the trough, depending of the type of the wave.

29 **Key words:** Fourth-order perturbation solution, constant vorticity, gravity-capillary waves,
30 Wilton singularity.

31

## 1. Introduction

32 Gravity-capillary (GC) waves can propagate at the surface of a fluid when both gravity and surface tension are important as restoring forces. Usually, they first appear during the initial

33

† Email address for correspondence: marc.francius@mio.osupytheas.fr
34

35

generation of waves by the wind, but they may also result from almost breaking gravity waves, rain or other millimetric-to-centimetric sized disturbances.

36

37

38

39

40

41

42

43

44

In view of their importance to wave energy transfer across the spectrum due to nonlinear wave-wave interactions, to air-sea momentum transfer as they contribute to some extent to the ocean surface stress, and to remote sensing of the ocean surface as their wavelengths (from 2 mm to 2 dm) contribute to a large extent to the sea-surface backscattering of electromagnetic microwaves in the X-band (centimeter range) and the Ku-band (subcentimeter range), they have been the subject of many studies. In fact, understanding their properties has been of utmost importance during the last decades with the advances of imaging airborne and spaceborn radars. The present investigation, however, is restricted to the problem of periodic GC waves steadily traveling on water of finite depth with a linear shear current.

45

46

47

48

49

50

51

52

53

54

55

56

57

Generally, methods used to solve problems related to periodic GC waves of arbitrary amplitude can be divided in two main types. One is based on perturbation expansion methods and the second category is purely numerical. The first to investigate GC waves was Harrison (1909) who obtained a third-order Stokes-type solution to the problem of periodic GC waves steadily traveling on water of infinite depth in irrotational flows. This work was extended to fifth-order by Wilton (1915), who also found that the classical Stokes expansion method breaks down for a countable set of critical wavenumbers. For the largest of these critical wavenumbers Wilton showed that two different solutions could exist, the so-called Wilton ripples. One of the profiles is gravity-like having a phase speed that increases with the amplitude, while the other is capillary-like having a phase speed that decreases as the amplitude increases. The phenomenon of Wilton's ripples was further analyzed by Pierson & Fife (1961) and Nayfeh (1970b), who found valid solutions when the wavenumber is near the first critical value and second critical value, respectively.

58

59

60

61

62

63

64

65

66

67

These works based on perturbation expansions culminated in the higher-order extensions by computer with the seminal reference works of Hogan (1980, 1981), using the methods pioneered by Schwartz (1974) for gravity waves. On extending Pierson & Fife (1961), the phenomenon of Wilton's ripples was resolved by Hogan (1981), who showed how to modify the perturbation scheme to avoid the non-uniformity in the ordering of the Fourier coefficients of the surface profile, which occurs when the wavenumber is both at and near singular values. In finite depth the classical perturbation method was carried out by Nayfeh (1970a) and Barakat & Houston (1968) to third- and fourth-order, respectively. Wilton ripples were also found to exist in finite depth, and Wilton-like solutions were derived at second-order by Barakat & Houston (1968) and at third-order by Nayfeh (1970a).

68

69

70

71

72

73

74

75

76

77

78

Though much attention has been given to GC waves at the free surface of irrotational flows, rather few studies have been made in the case when the flow is rotational. While irrotational flows are suitable for waves traveling into still water or over a uniform current, nonuniform currents usually give rise to water flows with vorticity. Quite apart from the difficult questions of existence, uniqueness and analyticity of small amplitude GC waves in the case of a rotational flow with finite depth (readers interested by these issues may read Martin (2013), Martin & Matioc (2013) and references therein), a brief, and certainly incomplete, review of the literature shows that the vast majority of existing works concerning GC waves on rotational currents assumes that the waves are two-dimensional and the background current is linearly varying with depth, as these assumptions simplify the analysis, and thus the flows have constant vorticity.

79

80

81

82

83

In contrast with the case when surface tension is neglected, studies devoted to the computations of traveling periodic GC waves in the presence of a linear sheared current are very limited. One can refer to the paper by Brantenberg & Brevik (1993) who used a two-layer model for the background current, the constant vorticity being confined in an upper layer on top of an infinite irrotational layer, and a third-order Stokes expansion for the GC
84 waves propagating on the surface of the upper layer. Subsequently, Hsu et al. (2016) derived
85 a third-order Stokes-like solution for periodic steady waves which includes the effects of
86 capillarity and constant vorticity, assuming a single layer model.

87 Using numerical methods, inclusion of surface tension effects in rotational flows with
88 constant vorticity has also been considered by Kang & Vanden-Broeck (2000), but the
89 focus was more on pure solitary waves and the so-called generalized solitary waves that are
90 characterized by oscillatory tails in the far field. Most recently, in the same vein, Guo et al.
91 (2014) have carried out numerical studies on the unsteady dynamics of GC solitary waves,
92 using the same time-dependent conformal mapping technique introduced by Choi (2009) for
93 pure gravity waves in a background linear sheared current.

94 Recently, there has been a renewed interest in several aspects of the dynamics of GC
95 waves in the case of contant vorticity and finite depth. Following Chabane & Choi (2019)
96 who analyzed, for irrotational flows, the interactions of three distinct GC modes that are
97 not necessarily collinear, Ivanov & Martin (2019) derived the amplitude equations to study
98 three-wave resonant interactions of one-directional propagating GC waves in the presence of
99 a linear sheared current. Moreover, from multiple scale expansion of the primitive equations
100 in finite depth and classical expansions series of the unknown variables, Hsu et al. (2018)
101 derived a third-order nonlinear Schrödinger equation to study the stability of weakly nonlinear
102 uniform GC wavetrain to one-dimensional disturbances. By analyzing the effect of vorticity
103 on the stability diagrams of the problem, Hsu et al. (2018) concluded that the properties
104 of modulational instabilities (side-band disturbances) of a uniform GC wavetrain depend
105 crucially on the existence of a linear shear current. At the same time, Curtis et al. (2018)
106 derived a higher-order nonlinear Schrödinger equation in the presence of surface tension,
107 to further investigate how a linearly sheared current affects the modulational instability
108 properties of weakly nonlinear GC waves in deep water, with a particular attention to the
109 motion and mean properties of particle paths. Later on, Dhar & Kirby (2023) extended the
110 results of Hsu et al. (2018) to fourth-order in the finite depth case, and they showed that their
111 results exhibit considerable deviations from those obtained with the third-order analysis.
112 Very recently, in order to get insights in the modulational instabilities Gao et al. (2021) have
113 developed weakly nonlinear theories for two particular resonant cases, second-harmonic
114 resonance and long-wave/short-wave interaction between capillary-gravity waves in finite
115 depth with a linear shear current. We should also mention that the modulational instability
116 has also been studied from heuristic nonlocal evolution equations with both surface tension
117 and constant vorticity, more precisely with an extension of the heuristic Whitham equation
118 by Hur & Johnson (2015) and with a corresponding set of shallow water equations by Hur
119 (2019) to permit bidirectional propagation.

120 Considering the results mentioned above and given that, to the best of the authors'
121 knowledge, all the existing Stokes-type solutions for GC waves with constant vorticity in
122 finite depth are only accurate up to the third-order, we shall extend in this paper the analysis
123 of Hsu et al. (2016) to include the fourth-order effects. The present paper also extends the
124 work of Barakat & Houston (1968) who obtained, without vorticity, the 4th-order Stokes-like
125 solutions for periodic gravity-capillary waves in finite depth. One source of motivation in
126 the present study has been the discovery by Fang et al. (2023) that experimental evidences
127 indicate that higher-order Stokes solutions are more accurate in describing the velocity
128 distributions, especially in strong following currents and positive vorticity conditions. These
129 authors have derived a new set of fifth-order Stokes solutions for periodic gravity-waves
130 interacting with a linear shear current.

131 To assess the accuracy of the third- and fourth-order theories, we shall obtain the reference
132 solution of the primitive equations using a numerical method based on a conformal mapping
133 method and Fourier series expansions of the unknowns. The differences between the exact
results and the approximations are examined in terms of surface profile, phase velocity and interior horizontal wave-induced velocities. It should be emphasized that, with constant vorticity and surface tension effects taken into account, previous studies based on perturbation methods have focused mainly on the wave profiles, but rarely on the wave-induced velocities beneath the free surface. As far as we know, only Curtis et al. (2018) have used a second-order approximation of the interior wave-induced velocities in infinite depth, in order to analyze the mean transport properties of weakly modulated plane gravity-capillary waves with constant vorticity. Our work complements thus the work of Curtis et al. (2018), and reveals that significant differences between the third-order and fourth-order predictions exist in the wave-induced velocities below the free surface, even though both approximations predict almost perfectly the phase velocity and the surface profiles.

The paper is organized as follows: The basic equations and the framework of the problem is presented in section 2. In section 3, we derive the fourth-order analytical solution for these periodic GC waves and discuss the effects of depth and vorticity on the occurrence of singularities in the classical perturbation theory. In contrast with the analysis presented by Hsu et al. (2016) the occurrence of singularities at any order is unambiguously elucidated as a function of depth and vorticity. In section 4, we first introduce a numerical method to compute exact steady GC wave solutions with constant vorticity in finite depth and, following, a comparison between the numerical results and the analytical results is presented. Final remarks are summarized in section 5.

## 2. Basic equations for steady waves

As is well known for two-dimensional incompressible flows of inviscid fluids, any perturbations of a background flow with constant vorticity, either small or of finite amplitude, are necessarily irrotational motions as a consequence of Kelvin's circulation theorem. Hence, letting $\Omega_0$ be the magnitude of the shear, it is possible to assume the existence of a velocity potential function $\Phi^*(x, y, t)$ such that the total velocity can be written as

$$ V = [\Omega_0(y + h_0)]\mathbf{i} + \nabla\Phi^*(x, y, t) \quad (2.1) $$

With this decomposition, the vorticity of the flow equals $-\Omega_0$. In this study, we consider the evolution of a steadily traveling periodic GC wave on the free surface of water of finite depth $h_0$. The undisturbed free surface is taken to the plane $y = 0$, where $y$ points vertically upwards, and the bed of the water is defined by $y = -h_0$.

Dimensionless variables are introduced by the following transformations

$$ \sqrt{\frac{k_0^3}{g}}\Phi^* \rightarrow \Phi, \quad k_0\mathbf{x} \rightarrow \mathbf{x}, \quad \sqrt{gk_0}t \rightarrow t, \quad \frac{k_0}{\rho g}P_0 \rightarrow P, \quad \sqrt{\frac{k_0}{g}}C_* \rightarrow c, $$

where $\mathbf{x} = (x, y)$ and $\rho, g, P_0, k_0, C_*$, represent respectively the fluid density, gravitational acceleration, dimensional pressure, wavenumber and wave phase velocity of a steadily traveling wave. With this scaling, the relevant dimensionless parameters are

$$ \Omega = \frac{\Omega_0}{\sqrt{gk_0}}, \quad \kappa = \frac{T k_0^2}{\rho g}, \quad \mu = k_0 h_0, $$

where $T$ represents the surface tension coefficient, and equation (2.1) can be expressed as

$$ V = [\Omega(y + \mu)]\mathbf{i} + \nabla\Phi(x, y, t). \quad (2.2) $$

Since the wave-induced motions are irrotational, the velocity potential $\Phi(x, y, t)$ satisfies the
Laplace equation
170

<figure><img src="image_1.png" /></figure>

$$
\Phi_{xx} + \Phi_{yy} = 0 \quad \text{for} \quad -\mu \le y \le \zeta(x,t) \qquad (2.3)
$$

where $y = \zeta(x, t)$ represents the dimensionless surface elevation, and the boundary conditions
173

174 are

<figure><img src="image_2.png" /></figure>

$$
\Phi_y = 0 \quad \text{at} \quad y = -h, \tag{2.4}
$$

<figure><img src="image_3.png" /></figure>

$$
\Phi_y = \zeta_t + [\Phi_x + \Omega(\zeta + \mu)]\zeta_x \quad \text{at} \quad y = \zeta(x, t), \quad (2.5)
$$

<figure><img src="image_4.png" /></figure>

$$
\Phi_t + \frac{1}{2} [\Phi_x + \Omega(\zeta + \mu)]^2 + \frac{1}{2} \Phi_y^2 + \zeta - \Omega\psi - \frac{\Omega^2}{2} (\zeta + \mu)^2
$$

<figure><img src="image_5.png" /></figure>

$$
- \kappa \zeta_{xx} (1 + \zeta_x^2)^{-3/2} = f(t) \quad \text{at} \quad y = \zeta(x, t). \qquad (2.6)
$$

Here $\psi$ is the harmonic conjugate function of $\Phi$, and $f(t)$ is an arbitrary function of time $t$,
180

which can be set to zero without loss of generality or absorbed into the definition of $\Phi$.
181

To seek steadily traveling waves of the above equations, it is also convenient to choose

a frame of reference moving with the dimensionless phase velocity $c$. Hence we make the
183

184 following change of variables

<figure><img src="image_6.png" /></figure>

$$
\xi = x - ct, \quad \eta = y, \quad \tau = t, \tag{2.7}
$$

so that the governing equations become
187

<figure><img src="image_7.png" /></figure>

$$
\Phi_{\xi\xi} + \Phi_{\eta\eta} = 0, \quad -\mu \le \eta \le E(\xi) \qquad (2.8)
$$

<figure><img src="image_8.png" /></figure>

$$
\Phi_{\eta} = (\Phi_{\xi} + \Omega E - C)E_{\xi}, \quad \eta = E(\xi) \tag{2.9}
$$

<figure><img src="image_9.png" /></figure>

$$
(\Phi_{\xi} + \Omega E - C)^2 - C^2 + \Phi_{\eta}^2 + 2E \\
- 2\kappa E_{\xi}\xi (1 + E_{\xi}^2)^{-3/2} = K, \quad \eta = E(\xi) \tag{2.10}
$$

<figure><img src="image_10.png" /></figure>

$$
\Phi_{\eta} = 0, \quad \eta = -\mu_{*}, \tag{2.11}
$$

<figure><img src="image_11.png" /></figure>

where $\Phi(\xi, \eta) = \Phi(x, y, t)$, $E(\xi) = \zeta(x, t)$, the Bernoulli's constant $K = 2\Omega\psi + \Omega^2 E^2 - 2\Omega\{E(c - \Omega\mu)\}$ and $C = c - \Omega\mu$.

195 where $\Phi(\xi, \eta) = \Phi(x, y, t)$, $E(\xi) = \zeta(x, t)$, the Bernoulli's constant $K = 2\Omega\psi + \Omega^2 E^2 - 2\Omega\{E(c - \Omega\mu)\}$ and $C = c - \Omega\mu$.

3. Asymptotic solutions
196

In order to solve (2.8-2.11), we use the classical Stokes expansion method and take
197

$\epsilon = k_0 a_1$ as a small parameter, where $a_1$ represents the amplitude of the first harmonic.
198

Like Hsu et al. (2016), we consider the following expansions for the unknown quantities
199

$\Phi(\xi, \eta)$, $\bar{E}(\xi)$, $C$ and $K$

$$
\Phi(\xi, \eta) = \sum_{m=1}^{\infty} \Phi_m(\xi, \eta), \quad E(\xi) = \sum_{m=1}^{\infty} E_m(\xi), \quad C = \sum_{m=0}^{\infty} \epsilon^m C_m, \quad K = \sum_{m=0}^{\infty} K_m, \quad (3.1)
$$

where the quantities $\Phi_m$, $E_m$ and $K_m$ are assumed to be of order $\epsilon^m$. By definition, the

quantities $C_m$ do not depend on $\epsilon$ and so the leading order term is $m = 0$, whereas each term

in the series of the velocity potential $\Phi(\xi, \eta)$ and the surface elevation amplitude $E(\xi)$ depend

on the amplitude parameter $\epsilon$, the subscript indicating its order. Substituting expansions (3.1)

in the system of equations (2.8-2.11) and then equating coefficients of different powers of $\epsilon$

on both sides, we obtain the approximations for consecutive orders of the nonlinear problem.

To the first-order in the expansion parameter $\epsilon$, we obtain the velocity potential and the
210 surface elevation

211
$$
\Phi_1(\xi, \eta) = \frac{C_0 \epsilon \cosh[(\eta + \mu)]}{\sigma \cosh(\mu)} \sin(\xi)
$$
(3.2)

212
$$
E_1(\xi) = \epsilon \cos(\xi), \tag{3.3}
$$

213
$$
K_1 = 0, \tag{3.4}
$$

215 as well as the linear dispersion relation connecting frequency and wavenumber of the primary
216 wave

217
$$
C_0^2 + \sigma C_0 \Omega = \sigma(1 + \kappa), \tag{3.5}
$$

219 where $\sigma = \tanh(k_0 h_0)$. Equation (3.5) is the same as the equation (3.7) of Hsu et al. (2016)
220 and the equation (38a) of Kang & Vanden-Broeck (2000) by returning to dimensional form.

221 Setting $\bar{\Omega} = \sigma\Omega/C_0$, the linear dispersion relation equation (3.5) can be rewritten like

222
$$
C_0^2 (1 + \bar{\Omega}) = \sigma (1 + \kappa), \tag{3.6}
$$

224 which shows that in any case $\bar{\Omega} > -1$. As we shall see, using $\bar{\Omega}$ and $\bar{\kappa} = \kappa/(\sigma C_0^2)$, introduced
225 for convenience, we have obtained compact formulas for the higher-order coefficients.

226 Calculation of the higher orders is quite standard and therefore omitted here to improve
227 readability. The analytical work was done with Mathematica. In the following we present
228 only the novel fourth-order equations and corresponding solutions. Since we have found few
229 misprints in some formulas of Hsu et al. (2016), both at second-order and third-order, the
230 intermediate results are given in Appendices A and B.

231 Next, we discuss the occurence of resonant cases taking into account vorticity and finite
232 depth effects. Notice that the misprints found in the formulas of Hsu et al. (2016) have no
233 consequences on the expressions of the first two critical values $\kappa_1$ and $\kappa_2$ (see appendices
234 A and B). Actually their expressions are identical to ours, when their (B, X) are replaced
235 by our ($\bar{\kappa}$, $\bar{\Omega}$). When $\kappa = \kappa_1$, the second-order solutions and the third-order solutions are
236 found to be singular, whereas when $\kappa = \kappa_2$ singularities occur only at third-order. Hsu et
237 al. (2016) have discussed the dependency of these critical values on depth and vorticity, but
238 their interpretation is ambiguous. We shall therefore clarify this discussion herein this study.

239
3.1. *Approximate solution of the fourth-order*

240 To the fourth-order in $\epsilon$, we obtain the sequence of equations given by

241
$$
\Phi_{4\xi\xi} + \Phi_{4\eta\eta} = 0, \quad -\mu \le \eta < 0 \tag{3.7}
$$

242
$$
\Phi_{4\eta} + C_0 E_{4\xi} = -E_1 \Phi_{3\eta\eta} - E_2 \Phi_{2\eta\eta} - E_3 \Phi_{1\eta\eta} - \frac{1}{2} E_1^2 \Phi_{2\eta\eta\eta} - E_1 E_2 \Phi_{1\eta\eta\eta} \\
\phantom{\Phi_{4\eta} + C_0 E_{4\xi} = -E_1 \Phi_{3\eta\eta} -} -\frac{1}{6} E_1^3 \Phi_{1\eta\eta\eta\eta} + (\Phi_{1\xi} + \Omega E_1) E_3 \xi + (\Phi_{2\xi} + E_1 \Phi_{1\xi\eta} + \Omega E_2 - C_2) E_2 \xi \\
\phantom{\Phi_{4\eta} + C_0 E_{4\xi} = -E_1 \Phi_{3\eta\eta} -} + (\Phi_{3\xi} + E_1 \Phi_{2\xi\eta} + E_2 \Phi_{1\xi\eta} + \frac{1}{2} E_1^2 \Phi_{1\xi\eta\eta} + \Omega E_3 - C_3) E_1 \xi, \quad \eta = 0 \tag{3.8}
$$
246

$$
\begin{align}
C_0 \Phi_{4\xi} - (1 - C_0 \Omega) E_4 + \kappa E_{4\xi\xi} &= \frac{1}{2} (\Phi_{2\xi}^2 + \Phi_{2\eta}^2) + \Phi_{1\xi} \Phi_{3\xi} + \Phi_{1\eta} \Phi_{3\eta} \nonumber \\
&\quad + E_1 (\Phi_{1\xi} \Phi_{2\xi\eta} + \Phi_{2\xi} \Phi_{1\xi\eta} + \Phi_{1\eta} \Phi_{2\eta\eta} + \Phi_{2\eta} \Phi_{1\eta\eta}) + \frac{1}{2} E_1^2 (\Phi_{1\xi}^2 \eta + \Phi_{1\xi} \Phi_{1\xi\eta\eta} \nonumber \\
&\quad + \Phi_{1\eta\eta}^2 + \Phi_{1\eta} \Phi_{1\eta\eta\eta}) + E_2 (\Phi_{1\xi} \Phi_{1\xi\eta} + \Phi_{1\eta} \Phi_{1\eta\eta}) + \frac{1}{2} \Omega^2 (E_2^2 + 2E_1 E_3) \nonumber \\
&\quad + \Omega (E_1 \Phi_{3\xi} + E_2 \Phi_{2\xi} + E_3 \Phi_{1\xi} + E_1^2 \Phi_{2\xi\eta} + 2E_1 E_2 \Phi_{1\xi\eta} + \frac{1}{2} E_1^3 \Phi_{1\xi\eta\eta}) \nonumber \\
&\quad - C_0 (E_1 \Phi_{3\xi\eta} + E_2 \Phi_{2\xi\eta} + E_3 \Phi_{1\xi\eta} + \frac{1}{2} E_1^2 \Phi_{2\xi\eta\eta} + E_1 E_2 \Phi_{1\xi\eta\eta} + \frac{1}{6} E_1^3 \Phi_{1\xi\eta\eta\eta}) \nonumber \\
&\quad - C_3 (\Phi_{1\xi} + \Omega E_1) - C_2 (\Phi_{2\xi} + \Omega E_2 + E_1 \Phi_{1\xi\eta}) + \frac{3}{2} \kappa E_1^2 \xi E_2 \xi \nonumber \\
&\quad + 3 \kappa E_1 \xi E_2 \xi E_1 \xi \xi - \frac{1}{2} K_4 = 0, \quad \eta = 0 \tag{3.9}
\end{align}
$$

254

$$
\Phi_{4\eta} = 0, \quad \eta = -\mu \tag{3.10}
$$

256
257 Classically, we seek the solutions for the fourth-order potential and elevation amplitude in the form

$$
258
\Phi_4(\xi, \eta) = A_{44} \frac{\cosh[4(\eta + \mu)]}{\cosh(4\mu)} \sin(4\xi) + A_{42} \frac{\cosh[2(\eta + \mu)]}{\cosh(2\mu)} \sin(2\xi) \quad (3.11)
$$

$$
\frac{268}{269} \qquad E_4(\xi) = a_{44} \cos(4\xi) + a_{42} \cos(2\xi) \tag{3.12}
$$

261
262 Then by inserting (3.11-3.12) into the boundary conditions (3.8-3.9) and using the first,
second and third-order solutions, after tedious algebra we obtain

$$
\begin{align*}
A_{44} &= \frac{3\epsilon}{8[(5+\sigma^2)-15\bar{\kappa}(1+\sigma^2)]} \frac{(1+6\sigma^2+\sigma^4)}{\sigma^3(1+3\sigma^2)} \left[ (1+3\sigma^2)\bar{\Omega} + 2(1-3\sigma^2-2\sigma^4) \right. \\
&\quad \left. + 15\bar{\kappa}\sigma^2(1+3\sigma^2) \right] A_{33} + \frac{C_0\epsilon}{8[(5+\sigma^2)-15\bar{\kappa}(1+\sigma^2)]} \frac{(1+6\sigma^2+\sigma^4)}{\sigma^4} \left[ \bar{\Omega}^2 + \left\{2+15\bar{\kappa}\sigma^2\right\}\bar{\Omega} \right. \\
&\quad \left. + (1-\sigma^2)+15\bar{\kappa}\sigma^2 \right] a_{33} + \frac{C_0\epsilon^4}{768(1-3\bar{\kappa})^2[(5+\sigma^2)-15\bar{\kappa}(1+\sigma^2)]} \frac{(1+6\sigma^2+\sigma^4)}{\sigma^{10}} \times \\
&\quad \left[ 3\bar{\Omega}^6 + \left\{3(9+4\sigma^2)+45\bar{\kappa}\sigma^2\right\}\bar{\Omega}^5 + \left\{12(9+3\sigma^2+\sigma^4)+18\bar{\kappa}\sigma^2(21+11\sigma^2)\right\}\bar{\Omega}^4 \right. \\
&\quad \left. + \left\{3(81-14\sigma^2+21\sigma^4)+9\bar{\kappa}\sigma^2(149+68\sigma^2+19\sigma^4)+270\bar{\kappa}^2\sigma^4(1+\sigma^2)\right\}\bar{\Omega}^3 \right. \\
&\quad \left. + \left\{12(27-29\sigma^2+21\sigma^4-4\sigma^6)+18\bar{\kappa}\sigma^2(139+18\sigma^2+41\sigma^4) + 27\bar{\kappa}^2\sigma^4(51-34\sigma^2+11\sigma^4) \right\}\bar{\Omega}^2 \right. \\
&\quad \left. + \left\{3(81-192\sigma^2+169\sigma^4-82\sigma^6) + 9\bar{\kappa}\sigma^2(273-110\sigma^2+253\sigma^4-12\sigma^6) + 54\bar{\kappa}^2\sigma^4(47-84\sigma^2+15\sigma^4) - 3240\bar{\kappa}^3\sigma^8 \right\}\bar{\Omega} \right. \\
&\quad \left. + (81-342\sigma^2+432\sigma^4-302\sigma^6+131\sigma^8) + 6\bar{\kappa}\sigma^2(171-177\sigma^2+359\sigma^4-283\sigma^6) + 144\bar{\kappa}^2\sigma^4(12-38\sigma^2+43\sigma^4) - 5940\bar{\kappa}^3\sigma^8 \right] \tag{3.13}
\end{align*}
$$
274

275

276

277

278

279

280

$$
\begin{align}
A_{42} ={}& \frac{3\epsilon}{4(1-3\bar{\kappa})} \frac{(1+\sigma^2)}{\sigma^3(1+3\sigma^2)} \left[ (1+3\sigma^2)\bar{\Omega} + 2(1-\sigma^4) + 3\bar{\kappa}\sigma^2(1+3\sigma^2) \right] A_{33} \nonumber \\
& + \frac{\epsilon}{4(1-3\bar{\kappa})} \frac{(1+\sigma^2)}{\sigma^3} \left[ \bar{\Omega} + 2(1-\sigma^2) + 3\sigma^2\bar{\kappa} \right] A_{31} - \frac{\epsilon^4}{8(1-3\bar{\kappa})^2} \frac{(1+\sigma^2)}{\sigma^6} \times \nonumber \\
& \left[ \bar{\Omega}^3 + \left\{ (5+2\sigma^2) + 3\bar{\kappa}\sigma^2 \right\} \bar{\Omega}^2 + \left\{ 3(3+\sigma^2) + 6\bar{\kappa}\sigma^2(2+\sigma^2) \right\} \bar{\Omega} \right. \nonumber \\
& \left. + (6-\sigma^2-\sigma^4) + 3\bar{\kappa}\sigma^2(5-\sigma^2) \right] C_2 + \frac{C_0\epsilon^4}{48(1-3\bar{\kappa})^2} \frac{(1+\sigma^2)}{\sigma^6} \left[ 18\bar{\Omega}^3 \right. \nonumber \\
& \left. + \left\{ 3(25+\sigma^2) + 45\bar{\kappa}\sigma^2 \right\} \bar{\Omega}^2 + \left\{ 39(3-\sigma^2) + 18\bar{\kappa}\sigma^2(7-\sigma^2) + 54\bar{\kappa}^2\sigma^4 \right\} \bar{\Omega} \right. \nonumber \\
& \left. + (63-52\sigma^2+13\sigma^4) + 3\bar{\kappa}\sigma^2(31-23\sigma^2) + 90\bar{\kappa}^2\sigma^4 \right], \tag{3.14}
\end{align}
$$

281

282

283

284

285

$$
\begin{equation}
\begin{aligned}
a_{44} ={}& \frac{\tanh(4\mu)}{C_0} A_{44} + \frac{3\epsilon}{2C_0} A_{33} + \frac{\epsilon}{2\sigma}(\bar{\Omega}+1)a_{33} + \frac{\epsilon^4}{192(1-3\bar{\kappa})^2} \frac{1}{\sigma^7} \left[ 3\bar{\Omega}^5 \right. \\
& \left. \qquad + 12(2+\sigma^2)\bar{\Omega}^4 + \left\{ 9(9+6\sigma^2+\sigma^4) + 18\bar{\kappa}\sigma^2(1+\sigma^2) \right\}\bar{\Omega}^3 + \left\{ 12(12+7\sigma^2+2\sigma^4) \right. \right. \\
& \left. \qquad \left. + 18\bar{\kappa}\sigma^2(5+\sigma^4) \right\}\bar{\Omega}^2 + \left\{ 3(45+12\sigma^2+13\sigma^4-2\sigma^6) + 18\bar{\kappa}\sigma^2(9-8\sigma^2+3\sigma^4) \right. \right. \\
& \left. \qquad \left. - 216\bar{\kappa}^2\sigma^6 \right\}\bar{\Omega} + 2(27-9\sigma^2+27\sigma^4-31\sigma^6) + 12\bar{\kappa}\sigma^2(9-21\sigma^2+28\sigma^4) - 396\bar{\kappa}^2\sigma^6 \right],
\end{aligned}
\tag{3.15}
\end{equation}
$$

286

and furthermore

287

288

289

$$
a_{42} = \frac{\tanh(2\mu)}{C_0} A_{42} + \frac{3\epsilon}{2C_0} A_{33} + \frac{\epsilon}{2C_0} A_{31} - \frac{\epsilon^2}{C_0} (C_2 a_2) + \frac{\epsilon^4}{24(1-3\bar{\kappa})} \frac{1}{\sigma^3} \times \\
\left[ 9\bar{\Omega}^2 + \left\{ 3(9+\sigma^2) + 18\bar{\kappa}\sigma^2 \right\} \bar{\Omega} + (27-19\sigma^2) + 30\bar{\kappa}\sigma^2 \right], \quad (3.16)
$$

291

292

293

294

295

296

297

$$
\begin{align}
K_4 ={}& \frac{C_0 \epsilon}{\sigma} (\bar{\Omega} + 1) A_{31} + \frac{C_0^2 \epsilon}{\sigma^2} (\bar{\Omega}^2 + \bar{\Omega} - \sigma^2) a_{31} - \epsilon^4 C_0 C_2 + \frac{C_0^2 \epsilon^4}{32(1 - 3\bar{\kappa})^2} \frac{1}{\sigma^8} \times \nonumber \\
& \left[ \bar{\Omega}^6 + 4(2 + \sigma^2)\bar{\Omega}^5 + \left\{ 4(7 + 4\sigma^2 + \sigma^4) + 6\bar{\kappa}\sigma^2(1 + \sigma^2) \right\}\bar{\Omega}^4 + \left\{ 2(27 + 8\sigma^2 + 9\sigma^4) \right. \right. \nonumber \\
& \left. \left. + 12\bar{\kappa}\sigma^2(3 - 2\sigma^2 + \sigma^4) \right\}\bar{\Omega}^3 + \left\{ 12(5 - 2\sigma^2 + 5\sigma^4) + 12\bar{\kappa}\sigma^2(7 - 14\sigma^2 + 2\sigma^4) \right. \right. \nonumber \\
& \left. \left. + 9\bar{\kappa}^2\sigma^4(1 - 10\sigma^2 + \sigma^4) \right\}\bar{\Omega}^2 + \left\{ 12(3 - 5\sigma^2 + 9\sigma^4 - 2\sigma^6) \right. \right. \nonumber \\
& \left. \left. + 6\bar{\kappa}\sigma^2(15 - 53\sigma^2 + 27\sigma^4 - \sigma^6) + 36\bar{\kappa}^2\sigma^4(1 - 6\sigma^2) \right\}\bar{\Omega} \right. \nonumber \\
& \left. + (9 - 36\sigma^2 + 78\sigma^4 - 40\sigma^6 + 9\sigma^8) + 12\bar{\kappa}\sigma^2(3 - 15\sigma^2 + 12\sigma^4 - 4\sigma^6) \right. \nonumber \\
& \left. + 36\bar{\kappa}^2\sigma^4(1 - 9\sigma^2 + 2\sigma^4) \right] \tag{3.17}
\end{align}
$$

398

$$
C_3 = 0. \tag{3.18}
$$
Equations (3.13-3.18) constitute the main analytical results of this study and extend the third-order analysis initiated by Hsu et al. (2016). In the next section we shall use these relations to describe the improvements gained with the fourth-order approximation by analyzing not only the surface profiles but also the internal flow velocities.

Not surprisingly, the fourth-order approximation is found to be invalid when $κ$ equals another critical value, $κ_3$, the root of the following equation

$$κ = \frac{(5 + σ^2)σ^2}{15(1 + \overline{Ω}) + 5(2 + 3\overline{Ω})σ^2 - σ^4} \quad (3.19)$$

Without vorticity, equation (3.19) becomes

$$κ_3 = \frac{(5 + σ^2)σ^2}{15 + 10σ^2 - σ^4} \quad (3.20)$$

in agreement with the results of Barakat & Houston (1968). This latter equation reduces in deep water ($σ = 1$) to the other well-known value $κ_3 = 1/4$.

## 3.2. Critical wavenumbers

Since the pioneering works of Harrison (1909) and Pierson & Fife (1961), it is well known that the application of the Stokes expansion method for the periodic steady irrotational GC wave problem, either in deep water or finite depth, is not legitimate when $κ$ is close to the critical values $κ_n$ ($n \ge 1$), which form, when they exist, an infinite discrete set of critical values. At these critical values, it is observed that there is a lack of uniqueness in the linear solution. Physically, it can be explained by a resonance mechanism due to the fact that infinitesimal waves of wavenumber $k_0$ and $nk_0$ have the same phase velocity when $κ = κ_n$.

For irrotational GC waves on finite depth, Barakat & Houston (1968) remarked that the nth approximation obtained with the classical Stokes expansion method breaks down when $κ = κ_{n-1}$, the roots of the following equations

$$κ = \frac{n \tanh μ - \tanh(nμ)}{n^2 \tanh(nμ) - n \tanh μ}, \quad n \ge 2 \quad (3.21)$$

In deep water these roots have the well-known values $κ_{n-1} = 1/n$ ($n \ge 2$). Solutions that bifurcate from these critical values are the so-called Wilton ripples, which correspond to resonances between the fundamental and the nth harmonic. For $n = 4$, with some algebra one can show that equation (3.21) corresponds to equation (3.20).

If the dependence between the dimensionless parameters $κ$ and $μ$ is not considered, graphs of $κ_n$ versus $μ$ may be misleading with respect to the interpretation of the influence of the depth of the fluid layer. It should be realized that the critical numbers are the roots $κ_{n-1}$ of the following equations for $n \ge 2$,

$$κ \left[ n^2 \tanh(nκ^{1/2}h) - n \tanh(nκ^{1/2}h) \right] - \left[ n \tanh(nκ^{1/2}h) - \tanh(nκ^{1/2}h) \right] = 0, \quad (3.22)$$

obtained from equation (3.21) and the relation $μ = κ^{1/2}h$, where $h = k_c h_0$ and $k_c = \sqrt{ρg/T}$. Numerical calculations were performed for water ($T = 74$ dynes/cm), by using a variant of the Newton-Raphson method to solve equation (3.22) for $κ$. The results for the first three critical values are shown in figure 1 as a function of $h$. The dotted lines represent the well-known asymptotic values of $κ_n$ ($n = 1, 2$ and $3$) for irrotational GC waves in deep water. The thin solid lines are drawn to show the dependency of the dimensionless number $κ$ on the depth layer $h$ for fixed values of $μ$. The thin solid lines of figure 1 are drawn for $μ = 10, 2$ and $1$.

As shown in figure 1 the critical values $κ_n$ decrease rapidly with decreasing $h$, and
<figure><img src="image_12.png" /><figcaption>Figure 1: Variations of $κ_1, κ_2, κ_3$ as a function of depth for water for irrotational GC waves. In any case $κ_1 > κ_2 > κ_3$. The dotted lines represent the corresponding deep water values. Using a separate y-axis, the relation $\mu = κ^{1/2}\tilde{h}$ is plotted for $\mu = 10, 2$ and $1$ (from the right to the left).</figcaption></figure>

eventually disappear when $\tilde{h} < \tilde{h}_c = \sqrt{3}$. It is noticed that these results match those of figure 2 in Barakat & Houston (1968). For water $\tilde{h}_c = \sqrt{3}$ corresponds to the particular water depth $h_0 \approx 4.76$ mm. In fact, when the surface tension is large, or equivalently, the fluid depth is sufficiently small so that $\tilde{h} < \tilde{h}_c$, the phase velocity is strictly a monotone function of the wavenumber, and thus no resonances are possible. In other words, when $\tilde{h} < \tilde{h}_c$ the asymptotic expansion of the solution is regular to any order, as demonstrated by Nayfeh (1970a) who considered the limit $\mu \ll 1$. In contrast, when $\tilde{h} > \tilde{h}_c = \sqrt{3}$ resonances between the fundamental and the nth harmonic are possible. Figure 1 also indicates that using the asymptotic theory with a constant value of $\mu$, the fundamental dimensional wavenumber $k_0$ being kept fixed, the first three singularities are found at the intersections between the corresponding thin solid line and the curves in thick solid lines. It should be emphasized that for a given $\mu$ each critical value of the wavenumber corresponds indeed to a different physical depth of the fluid layer. Nonetheless, for $\tilde{h} > 5$, the effects of the depth are negligible and the GC waves may be regarded as deep water phenomena.

It is possible to generalize the approach of Barakat & Houston (1968) to the case of GC waves with constant vorticity in finite depth. Namely, for steadily traveling periodic waves it is expected that the nth-order approximation obtained with the classical Stokes expansion method will break down when the wavenumber $k_0$ satisfies the following relation,

$$C_*(k_0) = C_*(nk_0)$$

or equivalently in dimensionless form, considering only forward propagating modes,

$$-\frac{\Omega}{2}\alpha_1 + \sqrt{\left(\frac{\Omega\alpha_1}{2}\right)^2 + (1 + \kappa)\alpha_1} = -\frac{\Omega}{2}\alpha_n + \sqrt{\left(\frac{\Omega\alpha_n}{2}\right)^2 + (1 + n^2\kappa)\alpha_n}, \quad (3.23)$$

where $n\alpha_n = \tanh(n\mu)$. With elementary relations for hyperbolic functions and some algebra, it can be shown that equation (3.23) yields the equations (A 11), (B 12) and (3.19) when
<figure><img src="image_13.png" /><figcaption>Figure 2: Variations of $κ_1, κ_2, κ_3$ as a function of depth for water for different values of the shear parameter: (a) $\tilde{\Omega} = \pm 0.5$; (b) $\tilde{\Omega} = \pm 0.9$. In any case $κ_1 > κ_2 > κ_3$; Dashed lines correspond to positive values of $\tilde{\Omega}$, dash-dot lines to negative ones and thin solid lines are drawn for $μ = 10, 2$ and $1$ (from the right to the left).</figcaption></figure>

$n = 2, 3$ and $4$ (respectively), and therefore represents the extension of equation (3.21) in the presence of current shear. Introducing $\tilde{\Omega} = \Omega_0/\sqrt{gk_c}$ and using the relation $\Omega = κ^{-1/4}\tilde{\Omega}$, one can see that equation (3.23) defines a nonlinear equation for $κ$ for fixed values of the two intrinsic dimensionless parameters, $\tilde{h}$ and $\tilde{\Omega}$. It should be emphasized here that in the analysis carried out by Hsu et al. (2016) the fully implicit definition of the critical wavenumbers was not properly recognized. More explicitly, the relation $\Omega = κ^{-1/4}\tilde{\Omega}$ was not taken into account when solving for the roots of equations (B11) and (B12). This gives rise to ambiguous results, as for instance the doubling of the critical wavenumbers at second and third order when $\Omega < 0$ for a given fixed "physical" dimensionless depth $\tilde{h}$ [see Figure 2 in Hsu et al. (2016)]. Though the plotted results in Figure 2 of Hsu et al. (2016) are correct, their interpretation may be misleading.

Numerical calculations were also performed for water ($T = 74$ dynes/cm) to solve equation (3.23) for $κ$. The results for the first three critical values are shown in figure 2 as a function of $\tilde{h}$ with four different values of the shear parameter $\tilde{\Omega}$. In any case it is seen that there is still a critical depth, $\tilde{h}_*$, below which resonances are not possible, although it depends on the value of the vorticity. When resonances are possible, namely when $\tilde{h} > \tilde{h}_*$, we always have $κ_1 > κ_2 > κ_3$ and their values are, of course, different from the irrotational case. Figure 2 shows that these critical values increase (decrease) with increasing positive (negative) vorticity. Interestingly, given that $κ_1$ decreases with increasing (negative) vorticity when $\tilde{h} > \tilde{h}_*$, the interval $[κ_1, ∞[$, over which the asymptotic theory is regular, also increases.

Actually the dependency of $\tilde{h}_*$ on the shear parameter $\tilde{\Omega}$ can be analyzed with the analytical results obtained by Martin (2013). As far as we are concerned with the forward propagating mode, we can use his Lemma 3 to write the appropriate relationship in dimensionless form, as

$$ \tilde{h}_*^2 = 3 + \tilde{\Omega}\tilde{h}_* \left[ \tilde{h}_* \sqrt{\tilde{h}_* + \frac{\tilde{\Omega}^2\tilde{h}_*^2}{4}} - \frac{\tilde{\Omega}\tilde{h}_*}{2} \right]. \quad (3.24) $$

Figure 3, in which the solution of equation (3.24) is plotted as a function of $\tilde{\Omega}$, shows that $\tilde{h}_* > \sqrt{3}$ when $\tilde{\Omega} > 0$ and $\tilde{h}_* < \sqrt{3}$ when $\tilde{\Omega} < 0$. In the region below the solid curve plotted in this figure, namely when $\tilde{h} < \tilde{h}_*$, resonances are not possible and thus the asymptotic theory is regular for any wavenumber. Notice that Gao et al. (2021) have also described the effects
<figure><img src="image_14.png" /><figcaption>Figure 3: Normalized critical depth $\tilde{h}_*$ as a function of $\tilde{\Omega}$</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 4: Variations of $\kappa_1, \kappa_2, \kappa_3$ for different values of $\tilde{h}$ for water as a function of the shear parameter $\tilde{\Omega}$: (a) $\tilde{h}/\sqrt{3} = 2$; (b) $\tilde{h}/\sqrt{3} = 0.9$. In any case $\kappa_1 > \kappa_2 > \kappa_3$. The dotted lines represent the corresponding well known values in deep water.</figcaption></figure>

of vorticity on the properties of the dispersion relation and obtained identical results for the forward propagating mode. Actually the first case in their inequalities (13) corresponds to the region in figure 3 where $\tilde{h} < \tilde{h}_*$, which they have named shallow-water regime. Figure 3 also reveals that for any given value of the physical depth, or equivalently $\tilde{h}$, there exist a threshold value of $\tilde{\Omega}$ above which there are no resonances. This particular value is positive when $\tilde{h} > \sqrt{3}$, but negative when $\tilde{h} < \sqrt{3}$. This is further illustrated in figure 4, where the critical values $\kappa_{n-1}$ ($n = 2, 3$ and 4) are plotted as a function of $\tilde{\Omega}$ for two particular values of the depth.

## 4. Numerical results

In this section we first present the numerical method used to compute steadily traveling wave solutions of the Euler equations with surface tension and constant vorticity in finite depth. Then, by comparison, we study the validity of the third and fourth-order solutions in two cases, namely, (a) GC waves without vorticity and (b) GC waves on linear sheared currents.

In the most general case, a family of steadily traveling wave solutions lies on a region of a four-parameter space, which involves a quadruplet of dimensionless parameters, either $(\kappa, \mu, \Omega, \epsilon)$ or $(\kappa, \tilde{h}, \tilde{\Omega}, \epsilon)$. We shall use the former quadruplet in this section, like in the recent
works of Hsu et al. (2018) and Dhar & Kirby (2023) concerning the stability analysis of GC waves in finite depth and constant vorticity within the framework of envelope equations. It is worth noting that with our choice of scaling the two quadruplets are equivalent only when $κ = 1$. As explained in section 3.2, when $κ ≠ 1$ the latter quadruplet (with $ϵ = 0$) is preferable to analyze unambiguously the effects of depth and vorticity on the occurrence of critical values of $κ$ for the perturbation scheme.

## 4.1. Fully nonlinear solutions

To compute steadily traveling periodic GC waves with constant vorticity in finite depth, we shall use a numerical method based on a conformal mapping method and Fourier series expansions of the unknowns, such as given by Guo et al. (2014) and Choi (2009) for two-dimensional GC solitary waves and gravity waves, respectively. It is emphasized here that, though these authors established the governing equations for GC waves in finite depth with constant vorticity, the numerical computations were carried out only in the deep water case. The formulation of the governing equations in the transformed plane (a canonical strip) is well documented by these authors, and therefore need not to be reproduced here.

To be concise and consistent, however, we present the transformed governing equations used in our study and adopt the same notations as Ribeiro et al. (2017) except for the coordinates in the canonical domain. Namely the conformal map is defined as $\tilde{Z}(\zeta, \vartheta) = \tilde{X}(\zeta, \vartheta) + i\tilde{Y}(\zeta, \vartheta)$ and the surface elevation as $z = \eta(X) = Y(\zeta)$ where $X = \tilde{X}(\zeta, 0)$ and $Y = \tilde{Y}(\zeta, 0)$. After performing similar calculations as those presented in Choi (2009), we can show that the free surface kinematic and dynamic boundary conditions can be combined into a single pseudo-differential equation for the free surface elevation $Y(\zeta)$,

$$ \begin{aligned} & -\frac{\tilde{c}^2}{2} - \frac{\tilde{c}^2}{2J} + Y + \gamma^2 \frac{(C[(Y+b)Y_\zeta]) - ((Y+b)Y_\zeta)^2}{2J} - \gamma\tilde{c}b \\ & + \frac{(\tilde{c} + \gamma C[(Y+b)Y_\zeta]) (\tilde{c} + \gamma (Y+b)(1-C[Y_\zeta]))}{J} + \tilde{\kappa} \frac{X_\zeta Y_\zeta \tilde{c} - Y_\zeta X_\zeta \tilde{c}}{J^{3/2}} = \tilde{B}, \quad (4.1) \end{aligned} $$

where $J = \tilde{X}_\zeta^2 + \tilde{Y}_\zeta^2$ is jacobian of the conformal map evaluated at the free surface. According to our choice of scaling (see section 2), we have $\tilde{c} = C/\sqrt{\mu}$, $\tilde{\kappa} = \kappa/\mu^2$, $\gamma = -\sqrt{\mu\Omega}$, and $\tilde{B}$ is the Bernoulli constant. As explained in Ribeiro et al. (2017), $b$ is a free parameter that determines the choice of the reference frame, or equivalently the depth at which the Eulerian mean horizontal velocity is zero in the limit of small amplitude waves. Since the governing equations have the property of Galilean invariance, the parameter $b$ affects neither the shape of the wave nor the streamlines in the frame of reference traveling with the wave. For the comparison with the analytical results derived in section 3 we shall take $b = 0$. Equation (4.1) contains three additional unknowns, the phase velocity $\tilde{c}$, the Bernoulli constant $\tilde{B}$ and the conformal depth $\tilde{d}$, the later being embedded in the pseudo-differential operator $C$, defined as $C[f] = \mathcal{F}^{-1}[i \coth(k\tilde{d})\hat{f}_k]$ where $\hat{f}_k = \mathcal{F}[f]$ and $\mathcal{F}$ represents the Fourier transform operator. Accordingly, we also have the relation

$$ X_{\zeta} = 1 - C[Y_{\zeta}]. \quad (4.2) $$

To obtain a closed system of equations for the unknowns, equation (4.1) must be supplemented by three additional equations. We choose to impose a zero mean level in the physical space,

$$ \int_{0}^{L/2} Y X_{\zeta} d\zeta = 0 \quad (4.3) $$
and, assuming the wave crest is at $\zeta = 0$, we fix the dimensionless wave height through

$$458 \qquad Y(0) - Y(L/2) = H, \qquad (4.4)$$

as well as the depth condition

$$461 \qquad \tilde{d} = 1 + \frac{1}{L} \int_{-L/2}^{L/2} Y d\zeta. \qquad (4.5)$$

A Fourier spectral discretization for the unknown function $Y(\zeta)$ allows one to obtain an efficient numerical method with fast Fourier transform, wherein all derivatives and pseudo-differential operators are calculated via Fourier multipliers while the nonlinear terms are computed pseudo-spectrally (with dealiasing). Hence we consider a truncated Fourier series with $N$ modes,

$$468 \qquad Y(\zeta) = \sum_{n=-N/2}^{N/2} \hat{Y}_n e^{in\frac{2\pi}{L}\zeta}, \qquad \hat{Y}_0 = 0, \qquad (4.6)$$

and assume that the solutions are periodic with a wavelength $L = 2\pi/\mu$ equal to the dimensionless wavelength in the physical space. Note that $\hat{Y}_{-n} = \hat{Y}_n^*$ since the surface elevation is a real function. Taking $N$ even and considering only symmetric waves, we discretize equation (4.1) at $N/2 + 1$ collocation points uniformly distributed over a half-period, $\zeta_j = (j-1)L/N$ with $j = 1, \dots, N/2 + 1$, to obtain $N/2 + 1$ equations for the unknown values $Y_j = Y(\zeta_j)$, $j = 1, \dots, N/2 + 1$. Performing the discretization of the additional equations (4.3-4.5), we can close the system for the $N/2+4$ unknowns $\tilde{c}, \tilde{B}, \tilde{d}$ and $Y_j$, $j = 1, \dots, N/2 + 1$.

To solve the resulting nonlinear system of algebraic equations, we use the Levenberg-Marquardt algorithm implemented in the fsolve function of the Optimisation toolbox of MATLAB. As an initial guess for this iterative procedure we use the first-order solution presented in section 3,

$$482 \qquad \tilde{c} = C_0/\mu, \qquad \tilde{B} = 0, \qquad \tilde{d} = 1, \qquad Y(\zeta) = H \cos\left(\frac{2\pi}{L}\zeta\right), \qquad (4.7)$$

where $H = \epsilon/\mu$ with our choice of scaling and $\epsilon$ is small. Unlike Ribeiro et al. (2017), who used the mean absolute residuals errors for the whole system of equations, our stopping criterion is based on the infinite norm of the residuals of equation (4.1). In our practice, with $N = 128$ or 256 both errors were found to be of the order $10^{-12}$.

The following numerical results were obtained by natural numerical continuation in the parameters space, using the prior converged solution as the initial guess to a new solution. Owing to the number of dimensionless parameters characterizing branches of solutions, it is obviously a difficult task to give a thorough description of the bifurcation diagrams. Even in the simpler case of irrotational flows in deep water, namely with $\Omega = 0$ and $\mu = \infty$, Shelton et al. (2021) have shown that the bifurcation space of the gravity-capillary problem is certainly non-trivial, rendering difficult to observe any clear structure. Crucially, this is due to the existence of multiple solutions, the so-called Wilton ripples or combined waves, which emerge from resonant interactions between two modes. On account of this difficulty, we take $\mu$ and $\kappa$ as constants on a branch of solutions and vary $\Omega$ and $\epsilon$ over their admissible range in most of the cases studied here.

## 4.2. Gravity-capillary waves in irrotational flows

We first study irrotational flows with $\mu = 10$ and several values of $\kappa$. The values $\kappa = 1$ and $0.8$ are taken to verify the correctness of the analytical solutions for GC waves, in
<figure><img src="image_16.png" /><figcaption>Figure 5: Plot of phase velocity $C$ against $\epsilon$ for $\mu = 10$, $\Omega = 0$ and some values of $\kappa = 1.0, 0.80, 0.58, 0.42, 0.30$ (from top to bottom); solid lines show the numerical results, dashed lines the analytical results and the boxed number indicates the type of the waves</figcaption></figure>

502 which the effects of both gravity and capillarity are equally important, whereas the values of $\kappa = 0.58, 0.42$ and $0.30$ are used to determine the behavior of the weakly nonlinear solutions
503 close to the critical wavenumbers.
504

505 The variations of phase velocity $C$ as a function of $\epsilon$ for $\Omega = 0$ and several values of $\kappa$ are shown in figure 5. The solid lines represent the numerical solutions of the Euler equations, and the dashed lines correspond to the approximate solutions. According to the asymptotic theory, the phase velocity decreases with the increase of $\epsilon$ for $\kappa = 1, 0.8$ and $0.58$ much like pure capillary waves, whereas it increases with $\epsilon$ for $\kappa = 0.42$ and $0.30$ much like pure gravity waves. For $\kappa = 1$ and $0.8$ the approximate solution agrees fairly well with the exact numerical solution over the range of $\epsilon$ values shown here. In contrast, with $\kappa = 0.58, 0.42$ and $0.30$ the approximate solution matches the exact solution only for small values of the amplitude parameter $\epsilon$. 511 512 513

514 To better understand the meaning of the curves in figure 5, it should be realized that in 515 general there is no unique branch of solutions for a fixed $\kappa$, in the neighbourhood or not 516 of the critical wavenumbers. This was demonstrated numerically by the seminal work of 517 Schwartz & Vanden-Broeck (1979) in the deep water case. With our numerical method and 518 continuation in either the amplitude parameter $\epsilon$ or the surface tension parameter $\kappa$, we find 519 that this is also the case in finite depth, as illustrated by figure 6 where the variations of the 520 squared phase velocity $C^2$ with $\kappa$ are shown. 521

521 In this figure three branches of numerical solutions with $\mu = 10$ and $\epsilon = 0.10$ are 522 plotted with solid lines, as well as the corresponding weakly nonlinear solution (dash-dotted 523 line obtained with equation (B 11)) and the linear solutions for infinitesimal waves with 524 (dimensional) wavenumber $nk_0$ for $n = 1, 2$ and $3$ (dashed lines obtained with equation 525 3.5). According to the classification of the different branches of solutions proposed in Schwartz & 526 Vanden-Broeck (1979), the types 1, 2, and 3 are represented in this figure, and correspond to 527 the number of observed 'dimples' or inflexion points on a (half) wave profile. Like in the deep 528 water case, for much of the figure the type number is seen to increase with $C^2$ for given $\kappa$. 529 For this value of the amplitude parameter, $\epsilon = 0.10$, the branches of types 1 and 3 follow the 530 linear solution closely as $\kappa$ increases from the critical values $\kappa_1$ and $\kappa_2$, respectively, whereas 531 the wave profiles become indistinguishable from a sinusoid with wavenumbers 1 and 3 (in
<figure><img src="image_17.png" /><figcaption>Figure 6: Plot of $C^2$ against $\kappa$ for $\mu = 10$, $\epsilon = 0.10$ and $\Omega = 0$. Solid line: numerical solutions; dash-dotted line: weakly nonlinear analytical solution; dashed lines: infinitesimal wave solutions.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 7: Comparison of surface profiles between third-order (dotted line), fourth-order (dashed line) and exact numerical solutions (solid line) for $\Omega = 0$, $\mu = 10$ and some values of $\epsilon = 0.3, 0.2, 0.1$ (from top to bottom); (a) $\kappa = 1$, (b) $\kappa = 0.8$.</figcaption></figure>

<footnote n="532">532 units of $k_0$), respectively. In contrast, the branch of type 2 reaches a limiting configuration, as $\kappa$ increases, through the trapping of a bubble (not shown). As $\kappa < \kappa_1$ decreases, all the branches reach a limiting configuration through the trapping of one or more bubbles (not shown). As expected the approximate solution for $C^2$, the dash-dotted line in figure 6, breaks down as $\kappa$ approaches the critical value $\kappa_1$. An interesting feature of this curve, however, is that for decreasing values of $\kappa < \kappa_1$ it approaches first the branch of type 2, and matches almost perfectly the branch of type 3 when $0.25 < \kappa < 0.38$.</footnote>

<footnote n="539">539 Knowing that initialization of the numerical continuation method with the linear solution with $\kappa$ in the open intervals $]\kappa_1, \infty[$ and $]\kappa_n, \kappa_{n-1}[$ (with $n \ge 2$) yields numerical solutions corresponding to the branch of type $n$, one can realize that, in figure 5, the branches with $\kappa = 1, 0.8$ and $0.58$ correspond to waves of type 1, while the branches with $\kappa = 0.42$ and $0.30$ correspond to waves of the type 2 and 3, respectively. This analysis is supported by the examination of typical surface profiles for several values of $\epsilon$ and $\kappa$, as shown in figures 7 and 8.</footnote>
<figure><img src="image_19.png" /><figcaption>Figure 8: Comparison of surface profiles between third-order (dotted line), fourth-order (dashed line) and exact numerical solutions (solid line) for $\Omega = 0$, $\mu = 10$ and some values of $\epsilon = 0.1, 0.05, 0.025$ (from top to bottom); (a) $\kappa = 0.58$, (b) $\kappa = 0.42$, (c) $\kappa = 0.30$.</figcaption></figure>

Figure 7 shows that the trough becomes deeper and narrower as $\kappa$ decreases from 1 to 0.8 and $\epsilon$ increases from 0.1 to 0.3. It is to be noted that it displays only wave profiles corresponding to the type 1. For these waves where surface tension balances gravity, we observe that the fourth-order analytical solution provides better results than the third-order solution for predictions of the surface elevation and is in good agreement with the exact numerical solution when $\epsilon \le 0.2$.

Figure 8 shows for smaller values of $\epsilon$ the surface profiles of near resonant waves with $\kappa = 0.58, 0.42$ and $0.30$. For waves with $\kappa = 0.58$, figure 8a shows that the troughs tend to sharpen with the increase of $\epsilon$, whereas the crests tend to flatten, which is typical of the waves of type 1. For the largest amplitude, $\epsilon = 0.1$, it is observed that the fourth-order solution makes a somewhat better job than the third-order solution, in the sense that it does not predict the formation of a crest dimple. Quantitatively, however, there is no significant improvement, neither near the crest nor near the trough. For waves with $\kappa = 0.42$ the numerical solution indicates the presence of a secondary maximum on the surface profile with a trough between the two maxima, which is typical of the waves of type 2 (see figure 8b). For the largest amplitude, $\epsilon = 0.1$, neither of the asymptotic solutions is in agreement with the numerical solution. In fact, the fourth-order solution incorrectly predicts the occurrence of a crest dimple that is absent in the profiles of the third-order and numerical solutions. With $\kappa = 0.30$, the plots of figure 8c show that when the amplitude increases, the numerical solution takes the appearance of a primary sine wave with a relatively small 3-cycle perturbation superimposed upon it, which is typical of the waves of type 3. Similarly, neither of the asymptotic solutions is in agreement with the numerical solution for the largest amplitude considered here, though they predict qualitatively waves of type 3. It is noteworthy that when $\kappa$ is close to $\kappa_1$, like in
<figure><img src="image_20.png" /><figcaption>Figure 9: Comparison of horizontal velocity profiles between third-order (dotted line), fourth-order (dashed line) and exact numerical solutions (solid line) for the waves plotted in figure 7a. The profiles beneath the crest(trough) are plotted on the left(right) and the steepness increases from top to bottom ($\epsilon = 0.1, 0.2$ and $0.3$).</figcaption></figure>

figures 8a and 8b, the fourth-order approximations fail more importantly than the third-order ones for the larger steepness, due to the fact that the associated small divisors in equations (3.13-3.17) are found to be squared in the fourth-order coefficients.

For further insight into the differences between the third and fourth order solutions, we have also compared the predictions of the interior flow with the numerical solutions. The method for the computation of the interior flow is explained in Ribeiro et al. (2017) and, therefore, need not be detailed here. To illustrate the differences in the interior flow, we have focused on the wave-induced horizontal velocities beneath the crest and the trough of the waves plotted in figure 7a, where $\kappa = 1$ is not close to $\kappa_1$ and $\epsilon = 0.1, 0.2$ and $0.3$. Figure 9 shows the comparison between the asymptotic solutions and the exact results over the dimensionless range $y = -1$ to the free surface. Note that the bottom is located at $y = -10$, where the wave-induced velocities vanish (not shown).

For the smallest steepness, it is found that the fourth-order predictions are closer to the exact results both beneath the crest and the trough. However, for the larger steepness the differences between the third-order and the fourth-order are much more important beneath the crest than beneath the trough. Actually, below the trough it is observed that those differences increase with the steepness and the fourth-order predictions are closer to the exact results. Beneath the crest, those differences also increase with the steepness, but neither such approximations are close to the exact results in the near surface region. Nonetheless, as the depth increases the fourth-order solutions are found to be better than the third-order ones. These results can be explained, noting that, firstly the local slopes at the trough are in closer agreement with the exact results than they are the crest, and secondly that a binomial expansion has been used in the dynamic boundary condition to approximate the nonlinear term related to the
<figure><img src="image_21.png" /><figcaption>Figure 10: Plot of phase velocity C against $\epsilon$ for $\mu = 10$ and some values of $\kappa = 1.0, 0.80, 0.58, 0.42, 0.30$ (from top to bottom): solid lines show the numerical results, dashed lines the analytical results and the boxed number indicates the type of the waves;</figcaption></figure>

(a) $\Omega = 0.5$, (b) $\Omega = 0.9$, (c) $\Omega = -0.5$ and (d) $\Omega = -0.9$.

surface tension. Only two terms in this binomial expansion are required for the third-order approximation, but three terms are necessary for the fourth-order approximation. Due to this approximation there is a mismatch between the approximate wave-induced horizontal velocities and the exact results, though the differences in the phase velocity are very small over the range of steepness considered here (see figure 5, the case $\kappa = 1$).

### 4.3. Gravity-capillary waves in flows with constant vorticity

Like in the previous subsection, we start by analyzing the nonlinear dependency of the phase velocity on the wave steepness for the same different values of $\kappa$ but with different values of the vorticity parameter $\Omega$. According to the numerical solutions the dependency of $C$ on $\epsilon$ is strongly affected by the vorticity whatever the value of $\kappa$, as shown in figure 10. In comparison with the case of irrotational flow (see figure 5), not only the values of $C$ in the small amplitude limit are modified, but also the type of the waves depending on the value of $\kappa$. In figure 10 the solid lines represent the numerical solutions of the Euler equations, and the dashed lines correspond to the approximate nonlinear dispersion relation presented in section 3. As expected from the linear dispersion relation (equation 3.5), it is observed that in the limit of small amplitude increasing negative vorticity ($\Omega > 0$) reduces the phase velocity, whereas increasing positive vorticity enhances it. In this limit the analytical results are in excellent agreement with the exact results whatever the value of $\kappa$.

When $\Omega > 0$, figures 10a and 10b show that for each value of $\kappa$, the phase velocity is a decreasing function of the steepness. Except for the two lowest values of $\kappa$, namely 0.30 and 0.42, the agreement between the asymptotic theory and the exact theory is excellent over
<table><caption>Table 1: Values of $ħ$ and $Ω̃$ for given values of $κ$ and $Ω$. Here $μ = 10$ for each wave. The deep water values of the first eight critical $κ_n$ are also reported, as well as the type of each wave.</caption><thead><tr><th>κ</th><th>Ω</th><th>ħ</th><th>Ω̃</th><th>κ₁</th><th>κ₂</th><th>κ₃</th><th>κ₄</th><th>κ₅</th><th>κ₆</th><th>κ₇</th><th>κ₈</th><th>type</th></tr></thead><tbody><tr><td rowspan="4">1.00</td><td>0.50</td><td rowspan="4">10.00</td><td>0.50</td><td>0.21</td><td>0.13</td><td>0.09</td><td>0.07</td><td>0.06</td><td>0.05</td><td>0.04</td><td>0.03</td><td>1</td></tr><tr><td>0.90</td><td>0.90</td><td>0.08</td><td>0.05</td><td>0.03</td><td>0.02</td><td>0.02</td><td>0.01</td><td>0.01</td><td>0.01</td><td>1</td></tr><tr><td>-0.50</td><td>-0.50</td><td>0.92</td><td>0.63</td><td>0.48</td><td>0.39</td><td>0.33</td><td>0.29</td><td>0.26</td><td>0.23</td><td>1</td></tr><tr><td>-0.90</td><td>-0.90</td><td>1.34</td><td>0.92</td><td>0.71</td><td>0.59</td><td>0.50</td><td>0.44</td><td>0.39</td><td>0.36</td><td>2</td></tr><tr><td rowspan="4">0.80</td><td>0.50</td><td rowspan="4">11.18</td><td>0.47</td><td>0.22</td><td>0.14</td><td>0.10</td><td>0.08</td><td>0.06</td><td>0.05</td><td>0.04</td><td>0.04</td><td>1</td></tr><tr><td>0.90</td><td>0.85</td><td>0.09</td><td>0.06</td><td>0.04</td><td>0.03</td><td>0.02</td><td>0.02</td><td>0.01</td><td>0.01</td><td>1</td></tr><tr><td>-0.50</td><td>-0.47</td><td>0.90</td><td>0.61</td><td>0.47</td><td>0.38</td><td>0.32</td><td>0.28</td><td>0.25</td><td>0.23</td><td>2</td></tr><tr><td>-0.90</td><td>-0.85</td><td>1.28</td><td>0.88</td><td>0.68</td><td>0.56</td><td>0.48</td><td>0.42</td><td>0.37</td><td>0.34</td><td>3</td></tr><tr><td rowspan="4">0.58</td><td>0.50</td><td rowspan="4">13.13</td><td>0.44</td><td>0.24</td><td>0.15</td><td>0.11</td><td>0.08</td><td>0.07</td><td>0.06</td><td>0.05</td><td>0.04</td><td>1</td></tr><tr><td>0.90</td><td>0.79</td><td>0.11</td><td>0.07</td><td>0.05</td><td>0.03</td><td>0.03</td><td>0.02</td><td>0.02</td><td>0.01</td><td>1</td></tr><tr><td>-0.50</td><td>-0.44</td><td>0.86</td><td>0.59</td><td>0.45</td><td>0.37</td><td>0.31</td><td>0.27</td><td>0.24</td><td>0.22</td><td>3</td></tr><tr><td>-0.90</td><td>-0.79</td><td>1.21</td><td>0.83</td><td>0.64</td><td>0.53</td><td>0.45</td><td>0.39</td><td>0.35</td><td>0.32</td><td>4</td></tr><tr><td rowspan="4">0.42</td><td>0.50</td><td rowspan="4">15.43</td><td>0.40</td><td>0.25</td><td>0.16</td><td>0.12</td><td>0.09</td><td>0.07</td><td>0.06</td><td>0.05</td><td>0.05</td><td>1</td></tr><tr><td>0.90</td><td>0.72</td><td>0.13</td><td>0.08</td><td>0.05</td><td>0.04</td><td>0.03</td><td>0.03</td><td>0.02</td><td>0.02</td><td>1</td></tr><tr><td>-0.50</td><td>-0.40</td><td>0.83</td><td>0.57</td><td>0.43</td><td>0.35</td><td>0.30</td><td>0.26</td><td>0.23</td><td>0.21</td><td>4</td></tr><tr><td>-0.90</td><td>-0.72</td><td>1.15</td><td>0.79</td><td>0.61</td><td>0.50</td><td>0.42</td><td>0.37</td><td>0.33</td><td>0.30</td><td>5</td></tr><tr><td rowspan="4">0.30</td><td>0.50</td><td rowspan="4">18.26</td><td>0.37</td><td>0.27</td><td>0.17</td><td>0.13</td><td>0.10</td><td>0.08</td><td>0.07</td><td>0.06</td><td>0.05</td><td>1</td></tr><tr><td>0.90</td><td>0.67</td><td>0.14</td><td>0.09</td><td>0.06</td><td>0.05</td><td>0.04</td><td>0.03</td><td>0.03</td><td>0.02</td><td>1</td></tr><tr><td>-0.50</td><td>-0.37</td><td>0.80</td><td>0.55</td><td>0.42</td><td>0.34</td><td>0.29</td><td>0.25</td><td>0.22</td><td>0.20</td><td>5</td></tr><tr><td>-0.90</td><td>-0.67</td><td>1.09</td><td>0.75</td><td>0.57</td><td>0.47</td><td>0.40</td><td>0.35</td><td>0.31</td><td>0.28</td><td>8</td></tr></tbody></table>

the range of steepness considered here, namely $0 < ϵ ≤ 0.25$. In any case, the exact results, as well as the corresponding approximations, show that the phase velocity is a decreasing function of steepness. This suggests that the waves are capillary-like. For the calculations with $κ = 0.30$ and $0.42$, it appears that the mismatch between the approximations and the exact results is reduced with the increase of $Ω$, the curves becoming closer to each other over a larger range of steepness. In contrast, when $Ω < 0$ the phase velocity is an increasing function of steepness for all $κ$ values, except for the branch of solutions with $κ = 1$ and $Ω = -0.5$ that is plotted in figure shown in 10c. In this figure, the exact branches of solutions with $κ = 1$ and $0.8$ agree with the approximations only when $ϵ < 0.1$. In figure 10d, $Ω = -0.9$ and the agreement between the approximations and the exact results is very satisfactory over the range of steepness shown here, namely $0 < ϵ ≤ 0.25$.

To identify the type of waves for each branch of solutions plotted in figure 10 we may adopt the reasoning used in section 4.2 for the GC waves in irrotational flows. This in turn requires to determine the critical values of $κ$ for given values of the intrinsic parameters $ħ$ and $Ω̃$. In our case $μ = 10$ and we have seen in section 3.2 that these critical values are mainly affected by the vorticity, the depth effects being negligible when $ħ ≫ ħ_*$ . For each value of $κ$ and $Ω$ used in figure 10, we have computed the first height $κ_n$. These values are reported in table 1. It should be emphasized that for given values of $μ$ and $Ω$, the intrinsic parameters $ħ$ and $Ω$ have different values that depend on the value of $κ$. The last column of table 1 indicates the type of the branch of solution for each triplet $(μ, κ, Ω)$ considered in figure 10. Solutions
<figure><img src="image_22.png" /><figcaption>Figure 11: Comparison of surface profiles and wave-induced horizontal velocity profiles under crest ($u_c$) and trough ($u_t$) for waves with $\mu = 10$, $\Omega = 0.5$, $\epsilon = 0.10$ and several values of $\kappa$. Third-order solution (dotted line), fourth-order solution (dashed line) and numerical solution (solid line). (a) $\kappa = 1$, (b) $\kappa = 0.80$, (c) $\kappa = 0.58$, (d) $\kappa = 0.42$ and (e) $\kappa = 0.30$.</figcaption></figure>

<footnote n="633">633 are of type $n + 1$ when the values of $\kappa$ falls in the interval $[\kappa_{n+1}, \kappa_n]$ for $n \ge 1$. Waves of type 634 1 have wavelengths such that $\kappa > \kappa_1$, where $\kappa_1$ is the largest critical wavenumber.</footnote>

<footnote n="635">635 For several values of $\kappa$ and a fixed wave steepness, $\epsilon = 0.10$, we have also compared 636 the surface profiles and the wave-induced horizontal velocities in the presence of a linear 637 shear current with $\mu = 10$. Figure 11 shows the comparison between the analytical solutions 638 and the exact results for waves of different wavelengths, namely $\kappa = 1, 0.80, 0.58, 0.42$ 639 and 0.30, and a vorticity parameter $\Omega = 0.5$. According to table 1 all the waves are of type 640 1, since their wavelength is such that $\kappa > \kappa_1$. Except for the longest wave with $\kappa = 0.30$, 641 there is a remarkable agreement between the analytical results and the exact ones, the</footnote>
<figure><img src="image_23.png" /><figcaption>Figure 12: Same legend as in figure 11 but with $\Omega = -0.5$.</figcaption></figure>

fourth-order results being more accurate than the third-order ones. As expected neither such approximations is valid for the waves with $\kappa = 0.30$, a value that is very close to $\kappa_1 \approx 0.27$.

The same comparison has been carried out with $\Omega = -0.5$, and the results are shown in figure 12. According to table 1, each wave is of a different type. The type number is seen to increase with increasing wavelength. The shortest wave ($\kappa = 1$) is the only capillary-like wave, in the sense that it belongs to a branch of solutions with decreasing phase velocity with increasing steepness (see figure 10c), and it is of type 1 as shown in figure 12a. Again, for this value of $\kappa$ that is close to the $\kappa_1 \approx 0.92$, both the third-order and fourth-order solutions differ significantly from the exact results, the disagreement being larger near the crest than near the trough. For the wave of type 2 shown in figure 12b we have $\kappa = 0.80$, which falls near the middle of the interval [$\kappa_2, \kappa_1$] = [0.61, 0.90]. Looking at the wave-induced horizontal velocities, it is found that the fourth-order theory does a better job than the third-order one beneath the trough, and that the situation is reversed beneath the crest. Nonetheless, both analytical solutions fail to match with the exact results in the near surface region. Actually,
<figure><img src="image_24.png" /><figcaption>Figure 13: Comparison of surface profiles and wave-induced horizontal velocity profiles under crest ($u_c$) and trough ($u_t$) for waves with $\mu = 10$, $\epsilon = 0.20$. Third-order solution (dotted line), fourth-order solution (dashed line) and numerical solution (solid line). (a) $\kappa = 1$, $\Omega = 0.5$, (b) $\kappa = 0.80$, $\Omega = 0.5$, (c) $\kappa = 0.42$, $\Omega = -0.5$, (d) $\kappa = 0.30$, $\Omega = -0.5$.</figcaption></figure>

these differences between the exact results and the approximations are expected in view of figure 10c, which shows that there is mismatch in the prediction of C for the waves with $\kappa = 1$ and 0.80.

In contrast, for the gravity-like waves with $\kappa = 0.58$, 0.42 and 0.30, figure 10c shows an excellent agreement between the exact phase velocity and the approximate one when $\epsilon = 0.1$. Not surprisingly, when the value of $\kappa$ is close to one of the first three critical values $\kappa_n$, both approximations fail to reproduce the near surface horizontal velocities (see figures 12c and 12d). However, when $\kappa = 0.42$ we have $|\kappa - \kappa_2| = 0.15$, $|\kappa - \kappa_3| = 0.01$, and although the approximate surface profiles are almost indistinguishable from the exact solution, it is found that there are important differences in the internal flow structure. Figure 12d shows that the fourth-order theory does a better job than the third-order one beneath the trough, the situation being reversed beneath the crest. As shown in figure 12e, the surface profiles and the wave-induced horizontal velocities are indistinguishable from the exact solution for the wave with $\kappa = 0.32$, which is in the interval $[\kappa_5, \kappa_4]$ and far enough the smallest critical value of either the third-order or fourth-order approximations.

Finally we present the results of a similar comparison for steeper waves with $\mu = 10$, $\Omega = \pm 0.5$ and $\epsilon = 0.2$. The four different waves have been chosen so that the approximation of the phase velocity is still in very good agreement with the exact solution. In view of the results of figures 10a and 10c, we consider $\kappa = 1$ and $\kappa = 0.8$ with $\Omega = 0.5$, whereas with $\Omega = -0.5$ we choose $\kappa = 0.42$ and $\kappa = 0.30$. For the waves with $\Omega = 0.5$, $\kappa$ is quite far from
the first critical value $\kappa_1$, and we observe that the third- and fourth-order surface profiles are almost indistinguishable from the exact solution. However, in each case the approximations fail to describe the wave-induced velocities close to the surface. Nonetheless it should be noticed that the fourth-order solution performs better than the third-order one, though this improvement is lessened beneath the crest. For the waves with $\Omega = -0.5$ and $\kappa = 0.42$, which falls in $[\kappa_4, \kappa_3] = [0.35, 0.43]$, figure 13c shows that nonlinearity has become sufficiently important to reveal the defect of the fourth-order theory in comparison with the third-order theory and the exact results (see also for comparison figure 12d). For the longest wave with $\kappa = 0.30$, shown in figure (13d), the approximate surface profiles are in very good agreement with the exact results. In this case the coefficients of the asymptotic series are well ordered, $\kappa$ being $O(\epsilon)$ far from $\kappa_1 = 0.80$, $\kappa_2 = 0.55$ and $\kappa_3 = 0.42$, though they fail to predict the wave-induced velocities close to the surface beneath both crest and trough. Actually there is no improvement using the fourth-order theory and its results are very similar to that obtained with the third-order theory.

## 5. Conclusion

In this paper we have derived a fourth-order asymptotic solution for weakly nonlinear GC waves with constant vorticity in finite depth and thereby, the works of Barakat & Houston (1968) and Hsu et al. (2016) have been extended.

We have shown that the determination of the critical values in the perturbation scheme, $\kappa_n$ (for any $n \ge 1$), should be done with two intrinsic dimensionless parameters $\tilde{h}$ and $\tilde{\Omega}$. It is emphasized here that, in the presence of surface tension, using instead $\mu$ and $\Omega$ may be misleading in the study of the effects of both vorticity and depth on the occurrence of these critical wavenumbers. We have found that vorticity changes the critical depth $\tilde{h}_*$, below which no resonances are possible and, thus, the perturbation scheme is valid for any wavelength specified with $\kappa$. This critical depth increases without bound with increasing positive values of $\tilde{\Omega}$ (negative vorticity). In contrast it decreases slowly with increasing positive vorticity.

To test the accuracy of the third-order and fourth-order solutions, we have used a numerical method of solution based on truncated Fourier series approximation in a conformal formulation of the problem. The fourth-order results are interesting as they reveal important differences not only in the surface profiles but also in the internal flow structure, when compared with the third-order predictions. With negative vorticity ($\Omega > 0$), the fourth-order solution improves on the third-order solution, inasmuch as the steepness is relatively small and the dimensionless squared wavenumber is far enough the first three critical values. With positive vorticity ($\Omega < 0$), the outcome of the comparison between the third and fourth order is rather surprising. Our analysis reveals that there are important differences in the wave-induced velocities below the free surface, even though the predicted surface profiles as well as phase velocities are in good agreement with the numerical solutions.

Finally, we point out that we have not tested the approximations in the shallow water limit, $\mu \ll 1$. We merely recall that without vorticity the classical Stokes expansion method become less accurate when both the amplitude $\epsilon$ and the Ursell number $\epsilon/\mu^3$ significantly increases. This usually implies a severe restriction on the steepness when moving to shallow water. How this restriction is affected, as vorticity effects become important in the presence of surface tension, remains to be explored.

## Acknowledgements

This work was started during a visit by Prof. A. K. Dhar to the University of Turin, Italy. During this work, Prof. Dhar was supported by a CPDA grant from the IIEST, Shibpur under Grant No. 1593/D(FW)/22. The award of the Senior Research Fellowship (Ref No:
724 276/ACAD/2023) by the IIEST, Shibpur to S. Halder is gratefully acknowledged.

725

**Declaration of interests.**

726

727 The authors report no conflict of interest.

REFERENCES

728 BARAKAT, R. & HOUSTON, A. 1968 Nonlinear periodic capillary-gravity waves on a fluid of finite depth. *J. Geophys. Res.* **73** (20), 6545–6554.

729 BRANTENBERG, C. & BREVIK, I. 1993 Higher order water waves in currents of uniform vorticity, in the presence of surface tension. *Phys. Scr.* **47** (3), 383.

730 CHABANE, M. & CHOI, W. 2019 On resonant interactions of gravity-capillary waves without energy exchange. *Stud. Appl. Math.* **142** (4), 528–550.

731 CHOI, W. 2009 Nonlinear surface waves interacting with a linear shear current. *Math. Comp. Simul.* **80** (1), 29–36.

732 CURTIS, C.W., CARTER, J.D. & KALISCH, H. 2018 Particle paths in nonlinear schrödinger models in the presence of linear shear currents. *J. Fluid Mech.* **855**, 322–350.

733 DHAR, A. K. & KIRBY, J. T. 2023 Fourth-order stability analysis for capillary-gravity waves on finite-depth currents with constant vorticity. *Phys. Fluids* **35** (2).

734 DJORDJEVIC, V. & REDEKOPP, L. G. 1977 On two-dimensional packets of capillary-gravity waves. *J. Fluid Mech.* **79** (4), 703–714.

735 FANG, H., LIU, P.L.-F., TANG, L. & LIN, P. 2023 The theory of fifth-order stokes waves in a linear shear current. *Proc. R. Soc. A* **479** (2280), 20230565.

736 GAO, T., MILEWSKI, P.A. & WANG, Z. 2021 Capillary-gravity solitary waves on water of finite depth interacting with a linear shear current. *Stud. Appl. Math.* **147** (3), 1–22.

737 GUO, D., TAO, B. & ZENG, X. 2014 On the dynamics of two-dimensional capillary-gravity solitary waves with a linear shear current. *Adv. Math. Phys.* **ID 480670**, 1–7.

738 HARRISON, W. J. 1909 The influence of viscosity and capillarity on waves of finite amplitude. *Proc. Lond. Math. Soc.* **2** (1), 107–121.

739 HOGAN, S. J. 1980 Some effects of surface tension on steep water waves. part 2. *J. Fluid Mech.* **96** (3), 417–445.

740 HOGAN, S. J. 1981 Some effects of surface tension on steep water waves. part 3. *J. Fluid Mech.* **110**, 381–410.

741 HSU, H. C., FRANCIUS, M., MONTALVO, P. & KHARIF, C. 2016 Gravity-capillary waves in finite depth on flows of constant vorticity. *Proc. R. Soc. Lond. A* **472**, 20160363.

742 HSU, H. C., KHARIF, C., ABID, M. & CHEN, Y. Y. 2018 A nonlinear schrödinger equation for gravity-capillary water waves on arbitrary depth with constant vorticity. part 1. *J. Fluid Mech.* **854**, 146–163.

743 HUR, V.M. 2019 Shallow water models with constant vorticity. *E. J. Mech.-B/Fluids* **73**, 170–179.

744 HUR, V.M. & JOHNSON, M.A. 2015 Modulational instability in the whitham surface with surface tension and vorticity. *Nonlinear Anal.* **129**, 104–118.

745 IVANOV, R. I. & MARTIN, C. I. 2019 On the time-evolution of resonant triads in rotational capillary-gravity water waves. *Phys. Fluids* **31** (11), 117103.

746 KANG, Y. & VANDEN-BROECK, J. M. 2000 Gravity-capillary waves in the presence of constant vorticity. *E. J. Mech.-B/Fluids* **19** (2), 253–268.

747 MARTIN, C.I. & MATIOC, B.-V. 2013 Existence of wilton ripples for water waves with constant vorticity and capillary effects. *SIAM J. Appl. Math.* **73** (4), 1582–1595.

748 MARTIN, C. I. 2013 Local bifurcation and regularity for steady periodic capillary-gravity water waves with constant vorticity. *Nonlinear Anal. Real World Appl.* **14**, 131–149.

749 NAYFEH, A. H. 1970a Finite amplitude surface waves in a liquid layer. *J. Fluid Mech.* **40(4)**, 671–684.

750 NAYFEH, A. H. 1970b Triple- and quintuple-dimped wave profiles in deep water. *Phys. Fluids* **13** (3), 545–550.

751 PIERSON, W. J. & FIFE, P. 1961 Some nonlinear properties of long-crested periodic waves with lengths near 2.44 centimeters. *J. Geophys. Res.* **66** (1), 163–179.

752 RIBEIRO, R., MILEWSKI, P. A. & NACHBIN, A. 2017 Flow structure beneath rotational water waves with stagnation points. *J. Fluid Mech.* **812**, 792–814.

753 SCHWARTZ, L. W. 1974 Computer extension and analytic continuation of stokes' expansion for gravity waves. *J. Fluid Mech.* **62** (3), 553–578.
777 SCHWARTZ, L. W. & VANDEN-BROECK, J. M. 1979 Numerical solution of the exact equations for capillary-gravity waves. *J. Fluid Mech.* **95** (1), 119–139.

778 SHELTON, J., MILEWSKI, P. A. & TRINH, P. H. 2021 On the structure of steady parasitic gravity-capillary waves in the small surface tension limit. *J. Fluid Mech.* **922**, A16.

781 WILTON, J. R. 1915 Lxxii. on ripples. *The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science* **29** (173), 688–700.

## 783 Appendix A. Approximate solution of the second-order

Equating the second power of $\epsilon$, we have the following sequence of equations

$$785 \qquad \Phi_{2\xi\xi} + \Phi_{2\eta\eta} = 0, \quad -\mu \le \eta < 0 \qquad (A1)$$

$$786 \qquad \Phi_{2\eta} + C_0 E_{2\xi} = -E_1 \Phi_{1\eta\eta} \\ 787 \qquad \qquad + (\Phi_{1\xi} + \Omega E_1 - C_1) E_{1\xi}, \quad \eta = 0 \qquad (A2)$$

$$788 \qquad C_0 \Phi_{2\xi} - (1 - C_0 \Omega) E_2 + \kappa E_{2\xi\xi} = \frac{1}{2}(\Phi_{1\xi}^2 + \Phi_{1\eta}^2) \\ 789 \qquad \qquad + \frac{1}{2}\Omega^2 E_1^2 + \Omega \Phi_{1\xi} E_1 - C_0 E_1 \Phi_{1\xi\eta}$$

$$790 \qquad - C_1 (\Phi_{1\xi} + \Omega E_1) - \frac{1}{2} K_2, \quad \eta = 0 \qquad (A3)$$

$$791 \qquad \Phi_{2\eta} = 0, \quad \eta = -\mu \qquad (A4)$$

Terms depending on $\eta$ appearing in equations for $\eta = 0$, as $E_1(\xi)\Phi_{1\eta\eta}$ and $E_1(\xi)\Phi_{1\xi\eta}$ are obtained by using Taylor series expansion given by

$$\Phi(\xi, \eta) \approx [\Phi(\xi, \eta)]_{\eta=E} + E(\xi)[\Phi(\xi, \eta)]_{\eta=E} + E^2(\xi)[\Phi(\xi, \eta)]_{\eta\eta}/2]_{\eta=E}$$

Using the second-order Laplace equation (A1) and the boundary condition (A4) we obtain the expressions for second-order potential and surface elevation in a similar way

$$795 \qquad \Phi_2(\xi, \eta) = A_2 \frac{\cosh[2(\eta + \mu)]}{\cosh(2\mu)} \sin(2\xi) \qquad (A5)$$

$$796 \qquad E_2(\xi) = a_2 \cos(2\xi) \qquad (A6)$$

Using (A5-A6) and the first-order solutions in the boundary conditions (A2-A3), we get the coefficients as follows

$$800 \qquad A_2 = \frac{C_0 \epsilon^2 (1 + \sigma^2)}{8(1 - 3\bar{\kappa})\sigma^4} [\bar{\Omega}^2 + 3(1 + \bar{\Omega}) + 3\sigma^2(\bar{\kappa}\bar{\Omega} + 2\bar{\kappa} - 1)], \qquad (A7)$$

$$801 \qquad a_2 = \frac{\epsilon^2}{4(1 - 3\bar{\kappa})\sigma^3} [\bar{\Omega}^2 + 3(1 + \bar{\Omega}) - \sigma^2(1 - \bar{\Omega})], \qquad (A8)$$

$$802 \qquad K_2 = \frac{C_0^2 \epsilon^2}{2\sigma^2} [(1 + \bar{\Omega})^2 - \sigma^2], \qquad (A9)$$

$$803 \qquad C_1 = 0, \qquad (A10)$$

Classically, the nonlinear correction of the phase velocity is unknown at this order and must be determined at the next order.

As found by Hsu et al. (2016), the expressions of $A_2$ and $a_2$ given by (A7) and (A8) are found to be singular when $\kappa$ equals $\kappa_1$, given implicitly as the root of the following equation

$$809 \qquad \kappa = \frac{\sigma^2}{(3 - \sigma^2) + 3\bar{\Omega}} \qquad (A11)$$
811 Without vorticity ($\overline{\Omega} = 0$), equation (A 11) reduces to the expressions obtained by Barakat & Houston (1968), Djordjevic & Redekopp (1977). In deep water ($\sigma = 1$) without vorticity, it further reduces to the well known value $\kappa_1 = 1/2$.

## 814 Appendix B. Approximate solution of the third-order

815 To derive the third-order approximate solution, we proceed in a similar manner and obtain
816 the following sequence of equations

$$817 \qquad \Phi_{3\xi\xi} + \Phi_{3\eta\eta} = 0, \quad -\mu \le \eta < 0 \qquad (B 1)$$

$$818 \qquad \Phi_{3\eta} + C_0 E_{3\xi} = -E_1 \Phi_{2\eta\eta} - E_2 \Phi_{1\eta\eta} - \frac{1}{2} E_1^2 \Phi_{1\eta\eta\eta} + (\Phi_{1\xi} + \Omega E_1) E_{2\xi} \\ 819 \qquad \qquad + (\Phi_{2\xi} + E_1 \Phi_{1\xi\eta} + \Omega E_2 - C_2) E_{1\xi}, \quad \eta = 0 \qquad (B 2)$$

$$820 \qquad C_0 \Phi_{3\xi} - (1 - C_0 \Omega) E_3 + \kappa E_{3\xi\xi} = -C_0 (E_1 \Phi_{2\xi\eta} + E_2 \Phi_{1\xi\eta} + \frac{1}{2} E_1^2 \Phi_{1\xi\eta\eta}) \\ 821 \qquad \qquad + (\Phi_{1\xi} + \Omega E_1)(\Phi_{2\xi} + E_1 \Phi_{1\xi\eta} + \Omega E_2 - C_2) + \Phi_{1\eta}(\Phi_{2\eta} + E_1 \Phi_{1\eta\eta}) \\ 822 \qquad \qquad + \frac{3}{2} \kappa E_{1\xi\xi} E_{1\xi}^2 - \frac{1}{2} K_3, \quad \eta = 0 \qquad (B 3)$$

$$823 \qquad \Phi_{3\eta} = 0, \quad \eta = -\mu \qquad (B 4)$$

825 Solving equations (B 1-B 4) we get the expressions for velocity potential and surface elevation amplitude as follows

$$827 \qquad \Phi_3(\xi, \eta) = A_{33} \frac{\cosh[3(\eta + \mu)]}{\cosh(3\mu)} \sin(3\xi) + A_{31} \frac{\cosh[(\eta + \mu)]}{\cosh(\mu)} \sin(\xi) \qquad (B 5)$$

$$828 \qquad E_3(\xi) = a_{33} \cos(3\xi) + a_{31} \cos(\xi) \qquad (B 6)$$

830 The derivation for obtaining the following results involves some algebra. On substituting (B 5-B 6) and using the first and second-order solutions into the boundary conditions (B 2-B 3),
we obtain

$$832 \qquad A_{33} = \frac{C_0 \epsilon^3}{64(1 - 3\bar{\kappa})[1 - \bar{\kappa}(3 + \sigma^2)]} \frac{(1 + 3\sigma^2)}{\sigma^7} \left[ \bar{\Omega}^4 + \left\{ 2(3 + \sigma^2) + 8\bar{\kappa}\sigma^2 \right\} \bar{\Omega}^3 + \left\{ (15 - \sigma^2) + \bar{\kappa}\sigma^2(43 + 19\sigma^2) \right\} \bar{\Omega}^2 + \left\{ 2(9 - 8\sigma^2 - \sigma^4) + 12\bar{\kappa}\sigma^2(7 + \sigma^2) + 24\bar{\kappa}^2\sigma^4(1 + \sigma^2) \right\} \bar{\Omega} + (9 - 22\sigma^2 + 13\sigma^4) + \bar{\kappa}\sigma^2(60 - 32\sigma^2 - 13\sigma^4) + 3\bar{\kappa}^2\sigma^4(16 + 5\sigma^2) \right], \quad (B 7)$$

$$836 \qquad A_{31} = \frac{C_2 \epsilon^3}{\sigma} + \frac{C_0}{\sigma} a_{31} - \frac{C_0 \epsilon^3}{8(1 - 3\bar{\kappa})} \frac{1}{\sigma^5} \left[ \bar{\Omega}^3 + (5 + 2\sigma^2)\bar{\Omega}^2 + \left\{ 3(3 + \sigma^2) + 3\bar{\kappa}\sigma^2(1 + \sigma^2) \right\} \bar{\Omega} + (6 - \sigma^2) + 3\bar{\kappa}\sigma^2(2 - \sigma^2) \right], \quad (B 8)$$

$$838 \qquad a_{33} = -a_{31} = \frac{\tanh(3\mu)}{C_0} A_{33} + \frac{\epsilon^3}{8(1-3\bar{\kappa})} \frac{1}{\sigma^4} \left[ \bar{\Omega}^3 + (5+2\sigma^2)\bar{\Omega}^2 + \left\{ 3(3+\sigma^2) + 3\bar{\kappa}\sigma^2(1+\sigma^2) \right\} \bar{\Omega} + (6-\sigma^2-2\sigma^4) + 3\bar{\kappa}\sigma^2(2+\sigma^2) \right], \quad (B 9)$$

$$840 \qquad K_3 = 0, \qquad (B 10)$$

$$841 \qquad C_2 = \frac{C_0}{8(\bar{\Omega} + 2)(1 - 3\bar{\kappa})} \frac{1}{\sigma^4} \left[ \bar{\Omega}^4 + 2(3 + \sigma^2)\bar{\Omega}^3 + \left\{ 3(5 + \sigma^2) + 3\bar{\kappa}\sigma^2(1 + \sigma^2) \right\} \bar{\Omega}^2 + \left\{ 2(9 - 2\sigma^2 + \sigma^4) + 12\bar{\kappa}\sigma^2(1 - \sigma^2) \right\} \bar{\Omega} + (9 - 10\sigma^2 + 9\sigma^4) + 3\bar{\kappa}\sigma^2(4 - 8\sigma^2 - \sigma^4) + 9\bar{\kappa}^2\sigma^6 \right] \quad (B 11)$$

At third-order approximation the first nonvanishing nonlinear correction of the phase velocity $C_2$ is given by equation (B 11) that matches equation (3.32) of Hsu et al. (2016), who first obtained this analytical result. The expressions for the third-order terms in equations (B 7) and (B 9) are found to be singular when $\kappa$ equals $\kappa_2$, given implicitly as the root of the following equation

$$850 \qquad \kappa = \frac{\sigma^2}{3(1 + \bar{\Omega}) + \sigma^2 \bar{\Omega}} \qquad (B 12)$$

Taking $\bar{\Omega} = 0$ in equation (B 12) yields $\kappa_2 = \sigma^2/3$, in agreement with the results of Barakat & Houston (1968). In deep water $\sigma = 1$ and we obtain the other well-known value $\kappa_2 = 1/3$.