2. Linear stability of a sheared two-fluid interface

Young & Wolfe (Reference Young and Wolfe2014) derived the eigenvalue problem associated with the linear stability of an inviscid parallel shear flow across a two-fluid interface, as sketched in figure 1. First, we outline the governing Rayleigh equation and boundary conditions, after which we describe the non-dimensionalization of the problem.

2.1. Eigenvalue problem

We consider a parallel shear flow, $U=U(z)$, monotonic in both air and water with a non-zero curvature. The air to water density ratio is $r\equiv \rho _{{a}}/\rho _{{w}}<1$, and the gravitational acceleration and surface tension are $g$ and $\sigma$, respectively. Incompressibility ensures that a perturbation of the flow is entirely determined by the streamfunction $\psi =\psi (x,z,t)$, where $t$ is the time. We use normal modes in the form

(2.1)\begin{equation} \psi(x,z,t)={\rm Re}\{ \hat{\psi}(z)\, {\rm e}^{{\rm i}k(x-ct)} \}, \end{equation}

where $k$ is a real wavenumber and $c$ a complex phase speed to be determined (the eigenvalue), conservation of vorticity yields the Rayleigh equation as

(2.2)\begin{equation} [U(z)-c][\hat{\psi}^{\prime\prime}(z)-k^2\hat{\psi}(z)] - U^{\prime\prime}(z) \hat{\psi}(z)= 0. \end{equation}

We require the function $U(z)$ to be continuous at $z=0$, which excludes a Kelvin–Helmholtz type of instability and ensures the continuity of $\hat {\psi }(z)$ (Drazin & Reid Reference Drazin and Reid1981). Note that the derivative $U^{\prime}(z)$ may have a finite jump at $z=0$ and the perturbation must decay in the far field. Provided that

(2.3)\begin{equation} \lim_{z\to\pm\infty} \frac{U^{\prime\prime}(z)}{U(z)-c} =0, \end{equation}

we can impose an exponential decay, viz.

(2.4)\begin{equation} \lim_{z\to\pm\infty} {\psi}^{\prime}(z)\pm k \psi(z) =0. \end{equation}

Finally, we impose continuity of the normal stress at $z=0$ and require the air–water interface to be a streamline, which yields the boundary condition

(2.5)\begin{align} &\frac{\sigma}{\rho_{{w}}} k^2 + (1-r) g + r \left[ (U_{{S}}- c)^2 \frac{\hat{\psi}^{\prime}(0^+)}{\hat{\psi}(0)} - (U_{{S}}- c) U^{\prime}(0^+)\right] \nonumber\\ &\quad - (U_{{S}}- c)^2 \frac{\hat{\psi}^{\prime}(0^-)}{\hat{\psi}(0)} + (U_{{S}}- c) U^{\prime}(0^-)=0, \end{align}

where $U_{{S}}= U(z=0)$ is the surface drift.

2.2. Non-dimensionalization

The length and velocity scales, $L$ and $V$, differ in the air and water, which we denote with the subscripts ${a}$ and ${w}$, respectively. Whereas, when there is no current in the water, we use the scales of the wind, $L_{{a}}$ and $V_{{a}}$, to non-dimensionalize the problem, because water is the primary medium of surface waves we use $L_{{w}}$ and $V_{{w}}$ when a current is present. Thus, the ratios

(2.6a,b)\begin{equation} \mathcal{R}_1 \equiv \frac{ V_{{a}}}{ V_{{w}}}\quad \text{and}\quad \mathcal{R}_2 \equiv \frac{L_{{a}}}{L_{{w}}} \end{equation}

act as additional control parameters. In a frame moving at speed $U_{{S}}$, we use the dimensionless profile

(2.7)\begin{equation} \mathcal{U}\equiv \frac{U-U_{{S}}}{V_{{w}}},\quad \text{where } \mathcal{U}(0)=0, \end{equation}

and the dimensionless intrinsic phase speed is

(2.8)\begin{equation} \mathcal{C}\equiv \frac{c-U_{{S}}}{V_{{w}}}. \end{equation}

The dimensionless wavenumber and vertical coordinate are $\mathcal {k} \equiv kL_{{w}}$ and $\mathcal {z}\equiv z/L_{{w}}$, respectively, and we define the dimensionless gravity and surface tension as

(2.9a,b)\begin{equation} \mathcal{G} \equiv \frac{gL_{{w}}}{V_{{w}}^2}\quad \text{and}\quad \mathcal{S} \equiv \frac{\sigma}{\rho_{{w}} V_{{w}}^2 L_{{w}}}. \end{equation}

In astrophysical contexts, $\sigma$ can represent a magnetic field in the lower fluid, whose direction is aligned with the flow (Alexakis, Young & Rosner Reference Alexakis, Young and Rosner2002). From (2.4), the far-field behaviour of the streamfunction has the form ${\rm e}^{\pm \mathcal {k}\mathcal {z}}$. Short waves are characterized by $\mathcal {k}\gg 1$, so the exponential decay is captured as follows:

(2.10)\begin{equation} \frac{\hat{\psi}(z)}{\hat{\psi}(0)} \equiv \begin{cases} {\rm e}^{-\mathcal{k} \mathcal{z}} f(\mathcal{z}) & \text{if } \mathcal{z}>0,\\ {\rm e}^{\mathcal{k} \mathcal{z}}\ h(\mathcal{z}) & \text{if } \mathcal{z}<0. \end{cases} \end{equation}

We introduce the small parameter $\epsilon \equiv 1/\mathcal {k}$ and use (2.10) in the Rayleigh equation (2.2), which gives

(2.11)$$\begin{gather} \epsilon f^{\prime\prime}(\mathcal{z})- 2 f^{\prime}(\mathcal{z})- \epsilon \frac{\mathcal{U}^{\prime\prime}(\mathcal{z})} { \mathcal{U}(\mathcal{z})-\mathcal{C} } f(\mathcal{z})=0, \quad f(0)=1, \end{gather}$$

(2.12)$$\begin{gather}\text{and}\quad \epsilon h^{\prime\prime}(\mathcal{z})+2 h^{\prime}(\mathcal{z})- \epsilon \frac{\mathcal{U}^{\prime\prime}(\mathcal{z})} { \mathcal{U}(\mathcal{z})-\mathcal{C} } h(\mathcal{z})=0, \quad h(0)=1. \end{gather}$$

We assume that

(2.13a,b)\begin{equation} f(\mathcal{z})=O(1)\text{ as } \mathcal{z}\to+\infty\quad \text{and}\quad h(\mathcal{z})=O(1)\text{ as } \mathcal{z}\to-\infty, \end{equation}

the veracity of which we check a posteriori. The dimensionless form of (2.5) is

(2.14)\begin{align} \frac{\mathcal{S}}{ \epsilon} + \mathcal{G}(1-r)\epsilon + [ r \mathcal{U}^{\prime}(0^+)-\mathcal{U}^{\prime}(0^-)] \epsilon \mathcal{C} - (1+r) \mathcal{C}^2 + [ r f^{\prime}(0^+)-h^{\prime}(0^-)]\epsilon\mathcal{C}^2=0.\end{align}

For a given profile $\mathcal {U}(\mathcal {z})$, our main task is to solve (2.11) and (2.12) subject to the boundary conditions (2.13a,b) and (2.14).

We consider the canonical situation where the wind blows over the water in which it induces a current (figure 1). Hence, $\mathcal {U}^{\prime}(\mathcal {z}\lessgtr 0)>0$, $\mathcal {U}(\mathcal {z}>0)>0$, and $\mathcal {U}(\mathcal {z}<0)<0$. We explore another situation in § 5.

2.3. Examples of profiles

A typical example of a wind and a wind-induced current is the double exponential profile

(2.15)\begin{equation} U(z)= \begin{cases} U_\infty + (U_{{S}}- U_\infty)\, {\rm e}^{{-}z/d_{{a}}} & \text{if } z>0,\\ U_{{S}}\, {\rm e}^{z/d_{{w}}} & \text{if } z<0, \end{cases}\end{equation}

where $U_\infty$ is the free-stream air velocity, and $d_{{a}}$ and $d_{{w}}$ denote the thicknesses of the air and water shear boundary layers, respectively. In the frame of the water surface, the dimensionless form of (2.15) is

(2.16)\begin{equation} \mathcal{U}(\mathcal{z})= \begin{cases} (\mathcal{R}_1- 1)(1- {\rm e}^{-\mathcal{z}/\mathcal{R}_2}) & \text{if } \mathcal{z}>0,\\ {\rm e}^{\mathcal{z}}-1 & \text{if } \mathcal{z}<0, \end{cases} \end{equation}

with $\mathcal {R}_1= U_\infty /U_{{S}}$ and $\mathcal {R}_2= d_{{a}}/d_{{w}}$. Young & Wolfe (Reference Young and Wolfe2014) showed that, for such a profile, the eigenvalue problem can be solved exactly in terms of hypergeometric functions.

In the context of physical oceanography, the double log profile

(2.17)\begin{equation} U(z)= \begin{cases} U_{{S}}+u_{{\star}{a}} \ln(1+z/z_{0{a}})/\kappa & \text{if } z>0,\\ U_{{S}} -u_{{\star}{w}} \ln(1-z/z_{0{w}})/\kappa & \text{if } z<0, \end{cases}\end{equation}

may be more realistic (Wu Reference Wu1975), where $u_{\star {a}}$ and $u_{\star {w}}$ are the friction velocities of air and water, respectively, $\kappa =0.4$ is the von Kármán constant and $z_{0{a}}$ and $z_{0{w}}$ are air and water roughness lengths, accounting for the presence of waves. Note that the velocity of the logarithmic current is negative for $z< z_{{min}}$, where

(2.18)\begin{equation} z_{{min}}\equiv{-}z_{0{w}} ({\rm e}^{\kappa U_{{S}}/u_{{\star}{w}}}-1). \end{equation}

Such a change of sign is unphysical, so we take $U(z)=0$ for $z< z_{{min}}$. The dimensionless form of (2.17) in the frame of the water surface is

(2.19)\begin{equation} \mathcal{U}(\mathcal{z})= \begin{cases} \mathcal{R}_1\ln(1+\mathcal{z}/\mathcal{R}_2)/\kappa & \text{if } \mathcal{z}>0,\\ -\ln(1-\mathcal{z})/\kappa & \text{if } \mathcal{z}<0, \end{cases} \end{equation}

with $\mathcal {R}_1= u_{\star {a}}/u_{\star {w}}$ and $\mathcal {R}_2= z_{0{a}}/z_{0{w}}$. In this case, exact analytical solutions are unknown.

3. Short wavelength expansions and exponential asymptotics

We treat the eigenvalue problem described in § 2.1 perturbatively, where the small parameter is the dimensionless inverse wavenumber, $\epsilon \ll 1$. We draw intuition for the approach from Miles (Reference Miles1957), who considered a simpler version of our problem, with the small parameter $r\ll 1$ and in the absence of a current. Hence, the leading-order eigenvalue corresponding to $r=0$ is real and equal to $c_s$, the phase speed of free surface waves. Moreover, due to the critical layer at height $z_c$, such that $U(z_c)=c_s$, the eigenvalue has an imaginary part of order $r$ as well as real corrections, also of order $r$.

In contrast to the treatment of Miles (Reference Miles1957), the leading-order eigenvalue is unknown. In fact, even the nature of the lowest-order behaviour in $\epsilon$ is murky and hence we only assume that

(3.1)\begin{equation} \mathcal{C} = \mathcal{C}_r(\epsilon) + {\rm i} \mathcal{C}_i(\epsilon), \quad \text{with } \mathcal{C}_i(\epsilon)\ll \mathcal{C}_r(\epsilon),\quad \epsilon\to 0. \end{equation}

We use the notation of Bender & Orszag (Reference Bender and Orszag1999) for ‘$\mathcal {C}_i(\epsilon )$ is much smaller than $\mathcal {C}_r(\epsilon )$ as $\epsilon$ tends to $0$’, namely that

(3.2)\begin{equation} \lim_{\epsilon\to 0} \frac{\mathcal{C}_i(\epsilon)}{\mathcal{C}_r(\epsilon)} =0. \end{equation}

The purpose of such an assumption is to replace $\mathcal {C}\in \mathbb {C}$ by $\mathcal {C}_r\in \mathbb {R}$ in (2.11) and (2.12). Recalling that $\mathcal {U}(\mathcal {z}>0)>0$ and $\mathcal {U}(\mathcal {z}<0)<0$, if $\mathcal {C}_r$ is positive (negative) then there is a critical layer in the air (water), associated with a singular point in (2.11) (2.12). However, at this stage in the development, we do not know the possible values for $\mathcal {C}_r$. Indeed, the term $rf^{\prime}(0^+)-h^{\prime}(0^-)$ depends on $\mathcal {C}$ in a non-trivial manner so that (2.14) is not necessarily quadratic in $\mathcal {C}$. Physically, $\mathcal {C}_r$ is the phase speed of sheared interfacial waves in the reference frame of the water surface. In that frame, a wave propagating in the direction of (against) the current has a positive (negative) $\mathcal {C}_r$ and Young & Wolfe (Reference Young and Wolfe2014) refer to these as prograde (retrograde) modes. In the case of a constant (zero shear) current $U$, there is one prograde and one retrograde mode, which are simply the Doppler-shifted forward and backward interfacial waves. We generalize this result to the case of arbitrary shear, which we check a posteriori. Thus we assume that there are two solutions, $\mathcal {C}_{r+}>0$ and $\mathcal {C}_{r-}<0$, that correspond to two critical levels, $\mathcal {z}_{{c}+}>0$ and $\mathcal {z}_{{c}-}<0$, such that

(3.3)\begin{equation} \mathcal{U}(\mathcal{z}_{{c}\pm}) = \mathcal{C}_{r\pm}. \end{equation}

Therefore, the prograde (retrograde) mode undergo the Miles (rippling) instability, and hence $\mathcal {C}_i \neq 0$ for both modes. We stress that, provided that the values of $\mathcal {C}_{r\pm }$ are within the bounds of the function $\mathcal {U}(\mathcal {z})$, critical layers actually exist. When $\mathcal {C}_{r\pm }$ are equal to these bounds, the system is marginally stable. In Appendix C, we derive a general asymptotic formula for the large wavenumber cutoff of the rippling instability.

Hence, although $\mathcal {C}_{r\pm }$ are unknown, we replace $\mathcal {C}$ by $\mathcal {C}_{r+}$ in (2.11) and by $\mathcal {C}_{r-}$ in (2.12). We then solve these equations using boundary layer theory in § 6, where we construct uniformly valid composite solutions. However, here we need only substitute the derivatives at $\mathcal {z}=0^\pm$ into (2.14). In Appendix B, we show that, for the prograde mode,

(3.4)$$\begin{gather} f^{\prime}(0^+)={-}{\rm i}{\rm \pi} \frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{{c}+})}{|\mathcal{U}^{\prime}(\mathcal{z}_{{c}+})|}\, {\rm e}^{{-}2\mathcal{z}_{{c}+}/\epsilon} + \epsilon \frac{\mathcal{U}^{\prime\prime}(0^+)}{2\mathcal{C}_{r+}(\epsilon)} -\frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{{c}+})}{\mathcal{U}^{\prime}(\mathcal{z}_{{c}+})} \frac{\epsilon}{2\mathcal{z}_{{c}+}}, \end{gather}$$

(3.5)$$\begin{gather}\text{and}\quad h^{\prime}(0^-) ={-} \epsilon \frac{\mathcal{U}^{\prime\prime}(0^-)}{2\mathcal{C}_{r+}(\epsilon)}, \end{gather}$$

whereas, for the retrograde mode,

(3.6)$$\begin{gather} f^{\prime}(0^+) = \epsilon \frac{\mathcal{U}^{\prime\prime}(0^+)}{2\mathcal{C}_{r-}(\epsilon)}, \end{gather}$$

(3.7)$$\begin{gather}\text{and}\quad h^{\prime}(0^-) = {\rm i}{\rm \pi} \frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{{c}-})}{|\mathcal{U}^{\prime}(\mathcal{z}_{{c}-})|} \, {\rm e}^{2\mathcal{z}_{{c}-}/\epsilon} - \epsilon \frac{\mathcal{U}^{\prime\prime}(r^-)}{2\mathcal{C}_{r-}(\epsilon)} +\frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{{c}-})}{\mathcal{U}^{\prime}(\mathcal{z}_{{c}-})} \frac{\epsilon}{2\mathcal{z}_{{c}-}}. \end{gather}$$

We can now solve (2.14). With the advent of results (3.4) and (3.7), our key assumption, (3.1), is made more precise by letting

(3.8)\begin{equation} \mathcal{C}_\pm= \mathcal{C}_{r\pm}(\epsilon) + {\rm i} A_\pm(\epsilon) \epsilon \,{\rm e}^{{-}2|\mathcal{z}_{{c}\pm}|/\epsilon}. \end{equation}

Therefore, in order to have solutions with a non-zero imaginary part, we must include exponentially small terms. When first discarding those terms, we obtain

(3.9)\begin{equation} \frac{\mathcal{S}}{ \epsilon} + \mathcal{G}(1-r)\epsilon + [ r \mathcal{U}^{\prime}(0^+)-\mathcal{U}^{\prime}(0^-)] \epsilon \mathcal{C}_r - (1+r + {\rm h.o.t.}) \mathcal{C}_r^2 =0,\quad \epsilon\ll1,\end{equation}

where h.o.t denotes higher-order terms, such as terms of order $\epsilon /\mathcal {C}_r$ and $\epsilon /z_{{c}}$.

We seek solutions of (3.9) as a series in powers of $\epsilon$. The nature of the leading-order term depends on the value of the dimensionless surface tension $\mathcal {S}$. When gravity is the only restoring force, we find

(3.10) \begin{equation} \mathcal{C}_{r\pm}^{{grav}}={\pm} \sqrt{\frac{1\,{-}\,r}{1\,{+}\,r} \mathcal{G}\epsilon} + \frac{r\mathcal{U}^{\prime}(0^+)\,{-}\,\mathcal{U}^{\prime}(0^-)}{2(1+r)} \epsilon \,{\pm}\, \frac{[r\mathcal{U}^{\prime}(0^+)\,{-}\,\mathcal{U}^{\prime}(0^-)]^2}{8\sqrt{\mathcal{G}(1-r)}} \left[\frac{\epsilon}{1\,{+}\,r}\right]^{{3}/{2}} \,{+}\, {\rm h.o.t}. \end{equation}

When capillary forces are included, we instead obtain

(3.11)\begin{equation} \mathcal{C}_{r\pm}^{{cap-grav}}={\pm} \sqrt{\frac{\mathcal{S}}{(1+r)\epsilon}} + \frac{r\mathcal{U}^{\prime}(0^+)-\mathcal{U}^{\prime}(0^-)}{2(1+r)} \epsilon \pm \frac{\mathcal{G}(1-r)}{2\sqrt{\mathcal{S}(1+r)}} \epsilon^{{3}/{2}} + {\rm h.o.t}.\end{equation}

Hence, apart from the effect of shear, for short waves gravity acts as a high-order correction to the effect of surface tension. Note that the third term in (3.11) can be derived by expanding the dispersion relation of interfacial capillary–gravity waves

(3.12)\begin{equation} \mathcal{C}(\epsilon)={\pm} \sqrt{\frac{\mathcal{S}}{(1+r)\epsilon} + \frac{1-r}{1+r} \mathcal{G}\epsilon} . \end{equation}

Armed with $\mathcal {C}_{r\pm }$, we can use assumption (3.8), and (3.4) and (3.7), in (2.14). We collect terms of order $\epsilon {\rm e}^{-2|\mathcal {z}_{{c}\pm }|/\epsilon }$ to find the amplitudes $A_\pm$, and infer the growth rates of the prograde ($+$) and retrograde ($-$) modes as

(3.13)$$\begin{gather} {\rm Im}\{\mathcal{k}\mathcal{C}_+\} ={-} \mathcal{C}_{r+}\ \frac{\rm \pi}{2} \frac{r}{1+r}\ \frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{{c}+})}{|\mathcal{U}^{\prime}(\mathcal{z}_{{c}+})|}\, {\rm e}^{{-}2\mathcal{k}\mathcal{z}_{{c}+}}, \end{gather}$$

(3.14)$$\begin{gather}\text{and}\quad {\rm Im}\{\mathcal{k}\mathcal{C}_-\} ={-} \mathcal{C}_{r-} \frac{\rm \pi}{2}\ \frac{1}{1+r} \frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{{c}-})}{|\mathcal{U}^{\prime}(\mathcal{z}_{{c}-})|}\, {\rm e}^{2\mathcal{k}\mathcal{z}_{{c}-}}, \end{gather}$$

where $\mathcal {k}=1/\epsilon$. We emphasize several aspects of the present results. Firstly, (3.13) for the prograde mode is a generalization of the growth rate obtained by Miles in an appendix of Morland & Saffman (Reference Morland and Saffman1992). That result was obtained from short wavelength asymptotic analysis of the exact solution of the Rayleigh equation for an exponential wind profile. Secondly, (3.14) for the retrograde mode is a generalization of the short wavelength limit of the growth rate of the rippling instability found by Young & Wolfe (Reference Young and Wolfe2014) for an exponential wind-induced current. Here, we have included the effect of the upper fluid on the rippling instability, showing the weakness of the effect for the air–water system because $r=O(10^{-3})$, consistent with the $r=0$ limit of Young & Wolfe (Reference Young and Wolfe2014). Finally, we find that the results of Shrira (Reference Shrira1993) are a special case of our own when $r=0$, but stress that his small parameter was the smallness of the deviation of the wave motion from potential flow, rather than the inverse wavenumber.

4. Interpretation

The behaviour of short waves in presence of a wind and a wind-induced current depends on the profile of the latter solely through the derivatives at $\mathcal {z}=0$ and at the critical levels, $\mathcal {z}=\mathcal {z}_{{c}\pm }$. In consequence, the results are qualitatively the same for the two profiles introduced in § 2.3. We note that, as shown by Young & Wolfe (Reference Young and Wolfe2014), the eigenvalue problem can be solved exactly in terms of hypergeometric functions in the case of double exponential profiles, whereas exact analytical solutions for double log profiles are unknown.

In figures 2 and 3, we show the phase speed and the growth rate of the prograde mode undergoing the Miles instability. We compare the results of our short wavelength asymptotic analysis with the exact solution for an exponential wind profile, in the case of gravity (a and d), capillary (b and e) and capillary–gravity waves (c and f), for different values of the density ratio $r$ and the control parameters $\mathcal {G}$ and $\mathcal {S}$. Figure 4 shows a similar comparison for the retrograde mode undergoing the rippling instability in the absence of air, that is $r=0$. The asymptotic results are in accord with the exact results.

Figure 2. Comparison of the short wavelength asymptotic results (dashed lines) for the prograde mode with the exact solution (solid lines) in the case of an exponential wind profile (no current), $U(z)= U_\infty (1-{\rm e}^{-z/d_{{a}}})$, a density ratio $r=0.001$ and different values (grey and black) of the dimensionless gravity, $\mathcal {G}=gd_{{a}}/U_\infty ^2$, and surface tension, $\mathcal {S}=\sigma /(\rho _{{w}} U_\infty ^2d_{{a}})$. The phase speed as a function of the wavenumber is shown for gravity (a), capillary (b) and capillary–gravity waves (c); the exact and the asymptotic solutions are indistinguishable. The growth rate as a function of the wavenumber for gravity (d), capillary (e) and capillary–gravity waves (f).

Figure 3. Same as figure 2 but with a density ratio of $r=0.9$. Whereas in figure 2(a–c) the asymptotic results (dashed lines) and the exact solution (solid lines) are indistinguishable, such is not the case here, where we also see the phase speed of free surface waves (dotted lines).

Figure 4. Comparison of the short wavelength asymptotic results (dashed lines) for the retrograde mode with the exact solution (solid lines) in the case of an exponential current (no air), $U(z)= U_{{S}}\, {\rm e}^{z/d_{{w}}}$ and different values (grey and black) of the dimensionless gravity, $\mathcal {G}=gd_{{w}}/U_{{S}}^2$, and surface tension, $\mathcal {S}=\sigma /(\rho _{{w}} U_{{S}}^2d_{{w}})$. The dotted lines depict the corresponding phase speed of Doppler-shifted surface waves. The growth rate as a function of the wavenumber for gravity (d), capillary (e) and capillary–gravity waves (f).

The phase speed of Doppler-shifted interfacial waves in presence of a constant, zero shear current is

(4.1)\begin{equation} c(k)=U_{{S}} \pm \sqrt{\frac{1-r}{1+r} \frac{g}{k}+\frac{\sigma}{\rho_{{w}}(1+r)} k},\end{equation}

which is the dimensional form of (3.12), where the plus (minus) sign corresponds to the prograde (retrograde) mode. Equations (3.10) and (3.11) predict that the non-uniformity of the current increases the phase speed by $r U^{\prime}(0^+)/2k$ and reduces it by $U^{\prime}(0^-)/2k$, respectively. Indeed, figure 3(a–c) shows that, when $r$ is close to $1$, and $\mathcal {k}$ approaches unity, the phase speed of sheared waves is significantly larger than the phase speed of free surface waves. Moreover, figure 4(a–c) shows the predicted reduction of the retrograde mode in the absence of air when $\mathcal {k}$ is of order unity.

Generally, for small values of $r$, the effect of the wind on dispersion is insignificant, as shown by the phase speed of free surface in figure 2(a–c). However, for large values of $r$, the wind has a major influence. For instance, as shown in figure 3(b), the wind is responsible for a minimum phase speed of capillary waves. Moreover, as $r$ increases, so does the growth rate of the prograde mode. Thus, we interpret the density ratio as a wind–wave coupling constant.

Finally, for both modes the growth rates increase when $\mathcal {G}$ and $\mathcal {S}$ are small, which is the case for large velocity scales (cf. figure 2 and 4). Thus, consistent with physical intuition, a small perturbation grows faster in the presence of strong winds and/or currents.

5. Prograde instability due to critical layers in both air and water

We now explore the case of a wind blowing in the air when there is also a current in the water, as shown in figure 5. Thus we have $\mathcal {U}(\mathcal {z}\lessgtr 0)>0$, $\mathcal {U}^{\prime}(\mathcal {z}>0)>0$, and $\mathcal {U}^{\prime}(\mathcal {z}<0)<0$, and again assume an infinite depth and proceed as in § 3. The key difference is that the prograde mode has a critical layer in both the air and the water, and the retrograde mode has no critical layer. Therefore, there are levels $\mathcal {z}_{{c}\pm }$ such that

(5.1a,b)\begin{equation} \mathcal{U}(\mathcal{z}_{{c}+}) = \mathcal{C}_{r+}\quad \text{and}\quad \mathcal{U}(\mathcal{z}_{{c}-}) = \mathcal{C}_{r+}, \end{equation}

and there is no value of $\mathcal {z}$ such that $\mathcal {C}_{r-}<0$ is equal to $\mathcal {U}(\mathcal {z})$. This implies that

(5.2)$$\begin{gather} f^{\prime}(0^+) ={-}{\rm i}{\rm \pi} \frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{{c}+})}{|\mathcal{U}^{\prime}(\mathcal{z}_{{c}+})|} \, {\rm e}^{{-}2 \mathcal{z}_{{c}+}/\epsilon} + \epsilon \frac{\mathcal{U}^{\prime\prime}(0^+)}{2\mathcal{C}_{r+}(\epsilon)} - \frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{{c}+})}{\mathcal{U}^{\prime}(\mathcal{z}_{{c}+})} \frac{\epsilon}{2\mathcal{z}_{{c}+}} \end{gather}$$

(5.3)$$\begin{gather}\text{and}\quad h^{\prime}(0^-) = {\rm i}{\rm \pi} \frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{{c}-})}{|\mathcal{U}^{\prime}(\mathcal{z}_{{c}-})|}\, {\rm e}^{2 \mathcal{z}_{{c}-}/\epsilon} - \epsilon \frac{\mathcal{U}^{\prime\prime}(0^-)}{2\mathcal{C}_{r+}(\epsilon)} + \frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{{c}-})}{\mathcal{U}^{\prime}(\mathcal{z}_{{c}-})} \frac{\epsilon}{2\mathcal{z}_{{c}-}}\end{gather}$$

for the prograde mode, while

(5.4)$$\begin{gather} f^{\prime}(0^+) = \epsilon \frac{\mathcal{U}^{\prime\prime}(0^+)}{2\mathcal{C}_{r-}(\epsilon)} \end{gather}$$

(5.5)$$\begin{gather}\text{and}\quad h^{\prime}(0^-) ={-} \epsilon \frac{\mathcal{U}^{\prime\prime}(0^-)}{2\mathcal{C}_{r-}(\epsilon)}, \end{gather}$$

for the retrograde mode. Hence, the retrograde mode is neutral but the prograde mode can undergo both the Miles and rippling instabilities, and has a growth rate of

(5.6)\begin{equation} {\rm Im}\{\mathcal{k}\mathcal{C}_+\} ={-} \frac{{\rm \pi} }{1+r} \frac{\mathcal{C}_{r+}}{2}\left[ r \frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{{c}+})}{|\mathcal{U}^{\prime}(\mathcal{z}_{{c}+})|}\, {\rm e}^{{-}2\mathcal{k}\mathcal{z}_{{c}+}}+\frac{\mathcal{U}^{\prime\prime}(\mathcal{z}_{{c}-})}{|\mathcal{U}^{\prime}(\mathcal{z}_{{c}-})|}\, {\rm e}^{2\mathcal{k}\mathcal{z}_{{c}-}} \right]. \end{equation}

The real parts $\mathcal {C}_{r\pm }$ are still given by (3.10) and (3.11). Such an enhancement of the Miles instability by the rippling instability may be difficult to observe in geophysical flows. However, it could be examined in a controlled laboratory setting through a refinement of the viscosity-stratified approach of Charles & Lilleleht (Reference Charles and Lilleleht1965) or the two-layer Couette flow approach of Barthelet, Charru & Fabre (Reference Barthelet, Charru and Fabre1995), as well as in the context of Holmboe wave experiments (e.g. Carpenter et al. Reference Carpenter, Tedford, Rahmani and Lawrence2010).

Figure 5. Schematic of a wind blowing in the air and a current in the water.

6. Asymptotic solution of the Rayleigh equation for short waves

Here, we solve (2.11) for $\mathcal {C}=\mathcal {C}_{r+}$ and note that a similar procedure is applicable to (2.12) when $\mathcal {C}=\mathcal {C}_{r-}$. For simplicity, we rewrite (2.11):

(6.1)\begin{equation} \epsilon f^{\prime\prime}(\mathcal{z})- 2 f^{\prime}(\mathcal{z})- \epsilon \frac{\mathcal{U}^{\prime\prime}(\mathcal{z})} { \mathcal{U}(\mathcal{z})-\mathcal{C}_r } f(\mathcal{z})=0, \quad f(0)=1, \end{equation}

and we drop the subscript $+$ for the rest of this section. There is a regular singularity at $\mathcal {z}=\mathcal {z}_{{c}}$ such that

(6.2)\begin{equation} \mathcal{U}(\mathcal{z}_{{c}}) = \mathcal{C}_r. \end{equation}

Because the small parameter $\epsilon$ multiplies the highest-order derivative in (2.11), we expect the solution to have a boundary layer somewhere in the domain, but do not know its location a priori. However, guidance is provided by the presence of the singularity. The Frobenius exponents are $0$ and $1$, so that the solution of (2.11) is finite at $\mathcal {z}=\mathcal {z}_{{c}}$ whereas its derivative has a logarithmic divergence (Drazin & Reid Reference Drazin and Reid1981). Therefore, we assume that an internal boundary layer emerges from the singularity. We assume that the point $\mathcal {z}=\mathcal {z}_{{c}}(\epsilon )$ is well separated from the lower boundary, $\mathcal {z} = 0$, as $\epsilon$ goes to zero. We check this a posteriori once the dependence of $\mathcal {C}_r$ on $\epsilon$ is known. In consequence, $\mathcal {C}_r$ can be treated as a constant in the following analysis.

The boundary layer is an inner region where the solution of (2.11) changes rapidly. We define two outer regions, where the solution changes slowly. One spans $\mathcal {z}=0$ to the inner region; the other spans the inner region to the far field. We call them ‘lower outer region’ and ‘upper outer region’, respectively (see figure 6).

Figure 6. Schematic of the domain of analysis of the Rayleigh equation for short waves; $\mathcal {k}\gg 1$. A boundary layer of thickness $\epsilon =1/\mathcal {k}$ emerges from the singularity at the critical level, $\mathcal {z}=\mathcal {z}_{{c}}$.

6.1. Outer solutions

Following Miles (Reference Miles1962), we seek outer solutions as a power series in $\epsilon$ as

(6.3)\begin{equation} f_{{out}}(\mathcal{z})=f_0(\mathcal{z}) + \epsilon f_1(\mathcal{z}) + O(\epsilon^2),\quad \epsilon\to0^+, \end{equation}

and find that

(6.4a,b)\begin{equation} f_0 = \text{const.}\quad \text{and}\quad f_1(\mathcal{z}) ={-}\frac{f_0}{2} \int^{\mathcal{z}}_{\mathcal{z}_1} {\rm d}\tilde{z} \frac{\mathcal{U}^{\prime\prime}(\tilde{z})}{\mathcal{U}(\tilde{z})-\mathcal{C}_r }.\end{equation}

The constant $f_0$ and the lower limit of the integral $\mathcal {z}_1$ are not a priori the same for the two outer solutions. In particular, they cannot be determined by the boundary conditions at $z=0$ and infinity, requiring us to find an inner solution.

6.2. Inner solution

Within the boundary layer, we introduce the stretched coordinate $Z\equiv (\mathcal {z}-\mathcal {z}_{{c}})/\delta$, where $0<\delta \ll 1$, and seek an inner solution, $F_{{in}}(Z)$, to (2.11). We approximate the coefficients by their Taylor series expansions about $\mathcal {z}=\mathcal {z}_{{c}}$, so that (2.11) becomes

(6.5)\begin{equation} \frac{\epsilon}{\delta^2 } F_{{in}}^{\prime\prime}(Z)-\frac{2}{\delta} F^{\prime}_{{in}}(Z) -\frac{\epsilon}{\delta} \frac{\mathcal{U}^{\prime\prime}_{{c}}}{\mathcal{U}^{\prime}_{c}Z} F_{{in}}(Z)=0;\quad Z=O(1), \end{equation}

where the subscript ‘$c$’ denotes evaluation at the critical level, $\mathcal {z}=\mathcal {z}_{{c}}$. By balancing the two first terms, we obtain the distinguished limit $\delta = \epsilon$ (Bender & Orszag Reference Bender and Orszag1999), and hence must solve

(6.6)\begin{equation} F_{{in}}^{\prime\prime}(Z) - 2 F^{\prime}_{{in}}(Z)=\epsilon \frac{\mathcal{U}^{\prime\prime}_{{c}}}{\mathcal{U}^{\prime}_{{c}}Z} F_{{in}}(Z),\quad \epsilon\to 0^+. \end{equation}

Because the boundary layer has an $O(\epsilon )$ thickness, the appropriate inner expansion is

(6.7)\begin{equation} F_{{in}}(Z)= F_0(Z) +\epsilon\ F_1(Z) + O(\epsilon^2), \quad \epsilon\to 0^+. \end{equation}

The general solutions are

(6.8)$$\begin{gather} F_0(Z) = A\, {\rm e}^{2Z} + B,\quad A,B \in \mathbb{C}, \end{gather}$$

(6.9)$$\begin{gather}\text{and}\quad F_1(Z) = \frac{\mathcal{U}^{\prime\prime}_{{c}}}{\mathcal{U}^{\prime}_{{c}}} \int^Z_b {{\rm d}\kern0.7pt x}\, {\rm e}^{2x} \int^x_a {\rm d}t \frac{{\rm e}^{{-}2t}}{t} F_0(t). \end{gather}$$

Next, we determine the integration constants, $A$ and $B$, and the bounds of integration, $a$ and $b$, by asymptotic matching.

6.3. Asymptotic matching and uniformly valid composite solutions

We must match the inner solution and the two outer solutions, as $\epsilon \to 0^+$. We use the superscripts ‘$\ell$’ and ‘u’ for lower and upper, respectively.

The lower outer solution, $f_{{out}}^{\ell }(\mathcal {z})$, is defined for $0\le \mathcal {z}\ll \mathcal {z}_{{c}}$, and must satisfy the lower boundary condition

(6.10a,b)\begin{equation} f_0^{\ell}(0)=1,\quad \text{and}\quad f_1^{\ell}(0)=0. \end{equation}

The upper outer solution, $f_{{out}}^{{u}}(\mathcal {z})$, is defined for $\mathcal {z}\gg \mathcal {z}_{{c}}$ and hence must satisfy condition (2.13a). We use the following matching conditions:

(6.11a,b)\begin{equation} \lim_{\mathcal{z}\to \mathcal{z}_{{c}}^- }f_{{out}}^{\ell}(\mathcal{z})=\lim_{Z\to-\infty}F_{{in}}(Z)\quad \text{and}\quad \lim_{\mathcal{z}\to \mathcal{z}_{{c}}^+}f_{{out}}^{{u}}(z)=\lim_{Z\to+\infty}F_{{in}}(Z),\end{equation}

and then apply the Van Dyke additive rule (Bender & Orszag Reference Bender and Orszag1999) to construct uniformly valid composite solutions

(6.12)\begin{equation} \text{uniform approx} = \text{inner}+\text{outer} - \text{common part}. \end{equation}

We stress that this is possible because the lower and upper outer solutions happen to have a common analytical expression.

6.3.1. Leading order

To leading order, the outer solutions are constant (see (6.4a)), and because of the boundary condition (6.10a) we have

(6.13)\begin{equation} f_0^{\ell}=1. \end{equation}

The leading-order inner solution is given by (6.8). To preempt divergence as $Z\to +\infty$, we impose $A=0$, and the matching condition (6.11a) implies that $B=1$. Therefore

(6.14)\begin{equation} F_0 =1, \end{equation}

which, upon imposition of the matching condition (6.11b), yields

(6.15)\begin{equation} f_0^{{u}}=1.\end{equation}

We combine the solutions (6.13), (6.14) and (6.15) using the additive rule (6.12), and thus obtain a uniformly valid composite solution at leading order

(6.16)\begin{equation} f_{{unif},0}= 1+O(\epsilon). \end{equation}

The effect of the boundary layer appears only at the next order.

6.3.2. Order $\epsilon$

Following (6.4b), the lower and upper outer solutions at order $\epsilon$ are

(6.17a,b)\begin{equation} f_1^{\ell}(\mathcal{z}) ={-}\frac{1}{2} \int^{\mathcal{z}}_{\mathcal{z}_1^{\ell}} {\rm d}\tilde{z} \frac{\mathcal{U}^{\prime\prime}(\tilde{z})}{\mathcal{U}(\tilde{z})-\mathcal{C}_r }\quad \text{and}\quad f_1^{{u}}(\mathcal{z}) ={-}\frac{1}{2} \int^{\mathcal{z}}_{\mathcal{z}_1^{{u}}} {\rm d}\tilde{z} \frac{\mathcal{U}^{\prime\prime}(\tilde{z})}{\mathcal{U}(\tilde{z})-\mathcal{C}_r }, \end{equation}

respectively. From the Laurent series expansion

(6.18)\begin{equation} \frac{\mathcal{U}^{\prime\prime}(\mathcal{z})} { \mathcal{U}(\mathcal{z})-\mathcal{C}_r }= \frac{\mathcal{U}^{\prime\prime}_{{c}}}{\mathcal{U}^{\prime}_{{c}}(\mathcal{z}-\mathcal{z}_{{c}})}+ \frac{1}{2}\left[\frac{\mathcal{U}^{\prime\prime}_{{c}}}{\mathcal{U}^{\prime}_{{c}}}\right]^2 + O(\mathcal{z} -\mathcal{z}_{{c}}), \end{equation}

we deduce that

(6.19)$$\begin{gather} f_1^{\ell}(\mathcal{z}) \sim{-} \frac{\mathcal{U}^{\prime\prime}_{{c}}}{2\mathcal{U}^{\prime}_{{c}}} \textrm{Log}\left(\frac{\mathcal{z}-\mathcal{z}_{{c}} }{\mathcal{z}_1^{\ell}-\mathcal{z}_{{c}}}\right),\quad \mathcal{z}\to\mathcal{z}_{{c}}^-, \end{gather}$$

(6.20)$$\begin{gather}\text{and}\quad f_1^{{u}}(\mathcal{z}) \sim{-}\frac{\mathcal{U}^{\prime\prime}_{{c}}}{2\mathcal{U}^{\prime}_{{c}}} \textrm{Log}\left(\frac{\mathcal{z}-\mathcal{z}_{{c}} }{\mathcal{z}_1^{{u}}-\mathcal{z}_{{c}}}\right),\quad \mathcal{z}\to\mathcal{z}_{{c}}^+, \end{gather}$$

where $\textrm {Log}$ denotes a continuation of the natural logarithm to the negative real numbers

(6.21)\begin{equation} \textrm{Log}(\mathcal{z}-\mathcal{z}_{{c}})= \begin{cases} \ln|\mathcal{z}-\mathcal{z}_{{c}}|- {\rm i}{\rm \pi} & \text{if }\mathcal{U}^{\prime}_{c}>0,\\ \ln|\mathcal{z}-\mathcal{z}_{{c}}|+ {\rm i}{\rm \pi} & \text{if }\mathcal{U}^{\prime}_{c}<0, \end{cases} \quad \text{for }\mathcal{z}<\mathcal{z}_{{c}}.\end{equation}

The choice of the branch cut – just above the negative real axis if $\mathcal{U}^{\prime}_{c} >0$, just below otherwise – follows from Lin (Reference Lin1955).

From (6.9), the inner solution at order $\epsilon$ is

(6.22)\begin{equation} F_1(Z) = \frac{\mathcal{U}^{\prime\prime}_{{c}}}{\mathcal{U}^{\prime}_{{c}}} \int^Z_b {{\rm d}\kern0.7pt x}\, {\rm e}^{2x} \int^x_a {\rm d}t \frac{{\rm e}^{{-}2t}}{t}. \end{equation}

We split the integral over $t$ as

(6.23)\begin{equation} \int^x_a = \int^x_{+\infty} + \int_a^{+\infty}. \end{equation}

Noting that the second integral is a constant, $C\in \mathbb {C}$, we obtain

(6.24)\begin{equation} F_1(Z)={-}\frac{\mathcal{U}^{\prime\prime}_{{c}}}{\mathcal{U}^{\prime}_{{c}}}\left\{ \int^Z_b {{\rm d}\kern0.7pt x}\, {\rm e}^{2x} E_1(2x)+ \frac{C}{2}( {\rm e}^{2Z}- {\rm e}^{2b})\right\}, \end{equation}

where

(6.25)\begin{equation} E_1(x) \equiv\int_x^{+\infty} {\rm d}t \frac{{\rm e}^{{-}t}}{t}, \quad |\arg(x)|<{\rm \pi}, \end{equation}

is the exponential integral. For large values of $|x|$, the divergent series (Bender & Orszag Reference Bender and Orszag1999)

(6.26)\begin{equation} E_1(x)=\frac{{\rm e}^{{-}x}}{x} \sum_{n=0}^{N-1} \frac{n!}{({-}x)^n}\ ,\quad x\to \infty;\quad |\arg(x)|<\frac{3{\rm \pi}}{2}, \end{equation}

yields

(6.27)\begin{equation} {\rm e}^{2x} E_1(2x)\sim\frac{1}{2x}, \quad x\to \pm\infty. \end{equation}

After integration, we find

(6.28)\begin{equation} F_1(Z) \sim{-}\frac{\mathcal{U}^{\prime\prime}_{{c}}}{2\mathcal{U}^{\prime}_{{c}}}\left\{ \textrm{Log}\left[ \frac{Z}{b} \right]+ C( {\rm e}^{2Z}- {\rm e}^{2b})\right\}, \quad Z\to \pm\infty.\end{equation}

Recalling that $Z=(\mathcal {z}-\mathcal {z}_{{c}})/\epsilon$, the matching of the inner limits, given by (6.19) and (6.20), and the outer limits, given by (6.28), yields

(6.29a,b)\begin{equation} C = 0 \quad \text{and}\quad \mathcal{z}_1^{\ell}-\mathcal{z}_{{c}} =\epsilon b = \mathcal{z}_1^{{u}}-\mathcal{z}_{{c}}. \end{equation}

Because the boundary condition (6.10b) requires that

(6.30)\begin{equation} \mathcal{z}_1^{\ell}=0, \end{equation}

the outer solution at order $\epsilon$ has the same expression in the lower and upper outer regions

(6.31)\begin{equation} f_1(\mathcal{z}) ={-}\frac{1}{2} \int_0^{\mathcal{z}} {\rm d}\tilde{z} \frac{\mathcal{U}^{\prime\prime}(\tilde{z})} { \mathcal{U}(\tilde{z})-\mathcal{C}_r }, \quad \mathcal{z}\ll\mathcal{z}_{{c}}\quad \text{and}\quad \mathcal{z}\gg\mathcal{z}_{{c}}.\end{equation}

The lower outer solution is real, because the path of integration along the real axis does not reach the branch point $\mathcal {z}=\mathcal {z}_{{c}}$. However, for the upper outer solution we have to make a detour in the complex plane. We deform the path of integration in the upper part if $\mathcal{U}^{\prime}_{{c}} >0$, and in the lower part otherwise. In Appendix B, we show that

(6.32)\begin{equation} {\rm Im}\{\kern1.5pt f_1^{{u}}\} ={-}\frac{\rm \pi}{2} \frac{\mathcal{U}^{\prime\prime}_{{c}}}{|\mathcal{U}^{\prime}_{{c}}|}. \end{equation}

Using the matching conditions (6.29a,b) together with (6.30), the inner solution (6.24) becomes

(6.33)\begin{equation} F_1(Z)={-}\frac{\mathcal{U}^{\prime\prime}_{{c}}}{\mathcal{U}^{\prime}_{{c}}} \int^Z_{-\mathcal{z}_{{c}}/\epsilon} {{\rm d}\kern0.7pt x}\, {\rm e}^{2x} E_1(2x). \end{equation}

Finally, we use the Van Dyke rule (6.12) to construct a uniformly valid composite solution at order $\epsilon$

(6.34)\begin{align} f_{{unif},1}(\mathcal{z}) &= 1+ \epsilon\left\{ -\frac{\mathcal{U}^{\prime\prime}_{{c}}}{\mathcal{U}^{\prime}_{{c}}} \int^Z_{-\mathcal{z}_{{c}}/\epsilon} {{\rm d}\kern0.7pt x}\, {\rm e}^{2x} E_1(2x)-\frac{1}{2} \int_0^\mathcal{z} {\rm d}\tilde{z} \frac{\mathcal{U}^{\prime\prime}(\tilde{z})} { \mathcal{U}(\tilde{z})-\mathcal{C}_r }\right. \nonumber\\ & \quad +\left.\frac{\mathcal{U}^{\prime\prime}_{{c}}}{2\mathcal{U}^{\prime}_{{c}} } \textrm{Log}\left( 1- \frac{\mathcal{z}}{\mathcal{z}_{{c}}} \right)\right\}+O(\epsilon^2). \end{align}

In Appendix B, we verify that $\chi (\mathcal {z}) = {\rm e}^{-\mathcal {k}\mathcal {z}}f_{{unif},1}(\mathcal {z})$ satisfies the global property:

(6.35)\begin{equation} {\rm Im}\{\chi^{\prime}(0^+)\} ={-}{\rm \pi} \frac{\mathcal{U}^{\prime\prime}_{{c}}}{|\mathcal{U}^{\prime}_{c}|}\ |\chi_{{c}}|^2\end{equation}

at leading order. In Appendix A we show that

(6.36)\begin{equation} f(\mathcal{z}_{{c}}) = 1 +\frac{\epsilon}{2} \frac{\mathcal{U}^{\prime\prime}_{{c}}}{\mathcal{U}^{\prime}_{{c}}} \left(\gamma_{{E}} +\ln\left[\frac{2\mathcal{z}_{{c}}}{\epsilon}\right] - i {\rm \pi}\right) +O(\epsilon^2). \end{equation}

Equation (6.36) generalizes a result that Miles derived in an appendix to Morland & Saffman (Reference Morland and Saffman1992) using an exact solution of the Rayleigh equation for an exponential wind profile.

6.4. Comparison with the numerical solution

We numerically solve (2.11) for the two standard wind profiles

(6.37a,b)\begin{equation} \mathcal{U}(\mathcal{z}) = 1-{\rm e}^{-\mathcal{z}} \quad \text{and}\quad \mathcal{U}(\mathcal{z}) = \ln(1+\mathcal{z})/\kappa, \end{equation}

where, as before, $\kappa =0.4$ is the von Kármán constant. We compare the numerical solution with our uniformly valid composite solution (6.34) for fixed values of $\mathcal {k}$ and $\mathcal {C}_r$ in figure 7. (Note that any dispersion relation can be retrieved with a proper choice of the control parameters.) For both profiles the composite solution and the numerical solution agree very well. We distinguish the lower (upper) outer solutions with their imaginary part being equal to zero (a positive constant). Consistent with the Frobenius solution (Drazin & Reid Reference Drazin and Reid1981), the solution within the inner layer depends on the wind profile and the dispersion relation solely through the scale factor $\mathcal{U}^{\prime\prime}_{{c}} /\mathcal{U}^{\prime}_{{c}}$ and the bound of integration $\mathcal {z}_{{c}}/\epsilon$. Since the phase of the solution of the Rayleigh equation changes only within the boundary layer, we conclude that the interaction of short waves with the wind principally occurs therein. In contrast, we showed that for $\mathcal {k}\ll 1$ the phase varies from $\mathcal {z}=0$ to $\mathcal {z}=\mathcal {z}_{{c}}$, so that long waves interact with the wind all the way from the mean water surface to the critical level (Bonfils et al. Reference Bonfils, Mitra, Moon and Wettlaufer2022).

Figure 7. Comparison of the uniformly valid composite solution (6.34) with the numerical solution of the Rayleigh equation for the exponential (a,b) and logarithmic (c,d) wind profiles, with $\mathcal {C}_r=0.1$ and different values of $\mathcal {k}$. The dots and the stars denote the real and imaginary parts of the numerical solution, respectively. The solid line shows the real part of (6.34) and the dashed line the imaginary part.

6.5. Similarity solution

We have non-dimensionalized the variables using external parameters characterizing the shear in the air. However, the general short-wave solution (6.34) can also be written in terms of $\xi \equiv k z$ and $\mathfrak {U}\equiv U/\tilde {c}_r$, where $\tilde {c}_r$ is the dimensional version of $\mathcal {C}_r$. The wind-generated ripples have a continuous wavenumber spectrum so that $k$, and subsequently $\tilde {c}_r$, are variables rather than parameters. In that sense, $\xi$ is a similarity variable introduced by Miles (Reference Miles1957), and the Rayleigh equation has a self-similar solution, $\phi =\phi (\xi )$, which for short waves is

(6.38) \begin{align} \phi(\xi) = {\rm e}^{-\xi} \left\{ 1-\frac{\mathfrak{U}^{\prime\prime}(\xi_{{c} })}{\mathfrak{U}^{\prime}(\xi_{{c} })} \int^{\xi-\xi_{{c} }}_{-\xi_{{c} }} {{\rm d}\kern0.7pt x}\, {\rm e}^{2x} E_1(2x)-\frac{1}{2} \int_0^\xi {\rm d}\tilde{z} \frac{\mathfrak{U}^{\prime\prime}(\tilde{z})} {\mathfrak{U}(\tilde{z})-1 }+ \frac{\mathfrak{U}^{\prime\prime}(\xi_{{c} })}{2\mathfrak{U}^{\prime}(\xi_{{c} })} \textrm{Log}\left( 1- \frac{\xi}{\xi_{{c} }} \right)\right\}. \end{align}

7. The WKB method in the short wavelength limit

Alexakis et al. (Reference Alexakis, Young and Rosner2004b) proposed a solution of the Rayleigh equation for short waves using the WKB method. However, because their solution does not satisfy the global property (6.35), their growth rate has an extra factor ${\rm \pi} /\mathcal {k}$; see their (A31). In figure 8, we show the growth rate of gravity waves for an exponential wind profile and compare the results of Alexakis et al. (Reference Alexakis, Young and Rosner2004b) with our short-wave asymptotic solution and the exact solution. Clearly the WKB method deviates from the others for all $k d_a \gtrsim 25$, and in what follows, we explain the origin of the deviation.

Figure 8. Comparison of the growth rate of gravity waves obtained by Alexakis et al. (Reference Alexakis, Young and Rosner2004b) with the predictions of our short wavelength asymptotic treatment (cf. (3.13)) and the growth rate calculated from the exact (hypergeometric) solution for a wind profile $U(z) = U_\infty (1-{\rm e}^{-z/d_{{a}}})$ and a density ratio $r=0.5$.

Setting $\epsilon =1/\mathcal {k}$, we write the Rayleigh equation in the form

(7.1)\begin{equation} \epsilon^2 \chi^{\prime\prime}(\mathcal{z})= Q(\mathcal{z},\epsilon) \chi(\mathcal{z}),\quad \text{with } Q(\mathcal{z},\epsilon)\equiv 1+ \epsilon^2 \frac{\mathcal{U}^{\prime\prime}(\mathcal{z})}{\mathcal{U}(\mathcal{z})-\mathcal{C}_r}. \end{equation}

A WKB expansion is (Bender & Orszag Reference Bender and Orszag1999)

(7.2)\begin{equation} \chi(\mathcal{z})\sim \exp\left[\frac{1}{\epsilon} \sum_{n=0}^{+\infty} \epsilon^n S_n(\mathcal{z})\right], \quad \epsilon\to0. \end{equation}

In the standard case, $Q(\mathcal {z},\epsilon )$ does not depend on the small parameter $\epsilon$. Then, we can take only two terms in the series (7.2); this physical optics approximation yields

(7.3)\begin{align} \chi(\mathcal{z})\sim C_1 [Q(\mathcal{z})]^{-{1}/{4}}\exp\left[\frac{1}{\epsilon}\int_a^\mathcal{z} {\rm d}\tilde{z} \sqrt{Q(\tilde{z})}\right]+ C_2 [Q(\mathcal{z})]^{-{1}/{4}}\exp\left[-\frac{1}{\epsilon}\int_a^\mathcal{z} {\rm d}\tilde{z} \sqrt{Q(\tilde{z})}\right],\end{align}

for $Q(\mathcal {z})>0$, and

(7.4)\begin{align} \chi(\mathcal{z})\sim C_3 [{-}Q(\mathcal{z})]^{-{1}/{4}}\exp\left[\frac{i}{\epsilon}\int_b^\mathcal{z} {\rm d}\tilde{z} \sqrt{-Q(\tilde{z})}\right]+ C_4 [{-}Q(\mathcal{z})]^{-{1}/{4}}\exp\left[-\frac{i}{\epsilon}\int_b^\mathcal{z} {\rm d}\tilde{z} \sqrt{-Q(\tilde{z})}\right],\end{align}

for $Q(\mathcal {z})<0$, where the bounds of integration, $a$ and $b$, are arbitrary. For $\mathcal {U}^{\prime\prime}(\mathcal {z})<0$ and $Q(\mathcal {z},\epsilon )$ introduced in (7.1), Alexakis et al. (Reference Alexakis, Young and Rosner2004b) used the physical optics approximation on three intervals, defined in figure 9, and obtained three solutions of the form of (7.3) or (7.4), in which they neglected the fact that $\epsilon \ll 1$ within $Q(\mathcal {z},\epsilon )$. Then, they had to match those solutions in order to determine the integration constants $C_1$, $C_2$, etc. The common matching procedure takes place at the simple turning point $\mathcal {z}_{\star }$ (e.g. Bender & Orszag Reference Bender and Orszag1999). However, they needed inner solutions near the critical level, $\mathcal {z}_{{c}}$, and it is at this juncture that the appeal of the WKB method is understood. Indeed, they found that the solution of the Rayleigh equation near the singularity can be represented in terms of Bessel functions for $\mathcal {z}>\mathcal {z}_{{c}}$, and in terms of modified Bessel functions for $\mathcal {z}<\mathcal {z}_{{c}}$. Moreover, the outer limits of those inner solutions formally match the inner limits of the physical optics approximations. Unfortunately, however, the distance between the critical point $\mathcal {z}_{{c}}$ and the turning point $\mathcal {z}_{\star }$ actually shrinks as $\epsilon$ tends to zero. For instance, in the case of the exponential profile $\mathcal {U}(\mathcal {z})=1-{\rm e}^{-\mathcal {z}}$, we have

(7.5)\begin{equation} \mathcal{z}_{{\star}}-\mathcal{z}_{{c}} = \ln(1+\epsilon). \end{equation}

Hence, the interval where $Q<0$ cannot be taken as an outer region. In fact, as seen in figure 7, the numerical solution of the Rayleigh equation does not exhibit oscillatory behaviour.

Figure 9. The three intervals in the WKB approach of Alexakis et al. (Reference Alexakis, Young and Rosner2004b) where $\mathcal {z}_{\star }$ is a simple turning point.

In the vicinity of the critical level, $\mathcal {z}_{{c}}$, we introduce the inner variable

(7.6)\begin{equation} {\unicode{x2C7F}}={-}\frac{\mathcal{U}^{\prime\prime}_{{c}}}{\mathcal{U}^{\prime}_{{c}}} (\mathcal{z}-\mathcal{z}_{{c}}), \end{equation}

in terms of which the inner solutions are

(7.7)$$\begin{gather} \chi_{{in+}}({\unicode{x2C7F}}) = \sqrt{{\unicode{x2C7F}}}[ A J_1(2\sqrt{{\unicode{x2C7F}}})+ B Y_1(2\sqrt{{\unicode{x2C7F}}})]\quad\text{for }{\unicode{x2C7F}}>0, \end{gather}$$

(7.8)$$\begin{gather}\text{and}\quad \chi_{{in-}}({\unicode{x2C7F}}) = \sqrt{-{\unicode{x2C7F}}}[ C I_1(2\sqrt{-{\unicode{x2C7F}}})+ D K_1(2\sqrt{-{\unicode{x2C7F}}})]\quad\text{for }{\unicode{x2C7F}}<0, \end{gather}$$

where $A$, $B$, $C$ and $D$ are complex constants, and $J_1$ and $Y_1$ ($I_1$ and $K_1$) are the Bessel functions (modified Bessel functions) of order $1$. For the solution to be continuous and its derivative to have the correct behaviour at $\mathcal {z}=\mathcal {z}_{{c}}$, we must establish some relationships between these constants. Although this was not done explicitly by Alexakis et al. (Reference Alexakis, Young and Rosner2004b), we could use their final matching formulae to determine such relations to find them to be incompatible with the following:

(7.9)\begin{align} I_1(\mathcal{z})={-}{\rm i} J_1({\rm i} \mathcal{z}) \quad\text{and}\quad K_1(\mathcal{z})={-}\frac{\rm \pi}{2}[J_1(-{\rm i}\mathcal{z})-{\rm i} Y_1(-{\rm i} \mathcal{z})],\quad -{\rm \pi}<\arg(\mathcal{z})\le \frac{\rm \pi}{2}. \end{align}

We find that

(7.10)\begin{equation} C ={-}A+{\rm i} B\quad \text{and}\quad D ={-}\frac{2}{\rm \pi}\ B. \end{equation}

To check the consistency of these results, we use (7.10) for the numerical integration of the Rayleigh equation and retrieve the growth rates calculated by Beji & Nadaoka (Reference Beji and Nadaoka2004).