Flow-driven interfacial waves: an inviscid asymptotic study

Abstract

Motivated by wind blowing over water, we use asymptotic methods to study the evolution of short wavelength interfacial waves driven by the combined action of these flows. We solve the Rayleigh equation for the stability of the shear flow, and construct a uniformly valid approximation for the perturbed streamfunction, or eigenfunction. We then expand the real part of the eigenvalue, the phase speed, in a power series of the inverse wavenumber and show that the imaginary part is exponentially small. We give expressions for the growth rates of the Miles (J. Fluid Mech., vol. 3, 1957, pp. 185–204) and rippling (e.g. Young & Wolfe, J. Fluid Mech., vol. 739, 2014, pp. 276–307) instabilities that are valid for an arbitrary shear flow. The accuracy of the results is demonstrated by a comparison with the exact solution of the eigenvalue problem in the case when both the wind and the current have an exponential profile.

1. Introduction

Waves at the interface between two fluids with different densities are ubiquitous in nature. A natural question concerns how a flow in either fluid affects these waves. Here, we consider fluid layers of infinite extent. A canonical example is the wind flowing over the ocean, both of which can be modelled as parallel flows of the form $\boldsymbol {U}=U(z)\ \boldsymbol {\hat {x}}$, with $z$ the vertical coordinate and $\boldsymbol {\hat {x}}$ a horizontal unit vector (figure 1). The stability of the arbitrary shear flow $\boldsymbol {U}$ under the influence of small two-dimensional perturbations has been studied extensively over the last seventy years. In the absence of a current in the water, Miles (1957) found an instability of the wind, leading to the growth of water waves. The theory of Miles is inviscid and assumes that the wind and the waves are weakly coupled because the air/water density ratio, $r$, is small. Hence, the wind transfers energy to the waves within a critical layer in the air, where the wind speed equals the phase speed of free surface waves. Recent laboratory measurements of the air flow over wind-generated waves provide evidence of this critical layer mechanism (Carpenter, Buckley & Veron 2022). Nonetheless, the growth rate can be calculated analytically in only a few cases (Bonfils et al. 2022).

Figure 1. Schematic of wind driving a current in the water. Both are modelled as a parallel shear flow $U=U(z)$ and perturbed by a wave with wavenumber $k$.

Miles (1957) argued that the critical layer should be above the viscous sublayer, and hence neglected viscosity. This rationale is now supported by the measurements of Carpenter et al. (2022). The experiments of Caulliez (2013) showed that viscous damping is the main dissipation mechanism for waves shorter than 4 cm, whereas, at larger wavelengths, the generation of capillary waves, micro-breaking and breaking also contribute to dissipation. Recent fully coupled direct numerical simulations further demonstrated that viscous stress plays a role in wave growth only in the case of strong wind forcing (Wu, Popinet & Deike 2022). Hence, the effect of viscosity is complex and has yet to be clarified (Zeisel, Stiassnie & Agnon 2008; Wu & Deike 2021). Thus, simplified inviscid models still provide valuable insights, as shown here.

When there is a laminar current in the water, the position of the critical layer responsible for the Miles instability is unknown, because the phase speed is itself unknown. This is associated with the fact that the surface waves are no longer free in the sense that the current modifies their dispersion relation in a non-trivial manner. Water currents are often wind induced and decay with depth. Stern & Adam (1974) were the first to suggest that, in such cases, sheared surface waves propagating slower than the water surface, which is dragged by the wind, undergo an instability. They considered a current with a broken-line velocity profile, a model further studied by Caponi et al. (1991), and extended to smooth velocity profiles by Morland, Saffman & Yuen (1992). Young & Wolfe (2014) showed that there is another critical layer in the water, where the unknown phase speed of the waves matches the speed of the current. They referred to this phenomenon as the ‘rippling instability’. Analytical progress on the rippling instability in deep water has only been made for piece-wise linear or exponential currents; see Young & Wolfe (2014) for a review. Finally, Kadam, Patibandla & Roy (2023) gave an exact analytical treatment of the stability of an exponential wind profile over a finite-depth water layer that is either quiescent or has linear or quadratic current profiles. For the quadratic current, they introduced spheroidal wave functions to assess the stability of a shear flow. Here, we use asymptotic methods to obtain general results on both the Miles and the rippling instabilities. The fully coupled numerical approach of Wu et al. (2022) resolved the development of the shear-induced drift layer beneath the water surface as well as the evolution of the air-side turbulent boundary layer. They showed that, in strongly forced cases, the flows are transient, in the sense that the waves feedback on the air flow while the water flow becomes turbulent. Although transient effects are beyond the scope of this work, we can explicitly account for the effect of turbulence using mean profiles.

In the absence of a current in the water and for $r$ less than unity, and yet not small, the Miles instability is difficult to study. Indeed, the waves are then strongly coupled to the wind so that the shear substantially changes the dispersion relation, yielding yet another situation in which the position of the critical layer is unknown. Not only is such a strong wind–wave coupling a central process on Earth, it may also play an important role in a number of astrophysical settings, including white dwarfs, one of the end products of stellar evolution, neutron stars and black holes (Shapiro & Teukolsky 2008). For example, although most recent high-resolution three-dimensional simulations (Casanova et al. 2011) suggest that Kelvin–Helmholtz instabilities driven by buoyant fingering may be able to explain the composition of white dwarfs, an alternative proposal by Rosner et al. (2001) and Alexakis et al. (2004a) is that the Miles instability is responsible.

Because most ocean waves have wavelengths much larger than the characteristic length scale of the wind profile, the Miles instability can be treated using asymptotic methods in the long-wave limit (Bonfils et al. 2022). Here, we focus on short waves with the goal of capturing the combined influences of an underlying current and a moderate density ratio. Whereas short waves may not be central in terrestrial oceanography, Alexakis, Young & Rosner (2004b) argued that they are at the core of the astrophysical setting described in the previous paragraph. White dwarfs are extremely dim and dense, with approximately one solar mass confined within an Earth scale radius, and hence possess large gravitational forces. Indeed, Alexakis et al. (2004b) showed that, for an exponential wind profile, the low wavenumber cutoff of the Miles instability is a growing function of gravity, and hence the low wavenumber cutoff of the Miles instability is in fact very large. Thus, the growing waves at the surface of white dwarfs are short with respect to astrophysical scales.

Therefore, rather than solving a particular geophysical or astrophysical problem, we treat a basic fluid mechanical question: the stability of a sheared two-fluid interface where the upper fluid is less dense than the lower fluid. For clarity of discussion, we refer to the upper and lower fluids as air and water, respectively, and refer to wind as the shear flow in the air, and current as the shear flow in the water.

In § 2, we describe the linear stability analysis of a sheared two-fluid interface. The eigenfunction satisfies the Rayleigh equation and the eigenvalue is a complex intrinsic phase speed. We focus on the case of a wind-induced current in § 3, and obtain the real part of the eigenvalue as a series in powers of the dimensionless inverse wavenumber. Moreover, we show that the imaginary part of the eigenvalue is small and obtain general formulae for the growth rates of the Miles and the rippling instabilities. In § 4, we treat the case of the wind and the current governed by an exponential profile and compare our asymptotic results with the exact eigenvalue. In § 5, we treat another type of current wherein a mode can have two critical layers, one in the air and one in the water. We then derive in § 6 the asymptotic solution of the Rayleigh equation that was used in § 3 for the calculation of the eigenvalue. In particular, we show that an internal boundary layer emerges from the singularity at the critical level. Finally, before concluding, in § 7 we demonstrate the limitations of the Wentzel–Kramers–Brillouin (WKB) approach to solving the Rayleigh equation.

2. Linear stability of a sheared two-fluid interface

Young & Wolfe (2014) 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 1981). 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 2002). 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 (2014) 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 1975), 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 (1957), 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 (1957), 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 (1999) 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 (2014) 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 (1992). 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 (2014) 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 (2014). Finally, we find that the results of Shrira (1993) 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 (2014), 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 (1965) or the two-layer Couette flow approach of Barthelet, Charru & Fabre (1995), as well as in the context of Holmboe wave experiments (e.g. Carpenter et al. 2010).

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 1981). 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 (1962), 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 1999), 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 1999) 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 (1955).

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 1999)

(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 (1992) 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 1981), 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. 2022).

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 (1957), 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. (2004b) 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. (2004b) 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. (2004b) 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 1999)

(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. (2004b) 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 1999). 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. (2004b) 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. (2004b), 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 (2004).

8. Conclusion

We have studied the effect of a shear flow on interfacial capillary–gravity waves when their wavelength is much smaller than the characteristic length scale of the flows in the fluids bounding the interface. Using the dimensionless inverse wavenumber as a small parameter, $\epsilon = 1/\mathcal {k}$, we asymptotically solved the eigenvalue problem for the stability of an arbitrary parallel flow $\boldsymbol {U}=U(z) \boldsymbol {\hat {x}}$ through a two-fluid interface with a density ratio that is not necessarily small. We constructed uniformly valid composite solutions of the governing Rayleigh equation, where the real part of the eigenvalue is a power series in $\epsilon$. We showed that including the effect of surface tension changes the nature of the leading-order solution, and that exponentially small terms must be considered in order to have a non-zero imaginary part of the solution. From a physical view point, there is a prograde mode whose phase speed is greater than the speed of the water surface, $U_{{S}}$, and a retrograde mode whose phase speed is smaller than $U_{{S}}$. When the velocity of the shear flow equals the phase speed of one these modes, there is a critical layer where the flow transfers energy to the wave. If the critical layer is in the air (water), this is called the Miles (rippling) instability. The only case considered in the literature thus far is when the prograde mode undergoes the Miles instability and the retrograde mode undergoes the rippling instability. In § 5 we studied the situation where the prograde mode can undergo both instabilities and the retrograde mode is neutral. In the short-wave limit, we found that (i) the effect of the shear on the phase speed of the two modes depends only on the derivatives of $U$ at $z=0^\pm$; and (ii) the Miles and rippling instabilities have a growth rate of the same form. Indeed, the interaction between the shear flow and the waves is mostly reduced to a narrow region around the critical level, an internal boundary layer of thickness $\epsilon$ where the solution of the Rayleigh equation has a self-similar structure. Heuristically speaking, the waves are barely influenced by the flow outside of this region. Nonetheless, we showed that there are significant effects on dispersion and growth rates when the characteristic velocity of the shear flow is large and the density ratio is close to $1$. Finally, we showed how the WKB approach to solving the Rayleigh equation breaks down.