Time-dependent nonlinear gravity–capillary surface waves with viscous dissipation and wind forcing

Abstract

We develop a time-dependent conformal method to study the effect of viscosity on steep surface waves. When the effect of surface tension is included, numerical solutions are found that contain highly oscillatory parasitic capillary ripples. These small-amplitude ripples are associated with the high curvature at the crest of the underlying viscous-gravity wave, and display asymmetry about the wave crest. Previous inviscid studies of steep surface waves have calculated intricate bifurcation structures that appear for small surface tension. We show numerically that viscosity suppresses these. While the discrete solution branches still appear, they collapse to form a single smooth branch in the limit of small surface tension. These solutions are shown to be temporally stable, both to small superharmonic perturbations in a linear stability analysis, and to some larger amplitude perturbations in different initial-value problems. Our work provides a convenient method for the numerical computation and analysis of water waves with viscosity, without evaluating the free-boundary problem for the full Navier–Stokes equations, which becomes increasingly challenging at larger Reynolds numbers.

1. Introduction

When surface tension is included in the classical formulation of a travelling nonlinear gravity wave, solutions are seen to exhibit highly oscillatory parasitic capillary ripples. These are small-scale ripples that travel with the same speed as the underlying gravity wave. In most previous numerical and asymptotic studies (Schwartz & Vanden-Broeck 1979; Shelton, Milewski & Trinh 2021; Shelton & Trinh 2022), this phenomenon has been studied with the use of potential flow theory, in which the effect of viscosity is neglected. This results in symmetric solutions, in which the capillary ripples are observed across the entire surface wave profile. This is in contrast to experimental observations (Ebuchi, Kawamura & Toba 1987), in which the parasitic ripples are focused on the forward face of the travelling wave. This disparity between the surface profiles viewed experimentally, and those from inviscid theory, is usually attributed to the effect of fluid viscosity (Dosaev, Troitskaya & Shrira 2021).

It is of significant interest to develop simplified models for gravity–capillary waves that include viscosity, but without resorting to solving the full Navier–Stokes equations. To this aim, the principal question is whether it is possible to combine the simplicity of the classical potential flow framework for an inviscid water wave, with effective conditions that govern and model viscous effects near key regions. Note that the viscous terms in the Navier–Stokes equations are identically zero for a velocity derived from a potential. In the limit of small viscosity, analytical progress is tractable through the study of a viscous boundary layer near to the fluid surface (Longuet-Higgins 1992). Previous authors, such as Ruvinsky, Feldstein & Freidman (1991) and Fedorov & Melville (1998), have applied these techniques to investigate the effect of viscosity on steep gravity–capillary waves. Both of these investigations produced solutions with asymmetric parasitic capillary ripples, which closely resembled previous experimental observations. In the work by Ruvinsky et al. (1991), kinematic and dynamic boundary conditions were presented, where viscosity feeds into the kinematic condition through a non-local condition involving time integration of the velocity potential – a step that was subsequently simplified into a local condition by Dias, Dyachenko & Zakharov (2008). Approximations valid in the limit of small amplitude were then made to determine explicit equations for the Fourier coefficients of the solution, which also yielded the damping rate as a function of the small-amplitude parameter. A steady formulation of this problem has also been developed by Fedorov & Melville (1998) using viscous boundary layer approximations. With the addition of a surface pressure forcing (modelling the effect of wind), steadily travelling solutions were calculated numerically. This model was subsequently used by Melville & Fedorov (2015) to demonstrate that in the ocean, parasitic capillary ripples, rather than the main gravity wave, can be responsible for the majority of the viscous damping required to offset the growth of the overall wave due to the effect of wind. We note that the earliest work investigating the effect of viscosity on steep waves with parasitic capillary ripples was by Longuet-Higgins (1963); however, these results compared poorly with the experimental observations by Perlin, Lin & Ting (1993). This theory would later be improved upon and updated by Longuet-Higgins (1995). While they proposed many ground-breaking ideas, the above works are often challenging to interpret mathematically on account of the number of approximations made.

In this work, we will study the formulation proposed by Dias et al. (2008), which incorporates viscosity in the surface boundary conditions of a potential flow boundary-value problem. In this model, the free-surface kinematic and dynamic boundary conditions are given by

(1.1)\begin{equation} \left. \begin{gathered} \phi_t +\frac{1}{2}\,(\phi_x^2+\phi_y^2) + \frac{\zeta}{F^2}-\frac{B}{F^2}\,\kappa +\frac{P}{F^2}\,\zeta_{x}+ \frac{2}{Re}\,\phi_{yy}=0,\\ \zeta_t = \phi_y - \phi_x \zeta_x + \frac{2}{Re}\,\zeta_{xx}, \end{gathered} \right\} \end{equation}

evaluated at the unknown free surface $y=\zeta (x,t)$. In boundary conditions (1.1), $\phi (x, t)$ is the velocity potential that satisfies Laplace's equation within the domain $-1/2 \leq x \leq 1/2$, $-\infty < y \leq \zeta (x,t)$. The formulation includes the effect of surface tension through the Bond number $B$ and curvature $\kappa =\zeta _{xx}/(1+\zeta _x^2)^{3/2}$, the effect of fluid viscosity through the Reynolds number $Re$, and surface wind forcing that depends on the wave slope $\zeta _{x}$. The wind forcing is required in order to obtain steadily travelling solutions in the presence of viscous dissipation. Full details of the mathematical formulation are given later, in § 2 and particularly in (2.1a)–(2.1d). Finding solutions to this formulation is difficult on account of the moving domain $-\infty \leq y \leq \zeta (x,t)$ and the dependence of the solution on three independent variables, $x$, $y$ and $t$.

Our work focuses on the development of a time-dependent conformal method that allows for the convenient numerical evaluation of a two-dimensional potential flow problem with viscosity. The main benefit of this conformal method is that the original two-dimensional domain, which is bounded by a moving surface wave, is reduced to a fixed one-dimensional domain for the time-dependent free surface. We apply the conformal mapping techniques of Dyachenko, Zakharov & Kuznetsov (1996) to system (1.1) to derive time-dependent equations that govern each variable on the free surface. These will depend on only one spatial coordinate, the conformal variable $\xi$, which parametrises the free surface. This conformal method preserves all features of the original boundary-value problem, and no further assumptions are made following the introduction of boundary conditions (1.1).

In previous investigations of inviscid surface waves, intricate bifurcation structures have emerged in the small surface tension limit (Champneys, Vanden-Broeck & Lord 2002; Shelton et al. 2021; Shelton, Milewski & Trinh 2023), an example of which is shown for steadily travelling gravity–capillary waves in figure 1(a). These consist of a countably infinite number of solution branches in the $(B,F)$ plane that manifest for fixed amplitude under the limit of small surface tension. Across each of the branches in figure 1(a), the wavenumber of the parasitic capillary ripples increases by one. Figure 1(b) shows an example of a travelling surface wave obtained in this inviscid framework, in which the capillary ripples are symmetric about the wave crest. The amplitude of these parasitic capillary ripples was measured by Shelton et al. (2021) to be exponentially small as $B \to 0$. Asymptotic solutions for these were later obtained by Shelton & Trinh (2022) using beyond-all-order asymptotics, in which the capillary ripples were produced by the Stokes phenomenon across Stokes lines associated with the high crest curvature of the leading-order nonlinear gravity wave. This also produced an asymptotic solvability condition for when perturbation solutions do and do not exist – values of non-existence are shown with black circles in figure 1(a). The exponential scaling of the capillary ripple amplitude is shown in the semi-log plot of figure 1(c).

Figure 1. Numerical results for steadily propagating inviscid gravity–capillary waves calculated by Shelton et al. (2021). (a) The location of solutions is shown in the $(B,F)$ plane, where the Bond number $B$ and Froude number $F$ are non-dimensional parameters defined in (2.2a–d). Black circles show the parameter values given from the failure of the nonlinear solvability condition derived by Shelton & Trinh (2022). (b) A typical solution containing oscillatory capillary ripples, which has $B=0.001648$ and $F=0.4245$. (c) Semi-log plot showing the exponentially-small amplitude of the capillary ripples for each solution branch in (a).

We use the nonlinear formulation (1.1) to investigate the effect of small capillarity on viscous gravity–capillary waves. First, we study steadily travelling solutions, for which surface wind forcing is introduced in the kinematic boundary condition to counteract the energy decay induced by viscous dissipation. We demonstrate numerically that in the presence of viscosity, the discrete branching structure obtained in the inviscid regime (figure 1a) does not persist below a certain value of the surface tension. While discrete branches are observed for larger values of the surface tension, they close off such that a single smooth branch exists in the limit of zero surface tension. This is observed numerically both when the viscosity is fixed, and when distinguished limits are chosen in which the viscosity decays alongside the surface tension in an algebraic manner.

Second, we study the temporal stability of our steady parasitic solutions. These are shown to be superharmonically stable through a linear stability analysis, in which an eigenvalue problem for small perturbations is studied and these are shown to decay in time. A nonlinear stability analysis is then considered with an initial-value problem that implements our time-dependent formulation. This allows us to comment on the global ‘attractiveness’ of the steadily travelling solutions, and convergence is observed when starting from both a large amplitude initial condition, and a small amplitude initial condition. In this latter case, the surface wind forcing initially dominates which causes the wave amplitude to increase, before eventually balancing with viscous dissipation as the target solution is converged upon.

1.1. Outline of our paper

We begin in § 2 by formulating the two-dimensional boundary-value problem that models nonlinear surface waves that travel upon a viscous fluid. Our conformal mapping of the unsteady problem, for which full details are presented in Appendix A, yields a one-dimensional formulation for the free surface. This is given in § 2.2 for unsteady solutions, and § 2.3 for steadily travelling solutions. Numerical solutions to the steady formulation are obtained in § 3, with an emphasis on detecting the bifurcation structure that emerges for small surface tension and viscosity. The unsteady formulation is implemented numerically in § 4, where we perform linear and nonlinear stability analysis to observe the temporal stability of the steady solutions when starting from different initial conditions.

2. Mathematical formulation

We consider the time evolution of a nonlinear surface wave subject to small viscous effects within a two-dimensional fluid extending to infinite depth. The travelling wave is assumed to be periodic in the direction of propagation, and is subject to the effects of both gravity and surface tension. The boundary-value problem describing this is specified on the fluid domain $-1/2 \leq x \leq 1/2$ and $-\infty < y \leq \zeta (x,t)$, with the velocity potential $\phi (x,y,t)$ and the free surface $y=\zeta (x,t)$ as solutions. After non-dimensionalisation, and moving into a co-moving frame of unit speed through a Galilean transformation, the boundary-value problem is given by

(2.1a)$$\begin{gather} \phi_t-\phi_{x} +\frac{1}{2}\,(\phi_x^2+\phi_y^2) + \frac{\zeta}{F^2}-\frac{B}{F^2}\,\kappa +\frac{P}{F^2}\,\zeta_x+ \frac{2}{Re}\,\phi_{yy}=0 \quad \text{at} \ y=\zeta(x,t), \end{gather}$$

(2.1b)$$\begin{gather}\zeta_t = \phi_y+\zeta_x(1-\phi_x) + \frac{2}{Re}\,\zeta_{xx}\quad \text{at}\ y=\zeta(x,t), \end{gather}$$

(2.1c)$$\begin{gather}\phi_{xx} +\phi_{yy}=0 \quad \text{for}\ y \leq \zeta(x,t), \end{gather}$$

(2.1d)$$\begin{gather}\phi_{x} \to 0, \quad \phi_{y} \to 0 \quad \text{as}\ y \to -\infty. \end{gather}$$

This system comprises of kinematic and dynamic boundary conditions (2.1a) and (2.1b) on the surface, Laplace's equation (2.1c) within the fluid, and the decay of motion in the deep-water limit (2.1d). The effect of viscosity appears in the boundary conditions (2.1a) and (2.1b); these model effects were proposed by Dias et al. (2008). Note also the inclusion of the surface wind forcing $P \zeta _{x}/F^2$ in (2.1a). We will refer to $P$ as the non-dimensional wind strength. Note that in the inviscid formulation with $1/Re=0$ and $P=0$, the wave energy, which measures the effects of kinetic, capillary and gravitational potential energies, is a time-conserved quantity of the unsteady formulation. However, when $1/Re \neq 0$ and $P=0$, dissipation renders the energy a decreasing function of time for solutions to formulation (2.1), which we demonstrate in Appendix B. Thus steadily travelling waves will not exist in this viscous formulation in the absence of wind forcing. This is the reason why we include the effects of both viscous dissipation and wind forcing in this work, in order to study and classify steadily travelling solutions, as well as to analyse their stability. This model for wind forcing is discussed in more detail in § 2.1.

Four non-dimensional constants appear in system (2.1). These are the Froude number $F$, the Bond number $B$, the Reynolds number $Re$, and the non-dimensional wind strength $P$, defined by

(2.2a–d)\begin{equation} F= \frac{c}{\sqrt{g \lambda}}, \quad B=\frac{\sigma}{\rho g \lambda^2}, \quad Re=\frac{c \lambda}{\nu}, \quad P=\frac{p}{g \lambda}. \end{equation}

Here, $c$ is the speed of the travelling frame (which for steady solutions is the wave speed), $g$ is the constant acceleration due to gravity, $\lambda$ is the wavelength, $\sigma$ is the constant coefficient of surface tension, $\rho$ is the fluid density, $\nu$ is the kinematic viscosity of the fluid, and $p$ is the amplitude of the physical wind forcing. Note that we have non-dimensionalised length scales and the free-surface elevation with respect to $\lambda$, the velocity potential with respect to $c \lambda$, and time with respect to $\lambda /c$.

Crucially, we note that when considering the full effects of both nonlinearity and viscosity, there is no known explicit formulation of the kinematic and dynamic boundary conditions that can be written in a form analogous to (2.1a) and (2.1b) for the velocity potential and surface height. In such cases, it is necessary to solve the full Navier–Stokes equations within the unknown fluid domain; this has been performed numerically by e.g. Hung & Tsai (2009) for the case of time-dependent nonlinear viscous gravity–capillary waves. However, explicit kinematic and dynamic boundary conditions for the velocity potential $\phi$ and free surface $\zeta$ can be derived in two cases.

1. (i) Inviscid flows, for which the kinematic and dynamic boundary conditions for gravity–capillary waves emerge, the former from the classical Bernoulli equation (cf. Vanden-Broeck 2010). These are the same as (2.1) but with the viscous terms, which contain the Reynolds number $Re$, removed.

2. (ii) Linear flows with non-zero viscosity: here, the kinematic and dynamic boundary conditions may be derived from the linearised two-dimensional Navier–Stokes equations. This is performed by decomposing the velocity field into irrotational and solenoidal components, for which the linearised forms of (2.1a) and (2.1b) were derived by Dias et al. (2008). These equations have subsequently been used to study viscous effects on a variety of free-surface formulations, such as for Faraday pilot waves in bouncing fluid droplets (Milewski et al. 2015; Blanchette 2016), three-dimensional solitary waves with forcing (Wang & Milewski 2012), and deriving dissipation rates for ocean swell (Henderson & Segur 2013).

Thus the equations that we use throughout this work are obtained by combining the viscous term found from (ii) above with the nonlinear inviscid boundary conditions from (i), yielding (2.1a) and (2.1b). We note that the ‘correct’ nonlinear generalisation of the term $\phi _{yy}$ in kinematic condition (2.1a) is $\phi _{nn}$, where $n$ is the normal vector to the free surface, which was originally derived by Ruvinsky et al. (1991). There is no known correct nonlinear generalisation of the term $\zeta _{xx}$ in (2.1b). These model equations were proposed by Dias et al. (2008), and one aim of our current work is to develop a time-dependent conformal map that efficiently solves the nonlinear version of this problem. This methodology is then used to study the effect that viscosity plays in the parasitic capillary ripples present on steep gravity waves, their associated bifurcation structure, and temporal stability. We are particularly interested in the small surface tension limit on account of the intricate bifurcation structures that exist in this singular regime. We note that these nonlinear equations have previously been applied to study dissipative effects on solitons by Brunetti et al. (2014) and Liao et al. (2023), both of whom derived forced nonlinear Schrodinger equations asymptotically.

2.1. Choice of wind forcing

To obtain steadily travelling solutions, we balance viscous dissipation with the addition of a wind forcing term in Bernoulli's equation (2.4b). The surface wind forcing that we consider is that developed by Jeffreys (1925) as a model for linear wave growth induced by wind forcing. This results in the addition of the term $p\hat {\zeta }_{\hat {x}}$ into the dimensional kinematic condition, which after non-dimensionalisation becomes $P \zeta _{x}/F^2$ in (2.1b). Here, $p$ is the dimensional wind strength, and $P=p/(g \lambda )$ is the non-dimensional wind strength. This is known as a sheltering model; for a wave travelling from left to right with $P>0$, it adds energy to the rear face of the wave, and removes energy from the forward face. In our numerical search for steadily travelling solutions in § 3, the constant $P$ is treated as an unknown of the formulation, and a unique value is obtained by solving an underdetermined system of equations.

In the study by Fedorov & Melville (1998), who included forcing in the kinematic condition, the term $P\cos (2 {\rm \pi}x)$ was used to mimic the effect of wind. Fedorov and Melville found that when all other parameters were fixed, a range of $P$ values was permitted, corresponding to different phase shifts between their fixed wind forcing and unknown solution profile. This complication occurred as a result of their choice of wind forcing breaking the translational invariance of the system. In contrast, the Jeffreys model for wind forcing that we use in this work both retains the translational invariance of the system and is a justifiable model for the shear stress induced by wind forcing. Since a unique value of $P$ is selected as an unknown to the steady solutions, the resultant bifurcation space is simpler than that investigated by Fedorov & Melville (1998), which allows for an easier study of solution branches.

We note that more accurate models exist to capture the influence of wind on surface water waves. The monograph by Janssen (2004) details elements incorporated in these, such as the transfer of energy between coupled air and water layers, and the consideration of turbulent effects in the air. For instance, the formulation from which Miles (1957) developed his theory of linear wave generation involved inviscid and incompressible air and water layers coupled by an interfacial stress condition. The water layer was irrotational, and the consideration of rotational effects in the air layer gave rise to a critical layer problem governed by the inviscid Orr–Sommerfeld equation for linear disturbances to the wave surface.

2.2. The time-dependent conformal mapping

Evolving the solutions $\zeta (x,t)$ and $\phi (x,y,t)$ to the boundary-value problem (2.1) in time is difficult. This is due to both the physical $(x,y)$ domain being two-dimensional, and the boundary conditions being imposed upon the free surface, which itself can change in time. In this subsection, we present a time-dependent mapping of system (2.1), under which the flow domain $-1/2 \leq x \leq 1/2$ and $-\infty < y \leq \zeta (x,t)$ is mapped to the lower-half $(\xi,\eta )$ plane. Upon evaluating the solutions at the fluid surface $\eta =0$, a one-dimensional surface formulation emerges, parametrised by $\xi$ and $t$. This conformal method was originally developed by Dyachenko et al. (1996) and Choi & Camassa (1999) for nonlinear gravity capillary waves.

In writing $x=x(\xi,\eta,t)$ and $y=y(\xi,\eta,t)$, the surface solutions $X$, $Y$, $\varPhi$ and $\varPsi$ are defined by

(2.3)\begin{equation} \left. \begin{array}{@{}ll@{}} X(\xi,t)=x(\xi,0,t), & \varPhi(\xi,t)=\phi(x(\xi,0,t),y(\xi,0,t),t), \\[2pt] Y(\xi,t)=\zeta(x(\xi,0,t),t), & \varPsi(\xi,t) = \psi(x(\xi,0,t),y(\xi,0,t),t). \end{array} \right\} \end{equation}

Here, the streamfunction $\psi$ is the harmonic conjugate of the velocity potential $\phi$. We then derive time-evolution equations for the surface solutions (2.3), for which full details are presented in Appendix A. We obtain evolution equations for the surface height $Y$,

(2.4a)\begin{align} Y_t=\frac{X_{\xi}(Y_{\xi}-\varPsi_{\xi})}{J}+\frac{2}{Re}\,\frac{X_{\xi} Y_{\xi \xi} - Y_{\xi} X_{\xi \xi}}{X_{\xi}J}-Y_{\xi}\,\mathscr{H} \left[\frac{Y_{\xi}-\varPsi_{\xi}}{J} +\frac{2}{Re}\,\frac{X_{\xi}Y_{\xi \xi}-Y_{\xi}X_{\xi\xi}}{X_{\xi}^2 J}\right], \end{align}

and the surface velocity potential $\varPhi$,

(2.4b) \begin{align} \varPhi_t&=\frac{\varPsi_{\xi}^2-\varPhi_{\xi}^2}{2J} - \frac{Y}{F^2}-\frac{P}{F^2}\,\frac{Y_{\xi}}{X_{\xi}} + \frac{B \kappa}{F^2}+\frac{X_{\xi}\varPhi_{\xi}}{J}\nonumber\\ &\quad -\varPhi_{\xi}\,\mathscr{H} \left[\frac{Y_{\xi}-\varPsi_{\xi}}{J} +\frac{2}{Re}\,\frac{X_{\xi}Y_{\xi \xi}-Y_{\xi}X_{\xi\xi}}{X_{\xi}^2 J}\right]\nonumber\\ &\quad +\frac{2}{Re}\left[\frac{Y_{\xi} X_{\xi \xi} - X_{\xi} Y_{\xi \xi}}{X^2_{\xi}J}\,\varPsi_{\xi} +\frac{Y_{\xi}X_{\xi \xi}(Y_{\xi}^2-3X_{\xi}^2)+X_{\xi}Y_{\xi\xi}(X_{\xi}^2-3Y_{\xi}^2)}{J^3}\,\varPsi_{\xi}\right.\nonumber\\ &\left.\quad {}+\frac{X_{\xi}X_{\xi \xi}(3Y_{\xi}^2-X_{\xi}^2)+Y_{\xi}Y_{\xi \xi}(Y_{\xi}^2-3X_{\xi}^2)}{J^3}\,\varPhi_{\xi} +\frac{X_{\xi}^2-Y_{\xi}^2}{J^2}\,\varPhi_{\xi \xi}+\frac{2X_{\xi}Y_{\xi}}{J^2}\,\varPsi_{\xi \xi}\right], \end{align}

where $X_{\xi }$ and $\varPsi$ are known from the harmonic relations

(2.4c,d)\begin{equation} X_{\xi}=1-\mathscr{H}[Y_{\xi}] \quad \text{and} \quad \varPsi = \mathscr{H}[\varPhi]. \end{equation}

In system (2.4), $J=X_{\xi }^2+Y_{\xi }^2$ is the Jacobian of the mapping, $\kappa =(X_{\xi }Y_{\xi \xi }-Y_{\xi }X_{\xi \xi })/J^{3/2}$ is the surface curvature, and $\mathscr {H}[Y]=\int _{-1/2}^{1/2}Y(\xi ') \cot {[{\rm \pi} (\xi '-\xi )]}\,\mathrm {d}\xi '$ is the periodic Hilbert transform. Given an initial condition for $Y(\xi,0)$, $\varPhi (\xi,0)$, and values of $B$, $F$, $Re$ and $P$, (2.4) may be used to evolve the solutions in time. In § 4, we use this scheme to study the stability of steadily travelling solutions with surface wind forcing, when viewed as convergence in a time-evolution problem.

2.3. Conformal mapping for steadily travelling waves

We now present integro-differential equations, depending only on the conformal domain $\xi$, that are satisfied by steady solutions of the mapped system (2.4). These are given by

(2.5a)\begin{align} &\frac{\varPhi_{\xi}^2+\varPsi_{\xi}^2}{2J}-\frac{X_{\xi}\varPhi_{\xi}+Y_{\xi} \varPsi_{\xi}}{J} + \frac{Y}{F^2} +\frac{P}{F^2}\,\frac{Y_{\xi}}{X_{\xi}} -\frac{B}{F^2}\,\frac{X_{\xi} Y_{\xi \xi} - Y_{\xi} X_{\xi \xi}}{J^{3/2}}\nonumber\\ &\qquad +\frac{2}{Re}\left[\frac{(Y_{\xi}^2-X_{\xi}^2)\varPhi_{\xi \xi}-2 X_{\xi}Y_{\xi}\varPsi_{\xi \xi}}{J^2} +\frac{\varPhi_{\xi}}{J^3}\left( X_{\xi \xi} X_{\xi}(X_{\xi}^2-3Y_{\xi}^2)+Y_{\xi \xi}Y_{\xi}(3X_{\xi}^2-Y_{\xi}^2)\right)\right.\nonumber\\ &\left.\qquad {}+\frac{\varPsi_{\xi}}{J^3}\left(X_{\xi \xi} Y_{\xi}(3X_{\xi}^2-Y_{\xi}^2)+Y_{\xi \xi} X_{\xi}(3Y_{\xi}^2-X_{\xi}^2)\right)\vphantom{\frac{(Y_{\xi}^2-X_{\xi}^2)\varPhi_{\xi \xi}-2 X_{\xi}Y_{\xi}\varPsi_{\xi \xi}}{J^2}}\right] =0, \end{align}

(2.5b)\begin{align} &\qquad\qquad\qquad\qquad\quad\qquad \varPsi_{\xi} = Y_{\xi} +\frac{2}{Re}\,\frac{X_{\xi} Y_{\xi \xi} - Y_{\xi} X_{\xi \xi}}{X_{\xi}^2}, \end{align}

together with harmonic relations

(2.5c,d)\begin{equation} X_{\xi}=1-\mathscr{H}[Y_{\xi}] \quad \text{and} \quad \varPhi_{\xi}={-}\mathscr{H}[\varPsi_{\xi}], \end{equation}

where $\mathscr {H}$ is the periodic Hilbert transform. This yields four equations for the four unknown functions $X$, $Y$, $\varPhi$ and $\varPsi$.

There are two ways to derive system (2.5). First, one may consider steady solutions of the original boundary-value problem (2.1), for which a steady conformal mapping analogous to that presented in Appendix A yields (2.5a) and (2.5b). Alternatively, one may consider steady solutions of the time-dependent conformal system (2.4), which after simplification also yields (2.5a) and (2.5b).

When solving (2.5) numerically, we will typically fix $B$, $Re$ and the wave energy $\mathscr {E}$, defined by

(2.6a) \begin{equation} \mathscr{E}=\frac{1}{E_{hw}}\int_{{-}1/2}^{1/2} \begin{bmatrix}\displaystyle \underbrace{\frac{F^2}{2}\,\varPsi_{\xi} \varPhi}_{\textit{kinetic}}+\underbrace{\vphantom{\frac{F^2}{2}}B(\sqrt{J}-X_{\xi})}_{\textit{capillary}}+\underbrace{\frac{1}{2}\,Y^2X_{\xi}}_{\textit{gravitational}}\end{bmatrix} \mathrm{d}\xi, \end{equation}

for which $F$ and $P$ are obtained as unknowns of the problem. This is performed in § 3 to calculate nonlinear solutions for small values of the surface tension. The numerical method used to solve system (2.5) is given in § 3.1. In (2.6a), we have normalised with respect to the energy of the highest inviscid Stokes wave, $E_{hw}=0.00184$, computed to three significant digits. The steady solutions calculated in this paper will have $\mathscr {E}=0.4$, which is chosen to compare with previous inviscid works (Shelton et al. 2021; Shelton & Trinh 2022) studying the classes of parasitic solutions that emerge in the small surface tension limit. Note that the energy is not a conserved quantity in this formulation, due to the presence of wind forcing and dissipation. An expression for the rate of change of energy in the absence of the surface wind forcing is derived in Appendix B, and this is shown to be inversely proportional to $Re$. Note also that the steady solutions provide information only about $\varPhi _{\xi }$, whereas evaluation of (2.6a) requires $\varPhi$, which we obtain by integration. The constant of integration for $\varPhi$ is determined explicitly from the condition

(2.6b)\begin{equation} \int_{{-}1/2}^{1/2} \varPhi X_{\xi} \,\mathrm{d}\xi =0. \end{equation}

Condition (2.6b) is equivalent to writing $\int _{-1/2}^{1/2}\phi (x,\zeta (x)) \,\mathrm {d}\kern0.06em x=0$ in conformal variables.

2.4. Linear solutions for the surface and the internal vorticity field

In the Dias et al. (2008) model, a Helmholtz decomposition was used to write the internal velocity field as $u(x,y,t)=\phi _{x}-A_{y}$ and $v(x,y,t)=\phi _{y}+A_{x}$, where $\phi$ is the velocity potential, and the function $A$ contributes to the vorticity field $\omega =v_{x}-u_{y}=A_{xx}+A_{yy}$. For our formulation of travelling solutions that are steady in a co-moving frame, $A(x,y)$ satisfies the boundary-value problem

(2.7a)$$\begin{gather} -\frac{\partial A}{\partial x} = \frac{1}{Re}\left( \frac{\partial^2 A}{\partial x^2}+\frac{\partial^2 A}{\partial y^2}\right) \quad \text{for}\ y\leq \zeta(x), \end{gather}$$

(2.7b)$$\begin{gather}\frac{\partial A}{\partial x}=\frac{2}{Re}\,\frac{\partial^2 \zeta}{\partial x^2} \quad \text{at}\ y=\zeta(x), \end{gather}$$

(2.7c)$$\begin{gather}A_{x} \to 0,A_{y} \to 0 \quad \text{as}\ y \to -\infty, \end{gather}$$

where $y=\zeta (x)$ is assumed known. In conformal variables, system (2.7) becomes

(2.8a)$$\begin{gather} y_{\xi}A_{\eta}-y_{\eta}A_{\xi}=\frac{1}{Re}\,(A_{\xi \xi}+A_{\eta \eta}) \quad \text{for}\ \eta \leq 0, \end{gather}$$

(2.8b)$$\begin{gather}X_{\xi}A_{\xi}-Y_{\xi}A_{\eta}=\frac{2}{Re}\,\frac{(X_{\xi}Y_{\xi \xi}-Y_{\xi}X_{\xi \xi})(X_{\xi}^2+Y_{\xi}^2)}{X_{\xi}^3} \quad \text{at}\ \eta=0, \end{gather}$$

(2.8c)$$\begin{gather}A_{\xi} \to 0,A_{\eta} \to 0 \quad \text{as}\ \eta \to -\infty, \end{gather}$$

and the vorticity field is given by $\omega =(A_{\xi \xi }+A_{\eta \eta })/(x_{\xi }^2+y_{\xi }^2)$. The functions $X(\xi )$ and $Y(\xi )$ in boundary condition (2.8b) parametrise the free surface $y=\zeta (x)$, and are assumed known from solving system (2.5). The functions $x(\xi,\eta )$ and $y(\xi,\eta )$ in partial differential equation (2.8a) and the expression for $\omega$ are the analytic continuations of these within the flow domain via Laplace's equation, and may be calculated by the Poisson integral formula (or in Fourier space with the multiplier $\exp ({2\, \lvert k \rvert \,{\rm \pi} \eta })$). As an example, if we know the coefficients of the Fourier series $Y(\xi )=\sum _{k=-\infty }^{\infty } c_k\,\mathrm {e}^{2 k {\rm \pi}\mathrm {i}\xi }$, then we may calculate $y(\xi,\eta )=\eta +\sum _{k=-\infty }^{\infty } c_k\,\mathrm {e}^{2\,\lvert k \rvert \,{\rm \pi} \eta }\, \mathrm {e}^{2 k {\rm \pi}\mathrm {i}\xi }$.

While both systems (2.7) and (2.8) are linear, they are difficult to solve – (2.7) on account of evaluation of boundary condition (2.7b) at the free surface $y=\zeta$, and (2.8) due to the coefficients being known nonlinear functions. We now proceed to solve for the vorticity field analytically, first by constructing linear solutions for the free surface in § 2.4.1, and then for the function $A$ and the vorticity field in § 2.4.2.

2.4.1. Linear theory for the free surface

We now analytically study small-amplitude solutions for the unknown free surface. These will be required for our small-amplitude study of the vorticity field, and will also provide a possible explanation for the behaviour of solution branches that we later observe numerically in the nonlinear regime in § 3.

We consider a Stokes expansion in powers of a small-amplitude parameter $\epsilon$. In substituting for solutions of the form $X \sim \xi + \epsilon X_1$, $Y \sim \epsilon Y_1$, $\varPsi \sim \epsilon \varPsi _1$ and $\varPhi \sim \epsilon \varPhi _1$ into (2.5a) and (2.5b), at $O(\epsilon )$, two equations are found for the first-order perturbation solutions. We eliminate $\varPhi _1$ from these two equations by using the harmonic relation $\varPhi _1=-\mathscr {H}[\varPsi _1]$, and then eliminate $\varPsi _1$ to find the single equation

(2.9)\begin{equation} X_{1 \xi} - \frac{Y_1}{F^2}-\frac{P}{F^2}\,Y_{1 \xi} +\frac{B}{F^2}\,Y_{1 \xi \xi} +\frac{4}{Re}\,X_{1 \xi \xi} + \frac{4}{Re^2}\,X_{1 \xi \xi \xi}=0. \end{equation}

In writing $Y(\xi )$ as a Fourier series of the form

(2.10)\begin{equation} Y(\xi)=a_0+\sum_{k=1}^{\infty} [a_k\cos(2k {\rm \pi}\xi) + b_k \sin(2k {\rm \pi}\xi)], \end{equation}

we use the harmonic relation $X_1=-\mathscr {H}[Y_1]$, and equate each coefficient of $\cos (2 k {\rm \pi}\xi )$ and $\sin (2 k {\rm \pi}\xi )$ to zero, to find

(2.11)\begin{equation} \left. \begin{gathered} \left(2k {\rm \pi}-\frac{1}{F^2} -\frac{4 k^2 {\rm \pi}^2 B}{F^2}-\frac{32 k^3{\rm \pi}^3}{Re^2}\right)a_k + \left(\frac{16 k^2 {\rm \pi}^2}{Re}-\frac{2 k {\rm \pi}P }{F^2}\right)b_k=0,\\ \left(2k {\rm \pi}-\frac{1}{F^2} -\frac{4 k^2 {\rm \pi}^2 B}{F^2}-\frac{32 k^3{\rm \pi}^3}{Re^2}\right)b_k - \left(\frac{16 k^2 {\rm \pi}^2}{Re}-\frac{2 k {\rm \pi}P}{F^2} \right)a_k=0. \end{gathered} \right\} \end{equation}

For non-degenerate solutions, there must exist a value of $k$ for which at least one of $a_k$ or $b_k$ is non-zero. This yields

(2.12a,b)\begin{equation} 2k {\rm \pi}-\frac{1}{F^2} -\frac{4 k^2 {\rm \pi}^2 B}{F^2}-\frac{32 k^3{\rm \pi}^3}{Re^2}=0 \quad \text{and} \quad P=\frac{8 k {\rm \pi}F^2}{Re}. \end{equation}

In the inviscid regime, for which $1/Re$ and $P$ are both zero, it is possible to find values of $B$ and $F$ for which (2.12a) is satisfied for two values of $k$. These are known as Wilton's ripples, after Wilton (1915). However, when $1/Re$ and $P$ are non-zero, this is not possible, due to condition (2.12b).

This seems to suggest that the discrete solution branch structure discovered in the inviscid problem, shown in figure 1(a), will not occur in our current viscous formulation. This is because in the beyond-all-orders asymptotic theory of Shelton & Trinh (2022), the values of the Bond number dividing adjacent solution branches (shown with black circles in figure 1a) were determined from the failure of a solvability condition. As the linear regime was approached ($\mathscr {E} \to 0$), these values of the Bond number tended towards those determined by Wilton (1915). Thus the lack of these linear Wilton ripples in our current viscous formulation suggests that the discrete branch structure of nonlinear solutions may not persist in the small surface tension limit, which was the asymptotic limit under which Shelton & Trinh (2022) obtained the solvability condition.

2.4.2. Linear theory for the vorticity field

We now consider linear solutions to the conformal formulation for $A(\xi,\eta )$ and the vorticity field $\omega =(A_{\xi \xi }+A_{\eta \eta })/(x_{\xi }^2+y_{\xi }^2)$. We substitute $x\sim \xi + \epsilon x_1$, $y \sim \eta + \epsilon y_1$ and $A \sim \epsilon A_1$ into the conformal system (2.8). At $O(\epsilon )$, this yields the partial differential equation $A_{1\xi }=-(A_{1 \xi \xi }+A_{1 \eta \eta })/Re$ for $\eta \leq 0$, boundary condition $A_{1 \xi }=2 Y_{1 \xi \xi } /Re$ at $\eta =0$, and decay conditions $A_{1 \xi } \to 0$, $A_{1 \eta } \to 0$ as $\eta \to -\infty$. The linear vorticity field may subsequently be determined as $\omega \sim \epsilon (A_{1 \xi \xi }+A_{1 \eta \eta })$. We solve the boundary-value problem for $A_1$ by separation of variables, which yields a Fourier series expansion in $\xi$, with coefficients depending on $\eta$. Two of the four coefficients must be zero to satisfy the decay conditions as $\eta \to -\infty$, and evaluation of the boundary condition at $\eta =0$ relates the remaining coefficients to those for the free-surface Fourier expansion from (2.10). Overall, this yields our solution for the vorticity field as

(2.13)\begin{align} \omega &\sim \epsilon \sum_{k=1}^{\infty} (2k{\rm \pi})^2\left[a_k \left\{ \exp\left({\sqrt{2k{\rm \pi}(2k{\rm \pi}+\mathrm{i}\, Re)}\,\eta}\right)+{\rm c.c.}\right\}\right.\nonumber\\ &\left.\quad{}+ \mathrm{i} b_k \left\{\exp\left({\sqrt{2k{\rm \pi}(2k{\rm \pi}+\mathrm{i}\,Re)}\,\eta}\right)-{\rm c.c.} \right\}\right] \cos(2 k {\rm \pi}\xi)\nonumber\\ &\quad -(2k{\rm \pi})^2\left[\mathrm{i} a_k \left\{ \exp\left({\sqrt{2k{\rm \pi}(2k{\rm \pi}+\mathrm{i}\,Re)}\,\eta}\right)-{\rm c.c.}\right\}\right.\nonumber\\ &\left.\quad {}-b_k \left\{ \exp\left({\sqrt{2k{\rm \pi}(2k{\rm \pi}+\mathrm{i}\,Re)}\,\eta}\right)+{\rm c.c.}\right\}\right]\sin(2 k {\rm \pi}\xi) \end{align}

(where $\textrm {c.c.}$ denotes complex conjugate), which decays with the behaviour $\exp [ ( k {\rm \pi}+[(2k{\rm \pi} )^2+Re^2]^{1/2}/2)^{1/2}\eta ]$ for $\eta \leq 0$.

3. Steadily travelling solutions

In this section, we numerically calculate steady solutions of our viscous formulation (2.5) for small values of the surface tension $B$. These correspond to travelling wave solutions that are steady in a co-moving frame.

3.1. The steady numerical method

We will use an iterative method to solve system (2.5) with Newton's method, in which each equation is evaluated at collocation points in the conformal variable $\xi$. Each derivative and Hilbert transform is evaluated efficiently in physical space by first utilising spectral relations for the Fourier multiplier of each operator in Fourier space.

We begin by assuming that an initial guess for $Y(\xi )$ is known at each of the $N$ collocation points $\xi _{l}=-1/2+l/N$. In practice, this is either a linear solution from § 2.4.1 or, for more energetic solutions, a previously computed numerical solution with different parameter values. When combined with the unknown constants, the Froude number $F$ and non-dimensional wind strength $P$, we have a total of $N+2$ unknown constants.

Each component of the governing equations is then evaluated efficiently by using spectral relations in Fourier space combined with the fast Fourier transform algorithm. Since the Fourier symbol for differentiation is $2 {\rm \pi}\mathrm {i} k$, and that for the Hilbert transform is $\mathrm {i} \times \textrm {sgn}(k)$ (where $\textrm {sgn}$ is the signum function), we have $Y_{\xi }=\mathscr {F}^{-1}[2 {\rm \pi}\mathrm {i} k\,\mathscr {F}[Y]]$ and $\mathscr {H}[Y]=\mathscr {F}^{-1}[\mathrm {i} \times \textrm {sgn}(k)\,\mathscr {F}[Y]]$. Here, $\mathscr {F}$ denotes the Fourier transform. Since an initial guess for $Y$ is known, we calculate first $X_{\xi }$ from the harmonic relation (2.5c), then $\varPsi _{\xi }$ from (2.5b), and finally $\varPhi _{\xi }$ from (2.5d).

Evaluation of our amplitude condition, the energy (2.6a), requires knowledge of $\varPhi$. This is determined spectrally with the Fourier multiplier of integration $1/(2 {\rm \pi}\mathrm {i} k)$ for $k \geq 1$, and is $0$ if $k=0$. Condition (2.6b), which correctly determines the constant of integration of $\varPhi$, is enforced in Fourier space by setting the constant level of $\varPhi X_{\xi }$ to zero.

We then evaluate Bernoulli's equation (2.5a) at the $N$ collocation points $\xi _{l}=-1/2+l/N$, and the energetic constraint (2.6a), which yields $N+1$ conditions, for the $N+2$ unknowns. We note that since solutions to the current formulation are translationally invariant, a further condition could be imposed to remove this and close the system. However, we have found that numerical convergence is faster when solving the underdetermined system with the Levenberg–Marquardt algorithm. We then seek to minimise the square of the $L^2$ norm such that it is smaller than $10^{-10}$.

For the numerical solutions presented in the following subsections, we have used $N=512$ in § 3.2 where $Re$ is fixed, and $N=1024$ in § 3.3 where $Re$ depends on $B$. Solutions were calculated on a desktop computer using fsolve in MATLAB, which usually took under a second to converge with approximately $10$ iterations. The residual was also typically minimised below $10^{-13}$. However, for some more difficult solutions, such as those near the end points of branches (solutions $a$ and $b$ in figure 3, for instance), up to $100$ iterations were required for the residual to fall below $10^{-10}$.

3.2. Viscous solutions with small surface tension

To demonstrate the effect that viscosity has on our solution profiles, we first take an inviscid gravity–capillary wave (determined in the current formulation as $Re \to \infty$), then use this as an initial guess to converge upon a numerical solution with a finite value of $Re$. The results of this are shown in figure 2, in which we start with a gravity–capillary wave with $B=0.002278$, $F=0.4307$ and $\mathscr {E}=0.4$. We then compute solutions with the same value of $B$ and $\mathscr {E}$, but different values for the Reynolds number, given by $Re=20\,000$, $10\,000$, $5000$ and $2500$. We see that the inviscid profile is symmetric about the wave crest, with parasitic capillary ripples distributed across both the forward and rear faces of the travelling wave. As the effect of viscosity increases (corresponding to decreasing $Re$), asymmetry is seen to develop, and the capillary ripples become less noticeable on the rear face of the wave. The surface profile of these viscous gravity capillary waves closely resembles the experimental observations by Ebuchi et al. (1987) and Perlin et al. (1993), although we note that this latter study produced capillary ripples whose amplitude fluctuated in time.

Figure 2. The effect of increasing viscosity, starting from an inviscid solution found as $Re \to \infty$, with $B=0.002278$, $F=0.4307$ and $\mathscr {E}=0.4$. As the viscosity increases, asymmetry develops in the parasitic capillary ripples, which are most noticeable near the forward face of the travelling wave. These solutions have $Re=\{20\,000,10\,000,5000,2500\}$, $F=\{0.4307,0.4308,0.4310,0.4310\}$ and $P=\{0.001424,0.001667,0.001770,0.002574\}$. For visibility, each profile has been shifted vertically by $0.015$.

The bifurcation structure associated with these viscous gravity–capillary waves is shown in figure 3 for fixed energy and viscosity. We have fixed $\mathscr {E}=0.4$, and have chosen the three values $Re=10\,000$, $7500$ and $5000$ to explore. Given one solution, the surrounding branch is explored in the $(B,F)$ plane by numerical continuation. We have focused on detecting the branches of solutions that exist for small values of the surface tension parameter, $B$. For inviscid solutions, a detailed bifurcation structure emerges in the small surface tension limit, shown in figure 1(a), in which a countably infinite number of adjacent solution branches pile up as $B \to 0$. The same phenomenon is not seen to occur in figure 3, in which the viscosity is held constant as the Bond number decreases. We see in figure 3(c) that while these adjacent branches do exist for $B>0.003$, they quickly disappear in the limit $B \to 0$. Furthermore, as the effect of viscosity increases, the fingering structure of the solution branches disappears at larger values of $B$. When $Re=5000$, the discrete branch structure occurs for $B>0.004$. The surface profiles of solutions labelled $a$ to $f$ in figure 3(c), with $Re=10\,000$, are shown in figure 4. The solutions shown in figures 4(c,d), which lie at the top of the solution branch in figure 3(c), resemble the parasitic ripples observed to form physically on steep travelling Stokes waves. They also have smaller values of the wind strength $P$ than the other displayed solutions, thus are expected to be more physically realisable.

Figure 3. Branches of solutions are shown in the $(B,F)$ plane for fixed energy $\mathscr {E}=0.4$. Each plot shows the solution branches computed for a different fixed value of the Reynolds number: (a) $Re=5000$, (b) $Re=7500$, and (c) $Re=10\,000$. The labelled points along the branch in (c) are the locations of the solutions plotted in figure 4.

Figure 4. The free surface $y=\zeta (x)$ for six numerical solutions across the same solution branch from figure 3(c). These solutions have $\mathscr {E}=0.4$ and $Re=10\,000$. Solutions (a) and ( f) correspond to where each side of the solution branch terminated, beyond which no further solutions could be obtained through numerical continuation.

3.3. Distinguished limit between viscosity and surface tension

In the previous subsection, we fixed the Reynolds number $Re$, and computed solution branches as the Bond number $B$ was varied. It was seen that the fingering branch structure, associated with inviscid solutions for small surface tension (figure 1a), does not emerge when $B$ is decreased and $Re$ held constant. The discrete branching structure, which occurs as $B \to 0$ in the inviscid regime, is associated with the exchange of energy between solutions dominated by gravitational energy (at the top of each solution branch) and solutions dominated by capillary energy (down the sides of each solution branch), as measured by each component of (2.6a). In this subsection, we explore the same bifurcation structure in the $(B,F)$ plane, but with a specified scaling for $Re$ given by

(3.1)\begin{equation} Re=\frac{\lambda_{\alpha}}{B^{\alpha}}. \end{equation}

The intention of this choice is to investigate possible distinguished limits of $Re$ for which viscosity has the same effect as capillarity on the parasitic capillary ripples, and also investigate whether a discrete branching structure can emerge as $B \to 0$.

We compute the bifurcation structure of solutions, between $B=0.001$ and $0.002$, for $\alpha =1$, $2$ and $3$. The constant $\lambda _{\alpha }$ in (3.1) is specified such that the distinguished curve passes through the point $B=0.005$ and $Re=5000$, which yields $\lambda _{\alpha }=5000 \times 0.005^{\alpha }$. These solution branches are shown in the $(B,F)$ bifurcation diagrams of figures 5(a–c), for $\alpha =1$, $2$ and $3$, respectively. The effect of viscosity is largest in the solutions along the branch in figure 5(a). We see that in all three cases, the discrete branching structure recedes as $B$ decreases, such that below a certain value of $B$, all solutions found by numerical continuation are dominated by the effect of gravity. The parasitic capillary ripples present in these solutions appear as a small perturbation to the base nonlinear viscous-gravity wave. Even with $\alpha =3$ in figure 5(c), the branch of solutions can be continued to $B=0$ in this manner.

Figure 5. Branches of solutions in the $(B,F)$ plane for fixed energy $\mathscr {E}=0.4$. The Reynolds number $Re$ is specified by $Re=\lambda _{\alpha }/B^{\alpha }$ from (3.1), thus the effect of viscosity decays proportionally to the surface tension. We have (a) $\alpha =1$, (b) $\alpha =2$, and (c) $\alpha =3$. Note that $\lambda _1=25$, $\lambda _2=0.125$ and $\lambda _3=0.000625$ are chosen such that the distinguished limit (3.1) intersects with $Re=5000$ and $B=0.005$. The marked points correspond to the solutions shown in figure 6.

We note that numerical verification of these trends for larger values of $\alpha$ would require computation of the bifurcation diagram to smaller values of $B$ to check if the discrete branching structure persists. However, due to the beyond-all-order (as $B \to 0$) nature of the parasitic capillary ripples, there exists a value of $B$ below which these are not captured by double precision accuracy. When this occurs, the vertical branches of the bifurcation diagram (in which energy is transferred into the oscillatory ripples) are unable to be computed through numerical continuation from the gravity-dominated solutions. In practice, for $\mathscr {E}=0.4$, this occurs when $B\leq 0.0008$. It is possible that the discrete branching structure of solutions is recovered in the small-surface-tension limit only when the effect of viscosity, as measured by $1/Re$, is exponentially small in comparison to $B$. However, due to the lack of any multiple-scales asymptotic theory for this regime, even in the inviscid framework, there is no analytical evidence to support this.

The parasitic capillary ripples present in the solution profile are shown in figure 6 for $\alpha =1$ and $2$. These were computed with $B=0.0015$, and correspond to the marked locations in figures 5(a,b). In this figure, we have subtracted the leading- and first-order solutions $Y_0(\xi )$ and $B\,Y_1(\xi )$, and to view the ripples in more detail, defined

(3.2)\begin{equation} y_{ripples} = Y(\xi)-Y_0(\xi)-B\,Y_1(\xi). \end{equation}

The leading-order solution $Y_0$ was computed numerically with $B=0$ and $\mathscr {E}=0.4$. For the first-order solution, rather than substituting asymptotic expansions into the governing system (2.5) and then solving the $O(B)$ equation numerically, we have estimated by computing $Y_1 \sim (Y-Y_0)/B$ with $B=0.0005$. We see in figure 6(a) that the ripple amplitude is highly asymmetric under this specification of $Re$, and that the ripple wavelength slightly decreases as we travel to the right away from the wave crest. When the effect of viscosity is weaker, as in figure 6(b) with $\alpha =2$, the asymmetry is less prominent.

Figure 6. The parasitic capillary ripples present in the free surface for two of the solutions with $B=0.0015$ in the bifurcation diagram of figure 5. These profiles have been estimated numerically by calculating $y_{ripples}=Y-Y_0-BY_1$, where $Y_0(\xi )$ is the solution found with $B=0$, and $Y_1(\xi )$ has been estimated from $(Y-Y_0)/B$ with $B=0.0005$. The effect of viscosity is stronger in solution (a), which produces substantial asymmetry in the parasitic ripple profile.

This behaviour is expected when analysing our previous exponential asymptotic theory for nonlinear gravity–capillary waves. This is because the inviscid parasitic ripples are generated by the Stokes phenomenon, which arises from branch points in the analytic continuation of the leading-order gravity wave solution, and this yields their exponential scaling

(3.3)\begin{equation} y_{ripples} \sim A(\xi)\,\mathrm{e}^{-\chi(\xi)/B} +{\rm c.c.}. \end{equation}

Here, $A$ is an amplitude function, and $\chi$ controls the exponential scaling as $B \to 0$. When $Re=O(1/B)$ as in figure 6(a), the viscous terms in the governing equations (2.5) will feed into the differential equations for both $A$ and $\chi$. However, when $Re=O(1/B^2)$ as in figure 6(b), the viscous terms will alter only the amplitude equation for $A$.

This is demonstrated in the semi-log plot of figure 7, in which the parasitic ripple amplitude is shown against $1/B$, for $\alpha =1$ and $2$. The slope of this behaviour corresponds to the exponential scaling in (3.3) evaluated at $\xi =1/2$, near to which the amplitude was measured, which is given by $-\chi (1/2)$. For inviscid capillary ripples, this gradient was determined analytically as $-0.0077$, as shown in figure 1(c), which is very close to the measurement $-0.0076$ for $Re=O(1/B^2)$ in figure 7. However, when $Re=O(1/B)$, this gradient is steeper, with value $-0.0097$. We also note that for $Re=O(1/B^2)$ in figure 7, all solutions along the branch have this asymptotic scaling as $B \to 0$, in contrast to that in the inviscid regime (see the spikes in figure 1c). It is therefore expected that the beyond-all-order solvability mechanism dividing adjacent branches of inviscid solutions (Shelton & Trinh 2022) will not manifest in the future exponential asymptotic analysis of our current viscous formulation.

Figure 7. The parasitic capillary ripple amplitude, measured near $\xi =0.5$, against the Bond number. This amplitude is shown for each solution from the branches in figures 5(a,b). The linear behaviour in a semi-log plot is indicative that this amplitude is exponentially small as $B \to 0$. For $Re=\lambda _1/B$ in (a), the gradient is approximately $-0.0097$, and for $Re=\lambda _2/B^2$ it is $-0.0076$. This is compared with the inviscid asymptotic theory of Shelton & Trinh (2022), which yields gradient $-0.0077$.

4. Time-dependent results

We now study the temporal stability of the steady solutions found in § 3. First, in § 4.1, we analyse the linear stability of these solutions. Growth rates are obtained for small superharmonic perturbations to the steady solutions by solving a linear eigenvalue problem for the perturbations. The eigenvalue $\sigma$ corresponds to the growth rate of the perturbations that contain a term of the form $\mathrm {e}^{\sigma t}$. We will find that the steady solutions are superharmonically stable since $\textrm {Re}[\sigma ] \leq 0$. Second, in § 4.2, we study the nonlinear stability of the solutions. This is achieved by considering the unsteady time-evolution problem, which was formulated in § 2.2 as a set of time-evolution equations for $Y(\xi,t)$ and $\varPhi (\xi,t)$ in terms of the conformal variable $\xi$. We will pick an initial condition at $t=0$, and for $t>0$ set all associated constants $B$, $F$, $Re$ and $P$ to the values of the steady solution whose nonlinear stability is to be studied. This allows us to comment on the ‘global’ stability and attractiveness of these solutions, as $t \to \infty$, with respect to different initial conditions, such as a steep gravity wave or an almost flat free surface.

4.1. Linear stability of the steady solutions

We now examine the linear stability of the steady solutions from § 3. We begin by detailing the methodology used to analyse the linear stability of these steady solutions. The method presented here is similar to that used by Tiron & Choi (2012) for pure capillary waves, and by Blyth & Părău (2022) for constant vorticity waves, with the exception that we study only superharmonic modes that fit with within the assumed period $-1/2 \leq \xi \leq 1/2$, in order to draw comparisons with the time-evolution problem studied in § 4.2.

4.1.1. Numerical implementation

We begin by linearising about a steady solution by writing $X \sim X_0(\xi ) + \epsilon \,X_1(\xi,t)$, $Y \sim Y_0(\xi ) + \epsilon \,Y_1(\xi,t)$, $\varPhi \sim \varPhi _0(\xi ) + \epsilon \,\varPhi _1(\xi,t)$ and $\varPsi \sim \varPsi _0(\xi ) + \epsilon \,\varPsi _1(\xi,t)$. Here, $X_0$ $Y_0$, $\varPsi _0$ and $\varPhi _0$ are solutions of equations (2.5a–d), whose temporal stability we will analyse by studying growth/decay rates of the perturbation quantities. Substitution of these expressions into the time-dependent equations (2.4a–d) yields at $O(\epsilon )$ four equations for $X_1$, $Y_1$, $\varPhi _1$ and $\varPsi _1$. Two of these equations are given by $X_{1 \xi }=-\mathscr {H}[Y_{1 \xi }]$ and $\varPsi _{1}=\mathscr {H}[\varPhi _{1}]$. The other two equations, obtained at $O(\epsilon )$ from (2.4a) and (2.4b), have coefficients involving the leading-order quantities, for which in practice we calculate expressions with a symbolic programming language. Next, for the time-dependent perturbations, we pose a Floquet ansatz of the form

(4.1)\begin{equation} \{ X_1,Y_1,\varPhi_1,\varPsi_1 \}=\mathrm{e}^{\sigma t} \sum_{m={-}M}^{M} \{a_m,b_m,c_m,d_m \}\, \mathrm{e}^{2 m {\rm \pi}\mathrm{i} \xi}, \end{equation}

where the real part of the complex-valued constant $\sigma$ will correspond to the growth/decay rate in time of the perturbation. Note that the ansatz (4.1) would differ if used to analyse the stability of time-dependent solutions, such as time-periodic standing waves whose stability was studied by Wilkening (2020). Since $\mathscr {H}[\mathrm {e}^{2 m {\rm \pi}\mathrm {i} \xi }]=\mathrm {i} \times \textrm {sgn}{(m)}\,\mathrm {e}^{2 m {\rm \pi}\mathrm {i} \xi }$, we immediately have from the $O(\epsilon )$ harmonic relations that $a_m=-\mathrm {i}\times \textrm {sgn}(m)\,b_m$ and $d_m=\mathrm {i}\times \textrm {sgn}(m)\,d_m$. We note that given a solution $\{\sigma,b_{-M},\ldots,b_{M},c_{-M},\ldots,c_{M}\}$ to the $O(\epsilon )$ equations, there exists another solution given by $\{ \sigma ^{*},b_{M}^{*},\ldots,b_{-M}^{*},c_{M}^{*},\ldots,c_{-M}^{*} \}$; when these are combined, they yield a real-valued solution to the $O(\epsilon )$ equations. Here, $^{*}$ denotes the complex conjugate.

To solve for the unknowns $\{\sigma,b_{-M},\ldots,b_{M},c_{-M},\ldots,c_{M}\}$, we collocate the two remaining $O(\epsilon )$ equations at the $N:=2M+1$ points $\xi =-1/2+l/N$, for $l=0,\ldots,N-1$. This yields $2N$ algebraic equations for $2N+1$ unknowns, which we write as the linear system $\sigma \boldsymbol {L} \boldsymbol {v} = \boldsymbol {R} \boldsymbol {v}$. Here, the eigenvalue $\sigma$ is the temporal coefficient from (4.1), $\boldsymbol {v}=[b_{-M},\ldots,b_{M},c_{-M},\ldots,c_{M}]^{T}$ is a $2N \times 1$ eigenvector consisting of the unknown Fourier coefficients, and $\boldsymbol {L}$ and $\boldsymbol {R}$ are $2N \times 2N$ matrices. The matrix $\boldsymbol {L}$ contains entries associated with the time derivative in the $O(\epsilon )$ equations, given by $\boldsymbol {L}_{ij}=\exp (2 (j-1-M) {\rm \pi}\mathrm {i} \xi _i)$ for $1 \leq i \leq N$ and $1 \leq j \leq N$, $\boldsymbol {L}_{ij}=\exp ({2 (j-N-1-M) {\rm \pi}\mathrm {i} \xi _{i-N}})$ for $N+1 \leq i \leq 2N$ and $N+1 \leq j \leq 2N$, and $\boldsymbol {L}_{ij}=0$ otherwise. For details on the construction of $\boldsymbol {R}$, we refer the reader to either the work of Blyth & Părău (2016) (cf. their equation (4.17)) for a similar boundary-integral problem, or the MATLAB code in our supplementary material available at https://doi.org/10.1017/jfm.2024.1227, which constructs $\boldsymbol {R}$ column by column.

We solve the generalised eigenvalue problem $\sigma \boldsymbol {L} \boldsymbol {v} = \boldsymbol {R} \boldsymbol {v}$ in MATLAB with the QZ algorithm. The code that implements this method, and produces the results shown next in § 4.1.2, is provided alongside this work as supplementary material. When $M=127$, solving this generalised eigenvalue problem on a standard desktop computer takes approximately $3.5$ s.

4.1.2. Linear stability results

We now calculate the growth rates $\sigma$ and the corresponding eigenvectors $Y_1(\xi,0)$, for instance, for perturbations to steady solutions calculated previously in § 3. The parameter values for the steady solutions, whose linear and nonlinear stability we will investigate in detail, are provided in table 1. We will investigate the effect of viscosity and surface wind forcing on the solution stability. Each of these solutions has the same Bond number $B=0.0026$, and different values of the Reynolds number, which are taken to be $Re=5000$, $7500$ and $10\,000$. Recall that the magnitude of the surface wind forcing $P$ was an unknown in the steady formulation. The solution with the largest viscous dissipation, $Re=5000$, has a larger value of $P$ to compensate.

Table 1. Parameter values used in the stability results of figure 8 and the time-dependent simulations of figures 11 and 12. These were obtained from the steady solutions of § 3, iterated upon such that the residual, defined by the square of the $L^{2}$-norm of (2.5), is smaller than $10^{-20}$.

The real and imaginary parts of the growth rates $\sigma$ for perturbations to each of these steady solutions are shown in figure 8. While only eigenvalues with $\textrm {Im}[\sigma ] \geq 0$ have been shown in figure 8, complex conjugate values are also solutions of the generalised eigenvalue problem. The real part of these eigenvalues corresponds to the growth rate of time-dependent perturbation. Note that our formulation does not possess time-reversal symmetry, therefore there is no fourfold symmetry in the eigenvalue spectrum (in which $\sigma$, $\sigma ^*$, $-\sigma$, $-\sigma ^*$ would be solutions). Such fourfold symmetry is associated with inviscid formulations that are Hamiltonian (Deconinck & Oliveras 2011). Aside from $\sigma =0$, corresponding to translational invariance of the steady formulation (a small perturbation of a steady solution can result in another steady solution), all eigenvalues have negative real part, $\textrm {Re}[\sigma ]<0$. The solutions are therefore superharmonically stable.

Figure 8. (a–c) The complex-valued eigenvalues $\sigma$ for small-amplitude perturbations to the steady solutions from § 3. The parameter values for these steady solutions are given in table 1, and these eigenvalues were calculated with the numerical methodology described in § 4.1.1. From (4.1), the real part of $\sigma$ is the temporal growth rate of the perturbation. The eigenvalue shown in black has the largest real part, with $\textrm {Re}[\sigma ]= \{ -0.059114 ,-0.045994, -0.0371606 \}$. The first eigenvalue on the real axis, with $\textrm {Im}[\sigma ]=0$, is shown in grey with values $\textrm {Re}[\sigma ]= \{ -0.230998,-0.232593,-0.199846 \}$. These values predict the growth rates observed later in the time-evolution simulations in figures 11 and 12. (d) Perturbations to the free-surface elevation predicted by the linear stability analysis, which correspond to the eigenvectors for the two labelled eigenvalues in the inset of (b). Here, (a) $Re = 5000$, (b) $Re = 7500$, (c) $Re = 10\,000$, (d) $Re = 7500$.

Note that under the inviscid limit $Re \to \infty$, the growth rates are observed to tend towards the imaginary axis. This is demonstrated in figure 9 for solutions without surface tension. Modulational and high-frequency instabilities (Creedon, Deconinck & Trichtchenko 2022) do not emerge here as we study only superharmonic perturbations with zero Floquet exponent in ansatz (4.1).

Figure 9. The real and imaginary components of the growth rates $\sigma$ from (4.1), for steady solutions with zero surface tension ($B=0$). The growth rates shown with white circles have $Re=5 \times 10^{3}$, $P=8.229 \times 10^{-4}$ and $F=0.4110$. The growth rates shown with grey circles have $Re=1 \times 10^{6}$, $P=4.115 \times 10^{-6}$, and $F=0.4110$. Black contours between these show the behaviour of the eigenvalues for $50$ intermediary values of $Re$. These tend towards the imaginary axis as $Re \to \infty$.

Typically, the eigenvalues with the largest real part will dominate the evolution of time-dependent simulations for solutions close to the steady solution. We investigate this in detail later, in § 4.2, where it is seen that different choices of initial conditions can result in different rates of convergence for the transient behaviour towards the stable steady state. The growth rates that dominate the time-dependent simulations of § 4.2 are: (i) the eigenvalue with the largest real part, shown with a black circle in figure 8; and (ii) the eigenvalue on the real axis with $\textrm {Im}[\sigma ]=0$, shown with a grey circle in figure 8. The steady solution with $Re=7500$, for which the growth rates are shown in figure 8(b), has $\textrm {Re}[\sigma ]=-0.045994$ and $\textrm {Re}[\sigma ]=-0.232593$ for these two values.

The perturbation solution $Y_1$ corresponding to each of these two growth rates is shown in figure 8(d). Here, we have computed the contribution to the free-surface elevation $Y_1(\xi,0)$ from the corresponding eigenvector and the complex conjugate eigenvector by the relation $Y_1(\xi,0)=\sum _{m=-M}^{M}[b_m\,\mathrm {e}^{2 m {\rm \pi}\mathrm {i} \xi }+b_m^{*}\,\mathrm {e}^{-2m {\rm \pi}\mathrm {i}\xi }]$ from (4.1). The eigenvector for the $\textrm {Re}[\sigma ]=-0.045994$ mode, shown in black, has more of an effect on the wave amplitude, whereas that for the $\textrm {Re}[\sigma ]=-0.232593$ mode, shown in grey, has a larger contribution to the oscillatory capillary ripples. We now turn to the nonlinear time-evolution problem to observe how these stability results emerge in practice.

We have also studied the linear stability of solutions for larger values of $Re$, to see if a critical Reynolds number exists beyond which instability emerges. However, all solutions that we have tested in the range $0 \leq Re \leq 10^{6}$ have been linearly stable.

4.2. Nonlinear stability of the steadily travelling solutions

We now perform time-evolution simulations to test the nonlinear stability of the steady solutions from § 3. To measure convergence towards the steady solution as $t \to \infty$, we will use the norm $\lvert \mathscr {E}(t)-\mathscr {E} \rvert$, which compares the wave energy $\mathscr {E}(t)$ in (2.6a) to that of the intended steady solution $\mathscr {E}$. The initial condition that we use at $t=0$ will be chosen as one of the following:

1. (i) an inviscid gravity wave that is a solution of the steady equations (2.5) solved for with $\mathscr {E}=0.4$, $B=0$, $P=0$ and $1/Re=0$;

2. (ii) a small-amplitude cosine profile $Y(\xi )=10^{-5}\cos (2 {\rm \pi}\xi )$.

For $t>0$, we then set the values for $B$, $F$, $Re$ and $P$ to those of a selected steady solution, whose nonlinear stability we wish to analyse. The parameter values used in this subsection are shown in table 1.

Surface profiles for both of these time-evolution simulations are shown in figure 10. We see in figure 10(a) that for the initial condition of a steep gravity wave, parasitic ripples form quickly within the interval $0 \leq t \leq 5$ due to the nonlinearity of the initial condition. When the initial condition is of smaller amplitude, as in figure 10(b), it takes much longer for the parasitic ripples to form. This is due to the surface wind forcing inducing a small growth rate in the solution, and the ripples then become noticeable within the plotted interval $200 \leq t \leq 300$. For the simulation with the large-amplitude initial condition, the solution tends towards the steady solution in an oscillatory manner, as shown in the inset of figure 10(a), which is associated with a complex-valued growth rate from the linear stability analysis. The simulation with a small-amplitude initial condition approaches the steady solution monotonically, which is associated with a real-valued growth rate in the linear stability analysis. These two dominant growth rates were highlighted in figure 8(a).

Figure 10. The time evolution of the free surface is shown from the initial conditions at $t=0$ of (a) a steep viscous-gravity wave with $\mathscr {E}=0.4$, and (b) a linear cosine profile with amplitude $10^{-5}$. The parameter values used for $t>0$ are given in table 1 for $Re=5000$, and the time interval between displayed solutions is $0.1$ in (a), and $2$ in (b). The insets show the difference in height of the wave crest of the steady solution, $y_s$, and that of the current numerical solution, $y$. For the simulation in (a), the solution approaches the steady solution in an oscillatory manner in time, while the solution in (b) approaches the same steady solution monotonically.

We begin by examining the stability of three different viscous solutions when starting from an inviscid gravity wave as an initial condition. The difference between the energy of the time-evolved solution and that of the target solution is shown in figure 11. In each case studied, the energy converges to within the tolerance of the computed steady solution, which occurs at approximately $t=250$ when $Re=5000$, $t=300$ when $Re=7500$, and $t=400$ when $Re=10\,000$. Convergence is seen to occur more quickly when the effect of viscosity, and the associated wind strength $P$, are larger. Note that in each of the three cases shown in figure 11, while the initial gravity wave at $t=0$ has $\mathscr {E}=0.4$, this value increases once the associated parameters are set to the values shown in table 1. For instance, at the first time step, we have $\mathscr {E}=0.4442$ in figure 11(a), $\mathscr {E}=0.4439$ in figure 11(b), and $\mathscr {E}=0.4437$ in figure 11(c). Two stages of convergence are observed in figure 11. First, there is the faster stage, during which the norm $\log ( \lvert \mathscr {E}-0.4 \rvert )$ decreases to $-16$, and then there is the slower stage, during which the energy oscillates about $0.4$, with the norm decreasing from $-16$ to $-27$. The gradients for each of these stages, predicted by the linear stability analysis of § 4.1, are shown within each plot. Both the fast and slow stages of convergence are predicted by the real part of the eigenvalues $\sigma$ obtained from this stability analysis. The fast growth rate in the initial stage of convergence is associated with the largest real-valued eigenvalue ($\textrm {Im}[\sigma ]=0$). For the $Re=7500$ simulation in figure 11(b), we have $\textrm {Re}[\sigma ]=-0.232593$. The slow growth rate in the last stage of convergence is associated with the complex-valued eigenvalue ($\textrm {Im}[\sigma ] \neq 0$) whose real part is largest. For the $Re=7500$ simulation, this is $\textrm {Re}[\sigma ]=-0.045994$. Note that in each of the three simulations shown in figure 11, the differing values at which $\mathscr {E}-0.4$ plateaus is due to the accuracy of the parameter values given in table 1.

Figure 11. The magnitude of the difference between the energy $\mathscr {E}(t)$ and the ‘target’ energy $0.4$ in a semi-log plot for the three cases (a) $Re=5000$, (b) $Re=7500$, and (c) $Re=10\,000$. The initial condition at $t=0$ was chosen to be a steep gravity wave with $\mathscr {E}=0.4$, $B=0$, $P=0$ and $1/Re=0$, and the parameter values used for $t>0$ in each of these simulations are given in table 1. The annotated gradients are predictions from the linear stability results of figure 8.

Convergence towards the same steady solutions is also studied in figure 12 from a small-amplitude initial condition at $t=0$. In figures 12(a i), 12(b i) and 12(c i), we plot the behaviour of $\log (\mathscr {E})$ in time to show the initial growth rate away from the linear solution. We then plot in figures 12(a ii), 12(b ii) and 12(c ii) the behaviour of $\log (\lvert \mathscr {E}-0.4 \rvert )$ in time to show the rate of convergence towards the steady solution. The total time taken for $\log (\lvert \mathscr {E}-0.4\rvert )$ to reach $-27$ is longer than that for the simulations in figure 11, which is due to the duration of the initial growth away from the linear profile at $t=0$. Once $\lvert \mathscr {E}-0.4 \rvert$ is smaller than $10^{-2}$, the steadily travelling solution is approached rapidly. The rate of convergence here is also predicted by the linear stability analysis of § 4.1. However, unlike the transient behaviour in figure 11, which contained two different growth rates, only one growth rate emerges for the simulations in figure 12. This growth rate is associated with not the complex-valued eigenvalue with largest real part (which emerged as a transient in figure 11), but rather the largest real-valued eigenvalue. For instance, the rate of convergence in the simulation with $Re=7500$ in figure 12 is controlled by the real-valued eigenvalue $\sigma =-0.232593$, and no transient behaviour associated with $\textrm {Re}[\sigma ]=-0.045994$ is present.

Figure 12. Convergence is shown for the initial condition at $t=0$ of a small-amplitude cosine profile $Y(\xi )=10^{-5}\cos (2 {\rm \pi}\xi )$. The parameter values used in each of these simulations are given in table 1, for (a) $Re=5000$, (b) $Re=7500$, and (c) $Re=10\,000$. (a i,b i,c i) Plots of $\mathscr {E}(t)$ show the initial growth rate from the initial condition. (a ii,b ii,c ii) Plots of $\lvert \mathscr {E}(t)-0.4 \rvert$ show the final convergence rate towards the steadily travelling solution. The annotated gradients are predictions from the linear stability results of figure 8.

In summary, the growth rates that emerge in this time-evolution problem can be determined from the formal linear stability analysis of § 4.1. However, it is not clear in advance which growth rate will emerge as the dominant transient behaviour; this depends on the choice of initial condition. The large-amplitude initial condition initially produces fast convergence, which subsequently is dominated by a slower mode (figure 11), whereas the convergence that emerges from the small-amplitude initial condition contains only the fast mode (figure 12). Further simulations are shown in a total energy versus capillary energy phase space in figure 13. The simulations in this figure have different values of the surface tension parameter $B$. Those that begin from a large-amplitude initial condition converge towards the steady solution in an oscillatory manner, as observed for other parameter values in figure 11. The radius of these oscillations decays rapidly as $B$ decreases, which is possibly related to the capillary modes of the steady solution being exponentially small as $B \to 0$ (figure 7). Conversely, the simulations with a small-amplitude initial condition tend towards the steady solution monotonically.

Figure 13. Phase space diagram for time-dependent simulations. We show the energy of the solution $\mathscr {E}$ from (2.6a) against the capillary energy defined by the middle term in the integrand on the right-hand side of (2.6a). The white circles indicate the positions of six steady solutions with $\mathscr {E}=0.4$ and $Re=7500$. These have $B=\{0.00222,0.0024,0.0026,0.00287,0.0033,0.00365\}$, in order from left to right in the figure. Two simulations have been performed for each set of parameter values associated with these steady solutions. Similar to the simulations from figures 10–12, one simulation (black lines/dashes) begins from a large-amplitude gravity wave with $\mathscr {E}=0.4$, and the other (grey lines/dashes) from a small-amplitude cosine profile.

Note that while the simulations that we show here all result in convergence, indicating nonlinear stability, this is not true in general. For instance, starting with a small-amplitude initial condition with multiple Fourier modes (in contrast to $Y(\xi,0)=10^{-5}\cos (2 {\rm \pi}\xi )$ used in figure 12) results in general unsteady motion persisting at large time. Similar time-evolution simulations have been performed by Hung & Tsai (2009) and Murashige & Choi (2017), who all used a steep gravity wave as an initial condition. However, these works studied the initial formation of capillary ripples, rather than characterising convergence towards steady solutions as $t \to \infty$ by selecting appropriate parameter values. Hung & Tsai (2009) observed that the amplitude of the parasitic capillary ripples fluctuated in time, much like that seen in the experiments of Perlin et al. (1993). We observe similar behaviour in time-evolution simulations when a different value of the wind strength $P$ is used.

5. Conclusion and discussion

While we have focused on viscous gravity–capillary waves travelling on water of infinite depth, the formulation with finite depth is also of significant interest. This is because viscosity requires a no-slip condition on the lower surface. A nonlinear potential flow model incorporating the no-slip condition at a lower boundary has previously been developed by Dutykh & Dias (2007), by first applying a Helmholtz decomposition to the linear Navier–Stokes equations, and then introducing a viscous boundary layer at the substrate. The dominant effect of viscosity in their model, of $O(Re^{-1/2})$, occurs from the lower boundary condition, rather than that of $O(Re^{-1})$ in the surface boundary conditions. We therefore expect that the discrete branching structure associated with capillarity in figure 4 would be further suppressed by viscosity in the presence of a finite-depth fluid.

It may also be possible to support steadily travelling surface wave solutions of a viscous fluid, by considering the effect of gravity on a flow down an inclined plane, rather than the surface wind forcing considered in this paper. The effect of capillarity and viscosity on the parasitic solutions that would emerge in this formulation remains unknown. More generally, it is also of interest to ask whether these viscous boundary layer techniques may be extended to derive potential flow formulations for water waves in the presence of more complicated lower topography and submerged obstacles. Such formulations have been studied extensively in the inviscid regime by e.g. Dagan & Tulin (1972), Dias & Vanden-Broeck (1989), Lustri, McCue & Binder (2012) and Ambrose et al. (2022), for which recirculation zones can emerge with the inclusion of viscosity (Biswas, Breuer & Durst 2004). The analytical insight available in these viscous potential flow models would provide a distinct advantage to that from the Navier–Stokes equations with a free boundary.

Due to the assumed periodicity of our solutions, the linear and nonlinear temporal stability results in § 4 should be regarded only as evidence for their superharmonic stability with respect to perturbations that lie within the periodic window. Their modulational stability (with respect to subharmonic perturbations wider than the periodic domain) has not been investigated in this work. While nonlinear gravity waves are known to exhibit the modulational instability (Longuet-Higgins 1978), there is the possibility that this can be suppressed by the effect of viscosity, which has been studied in a damped nonlinear Schrödinger equation by Segur et al. (2005).

Previously, Longuet-Higgins (1995) had obtained approximate solutions for the asymmetric parasitic capillary waves that form upon steep gravity waves with viscosity, which built upon an earlier theory by Longuet-Higgins (1963). However, the methodology used there involves a number of modelling assumptions, rather than a systematic asymptotic study of the governing equations. The direct asymptotic analysis of these, in the limit $B \to 0$, is therefore of interest. In the inviscid study by Shelton & Trinh (2022, 2024), the capillary ripples observed on steep gravity–capillary waves were determined asymptotically. These were obtained by resolving the Stokes phenomenon generated from branch points in the analytic continuation of the leading-order gravity wave. We note also that the solutions upon which we have focused in this work are steady in a co-moving frame, which may not be the precise manner in which they occur physically. In experiments performed by Perlin et al. (1993) in a wave tank with a wave maker, parasitic ripples were seen to travel at the same speed as the propagating gravity wave, but the ripples had an amplitude that varied in time. Beyond-all-order asymptotic techniques are well suited to investigate if such time-dependent behaviour can emerge in this formulation for small surface tension, and the application of our viscous formulation to study this is an exciting area of future research. Here, the leading-order gravity wave solution would be time-independent, and the exponentially subdominant order (the parasitic capillary ripples) and associated Stokes lines would be time-dependent. The future extension of this methodology to our viscous formulation is also expected to yield insight into the competing effects of viscosity and surface tension on the observed ripple asymmetry.