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.

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

I. 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 modeled as parallel flows of the form U = U (z) x, with z the vertical coordinate and x a horizontal unit vector (Figure 1).The stability of 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 [1] 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 [3].Nonetheless, the growth rate can be calculated analytically in only a few cases [4].
Miles [1] 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. [3].The experiments of Caulliez [5] showed that viscous damping is the main dissipation mechanism for waves shorter than 4 cm whereas at larger wavelengths the generation of capillary waves, microbreaking and breaking also contribute to dissipation.Recent fully coupled direct numerical simulations further demonstrated that the viscous stress plays a role on wave growth only in the case of strong wind forcing [6].Hence, the effect of viscosity is complex and has yet to be clarified [7,8].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 and Adam [9] 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. [10], and extended to smooth velocity profiles by Morland et al. [11].Young and Wolfe [2] 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 and Wolfe [2] for a review.Finally, Kadam et al. [12] gave an exact arXiv:2211.02942v3[physics.flu-dyn]25 Oct 2023 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, however, we obtain general results on both the Miles and rippling instabilities using asymptotic analysis.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. [6] 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, but it may 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 [13].For example, although most recent high-resolution three-dimensional simulations [14] 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. [15] and Alexakis et al. [16], 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 [4].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 et al. [17] argued that they are at the core of the astrophysical setting described in the previous paragraph.White dwarfs are extremely dense and dim objects, with many solar masses confined within Earth scale radii, and hence possess large gravitational forces.Indeed, Alexakis et al. [17] showed that for an exponential wind profile, the low wavenumber cut-off of the Miles instability is a growing function of gravity, and hence the low wavenumber cut-off 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 §II, 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 windinduced current in §III, 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 §IV, we treat the case of the wind and the current governed by an exponential profile and compare our asymptotic results to the exact eigenvalue.In §V, 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 §VI the asymptotic solution of the Rayleigh equation that was used in §III 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 §VII we demonstrate the limitations of the Wentzel-Kramers-Brillouin (WKB) approach to solving the Rayleigh equation.

II. LINEAR STABILITY OF A SHEARED TWO-FLUID INTERFACE
Young and Wolfe [2] derived the eigenvalue problem associated with the linear stability of an inviscid parallel shear flow across a two-fluid interface as sketched in Fig. 1.First we outline the governing Rayleigh equation and boundary conditions, after which we describe the non-dimensionalization of the problem.

A. 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 ≡ ρ a /ρ w < 1, and the gravitational acceleration and surface tension are g and σ respectively.Incompressibility ensures that a perturbation of the flow is entirely determined by the streamfunction ψ = ψ(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 ψ(z) [18].Note that the derivative U ′ (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.

B. 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 k ≡ kL w and z ≡ z/L w , respectively, and we define the dimensionless gravity and surface tension as From equation ( 4), the far field behavior of the streamfunction has the form e ±kz .Short waves are characterized by k ≫ 1, so the exponential decay is captured as follows, We introduce the small parameter ϵ ≡ 1/k and use (10) in the Rayleigh equation ( 2), which gives and We assume that ) the veracity of which we check a posteriori.The dimensionless form of equation ( 5) is For a given profile U(z), our main task is to solve equations (11) and ( 12) subject to the boundary conditions ( 13) and ( 14).
We consider the canonical situation where the wind blows over the water in which it induces a current (fig.1).Hence, U ′ z ≶ 0 > 0, U(z > 0) > 0, and U(z < 0) < 0. We explore another situation in §V.

C. Examples of profiles
A typical example of a wind and a wind-induced current is the double exponential profile, where U ∞ 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 Eq. ( 15) is Young and Wolfe [2] 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 [19], where u ⋆a and u ⋆w are the friction velocities of air and water, respectively, κ = 0.4 is the von Kármán constant, and z 0a and z 0w 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 (17) in the frame of the water surface is with R 1 = u ⋆a /u ⋆w and R 2 = z 0a /z 0w .In this case, exact analytical solutions are unknown.

III. SHORT WAVELENGTH EXPANSIONS AND EXPONENTIAL ASYMPTOTICS
We treat the eigenvalue problem described in §II A perturbatively, where the small parameter is the dimensionless inverse wavenumber, ϵ ≪ 1.We draw intuition for the approach from Miles [1], who considered a simpler version of our problem, with the small parameter r ≪ 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 [1], the leading order eigenvalue is unknown.In fact, even the nature of the lowest order behavior in ϵ is murky and hence we only assume that We use the notation of Bender and Orszag [20] for "C i (ϵ) is much smaller than C r (ϵ) as ϵ tends to 0", namely that The purpose of such an assumption is to replace C ∈ C by C r ∈ R in equations ( 11) and (12).Recalling that then there is a critical layer in the air (water), associated with a singular point in Eq. 11 (Eq.12).However, at this stage in the development, we do not know the possible values for C r .Indeed, the term rf ′ (0 + ) − h ′ (0 − ) depends on C in a non-trivial manner so that Eq. ( 14) is not necessarily quadratic in C. Physically, 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) C r and Young and Wolfe [2] 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, C r+ > 0 and C r− < 0 that correspond to two critical levels, z c+ > 0 and z c− < 0, such that Therefore, the prograde (retrograde) mode undergo the Miles (rippling) instability, and hence C i ̸ = 0 for both modes.We stress that, provided that the values of C r± are within the bounds of the function U(z), critical layers actually exist.When C r± are equal to these bounds, the system is marginally stable.In Appendix C, we derive a general asymptotic formula for the large wavenumber cut-off of the rippling instability.
Hence, although C r± are unknown, we replace C by C r+ in Eq. ( 11) and by C r− in Eq. (12).We then solve these equations using boundary layer theory in §VI, where we construct uniformly valid composite solutions.However, here we need only substitute the derivatives at z = 0 ± into Eq.( 14).In Appendix B, we show that for the prograde mode, and whereas for the retrograde mode, and We can now solve Eq. ( 14).With the advent of results ( 23) and ( 26), our key assumption, Eq. ( 20), 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 ϵ/C r and ϵ/z c .We seek solutions of Eq. ( 28) as a series in powers of ϵ.The nature of the leading order term depends on the value of the dimensionless surface tension 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 Eq. ( 30) can be derived by expanding the dispersion relation of interfacial capillary-gravity waves, Armed with C r± , we can use assumption (27), and Eqs. ( 23) and ( 26), in Eq. ( 14).We collect terms of order ϵ e −2|zc±|/ϵ to find the amplitudes A ± , and infer the growth rates of the prograde (+) and retrograde (-) modes as and where k = 1/ϵ.We emphasize several aspects of the present results.Firstly, equation (32) for the prograde mode is a generalization of the growth rate obtained by Miles in an Appendix of Morland and Saffman [21].That result was obtained from short wavelength asymptotic analysis of the exact solution of the Rayleigh equation for an exponential wind profile.Secondly, equation (33) for the retrograde mode is a generalization of the short wavelength limit of the growth rate of the rippling instability found by Young and Wolfe [2] 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 and Wolfe [2].Finally, we find that the results of Shrira [22] 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.

IV. INTERPRETATION
The behavior 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 z = 0 and at the critical levels, z = z c± .In consequence, the results are qualitatively the same for the two profiles introduced in §II C. We note that, as shown by Young and Wolfe [2], 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 G and 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 Eq. ( 31), where the plus (minus) sign corresponds to the prograde (retrograde) mode.Equations ( 29) and (30) predict that the non-uniformity of the current increases the phase speed by rU ′ (0 + )/2k and reduces it by U ′ (0 − )/2k respectively.Indeed, Figures 3(a)-(c) show that when r is close to 1, and k approaches unity, the phase speed of sheared waves is significantly larger than the phase speed of free surface waves.Moreover, Figures 4(a)-(c) show the predicted reduction of the retrograde mode in the absence of air when 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 Figures 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 G and S are small, which is the case for large velocity scales (c.f., Figs 2 and 4).Thus, consistent with physical intuition, a small perturbation grows faster in the presence of strong winds and/or currents.

V. 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 U z ≶ 0 > 0, U ′ (z > 0) > 0, and U ′ (z < 0) < 0, and again assume an infinite depth and proceed as in §III.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 z c± such that and there is no value of z such that C r− < 0 is equal to U(z).This implies that and  for the prograde mode, while and 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 C r± are still given by equations ( 29) and (30).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 refinements of the viscosity-stratified approach of Charles and Lilleleht [23] or the two-layer Couette flow approach of Barthelet et al.

VI. ASYMPTOTIC SOLUTION OF THE RAYLEIGH EQUATION FOR SHORT WAVES
Here we solve equation (11) for C = C r+ and note that a similar procedure is applicable to equation ( 12) when C = C r− .For simplicity, we rewrite equation ( 11) and we drop the subscript + for the rest of this section.
There is a regular singularity at z = z c such that Because the small parameter ϵ multiplies the highest order derivative in equation ( 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 equation ( 11) is finite at z = z c whereas its derivative has a logarithmic divergence [18].Therefore, we assume that an internal boundary layer emerges from the singularity.We assume that the point z = z c (ϵ) is well separated from the lower boundary, z = 0, as ϵ goes to zero.We check this a posteriori once the dependence of C r on ϵ is known.
In consequence, C r can be treated as a constant in the following analysis.
The boundary layer is an inner region where the solution of equation ( 11) changes rapidly.We define two outer regions, where the solution changes slowly.One spans 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).

A. Outer solutions
Following Miles [26], we seek outer solutions as a power series in ϵ as and find that (42a, b) The constant f 0 and the lower limit of the integral 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.

B. Inner solution
Within the boundary layer, we introduce the stretched coordinate Z ≡ (z − z c )/δ, where 0 < δ ≪ 1, and seek an inner solution, F in (Z), to equation (11).We approximate the coefficients by their Taylor series expansions about z = z c , so that equation (11) becomes ) where the subscript 'c' denotes evaluation at the critical level, z = z c .By balancing the two first terms, we obtain the distinguished limit δ = ϵ [20], and hence must solve Because the boundary layer has an O(ϵ) thickness, the appropriate inner expansion is The general solutions are and Next we determine the integration constants, A and B, and the bounds of integration, a and b, by asymptotic matching.

C. Asymptotic matching and uniformly valid composite solutions
We must match the inner solution and the two outer solutions, as ϵ → 0 + .We use the superscripts 'ℓ' and 'u' for lower and upper, respectively.
The lower outer solution, f ℓ out (z), is defined for 0 ≤ z ≪ z c , and must satisfy the lower boundary condition; The upper outer solution, f u out (z), is defined for z ≫ z c and hence must satisfy condition (13a).We use the following matching conditions, and lim and then apply the Van Dyke additive rule [20] 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.

Leading order
To leading order, the outer solutions are constant (see Eq. 42a), and because of the boundary condition (48a) we have The leading order inner solution is given by equation (46).To preempt divergence as Z → +∞, we impose A = 0, and the matching condition (49) implies that B = 1.Therefore, which, upon imposition of the matching condition (50), yields We combine the solutions (52), ( 53) and (54) using the additive rule (51), and thus obtain a uniformly valid composite solution at leading order; The effect of the boundary layer appears only at the next order.

Order ϵ
Following equation (42b), the lower and upper outer solutions at order ϵ are and respectively.From the Laurent series expansion we deduce that and where Log denotes a continuation of the natural logarithm to the negative real numbers, (61) The choice of the branch cut -just above the negative real axis if U ′ c > 0, just below otherwise -follows from Lin [27].
From equation (47), the inner solution at order ϵ is We split the integral over t as Noting that the second integral is a constant, C ∈ C, we obtain (64) where is the exponential integral.For large values of |x|, the divergent series [20] After integration, we find (68) Recalling that Z = (z − z c )/ϵ, the matching of the inner limits, given by equations ( 59) and (60), and the outer limits, given by equation (68), yields Because the boundary condition (48b) requires that the outer solution at order ϵ has the same expression in the lower and upper outer regions; (71) The lower outer solution is real, because the path of integration along the real axis does not reach the branch point z = 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 U ′ c > 0, and in the lower part otherwise.In Appendix B, we show that Using the matching conditions (69) together with equation (70), the inner solution (64) becomes Finally, we use the Van Dyke rule (51) to construct a uniformly valid composite solution at order ϵ; In Appendix B, we verify that χ(z) = e −kz f unif,1 (z) satisfies the global property; at leading order.In Appendix A we show that Equation ( 76) generalizes a result that Miles derived in an Appendix to Morland and Saffman [21] using an exact solution of the Rayleigh equation for an exponential wind profile.

D. Comparison with the numerical solution
We numerically solve equation (11) for the two standard wind profiles where, as before, κ = 0.4 is the von Kármán constant.We compare the numerical solution with our uniformly valid composite solution (74) for fixed values of k and 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 [18], the solution within the inner layer depends on the wind profile and the dispersion relation solely through the scale factor U ′′ c /U ′ c and the bound of integration z c /ϵ.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 k ≪ 1 the phase varies from z = 0 to z = z c , so that long waves interact with the wind all the way from the mean water surface to the critical level [4].

E. Similarity solution
We have non-dimensionalized the variables using external parameters characterizing the shear in the air.However, the general short wave solution (74) can also be written in terms of ξ ≡ kz and U ≡ U/c r , where cr is the dimensional version of C r .The wind-generated ripples have a continuous wavenumber spectrum so that k, and subsequently cr , are variables rather than parameters.In that sense, ξ is a similarity variable introduced by Miles [1], and the Rayleigh equation has a self-similar solution, ϕ = ϕ(ξ), which for short waves is

VII. THE WKB METHOD IN THE SHORTWAVELENGTH LIMIT
Alexakis et al. [17] proposed a solution of the Rayleigh equation for short waves using the WKB method.However, because their solution does not satisfy the global FIG.8: Comparison of the growth rate of gravity waves obtained by et al. [17] with the predictions of our short wavelength asymptotic treatment (cf.Eq. 32) and the growth rate calculated from the exact (hypergeometric) solution for a wind profile U (z) = U ∞ (1 − e −z/da ) and a density ratio r = 0.5.
property (75), their growth rate has an extra factor π/k; see their equation (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. [17] with our short wave asymptotic solution and the exact solution.Clearly the WKB method deviates from the others for all kd a ≳ 25, and in what follows we explain the origin of the deviation.Setting ϵ = 1/k, we write the Rayleigh equation in the form A WKB expansion is [20] In the standard case, Q(z, ϵ) does not depend on the small parameter ϵ.Then, we can take only two terms in the series (80): this, physical optics approximation, yields for Q(z) > 0, and for Q(z) < 0, where the bounds of integration, a and b, are arbitrary.
For U ′′ (z) < 0 and Q(z, ϵ) introduced in equation (79), Alexakis et al. [17] used the physical optics approximation on three intervals, defined in Figure 9, and obtained three solutions of the form of (81) or (82), in which they neglected the fact that ϵ ≪ 1 within Q(z, ϵ).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 z ⋆ [e.g.20].However, they needed inner solutions near the critical level, 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 z > z c , and in terms of modified Bessel functions for z < 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 z c and the turning point FIG.9: The three intervals in the WKB of Alexakis et al. [17] where z ⋆ is a simple turning point.
z ⋆ actually shrinks as ϵ tends to zero.For instance, in the case of the exponential profile U(z) = 1 − e −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 behavior.
In the vicinity of the critical level, 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 behavior at z = z c , we must establish some relationships between these constants.Although this was not done explicitly by Alexakis et al. [17], we could use their final matching formulae to determine such relations to find them to be incompatible with the following; ) for −π < arg(z) ≤ π/2.We find that To check the consistency of these results, we use equation (88) for the numerical integration of the Rayleigh equation and retrieve the growth rates calculated by Beji and Nadaoka [28].

VIII. 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, ϵ = 1/k, we asymptotically solved the eigenvalue problem for the stability of an arbitrary parallel flow U = U (z) 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 ϵ.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 § V 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 ± ; 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 ϵ 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 signifiant 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 of solving the Rayleigh equation breaks down.

FIG. 2 :FIG. 3 :
FIG. 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 ∞ (1 − e −z/da ), a density ratio r = 0.001, and different values (grey and black) of the dimensionless gravity, G = gd a /U 2 ∞ , and surface tension, S = σ/(ρ w U 2 ∞ d 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).

FIG. 4 :FIG. 5 :
FIG. 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 e z/dw , and different values (grey and black) of the dimensionless gravity, G = gd w /U 2 S , and surface tension, S = σ/(ρ w U 2 S d 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).

FIG. 6 :
FIG. 6: Schematic of the domain of analysis of the Rayleigh equation for short waves; k ≫ 1.A boundary layer of thickness ϵ = 1/k emerges from the singularity at the critical level, z = z c .

FIG. 7 :
FIG. 7: of the uniformly valid composite solution (74) with the numerical solution of the Rayleigh equation for the exponential (a,b) and logarithmic (c,d) wind profiles, with C r = 0.1 and different values of 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 (74) and the dashed line the imaginary part.