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 (Reference Jeffreys1925) 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 (Reference Fedorov and Melville1998), 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 (Reference Fedorov and Melville1998), 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 (Reference Janssen2004) 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 (Reference Miles1957) 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. (Reference Dyachenko, Zakharov and Kuznetsov1996) and Choi & Camassa (Reference Choi and Camassa1999) 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. Reference Shelton, Milewski and Trinh2021; Shelton & Trinh Reference Shelton and Trinh2022) 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. (Reference Dias, Dyachenko and Zakharov2008) 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 (Reference Wilton1915). 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 (Reference Shelton and Trinh2022), 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 (Reference Wilton1915). 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 (Reference Shelton and Trinh2022) 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.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. (Reference Ebuchi, Kawamura and Toba1987) and Perlin et al. (Reference Perlin, Lin and Ting1993), 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 Reference Shelton and Trinh2022) 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 (Reference Shelton and Trinh2022), which yields gradient $-0.0077$.

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 (Reference Tiron and Choi2012) for pure capillary waves, and by Blyth & Părău (Reference Blyth and Părău2022) 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 (Reference Wilkening2020). 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 (Reference Blyth and Părău2016) (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 Reference Deconinck and Oliveras2011). 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 Reference Creedon, Deconinck and Trichtchenko2022) 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 (Reference Hung and Tsai2009) and Murashige & Choi (Reference Murashige and Choi2017), 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 (Reference Hung and Tsai2009) observed that the amplitude of the parasitic capillary ripples fluctuated in time, much like that seen in the experiments of Perlin et al. (Reference Perlin, Lin and Ting1993). We observe similar behaviour in time-evolution simulations when a different value of the wind strength $P$ is used.