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 (Reference Miles1957) 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 Reference Carpenter, Buckley and Veron2022). Nonetheless, the growth rate can be calculated analytically in only a few cases (Bonfils et al. Reference Bonfils, Mitra, Moon and Wettlaufer2022).
Miles (Reference Miles1957) 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. (Reference Carpenter, Buckley and Veron2022). The experiments of Caulliez (Reference Caulliez2013) 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 Reference Wu, Popinet and Deike2022). Hence, the effect of viscosity is complex and has yet to be clarified (Zeisel, Stiassnie & Agnon Reference Zeisel, Stiassnie and Agnon2008; Wu & Deike Reference Wu and Deike2021). 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 (Reference Stern and Adam1974) 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. (Reference Caponi, Yuen, Milinazzo and Saffman1991), and extended to smooth velocity profiles by Morland, Saffman & Yuen (Reference Morland, Saffman and Yuen1992). Young & Wolfe (Reference Young and Wolfe2014) 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 (Reference Young and Wolfe2014) for a review. Finally, Kadam, Patibandla & Roy (Reference Kadam, Patibandla and Roy2023) 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. (Reference Wu, Popinet and Deike2022) 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 Reference Shapiro and Teukolsky2008). For example, although most recent high-resolution three-dimensional simulations (Casanova et al. Reference Casanova, José, Garcia-Berro, Shore and Calder2011) 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. (Reference Rosner, Alexakis, Young, Truran and Hillebrandt2001) and Alexakis et al. (Reference Alexakis2004a) 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. Reference Bonfils, Mitra, Moon and Wettlaufer2022). 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 (Reference Alexakis, Young and Rosner2004b) 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. (Reference Alexakis, Young and Rosner2004b) 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 (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
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
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
we can impose an exponential decay, viz.
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
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
act as additional control parameters. In a frame moving at speed $U_{{S}}$, we use the dimensionless profile
and the dimensionless intrinsic phase speed is
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
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:
We introduce the small parameter $\epsilon \equiv 1/\mathcal {k}$ and use (2.10) in the Rayleigh equation (2.2), which gives
We assume that
the veracity of which we check a posteriori. The dimensionless form of (2.5) is
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
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
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
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
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
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
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
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
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,
whereas, for the retrograde mode,
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
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
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
When capillary forces are included, we instead obtain
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
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
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.
The phase speed of Doppler-shifted interfacial waves in presence of a constant, zero shear current is
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
and there is no value of $\mathcal {z}$ such that $\mathcal {C}_{r-}<0$ is equal to $\mathcal {U}(\mathcal {z})$. This implies that
for the prograde mode, while
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
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).
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):
and we drop the subscript $+$ for the rest of this section. There is a regular singularity at $\mathcal {z}=\mathcal {z}_{{c}}$ such that
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).
6.1. Outer solutions
Following Miles (Reference Miles1962), we seek outer solutions as a power series in $\epsilon$ as
and find that
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
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
Because the boundary layer has an $O(\epsilon )$ thickness, the appropriate inner expansion is
The general solutions are
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
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:
and then apply the Van Dyke additive rule (Bender & Orszag Reference Bender and Orszag1999) to construct uniformly valid composite solutions
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
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
which, upon imposition of the matching condition (6.11b), yields
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
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
respectively. From the Laurent series expansion
we deduce that
where $\textrm {Log}$ denotes a continuation of the natural logarithm to the negative real numbers
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
We split the integral over $t$ as
Noting that the second integral is a constant, $C\in \mathbb {C}$, we obtain
where
is the exponential integral. For large values of $|x|$, the divergent series (Bender & Orszag Reference Bender and Orszag1999)
yields
After integration, we find
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
Because the boundary condition (6.10b) requires that
the outer solution at order $\epsilon$ has the same expression in the lower and upper outer regions
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
Using the matching conditions (6.29a,b) together with (6.30), the inner solution (6.24) becomes
Finally, we use the Van Dyke rule (6.12) to construct a uniformly valid composite solution at order $\epsilon$
In Appendix B, we verify that $\chi (\mathcal {z}) = {\rm e}^{-\mathcal {k}\mathcal {z}}f_{{unif},1}(\mathcal {z})$ satisfies the global property:
at leading order. In Appendix A we show that
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
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).
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
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.
Setting $\epsilon =1/\mathcal {k}$, we write the Rayleigh equation in the form
A WKB expansion is (Bender & Orszag Reference Bender and Orszag1999)
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
for $Q(\mathcal {z})>0$, and
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
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.
In the vicinity of the critical level, $\mathcal {z}_{{c}}$, we introduce the inner variable
in terms of which the inner solutions are
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:
We find that
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).
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.
Acknowledgements
We thank C. Arratia for suggesting the possibility of a prograde mode with two critical layers, developed in § 5.
Funding
A.F.B. thanks N. Balmforth and R. Rosner for discussions of the WKB approach and the 2022 Geophysical Fluid Dynamics Summer Study Program at the Woods Hole Oceanographic Institution, which is supported by the National Science Foundation and the Office of Naval Research under no. OCE-1332750. All authors acknowledge Swedish Research Council grant no. 638-2013-9243. Nordita is supported in part by NordForsk.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Solution of (2.11) at the singular point $\mathcal {z}=\mathcal {z}_{{c}}$
The regular singular point $\mathcal {z}=\mathcal {z}_{{c}}$ is located in the inner region where the solution at order $\epsilon$ is given by (6.33), which we reproduce here:
is the inner variable.
A.1. Imaginary part
From the series representation of the exponential integral
where $\gamma _{{E}}=0.577$ is the Euler constant, we see that $x=0$ is a logarithmic branch point for the integral in (A1). The only contribution to the imaginary part of $F_1$ arises from the integration of $\textrm {Log} (x)$ along the branch cut. Namely, if $\varTheta$ is the Heaviside step function, then
where we select the plus sign when $\mathcal{U}^{\prime}_{{c}} >0$, and the branch cut is just above the negative real axis, and the minus sign otherwise. From (A3), we find
Then, upon taking the limit $Z\to 0$, we obtain
It is interesting to take the outer limits $Z\to \pm \infty$ of the expression (A4). In particular, if the matching was done properly, these should be equal to the inner limits of the imaginary parts of the lower and upper outer solutions, viz.
In each case, we shall check that
which is the condition of separation of the boundary layer of thickness $\epsilon$ from the lower boundary. Thus, we consistently retrieve a zero imaginary part in the lower outer region. We recall that the upper outer solution is given by (6.31)
According to the Sokhotski–Plemelj theorem
where ${\unicode{x01A4}}$ denotes the Cauchy principal value. Hence, the imaginary part in the upper outer region is constant, equal to the upper outer limit of ${\rm Im}\{F_1\}$ up to exponentially small terms.
A.2. Real part
At the leading order, we have
Here, we calculate the next-order contribution to the real part. For convenience, we define
and we express the exponential integral using the Ramanujan series as
where $\lfloor x \rfloor$ denotes the floor function. This series converges faster than (A2). We perform a direct integration as
and after some algebra, we obtain
with
in which
is the upper incomplete gamma function. Note that
and that
Therefore
We define the ‘remainder’ as
which is a negligible quantity for $\zeta \gg 1$, as shown in figure 10. Discarding the exponentially small terms, we arrive at
Appendix B. Global property of the solution of the Rayleigh equation
The dimensionless Rayleigh equation is
where $\mathcal {C}_r=\mathcal {C}_r(\mathcal {k})$ is a known function. Assuming that there is a unique level $\mathcal {z}_{{c}}>0$ such that $\mathcal {U}(\mathcal {z}_{{c}})=\mathcal {C}_r$, Miles (Reference Miles1957) showed that
where the subscript ‘$c$’ denotes evaluation at $\mathcal {z}=\mathcal {z}_{{c}}$. The derivation of the global property (B2) is analogous to the canonical derivation of Rayleigh's inflection point theorem as follows. Multiply the Rayleigh equation (B1) by the complex conjugate of $\chi (\mathcal {z})$, integrate by parts, use the boundary conditions to evaluate the integrated term and take the imaginary part; the result follows from the Sokhotski–Plemelj theorem.
Here, we show that our asymptotic solution for short waves satisfies (B2). For $\mathcal {z}>0$, we let
We take the derivative of the composite solution, $f_{{unif},1}(\mathcal {z})$, obtained in (6.34) and obtain
where we select $+ {\rm i}{\rm \pi}$ if $\mathcal{U}^{\prime}_{{c}} >0$, and $- {\rm i}{\rm \pi}$ otherwise. The contribution to the real part of the exponential integral multiplied by the decaying exponential is negligible, and hence is discarded in (3.4) and (3.6).
We use (B4) to calculate the left-hand side of the identity (B2)
Using the results of Appendix A, we have
Using (B6) on the right-hand side of (B2), we recover the left-hand side of (B5).
Appendix C. Short wavelength cutoff of the rippling instability
As shown in figure 4, the rippling instability has a short wavelength or high wavenumber cutoff when the effect of surface tension is taken into account. Following Young & Wolfe (Reference Young and Wolfe2014), we denote this cutoff $k^+_{{neut}}$, where the subscript ‘neut’ denotes neutral. In other words, the wave with wavenumber $k^+_{{neut}}$ is marginally stable because its phase speed equals the lower bound of the velocity profile in the water, namely
where all variables are dimensional. Hence, using (3.11) for the phase speed of short waves, $k^+_{{neut}}L\equiv \mathcal {k}^+_{{neut}}$ is solution of
For convenience, we introduce the following notation:
We assume a wind-induced current in the water so that $\mathcal {B}>0$. Following Young & Wolfe (Reference Young and Wolfe2014), we let the inverse Weber number be a small parameter, $\tilde {\mathcal {S}}\ll 1$, and solve (C2) in the form
Taking the square of this result and letting $X\equiv \tilde {\mathcal {S}}\mathcal {k}$ we obtain
Equation (C5) can be readily solved using regular perturbation theory. We seek solutions in the form
We find
and infer
For the double exponential profile introduced in § 2.3 we find
which is a generalization of (5.10) in Young & Wolfe (Reference Young and Wolfe2014).