2.1. Lagrangian governing equations

In the Lagrangian frame, fluid particles (or parcels) are tracked using fixed labels $(a,b)$ and time $t$ as independent variables. Time derivatives are taken with respect to fixed particle labels $(a,b)$ . Trajectories are then described by $\{(x(a,b,t), z(a,b,t))\}_{t \in \mathbb{R}^+}$ . The mapping between label space $(a,b)$ and physical space $(x,z)$ is invertible, as no two particles can occupy the same physical location simultaneously. This implies a non-zero Jacobian:

(2.1) \begin{equation} \mathcal{J} \overset {\textrm{def}}{=} \frac {\partial (x,z)}{\partial (a,b)}= x_az_b-x_bz_a\neq 0, \end{equation}

where $x_a$ denotes the partial derivative of the $x$ -coordinate with respect to particle label $a$ , and similarly for the other derivatives. While this mapping is generally time-dependent, our choice of labels based on particle locations at time $t=0$ ensures that for an incompressible fluid, the Jacobian is constant along a trajectory, i.e. $\dot {\mathcal{J}}=0$ . Beyond these constraints, there remains considerable flexibility in particle label assignment.

We choose particle labels such that the free surface corresponds to $b=0$ . We assign labels based on particle positions at time $t=0$ , with label $a$ following the $x$ -axis. Our two-fluid system consists of water ( $b\lt 0$ ) and air $(b\gt 0)$ . This Lagrangian formulation offers a significant geometric advantage: the time-varying free surface can be described as a simple, fixed plane in label space rather than an additional function $z=\eta (x,t)$ . The conservation of particle labels implies that particles at the surface remain at the surface for all time. This approach directly circumvents the well-known difficulties of defining physically meaningful averages in the Eulerian frame near an undulating interface, a challenge that has motivated the development of analogous coordinate transformations in experimental and numerical work (e.g. Sullivan, McWilliams & Patton Reference Sullivan, McWilliams and Patton2014; Buckley & Veron Reference Buckley and Veron2016). Furthermore, this formulation has a significant analytical advantage in that Lagrangian approximations are generally more accurate and converge faster than Eulerian ones of the same order. For moderately steep waves, for instance, an $n$ th-order Lagrangian approximation of the surface can match the accuracy of an Eulerian approximation of order $n+2$ (Clamond Reference Clamond2007).

The Euler equations in Lagrangian coordinates are given by (Lamb Reference Lamb1924, Art.15)

(2.2a) \begin{align} \mathcal{J}\ddot {x} + \frac {1}{\rho }(p_az_b - p_bz_a)&=0, \end{align}

(2.2b) \begin{align} \mathcal{J}\ddot {z} + \frac {1}{\rho }(p_bx_a-p_ax_b + \mathcal{J} g)&=0, \end{align}

where $p$ is pressure, $\rho$ is density and $g$ is the gravitational constant. Since particle density remains constant along trajectories, we take $\rho =\rho _0$ as constant.

Lastly, the vorticity is defined as

(2.3) \begin{equation} \varGamma \overset {\textrm{def}}{=} v_x-u_y = \frac {1}{\mathcal{J}}(\dot {x}_a x_b - \dot {x}_bx_a + \dot {z}_az_b - \dot {z}_bz_a). \end{equation}

Vorticity is always conserved, so that $\dot {\varGamma }=0$ . For irrotational flows, $\varGamma =0$ .

As in the classical stability analysis (Miles Reference Miles1957), we neglect viscosity and turbulent stresses, and model the evolving perturbation as a normal mode on a prescribed, smooth background shear $U^{(0)}$ . We thus interpret the theory as most applicable to non-breaking conditions in which the airflow remains predominantly attached and the wave steepness is small, $\epsilon \ll 1$ . Laboratory observations indicate that airflow separation is uncommon over low-wind, small-slope waves, but becomes increasingly prevalent as wind forcing and wave slope increase (Buckley et al. Reference Buckley, Veron and Yousefi2020). Within this attached-flow regime, we isolate the wave-induced mean flow as the leading higher-order correction, yielding the order $\epsilon ^2$ modification to the Miles growth rate derived in $\S$ 3.

2.2. Asymptotic expansions

We consider permanent, progressive, monochromatic, spatially periodic, two-dimensional waves progressing with velocity $c$ (Clamond Reference Clamond2007; Pizzo et al. Reference Pizzo, Lenain, Rømcke, Ellingsen and Smeltzer2023). The coordinates $(x,z)$ are expanded as

(2.4) \begin{equation} x= a +U(b)t + \sum _{n=1}^\infty x_n(b)\sin (\theta _n), \quad z = b +z_0(b) + \sum _{n=1}^\infty z_n(b)\cos (\theta _n), \end{equation}

where $\theta _n\overset {\textrm{def}}{=} nk(a-(c-U(b))t)$ , $k$ is the wavenumber and $U$ is the Lagrangian mean flow (and may be identified with $\overline {u}_L$ evaluated in Lagrangian coordinates). The mean surface level is $z_0(0)\overset {\textrm{def}}{=} \langle z\rangle |_{b=0}$ , where $\langle \boldsymbol{\cdot }\rangle = ({k}/{2\pi })\int _0^{2\pi /k} \boldsymbol{\cdot }\;\mathrm{d}a$ denotes phase averaging. Pressure is expanded as

(2.5) \begin{equation} p = p_0 - \rho _0gb + \sum _{n=1}^\infty p^{(n)}, \end{equation}

where $p_0$ denotes the constant background pressure, so that $p_0=p_{\textit{air}}$ in the air and $p_0=p_w$ in the water. As done by Pizzo et al. (Reference Pizzo, Lenain, Rømcke, Ellingsen and Smeltzer2023), the pressure can be solved as

(2.6) \begin{equation} p=\rho _0\left (-gz-\frac {(\dot {x}-c)^2 + \dot {z}^2}{2} + \int _{-\infty }^b (c-U)\varGamma\! \mathcal{J}\text{d}\beta + f(a,t)\right )\!, \end{equation}

where $f(a,t)$ is a constant of integration.

Imposing ordering in (2.4), through third order in the wave slope $\epsilon$ ,

(2.7a) \begin{align} x&= a +U^{(0)}(b)t + \epsilon x_1(c,k;b)\sin (\theta _1)+ \epsilon ^2 U^{(2)}(c,k;b)t+\epsilon ^2 x_2(c,k;b) \sin (\theta _2) \nonumber \\ &\quad + \epsilon ^3 U^{(3)}(c,k;b)t + \epsilon ^3 x_3(b)\sin (\theta _3) + O(\epsilon ^4), \end{align}

(2.7b) \begin{align} z &= b + \epsilon z_1(c,k;b) \cos (\theta _1) + \epsilon ^2z_0^{(2)} + \epsilon ^2z_2(b)\cos (\theta _2)+ \epsilon ^3z_0^{(3)} \nonumber \\ &\quad + \epsilon ^3 z_3(b)\cos (\theta _3) +O(\epsilon ^4), \end{align}

(2.7c) \begin{align} p &= p_0 - \rho _0gb + \epsilon p^{(1)} + \epsilon ^2 p^{(2)} + \epsilon ^3 p^{(3)} + O(\epsilon ^4), \end{align}

where the total Lagrangian mean flow has been written as $U(b) = \sum _{n=0}^\infty \epsilon ^n U^{(n)}(b)$ , so $U^{(0)}(b)$ denotes the leading-order portion and $U^{(2)}(b)$ is associated with the Stokes drift. The terms $z_0^{(2)}$ and $z_0^{(3)}$ in $z$ are to ensure $\langle zx_a|_{b=0}\rangle =0$ (note that $z_0^{(0)}, z_0^{(1)}=0$ ). Additionally, the phase speed is expanded as $c=c^{(0)}+\epsilon c^{(1)}+\epsilon ^2c^{(2)}+O(\epsilon ^3)$ . The phase speed is also expanded in this manner by Clamond (Reference Clamond2007), and this expansion makes it possible to find the modified growth rate in § 3.

2.3. Instability growth rate and critical level in Lagrangian coordinates

The linear stability analysis, conducted in the Lagrangian frame, recovers the established Miles instability growth rate. This approach, however, offers a more intuitive understanding of the critical level’s geometry, showing it is a wavy surface that follows the interface, which is less apparent in a traditional Eulerian analysis. To calculate the instability growth rate, we first formulate the Rayleigh equation in Lagrangian coordinates, which requires a change of variables:

(2.8a) \begin{align} x^{(1)}(a,b,t) &= A(a+U^{(0)}(b)t, b, t), \;\; A(a,b,t) = \xi (b)\sin (k(a-c^{(0)}t)) , \end{align}

(2.8b) \begin{align} z^{(1)}(a,b,t) &= B(a+U^{(0)}(b)t, b, t), \;\; B(a,b,t)= \frac {\varphi (b)}{k\big(c^{(0)}-U^{(0)}\big)}\cos (k(a-c^{(0)}t)). \end{align}

The Rayleigh equation depends only on $b$ and is derived from the linear approximations of (2.2a)–(2.2b) and mass continuity (see Appendix A) as

(2.9) \begin{equation} \varphi _{bb}- \left (k^2 - \frac {\frac {{\rm d}^2}{{\rm d}b^2}U^{(0)}}{c^{(0)}-U^{(0)}}\right )\varphi = 0. \end{equation}

We will use the Rayleigh equation after deriving a dispersion relation from the linearised dynamic boundary condition. In Lagrangian coordinates, the linearised dynamic boundary condition at the air–sea interface is given by

(2.10) \begin{equation} p^{(1)}(a,0^+,t)-p^{(1)}(a,0^-,t)=T z_{aa}^{(1)}(a,0,t), \end{equation}

where $T$ is the coefficient of surface tension.

Denote $\tilde {a} \overset {\textrm{def}}{=} a+U_s^{(0)}t$ , where $U_s^{(0)}\overset {\textrm{def}}{=} U^{(0)}(b=0)=U^{(0)}(b=0^+,0^-)$ . In the new variables (2.8b ), using (2.6), (2.10) reduces to

(2.11) \begin{align} p^{(1)}(\tilde {a},0^+,t)-p^{(1)}(\tilde {a},0^-,t)&= -\rho _{\textit{air}}gB(\tilde {a},0^+,t)+\rho _{w}gB(\tilde {a},0^-,t) \nonumber\\ &\quad + \big(c^{(0)}-U^{(0)}_s\big)(\rho _{\textit{air}}B_b(\tilde {a},0^+,t) - \rho _{w}B_b(\tilde {a},0^-,t)). \end{align}

Define the constants

(2.12) \begin{equation} \epsilon _\rho \overset {\textrm{def}}{=} \frac {\rho _{\textit{air}}}{\rho _{\textit{air}}+\rho _w}, \quad \gamma \overset {\textrm{def}}{=} \frac {T}{\rho _{\textit{air}}+\rho _w} \quad \text{and} \quad g'\overset {\textrm{def}}{=} g\frac {\rho _w-\rho _{\textit{air}}}{\rho _{\textit{air}}+\rho _w}, \end{equation}

and the functions

(2.13a) \begin{align} \varXi _{\textit{air}}\big(c^{(0)},k\big) &\overset {\textrm{def}}{=} -\frac {\varphi _b\left(c^{(0)},k;0^+\right)}{\varphi \left(c^{(0)},k;0\right)}, \quad \varXi _{w}\big(c^{(0)},k\big) \overset {\textrm{def}}{=} \frac {\varphi _b\left(c^{(0)},k;0^-\right)}{\varphi \left(c^{(0)},k;0\right)}\quad \text{ and} \end{align}

(2.13b) \begin{align} S&\overset {\textrm{def}}{=} (1-\epsilon _\rho )\frac {\rm d}{{\rm d}b}U^{(0)}(0^-)-\epsilon _\rho \frac {\rm d}{{\rm d}b}U^{(0)}(0^+). \end{align}

After substituting (2.10) into (2.11) and some algebra, we obtain the dispersion relation

(2.14) \begin{equation} (\epsilon _\rho \varXi _{\textit{air}} + (1-\epsilon _\rho )\varXi _w)\big(c^{(0)}-U_s^{(0)}\big)^2+S\big(c^{(0)}-U_s^{(0)}\big)-g'-\gamma k^2 =0. \end{equation}

As done by Young & Wolfe (Reference Young and Wolfe2014), we rewrite (2.14) as

(2.15) \begin{equation} \mathcal{D}_0(c^{(0)},k) + \epsilon _\rho \mathcal{D}_1\big(c^{(0)},k\big) = 0, \end{equation}

where

(2.16) \begin{align} \mathcal{D}_0(c^{(0)},k) &\overset {\textrm{def}}{=} \varXi _w\big(c^{(0)},k\big)\big(c^{(0)}-U_s^{(0)}\big)^2 + S_0\big(c^{(0)}-U_s^{(0)}\big)-g-\gamma k^2, \nonumber\\ \mathcal{D}_1\big(c^{(0)},k\big)&\overset {\textrm{def}}{=} \left(\varXi _{\textit{air}}\big(c^{(0)},k\big)-\varXi _w\big(c^{(0)},k\big)\right)\big(c^{(0)}-U_s^{(0)}\big)^2 + S_1\big(c^{(0)}-U_s^{(0)}\big) + 2g, \nonumber\\ S_0 &\overset {\textrm{def}}{=} \frac {\rm d}{{\rm d}b}U^{(0)}(0^-), \quad S_1\overset {\textrm{def}}{=}-\frac {\rm d}{{\rm d}b}U^{(0)}(0^-)-\frac {\rm d}{{\rm d}b}U^{(0)}(0^+). \end{align}

The solution to (2.15) will be of the form $c^{(0)}(k,\epsilon _\rho )$ , so we expand

(2.17) \begin{equation} c^{(0)}(k,\epsilon _\rho )=c_0^{(0)}(k) + \epsilon _\rho c_1^{(0)}(k) +O\left(\epsilon _\rho ^2\right)\!. \end{equation}

We then obtain the leading-order balance

(2.18) \begin{equation} \mathcal{D}_0\big(c_0^{(0)}(k),k\big)=0, \end{equation}

which only involves flow in the water, so does not itself affect the Miles instability. The Miles instability results from a critical level in the air, $b_c\gt 0$ , such that $c_0^{(0)}(k) = U^{(0)}(b_c).$ At next order in $\epsilon _\rho$ , using a Taylor expansion in terms of $c^{(0)}$ , we can solve

(2.19) \begin{equation} c_1^{(0)} = -\frac {\mathcal{D}_1\big(c_0^{(0)},k\big)}{\partial _{c^{(0)}}\mathcal{D}_0\big(c^{(0)},k\big)}. \end{equation}

The instability growth rate is that of the unstable solutions to the eigenvalue problem. Notice the Rayleigh equation (2.9) becomes singular when $U^{(0)}=c^{(0)}$ ; to avoid this singularity, we introduce a complex phase speed

(2.20) \begin{equation} c^{(0)}=c_r^{(0)}+ic_i^{(0)}, \end{equation}

where $U^{(0)}(b_c)=c_r^{(0)}$ and $c_i^{(0)} \ll c_r^{(0)}$ . The form of the perturbation in (2.8b ) grows in time according to $kc_i^{(0)}$ , so that the instability growth rate (in this linear analysis) is $\omega = k \text{Im } c^{(0)}$ . Noting that $c_0^{(0)}\in \mathbb{R}$ , we approximate the growth rate $\omega$ as

(2.21) \begin{equation} \omega = \epsilon _\rho k\text{Im } c_1^{(0)}. \end{equation}

We then compute that

(2.22) \begin{equation} \text{Im } c_1^{(0)} = - \frac {\text{Im }\mathcal{D}_1}{\text{Re } \partial _{c^{(0)}}\mathcal{D}_0}, \end{equation}

where (see Appendix A)

(2.23) \begin{equation} \text{Im } \mathcal{D}_1 = (c_0^{(0)}-U_s^{(0)})^2\pi \frac {\frac {{\rm d}^2}{{\rm d}b^2}U^{(0)}_c}{\Big |\frac {\rm d}{{\rm d}b}U^{(0)}_c\Big |}\frac {|\varphi _c|^2}{|\varphi _s|^2}. \end{equation}

Here, $U_c^{(0)}\overset {\textrm{def}}{=} U^{(0)}(b_c)$ , $\varphi _c\overset {\textrm{def}}{=} \varphi (b_c)$ and $\varphi _s\overset {\textrm{def}}{=} \varphi (0)$ .

By inverting $x=x(a,b)$ and $z=z(a,b)$ iteratively at each order of $\epsilon$ , using (2.7a ) and (2.7b ), we can write $a$ and $b$ in terms of $x$ and $z$ :

(2.24a) \begin{align} a(x,z) &= x-U^{(0)}(z)t-\epsilon \xi (z)\sin \big(k\big(x-c^{(0)}t\big)\big), \end{align}

(2.24b) \begin{align} b(x,z) &= z-\epsilon \frac {\varphi (z)}{k(c^{(0)}-U^{(0)}(z))}\cos \big(k\big(x-c^{(0)}t\big)\big). \end{align}

We thus recover the known growth rate for the Miles instability, but in the Lagrangian frame; $U_c^{(0)}$ and $\varphi _c$ are functions of $b$ rather than $z$ , but to leading order $U^{(0)}(b)=U^{(0)}(z)$ and $\varphi (b)=\varphi (z)$ , per (2.24a ) and (2.24b ). However, the representation of the growth rate (2.21) given by (2.22) and (2.23) differs from the Eulerian representation because the arguments of $U^{(0)}$ and $\varphi$ should be evaluated at the perturbed $z$ -coordinate described by the right-hand side of (2.24b ), rather than simply $z$ . The critical level, defined by a constant $b$ (or equivalently, perturbed $z$ ), is therefore ‘wavy’, following the shape of the interface. Indeed, (2.24b ) shows that a line of constant $b$ is not one of constant physical height $z$ , but rather a wavy surface oscillating in $x$ (figure 2). Note that at the linear order considered here, the waviness of the critical level is a geometric consequence of the mapping (2.24b ) and does not introduce an additional phase offset beyond the classical Miles mechanism, as reflected in the recovery of the standard linear growth rate. The linear growth rate varies across the interface, since $\varphi _c$ and $U_c^{(0)}$ in (2.23) occur at a constant level $b_c$ , which is not a constant height $z_c$ per (2.24b ). The waviness of the critical layer, stated by Lighthill (Reference Lighthill1962) and apparent through the use of orthogonal coordinates by Benjamin (Reference Benjamin1959), emphasises that the critical height $z_c$ is better regarded as height from the actual, rather than mean, water surface.

This explicit description of the geometry of the critical level is one advantage of the Lagrangian analysis. A more significant outcome of this analysis, however, is the Lagrangian frame’s capacity to show how the wave-induced mean flow modifies the instability growth rate, a task that requires extending the stability analysis to higher order.

3.1. Nonlinear stability analysis

Computing the modification to the growth rate by the wave-induced current requires considering the system through third order. This is similar to other nonlinear stability analyses, e.g. Nayfeh & Saric (Reference Nayfeh and Saric1972), but in the Lagrangian frame, we expand in terms of monochromatic waves (2.4) and expand the phase speed and mean velocity (e.g. Clamond Reference Clamond2007) and not the independent variables $a$ , $b$ or $t$ . In the higher-order analysis, we first consider the special case in which $\gamma =0$ in (2.12), as in Miles’ original analysis. The full details of the higher-order analysis are provided in Appendix B. While it turns out that $c^{(1)}=0$ due to the second-order conservation of vorticity, $c^{(2)}$ is described by the dispersion relation

(3.1) \begin{equation} \rho _{\textit{air}}\big(F(0^+)+c^{(2)}G(0^+)\big)-\rho _w\big(F(0^-)+c^{(2)}G(0^-)\big)=0, \end{equation}

where

(3.2a) \begin{align} F(b) &\overset {\textrm{def}}{=} \varphi '(b) \left (\frac {2}{k}\big (U^{(2)}(b) + k\big(c^{(0)}-U_s^{(0)}\big)x_2(b)\big ) \right )- 2k\big(c^{(0)}-U_s^{(0)}\big)\varphi _s z_2(b) \nonumber \\ &\quad + \varphi _s\left ( \frac {1}{k}\left (\frac {\rm d}{{\rm d}b}U^{(0)}(b)\left ( \frac {3U^{(2)}(b)}{c^{(0)}-U_s^{(0)}}+ 2kx_2(b)\right ) + \frac {\rm d}{{\rm d}b}U^{(2)}(b)\right )\right )\!, \end{align}

(3.2b) \begin{align} G(b) &\overset {\textrm{def}}{=} -\frac {1}{k\big(c^{(0)}-U_s^{(0)}\big)} \left ( 2\varphi '(b)\big(c^{(0)}-U_s^{(0)}\big) + 3\varphi _s\frac {\rm d}{{\rm d}b}U^{(0)}(b) \right )\!. \end{align}

In contrast to our approach to obtain the linear growth rate, since the dispersion relation (3.1) is only linear in $c^{(2)}$ , we obtain a relatively simple analytical solution without using an asymptotic expansion in $\epsilon _\rho$ :

(3.3) \begin{equation} c^{(2)} = -\frac {\rho _{\textit{air}}F(0^+)-\rho _wF(0^-)}{\rho _{\textit{air}}G(0^+)-\rho _wG(0^-)}. \end{equation}

Additionally, under the mild assumption that derivatives of $x_n, z_n, n\geq 1$ introduce a constant coefficient of $n$ (as in a Gerstner or Stokes wave), an analytical formula for $x_2$ and $z_2$ in terms of the known lower order quantities $c^{(0)}$ , $U^{(0)}$ and $z_1$ is possible,

(3.4) \begin{equation} kx_2=z_2^{\prime} \quad \text{ and }\quad 4z_2^{\prime}+z_1^{\prime\prime}z_1-(z_1^{\prime})^2=0. \end{equation}

The leading-order wave-induced mean flow $U^{(2)}$ can also be expressed in terms of known lower-order quantities (see Appendix E),

(3.5) \begin{equation} U^{(2)}=\frac {1}{4}\big(c^{(0)}-U^{(0)}\big)\frac {{\rm d}^2}{{\rm d}b^2}\left (\frac {\varphi ^2}{k^2\big(c^{(0)}-U^{(0)}\big)^2}\right )\!. \end{equation}

Thus, (3.3) is entirely in terms of known quantities, since we can also use the linear growth rate when considering $c^{(0)}$ , (2.21). However, it may not be computable in general, as $x_2$ may be an indefinite integral that cannot be computed. Nonetheless, (3.3) readily gives an exact solution for the modified growth rate for many relevant background wind profiles $U^{(0)}$ . Regardless of the ability to compute an exact solution, the dispersion relation clearly indicates that the leading-order growth-rate modification will be directly affected by the wave-induced mean flow, due to the appearance of $U^{(2)}$ throughout, as well as its shear.

3.2. Growth rates with respect to the double exponential profile

To illustrate the impact of the modified growth rate, we consider an explicit base-state velocity profile

(3.6) \begin{equation} U^{(0)}(b) = \begin{cases} U_\infty ^{(0)} - (U_\infty ^{(0)} - U_s^{(0)})e^{-b/h_{\textit{air}}} &\text{ if } b\gt 0, \\[5pt] U_s^{(0)}e^{b/h_w} &\text{ if } b\lt 0.\end{cases} \end{equation}

This double exponential profile is defined by the scale heights in the air and water ( $h_{\textit{air}}$ and $h_w$ ), and the velocities at the top of the boundary layer and at the surface ( $U_\infty ^{(0)}$ and $U_s^{(0)}$ ). For a given profile, the Rayleigh equation has an exact solution

(3.7) \begin{equation} \varphi (c,k;b) = e^{-k|b|} \frac {F\left(\alpha _j, \beta _j, 2\kappa _j+1, X_j e^{-|b|/h_j}\right)}{F(\alpha _j, \beta _j, 2\kappa _j+1, X_j)}, \end{equation}

where the parameters to the Gaussian hypergeometric function $F(a,b,c,\xi )$ are defined as

(3.8) \begin{equation} (\alpha _j, \beta _j, \kappa _j, X_j) = \begin{cases} \left (\kappa _{\textit{air}}-\sqrt {1+\kappa _{\textit{air}}^2}, \kappa _{\textit{air}}+\sqrt {1+\kappa _{\textit{air}}^2}, kh_{\textit{air}}, \frac {U_\infty ^{(0)}-U_s^{(0)}}{U_\infty ^{(0)} -c^{(0)}}\right ) & \text{for } j=air, \\[1.2em] \left (\kappa _w-\sqrt {1+\kappa _w^2}, \kappa _w+\sqrt {1+\kappa _w^2}, kh_w, \frac {U_s^{(0)}}{c^{(0)}}\right ) & \text{for } j=w. \end{cases} \end{equation}

The amplitude corresponding to a given background shear profile therefore depends on both the wavenumber $k$ and the parameters in the definition of the background profile.

The Gaussian hypergeometric function in (3.7) is given by

(3.9) \begin{equation} F(y_1, y_2, y_3; \xi )\overset {\textrm{def}}{=} 1 + \frac {y_1y_2}{y_3}\frac {\xi }{1!} + \frac {y_1(y_1+1)y_2(y_2+1)}{y_3(y_3+1)}\frac {\xi ^2}{2!}+\ldots \end{equation}

In the case of this double exponential profile (3.6) (figure 1), $\varphi _c, \varphi _s$ , $\varphi '(0^+)$ and $\varphi ''(0^+)$ may be computed exactly using (3.7), the differentiation identity for hypergeometric functions, and the Rayleigh equation, respectively. In addition, if $U_s^{(0)}=0$ , as in the original case considered by Miles (Reference Miles1957), (3.3) may be explicitly, exactly computed, whereas more complex background profiles may necessitate approximation via numerical methods. Assuming $\gamma =0$ ((2.12), no surface tension) also as in the original work, the growth rate (2.21) is evaluated as

(3.10) \begin{equation} \omega ^{(0)}(k) =\epsilon _\rho k \text{Im }c^{(0)}_1=-\epsilon _\rho \left (\frac {g}{k}\right )^{\frac {1}{2}}\frac {\pi }{2} \frac {\frac {{\rm d}^2}{{\rm d}b^2}U_c^{(0)}}{\Big |\frac {\rm d}{{\rm d}b} U_c^{(0)}\Big |}\frac {|\varphi _c|^2}{|\varphi _s|^2}. \end{equation}

Figure 1. The set-up of the Miles instability in Lagrangian coordinates, adapted from Young & Wolfe (Reference Young and Wolfe2014). The Lagrangian mean velocity profile, $\overline {u}_L(z)$ (identified as $U$ in (2.4)), is shown alongside the Eulerian mean, $\overline {u}_E(z)$ . While the two profiles are identical at leading order, they differ at higher order due to the surface wave; note that $\overline {u}_L$ is more naturally a function of the particle label $b$ , but is shown in terms of $z$ for comparison. Example trajectories (dashed), with the same initial $x$ -coordinate, are in the water. As the Lagrangian mean velocity approaches zero, the orbits of particles become closed circles.

Figure 2. The critical layer defined by $b=b_c$ , the surface defined by $b=0$ and intermediate lines of constant $b$ , shown in Eulerian coordinates (the $x{-}z$ plane) by using the transformations (2.24a ) and (2.24b ). These levels of constant $b$ are also shown in the Lagrangian frame (the $a{-}b$ plane) in the inset, illustrating the simple representation of the wavy critical level when considering the instability in Lagrangian coordinates. The relative wind is shown on the left. Here, the background wind is the double exponential profile described in § 3.2 with $U_s^{(0)}=0$ (the original case Miles considered, together with $\gamma =0$ ), $\rho _{\textit{air}} = 1.25 \ \mathrm{kg \ m^{-3}}$ , $\rho _w = 1025 \ \mathrm{kg \ m^{-3}}$ , $h_{\textit{air}} = 1 \ \mathrm{m}$ and $h_w = 0.54 \ \mathrm{cm}$ . In this case, $c^{(0)}_r =(g/k)^{{1}/{2}}$ , where a wavelength of $k=2\pi / 20$ and a wave slope of $\epsilon =0.1$ were selected (arbitrarily) for visualisation. This picture would evolve in time according to (2.24b ); here, $t=0$ is shown.

In total, the modified growth rate is given by (3.10) and (3.3) via

(3.11) \begin{equation} \omega = k\big(\epsilon _\rho \text{Im }c^{(0)}_1 + \epsilon ^2\text{Im }c^{(2)}\big), \end{equation}

which is plotted for different values of $\epsilon$ in figure 3. Note that the $\epsilon _\rho$ and the subscript 1 are due to the asymptotic expansion done to find $c^{(0)}$ , which is unnecessary when computing $c^{(2)}$ . The formula for $c^{(2)}$ (3.3) also incorporates the relative densities. Consistent with linear theory (figure 3, solid curves), the growth rate peaks at an intermediate wavenumber and increases with background shear (modified via increasing $U_\infty ^{(0)}$ ). For a finite but non-zero wave slope $\epsilon$ , this growth is systematically suppressed, an effect that intensifies with increasing $\epsilon$ . For moderate values of $U_\infty ^{(0)}$ , this suppression is largely confined to high wavenumbers, but at larger $U_\infty ^{(0)}$ values, it extends even to intermediate wavenumbers, reducing the peak growth rate by as much as a factor of four for the largest wave slope. The aggregate modulation, obtained by trapezoidal integration of the growth curves shown in figure 3, exhibits a maximum growth rate decrease (among the parameters shown) of $39.85\,\%$ for the $\epsilon =0.3$ , $U_\infty =10$ case and a minimum decrease of $0.30\%$ for the $\epsilon =0.1$ , $U_\infty =5.0$ case.

Figure 3. Modified growth rate as a function of wavenumber $k$ , for four wave slopes $\epsilon$ and three values of $U_\infty ^{(0)}$ (which increases the background shear). The linear growth rate ( $\epsilon =0$ , solid) computed in the Lagrangian frame equals the known result computed in the Eulerian frame, while the modified growth rates, dependent on the leading-order wave-induced mean flow ( $\epsilon \neq 0$ ), show that increased wave slope combined with increased background shear can lead to a significant suppression of the instability at high wavenumbers. For larger values of $U_\infty ^{(0)}$ , the modification can also significantly reduce the peak growth rate. Growth rates were computed for the double-exponential profile with parameters: $U^{(0)}_s=0$ , $\rho _{\textit{air}} = 1.25 \ \mathrm{kg \ m^{-3}}$ , $\rho _w = 1025 \ \mathrm{kg \ m^{-3}}$ , $h_{\textit{air}} = 1 \ \mathrm{m}$ and $h_w = 0.54 \ \mathrm{cm}$ (as in Young & Wolfe Reference Young and Wolfe2014). The values of $\epsilon$ chosen are based on observational data, e.g. those of Peirson & Garcia (Reference Peirson and Garcia2008). The linear growth rate shown here is the asymptotic solution, but it exactly aligns with the numerical solution to (2.14) (cf. Young & Wolfe Reference Young and Wolfe2014).

Up to this point, we have considered the same system as Miles (Reference Miles1957) and neglected the effects of surface tension. Since the predicted impacts of the higher-order theory are largest for shorter waves (figure 3), which are likely impacted by surface tension, we also compute the impacts of surface tension on both $c^{(0)}$ and $c^{(2)}$ . The inclusion of surface tension does not change that $c^{(1)}=0$ . The full calculation of these effects involves computing the curvature at points on the surface in Lagrangian coordinates by considering it as a parametric curve $\boldsymbol{r}(a)=(x(a;b=0,t),z(a;b=0,t))$ and expanding asymptotically, and is included in Appendix C. A typical value of the kinematic surface tension $\gamma$ is $7.2\times 10^{-5} \;\text{m}^3\;\text{s}^{-2}$ and the modified growth rate computed with this value is shown in figure 4. While the intermediate wavenumbers affected by the modified growth rate (i.e. the inclusion of the wave-induced current) have the same growth whether surface tension is included or not, the highest wavenumbers ( $k\gt 30\; \text{m}^{-1}$ ) are significantly impacted by capillary effects. The impacts of surface tension, like those of the inclusion of the wave-induced mean flow, increase with wave slope. In all, the stabilising effect of surface tension for short waves is amplified by the nonlinear modification of the growth rate shown in figure 3.

Figure 4. Modified growth rate with capillarity as a function of wavenumber $k$ , for four wave slopes $\epsilon$ . The grey curves are as in figure 3 ( $\gamma =0$ ) whereas the black curves show the impact of including surface tension effects ( $\gamma \neq 0$ , here $\gamma =7.2\times 10^{-5} \;\text{m}^3\;\text{s}^{-2}$ ). The parameters are as in figure 3, but all growth rate curves are with respect to $U_\infty ^{(0)}=10.0$ m s $^{-1}$ . In addition to the modification from the wave-induced current, surface tension further suppresses growth at high wavenumbers. However, surface tension does not alter the effect of the modification to growth rate shown in figure 3 at intermediate wavenumbers.

3.3. Growth rates with respect to a logarithmic profile

While the double exponential profile conveniently permits analytical solutions, and thus more accurate solutions for the growth rate, it is representative of laminar flow; in contrast, a background logarithmic profile approximates a (mean) turbulent wind velocity, providing a realistic model better suited for comparing growth rates against empirical data (Morland & Saffman Reference Morland and Saffman1993). To that end, we extend Miles’ (Reference Miles1957) method for approximating the growth rate given a background logarithmic profile

(3.12) \begin{equation} U^{(0)}(b) = \begin{cases} U_0 + U_1 \log \left (\dfrac {b}{b_0}\right ) &\text{if } b\geq b_0, \\ 0 &\text{if } b\lt 0, \end{cases} \end{equation}

to our Lagrangian, higher-order analysis. In (3.12), $U_1$ is a reference wind speed, $b_0$ is a label near $b=0$ (analogous to the near-surface streamline in Miles’ analysis) and $U_0\overset {\textrm{def}}{=} U^{(0)}(b_0)$ (see Appendix F). We explicitly leave the profile undefined in the viscous/transition layer $0\leq b \lt b_0$ ; as done by Miles (Reference Miles1957), we assume this layer ‘moves with the surface wave’, and treat $b= b_0$ as the effective air-side interface. We then evaluate air-side quantities at $b=b_0$ and water-side quantities at $b=0^-$ , in the growth rate approximation.

It is often preferable, especially when comparing against empirical data, to rewrite (3.12) in terms of an aerodynamic roughness length $\tilde {b}_0$ ,

(3.13) \begin{equation} U^{(0)}(b)= \begin{cases} U_1\log \left (\dfrac {b}{\tilde {b}_0}\right ) &\text{ if } b\geq \tilde b_0, \\ 0 &\text{ if } b\lt 0.\end{cases} \end{equation}

The representations (3.12) and (3.13) are equivalent provided $U_0 = U_1\log ( b_0/\tilde {b}_0).$ The near-surface label $b_0$ may then be chosen in relation to the aerodynamic roughness $\tilde {b}_0$ , e.g. $b_0=30\tilde {b}_0$ as in the classical case (Miles Reference Miles1957, following (7.5a,b)).

To obtain the leading-order growth rate (in Eulerian coordinates), Miles (Reference Miles1957) assumes a particular form of the pressure that is proportional to a dimensionless coefficient, i.e. $p\propto \alpha + i\beta$ , where $\alpha =\alpha (k,c)$ and $\beta =\beta (k,c)$ . For consistency with that presentation, we will then find the higher-order expression for $\beta$ , which is a normalised growth rate proportional, but not equal, to $\omega$ . For instance, $\beta$ equals $\text{Im } c_1^{(0)}$ in (3.10) (Miles Reference Miles1957, (3.3)). It is also assumed that the real part of the phase speed is not significantly altered by the process of interest, so that

(3.14) \begin{equation} c \approx c_r^{(0)}\left (1+\frac {\rho _{\textit{air}}}{\rho _w}(\alpha +i\beta )\left (\frac {U_1}{c_r^{(0)}}\right )^2\right )\!, \end{equation}

with $c_i^{(0)}$ and $c_r^{(0)}$ again as in (2.20). Miles then finds that

(3.15) \begin{equation} \beta (c,U_1)\approx \frac {2c_r^{(0)}}{U_1^2}\frac {\rho _w}{\rho _{\textit{air}}}\text{Im }c. \end{equation}

When extending (3.15) to account for not only $c^{(0)}$ but also $c^{(2)}$ , to remain consistent with the original approximation, one must still assume that

(3.16) \begin{equation} \text{Re } c\approx c_r^{(0)} = \left (\frac {g}{k}\right )^{\frac {1}{2}}, \end{equation}

and thus simply substitute in $\text{Im}(c^{(0)}+\epsilon ^2 c^{(2)})$ for $\text{Im }c$ in (3.15). Substituting (3.12) into all relevant quantities in the definition of $c^{(2)}$ (3.3) yields

(3.17) \begin{equation} c^{(2)} = \frac {2c^{(0)}-4U_0(1-\epsilon _\rho )}{c^{(0)}\left (2\big(c^{(0)}-U^{(0)}\big)-\epsilon _\rho \frac {U_1}{kb_0}\right )} \end{equation}

(see Appendix F). To find the imaginary part of (3.17), we expand $c^{(2)}$ in terms of $c^{(0)}_i$ , akin to (2.17), which yields a simple form for $\text{Im } c^{(2)}$ ,

(3.18) \begin{equation} \text{Im }c^{(2)}=c_i^{(0)} \left (\frac {\partial c^{(2)}}{\partial c^{(0)}}\right )\Bigg |_{c^{(0)}=c_r^{(0)}}+ O\left ((c_i^{(0)})^2\right )\!. \end{equation}

Evaluating (3.18) per (3.17) gives the approximate modified growth rate

(3.19) \begin{equation} \beta \approx \beta _0\left (\!1-\epsilon ^2\left(c^{(0)}_r-U_0\right)^2\frac { 4\left ( 2\epsilon _\rho \left(c^{(0)}_r\right)^2 +(1-\epsilon _\rho )\left(2\left(c^{(0)}_r-U_0\right)^2 +\epsilon _\rho U_0\frac {U_1}{k b_0}\right)\! \right ) }{ \left(c^{(0)}_r\right)^2 \left (2\left(c^{(0)}_r-U_0\right)-\epsilon _\rho \frac {U_1}{k b_0}\right )^2} \!\right )\!, \end{equation}

with $\epsilon _\rho$ as in (2.12). In (3.19), $\beta _0$ is the leading-order approximate growth rate, obtained by translating the calculation of $\beta (c^{(0)}, U_1)$ from Miles (Reference Miles1957) into Lagrangian coordinates,

(3.20) \begin{equation} \beta _0 = \pi (kb_c)^{3}\left ( \int _{1}^{\infty }e^{-kb_c r}(\log r)^2\,{\rm d}r\right )^2. \end{equation}

In (3.20), one may calculate from (3.12) that the critical level is $b_c=b_0e^{c^{(0)}_r/U_1}$ .

Figure 5 shows $\beta$ according to (3.19) for a range of wave steepness values. Note that, to reproduce figure 1 of Miles (Reference Miles1957), the $y$ -axis is not on a logarithmic scale (unlike figures 3 and 4) and recall that $\beta$ is proportional, rather than equal, to $kc_i$ . Nonetheless, $\beta$ is clearly suppressed as $\epsilon$ increases, confirming that the growth rate reduction with steepness also holds for the logarithmic profile. Other features also persist: increasing $\epsilon$ lowers the peak value of $\beta$ , shifts the wavenumber corresponding to that maximum and narrows the instability window such that it terminates at progressively smaller $\xi$ . We describe these shifts in terms of $\xi$ rather than $k$ because $\xi$ is not necessarily monotone in $k$ ; it depends on $b_c$ , which itself is a function of $k$ through $U^{(0)}(b_c)=c^{(0)}_r(k)$ . Although the general suppressive trend is the same for both background profiles, the relative deformation of the curves (peak shifts and cutoff shifts) differs, as expected due to the different near-surface structure of the profiles. In the logarithmic case, shear and curvature are concentrated near the cutoff at $b_0$ (scaling with $1/b$ and $1/b^2$ , respectively), so the critical-layer contribution is sensitive to how close $b_c$ is to $b_0$ . Since the corrections due to the wave-induced mean flow are intensified near the surface, waves whose corresponding critical layer lies closest to $b_0$ are preferentially suppressed. In contrast, the double exponential profile has bounded near-surface shear and curvature, which decay exponentially away from the interface on the scale $h_{\textit{air}}$ ; consequently, this decay scale controls the steepness-dependent shifts.

Figure 5. Normalised growth rate $\beta$ as a function of the dimensionless critical level $kb_c$ . All curves are calculated per (3.19). The solid line ( $\epsilon =0$ ) is an exact reproduction of figure 1 of Miles (Reference Miles1957). The dashed, dash-dotted and dotted curves show the (modified) $\beta$ value corresponding to five non-zero wave steepness values $\epsilon$ . These calculations use a representative wind speed of $U_1=1\,\text{m s}^{-1}$ and a relatively smooth surface roughness of $\tilde {b}_0=10^{-4} \text{ m}$ , consistent with the initial generation of waves on a quiescent sea. As done by Miles (Reference Miles1957), the near-surface label $b_0$ is related to the aerodynamic roughness $\tilde {b}_0$ via $b_0=30\tilde {b}_0$ .

3.4. Comparison to experimental data

Importantly, (3.19) enables a direct comparison of the theoretical correction to empirical data, offering a new explanation for the observed steepness dependence of growth rates. Peirson & Garcia (Reference Peirson and Garcia2008) systematically investigate the significance of wave steepness for wave growth by combining new wave-tank experiments with several reanalysed datasets (1967–1998). Their formula for $\beta$ ((7) therein) is given by

(3.21) \begin{equation} \beta = \frac {\rho _w}{\rho _{\textit{air}}}\frac {S_{in}}{\omega E}\left (\frac {c}{u_*}\right )^2, \end{equation}

where $S_{in}$ is the energy input from the wind, $E$ is the local total energy density and $u_*$ is the friction velocity. While Miles has $U_1$ in place of $u_*$ , the logarithmic profile implies $U_1=u_*/\kappa$ (Miles Reference Miles1957, (5.3b) therein), where $\kappa =0.41$ is the von Kármán constant; we therefore rescale $\beta$ calculated via (3.19) by $\kappa ^{-2}$ .

Because the experiments analysed by Peirson & Garcia (Reference Peirson and Garcia2008) originate from distinct studies, the background parameters span a broad range: friction velocities $u_*$ vary from roughly 0.2 m s $^{-1}$ to over 1.0 m s $^{-1}$ and the intrinsic frequencies $f$ of the mechanically generated waves vary from approximately 0.5 to 6.0 Hz. Consequently, the set of empirical data could be considered as lying on a family of curves rather than a single curve, as (3.19) depends on these different parameters. To assess the general steepness trend rather than predicting individual data points (which can also be subject to large error bars), figure 6 evaluates (3.19) using fixed background parameters within the experimental range: aerodynamic roughness $\tilde {b}_0=10^{-4} \text{ m}$ , friction velocity $u_*=0.3 \text{ m s}^{-1}$ and $f=2.0 \text{ Hz}$ . The sensitivity of the trend to these specific choices is examined in figure 7.

Figure 6. Comparison of the original Miles $\beta$ , the modified $\beta$ (3.19) and the experimental compilation of Peirson & Garcia (Reference Peirson and Garcia2008) (their figure 6), against wave steepness $\epsilon$ . Symbols and error bars are digitised from Peirson & Garcia (Reference Peirson and Garcia2008), with each marker type denoting a different set of experiments. The Miles baseline (solid curve) and the modified growth parameter (thick dashed curve) are evaluated at fixed background parameters $f=2.0 \text{ Hz}, u_* = 0.3 \text{ m s}^{-1}$ and $\tilde {b}_0=10^{-4}\text{ m}$ (as in figure 5, $b_0=30\tilde {b}_0$ ). Dash-dotted and dotted curves show the reference trends of Belcher (Reference Belcher1999) and Longuet-Higgins (Reference Longuet-Higgins1969), respectively, as discussed by Peirson & Garcia (Reference Peirson and Garcia2008).

A comparison of the theoretical curve to experimental data demonstrates that the correction due to the wave-induced mean flow (3.19) again yields decreasing wave growth with wave steepness, which captures the qualitative steepness dependence of the data and closely matches the measurements in magnitude over much of the steepness range (figure 6). Although the growth rate is often within the reported error bars at both high and low wave steepness values in figure 6, the fit between theory and experimental data could be more robust if the growth rate were calculated with respect to the background parameters used in each individual experiment, rather than a fixed set. However, not all experiments were reported with specific values of each parameter (as described by Peirson & Garcia (Reference Peirson and Garcia2008) and references therein). Even though the prediction shown in figure 6 is thus not fitted to individual data points, the growth rate modification shows a remarkable correspondence with the observed trend of growth suppression.

Figure 7. Sensitivity of the $\beta (\epsilon )$ trend calculated according to (3.19) and displayed in figure 6 to the background parameters in the logarithmic profile, with the corresponding Miles $\beta$ for each set of parameters and a linear fit to the Peirson & Garcia (Reference Peirson and Garcia2008) data plotted for reference. In each panel, the Miles (Reference Miles1957) baseline (solid, light grey) and corresponding modified prediction (black dashed) are compared with the linear fit to the mean trend reported by Peirson & Garcia (Reference Peirson and Garcia2008) (long dashes; shaded bands show one standard deviation of variation calculated from the data, on either side). (a) Dependence on aerodynamic roughness parameter $\tilde {b}_0$ for fixed $f=2.0 \text{ Hz}$ and $u_* = 0.30 \text{ m s}^{-1}$ ( $\tilde {b}_0$ values are as marked in the legend). (b) Dependence on friction velocity $u_*$ for fixed $f=2.0 \text{ Hz}$ and $\tilde {b}_0=10^{-4}\text{ m}$ ( $u_*$ values are as marked in the legend). (c) Dependence on wave frequency $f$ for fixed $u_*=0.3 \text{ m s}^{-1}$ and $\tilde {b}_0=10^{-4}\text{ m}$ ( $f$ values are as marked in the legend).

Moreover, a calculation of the (approximate) modified growth rate across a range of realistic parameter values (figure 7) demonstrates that in all cases, the modification yields a strict suppression relative to the original Miles growth rate evaluated at the same background parameters. While the magnitude of the suppression depends on parameters – and can even result in no growth at very high steepness given certain parameter values (e.g. high $u_*$ or high $f$ ) – the reduction in growth with increasing wave steepness is seemingly universal. That the modified prediction can reproduce not only the trend, but also the magnitude of experimental results reasonably well may be somewhat surprising given that the present calculation does not include any effects of viscous or turbulent losses, unlike the data from Peirson & Garcia (Reference Peirson and Garcia2008). We also note that the sensitivity to background parameters shown in figure 7 is consistent with the mixed performance of the original Miles-type predictions across experiments (e.g. van Gastel et al. Reference van, Klaartje, Peter A.E. and Komen1985 or Wilson et al. Reference Wilson, Banner, Flower, Michael and Wilson1973). Since the modified growth rate equals the original Miles growth rate at $\epsilon =0$ , its correspondence with experimental data is naturally limited by that of the unmodified growth rate. The Miles (Reference Miles1957) prediction itself of course has no explicit dependence on wave steepness and therefore does not reflect the decreasing trend, though it can agree with the magnitude at low steepness.

Despite the obvious dependence on background parameters (figure 7), it is possible to formulate a very simple estimate of the relative change in the growth rate, $|\beta -\beta _0|/\beta _0$ , to assess the practical significance of the wave-slope dependent correction in experimental or observational contexts. This deviation is given by the $\epsilon ^2$ term in (3.19); simplifying with $U_0=0$ and $\epsilon _\rho \to 0$ , its magnitude approaches $2\epsilon ^2$ . Wave steepness values in developing seas typically range from $ak\sim 0.05$ – $0.15$ (e.g. Donelan, Hamilton & Hui Reference Donelan, Hamilton and Hui1985), but can reach $ak\sim 0.3$ – $0.4$ in strong forcing or near-breaking regimes (Perlin, Choi & Tian Reference Perlin, Choi and Tian2013). Applying this rough estimate, a wave steepness of 0.15 yields a $4.5\,\%$ effect and a wave steepness of 0.3 yields an $18\,\%$ effect, which is further amplified by the exponential nature of wave growth.

Peirson & Garcia (Reference Peirson and Garcia2008) interpret the $\beta (\epsilon )$ trend by appealing to different mechanisms in different steepness ranges, as marked in figure 6. For the higher-steepness regime, they compare against Belcher’s framework, though this requires tuning the tangential-stress contribution to match observations. While this representation consequently provides a reasonable match at high steepness, it is a poor representation at low steepness, so they instead attribute strong growth in that regime to the Longuet-Higgins ‘maser-like’ mechanism. In contrast, figures 6 and 7 show that the observed trend can alternatively be accounted for across the steepness range with a single Miles-type mechanism. Although this monotonic suppression is clearly evident across the many experiments compiled by Peirson & Garcia (Reference Peirson and Garcia2008), the underlying trend could be obscured by other factors in other empirical datasets. In studies such as that of Peirson & Garcia (Reference Peirson and Garcia2008), ‘growth’ is inferred from net wave-energy evolution. However, Grare et al. (Reference Grare, Peirson, Branger, Walker, Giovanangeli and Makin2013b ) find wind-input levels that exceed the net growth rates assembled by Peirson & Garcia (Reference Peirson and Garcia2008), implying that substantial turbulent effects cause this residual metric to underestimate true wind input, especially when breaking is frequent. Extrapolation biases in drag estimates, linked to viscous-stress decreases near the interface (Buckley et al. Reference Buckley, Veron and Yousefi2020), could also make this relationship less apparent. Irrespective of the various observational conditions that may dilute the steepness dependence, a specific physical mechanism driving the suppression can be identified, as shown through the integral momentum budget formulated in the next section.

3.5. Momentum budget

Having demonstrated that the growth rate modification due to the wave-induced mean flow consistently acts to suppress the instability across different background profiles, we now seek to physically interpret this effect by formulating a momentum budget. Expressed in the unapproximated Lagrangian coordinates (2.4), this budget emphasises that the critical-layer interaction is not governed by the Eulerian background shear alone, but by the coupling between the total phase speed $c$ and the total Lagrangian mean flow $U$ :

(3.22) \begin{equation} \frac {\partial }{\partial t}\!\int _0^{b_c}\!\!\mathcal{J}\dot {x} \; \text{d}b + \frac {\partial }{\partial a}\!\int _0^{b_c}\!\! pz_b \, \text{d}b = \rho _0 g z \frac {\dot {z}}{c-U}\Bigg \vert _{b=0}^{b=b_c}\! - \frac {1}{2}\rho _0 \big(\dot {x}^2\!+\!\dot {z}^2\big)\dot {z}\Bigg \vert _{b=0}^{b=b_c}\!\!+\!\rho _0 f(b,t)\Big \vert ^{b=b_c}_{b=0}. \end{equation}

In (3.22), $f(b,t)$ is a constant of integration determined by the phase-averaged pressure at the surface. This budget results from integrating the momentum equation and using an alternative form of pressure which can be derived from the momentum equations (see Appendix E),

(3.23) \begin{equation} p = \rho _0 \left (- g z + \frac {1}{2}(c-U)(\dot {x}^2+\dot {z}^2)+f(b,t)\right )\!. \end{equation}

The first term on the left-hand side of (3.22) represents the change in momentum density over time, while the second term on the left-hand side is momentum flux. Physically, the second left-hand term accounts for redistribution of momentum laterally within the layer between the surface and the critical level. The right-hand side accounts for the net momentum input across the boundaries at the surface ( $b=0$ ) and the critical level computed in § 2 $(b=b_c)$ , which is the physical manifestation of form drag. This form drag is split into two components, originating from the hydrostatic and non-hydrostatic components of the pressure (3.23), respectively. The strength of the forcing associated with the first component is proportional to the wave amplitude via $z(b)$ ; however, this is modulated by the key phase shift between wave motion and the pressure field. The second component of the form drag takes the form of a vertical flux of kinetic energy, representing the work done by that pressure component at the boundaries.

Crucially, (3.22) reveals that the critical-layer interaction is governed by the coupling between the total phase speed $c$ and the total Lagrangian mean flow $U$ . The suppression of wave growth with increasing wave slope, reflected in both the modified growth rate (figures 3, 4 and 5) and the experimental comparisons (figures 6 and 7), can likewise be seen as a direct consequence of the wave’s nonlinear feedback on its own growth via (3.22). The dependence of the growth rate on wave steepness emerges from the asymptotic expansion; the well-known $O(\epsilon ^2)$ scaling of the wave-induced mean flow results in a growth rate modification of the same order, providing a causal link between steepness and the suppression. While at linear order the critical layer (defined as $b_c: c^{(0)}(b_c)=U^{(0)}(b_c)$ , as it is not necessarily the case that there is any $b$ such that $c(b)=U(b)$ ) resonance is perfect and the momentum transfer is maximally efficient, the higher-order induced mean flow creates a mismatch between the total phase speed and Lagrangian mean flow ( $c\neq U$ in general). This detunes the resonance and ultimately results in less efficient momentum transfer. Note that this mean-flow mechanism is distinct from a simple alteration of the surface current (as considered in prior studies); it modifies the Lagrangian flow profile without changing the overall shear of the background Eulerian flow. Thus, while the observed growth suppression (Peirson & Garcia Reference Peirson and Garcia2008) has been explained as a shift from a highly efficient, extrinsic ‘maser-like’ mechanism (Longuet-Higgins Reference Longuet-Higgins1969) dominant at low steepness, our analysis shows it can be understood as an inherent consequence of the wave’s own dynamics.