On the intermediate long wave propagation for internal waves in the presence of currents

A model for the wave motion of an internal wave in the presence of current in the case of intermediate long wave approximation is studied. The lower layer is considerably deeper, with a higher density than the upper layer. The flat surface approximation is assumed. The fluids are incompressible and inviscid. The model equations are obtained from the Hamiltonian formulation of the dynamics in the presence of a depth-varying current. It is shown that an appropriate scaling leads to the integrable Intermediate Long Wave Equation (ILWE). Two limits of the ILWE leading to the integrable Benjamin-Ono and KdV equations are presented as well.


Introduction
The equatorial region in the Pacific is characterised by a shallow layer of warm and less dense water over a much deeper layer of cold denser water.The two layers are separated by a sharp thermocline (where the temperature gradient has a maximum, it is very close to the pycnocline, where the pressure gradient has a maximum) at a depth, depending on the location, but usually at 100 -200 m beneath the surface.Both layers are supposed to be homogeneous and their sharp boundary is the thermocline/pycnocline.The strong stratification confines the wind-driven currents to a shallow nearsurface region.In the Atlantic and Pacific, the westward trade winds induce a westward surface flow at speeds of 25-75 cm/s, while a jet-like current, the Equatorial Undercurrent (EUC), flows below it toward the East (counter to the surface current), attaining speeds of more than 1 m/s at a depth of nearly 100 m.Since the wind usually moves the top 20 m of the sea surface, and that the wave speed is usually less than 4 m/s, the vorticity values are not exceeding by magnitude 0.2 rad/s.The values for the lower layer are typically lower.
The EUC flows in a region that is in a very narrow region close to the Equator, it is symmetric about the Equator and extends nearly across the whole length (more than 12000 km) of the Pacific Ocean basin [32,26].With speeds in excess of 1 m/s, the EUC is one of the fastest permanent currents in the world.The flow has nearly two-dimensional character, with small meridional variations, being combinations of longitudinal non-uniform currents and waves, and presenting a significant fluid stratification that results in a pycnocline/thermocline separating two internal layers of practically constant density (see [36]).While at depths in excess of about 240 m there is, essentially, an abyssal layer of still water, the ocean dynamics near the surface is quite complex.In this region the wave motion typically comprises surface gravity waves with amplitudes of 1-2 m and oscillations with an amplitude of 10-20 m at the thermocline (of mean depth between 50 m and 150 m).These internal waves interact with the underlying current (EUC).The topic of wave-current interactions is a large and important research area within the fluid dynamics and we refer to [14,15,16,17,22,29,30,35,48,53,54] and the references therein.
Our goal is to model the wave propagation at the thermocline in the case of intermediate long wave approximation, when the wavelength is comparable to the depth of the lower layer, which in turn is much deeper than the upper layer.We obtain several approximate models and one of them is an equation, integrable by the inverse scattering method.
The Hamiltonian approach is central to our modelling of internal waves in the presence of current.It originates from Zakharov's paper [57] for irrotational surface waves over infinitely deep water.The Hamiltonian formalism is often utilized in the study of nonlinear waves in continuous media, see for example the review article [58].The Lagrangian and the Hamiltonian approaches for surface and internal waves have been developed further in many papers, from which here we mention only [3,4,5,23,24,25,44,45,47]. Currents with linear profiles and constant vorticity could be included in the Hamiltonian formulation of the wave-current dynamics, recent studies are published in [10,11,12,13,18,19,20,31,33].
The setup and the governing equations for internal waves are described briefly in Sections 2 and 3.The Hamiltonian Formulation of the problem is presented in Section 4. The scales leading to the intermediate long wave propagation regime are introduced in Section 5 where the model equations are derived as well.The integrable Intermediate Long Wave Equation (ILWE) is derived in Section 6. Mathematical details about some integro-differential operators and the integrability of ILWE are given in the Appendix.

System setup for internal equatorial waves
The internal wave along the equator is modelled by two layers of water in two dimensions, since the meridional motion is neglected.The layers are separated by a free common interface (the thermocline/pycnocline) as per Figure 1.
The system is bounded at the bottom by an impermeable flatbed and is considered as being bounded on the surface by an assumption of absence of surface motion, i.e. rigid lid approximation.The two fluid domains are The functions and parameters associated with the upper layer will be marked with subscript 1.Also, subscript c (implying common interface) will be used to denote evaluation on the internal wave.Propagation of the internal wave is assumed to be in the positive x-direction which is considered to be 'eastward'.The direction of the gravity force is in the negative y-axis.
The function η(x, t) describes the elevation of the internal wave with the mean of η assumed to be zero, R η(x, t)dx = 0. Actually it is sufficient R η(x, t)dx = p < ∞, because the mean value involves division of the length of the corresponding interval, which in this case is R dx = ∞.The fluids are incompressible with densities ρ and ρ 1 .The stable stratification is given by the immiscibility condition ρ > ρ 1 .
The stream functions, ψ and ψ 1 , are related to the velocity fields u = (u, v) and u 1 = (u 1 , v 1 ) via the relations due to the incompressibility assumption ∇ • u = 0, ∇ • u 1 = 0.The velocity potentials, ϕ and ϕ 1 , are introduced such that where γ = −v x + u y and γ 1 = −v 1,x + u 1,y are the constant vorticities.This setup allows for modelling of an undercurrent, such as the Equatorial Undercurrent.A piecewise linear current profile can be represented by the velocity fields of the form (2), [18] by writing where κ and κ 1 are constants representing the current horizontal velocities at y = 0.The wave-only components have been separated out by introducing a tilde notation.Effectively the current profiles in Ω and Ω 1 are U(y) = γy + κ and U 1 (y) = γ 1 y + κ 1 , correspondingly.Here κ and κ 1 are constants, representing the average flow in each respective domain, see [20,19], where the situation with a free surface is studied.We assume that the functions η(x, t), ϕ(x, y, t) and ϕ 1 (x, y, t) belong to the Schwartz class S(R) with respect to x (for any y and t).The assumption of course implies that for large absolute values of x the internal wave attenuates, and is vanishing at infinity, and therefore

Governing equations
The Euler equations for the two layers are where are the Coriolis forces per unit mass at the equator with ω being the rotational speed of the Earth, g = (0, 0, −g) is the Earth acceleration, i.e. g is the magnitude of the acceleration due to gravity, ρ and ρ 1 (due to the assumption of incompressibility) are the constant densities and p and p 1 are the corresponding pressures.The pressure gradients can be expressed as The dynamic boundary condition at the interface p = p 1 gives the Bernoulli equation where the subscript c is a notation for an evaluation at the common interface y = η(x, t), χ = ψ(x, η, t) and χ 1 = ψ 1 (x, η, t).The equation (7) gives the time evolution of ξ := ρ( ϕ) c − ρ 1 ( ϕ 1 ) c which is the momentum type variable in the Hamiltonian formulation of the problem.The coordinate variable is η(x, t) and it evolves according to the kinematic boundary condition at the interface [34] η This can be expressed in terms of the velocity potentials, using (3), as The kinematic boundary condition at the bottom, requiring that there is no velocity component in the y-direction on the flat bed, is given by and, additionally, there is a kinematic boundary condition at the top, requiring that there is no velocity component in the y-direction on the surface, given by ϕ 1 (x, h 1 , t) y = 0 and

Hamiltonian formulation
Similar to the case of a single layer with a constant vorticity [21], the wavemotion on the interface (thermocline/pycnocline) admits a (nearly) Hamiltonian representation.The details could be found for example in [8,9,11].Let us introduce the following notations for the velocity potentials values at the interface The variable introduced in [3,4] plays the role of a generalised momentum.We assume that ξ ∈ S(R) for all t.The (non-canonical) equations of motion are where H = H(η, ξ) is the total energy of the system, and is the stream function, evaluated at y = η(x, t).
The dynamics can be formally written in the canonical form under the transformation (cf.[55,56]) The equations ( 14) could be further expressed in terms of the variable u = ξ x for a Hamiltonian written in terms of u and η : In order to obtain the Hamiltonian in terms of variables defined at the interface we need to introduce the so-called Dirichlet-Neumann (DN) operators [24,25] G(η)φ = ( φn where ϕ n and ϕ 1,n 1 are the normal derivatives of the velocity potentials, at the interface, and n and n 1 are the outward normals to the corresponding domains, with n = −n 1 .We need also the operator Using the boundary conditions we obtain With (13) we obtain and we can solve for φ and φ 1 : The expression for the Hamiltonian in terms of ξ, η ∈ S(R) has the form [12] In the special case γ 1 = γ, κ 1 = κ and µ = 0 the Hamiltonian acquires the form thus recovering the Hamiltonian determined in [11].When γ 1 = γ and κ 1 = κ = 0 then µ = (γ − γ 1 )ηη x and the Hamiltonian becomes which, given the definition of B from ( 21), recovers the result in [9] (noting that there is a sign difference is the second and fourth terms due to the different stream function convention used in [9]).
In the situation with κ 1 = κ which is physically realistic, because the unperturbed currents at the two layers have the same speed at the interface (i.e.absence of a vortex sheet), we have also µ = (γ − γ 1 )ηη x and The stability of the flow under disturbances of arbitrary wavenumber is assumed.For a detailed study of the stability (which is beyond the scope of our study) we refer to [1] and the references therein.
In what follows appropriate scales will be introduced and the corresponding expansions of the DN operators will be provided in terms of the smallorder scale parameters.

Scales and approximations
Let us introduce a scale ε = a h 1 (30) where the constant a represents the average amplitude of the waves η(x, t) under consideration and ε ≪ 1 is a small parameter which will be used to separate the order of the terms in the model.The Dirichlet-Neumann (DN) operators have the following structure where The corresponding expansions are [24,25] where D = −i∂/∂x.
We will study the equations of motion under the additional approximation that the wavelengths L are much bigger than h 1 , i.e.
Noting that the wave number k = 2π/L is an eigenvalue or a Fourier multiplier for the operator D (when acting on waves of the form e ikx ) we make the following further assumptions about the scales we can determine explicitly the scale factors where the barred quantities and operators together with h and h 1 are assumed to be of order 1.Introducing t h := tanh(hD), D 1 := h 1 D and omitting the bars for convenience we truncate the DN expansions as follows: Note that h appears only in the definition of the operator t h , which is of order 1.Since h 1 is assumed of order 1, then formally the order of the differentiation ∂ x is δ.Hence the order of the integration measure dx is 1/δ.We see now that the leading order terms in G are O(δ) and the leading order terms in G 1 are δ 2 hence G 1 G −1 ∼ δ ≪ 1 and hence we can expand as follows: Since both G and G 1 are self-adjoint, it is now evident that GB −1 G 1 is self-adjoint too.The substitution of ( 38) and ( 39) in (40) gives where the notation T h := −i coth(hD) = (it h ) −1 is introduced, more details are given in the Appendix Appendix A.1.We need also The quantity ηη x = ih 3 1 (η/h 1 )D 1 (η/h 1 ) ∼ δ 3 .The contribution of the integral density ηη x B −1 ηη x in the Hamiltonian is therefore of order δ 5 .Recall that dx ∼ 1/δ.Hence, by keeping terms up to δ 3 in (29) we have where the constant A = g(ρ − ρ 1 ) + κ(ργ − ρ 1 γ 1 ).We notice that H is of order δ.This gives the proper scaling of ∂ t which should be also of order δ, same as the order of ∂ x .The variation δu bears a scale factor δ as well.The equations (19) with the scaling written explicitly therefore are producing the coupled system These equations can be viewed as a generalisation of the irrotational case (Γ = γ 1 = γ = 0, κ = 0) derived in [24].

The intermediate long wave equation (ILWE)
Neglecting the terms of order δ 2 in ( 46), ( 47) we obtain We can perform a Galilean transformation of coordinates and taking into account that for the typical values of κ of several m/s, g ≫ 2ωκ the equations of motion can be written as A two-dimensional version of ( 53) -( 54) in the irrotational case (Γ = γ 1 = γ = 0) is obtained in [6].The leading order terms (i.e.neglecting the terms with δ above) produce a system of linear equations with constant coefficients from where the speed(s) of the travelling waves (in the leading order) could be determined: The plus sign is for the right-running waves and the minus sign is for the left-running waves.These speeds coincide with the speeds in the case of infinitely deep lower layer [13], indeed, the h-dependence comes only from the term T h which is of order δ.
In what follows, c could be either of the two solutions of the dispersion equation, the choice of solution determines if the results are relevant for the left or the right-running waves.
For the travelling wave, which depends on the characteristic variable X − cT, we also have u = ρ 1 h 1 cη and in order to obtain a single nonlinear equation for η we expect a relation which involves terms of order δ as well.In other words, we consider an expansion of the form for some yet undetermined constants α and β.This type of relation is known also as the Johnson transformation.The substitution of u from ( 56) in ( 51) and ( 52) when keeping only the terms up to order δ leads to two equations for η and therefore these two equations must coincide.This leads to equality of the coefficients in front of the terms of the same type, which further allows to determine the previously unknown and The equation for η is (59) The obtained equation is known as the Intermediate Long Wave Equation (ILWE) introduced in [37,43].It is an integrable equation.The soliton theory for ILWE has been developed in a number of works of which we mention [7,38,40,42,52].
The ILWE in the irrotational case (γ = where, from ( 55) Let us write the ILWE in the form where Details about the operator T h and the integrability of the ILWE are given in the Appendix Appendix A.1 and Appendix A.2 respectively.Here we only mention that the one-soliton solution of (61) has the form (cf. (A.10), Appendix Appendix A.2 ) 0 < k 0 < 2π h .
In the above formula X 0 and k 0 are the soliton parameters, i.e. arbitrary constants within their range of allowed values.X 0 is the initial position of the crest of the soliton and k 0 is related to its amplitude.The wavespeed of the soliton is c − δBk 0 cot(k 0 h) and the correction of order δ depends on the coefficient B and the dispersion law related to the dispersive term and also on the parameter k 0 which is related to the discrete eigenvalue of the spectral problem (A.2); in the 1-soliton case there is only one discrete eigenvalue.We observe that both the amplitude of the soliton and its speed are related through k 0 .Another feature of this solution is the fact that the function cot is unbounded.The physical relevance of the solution however requires that the choice of k 0 should be such that the quantity k 0 cot(k 0 h) is of order 1.As it will be shown in the next section, there is no such anomaly for the related Benjamin-Ono equation, and the limiting procedure requires a special choice of the parameter k 0 .

Connection to the Benjamin-Ono equation
The BO model of waves in the presence of uniformly-sheared currents has been derived previously in [13].The mathematical facts about the BO equation are given in Appendix Appendix A.2 for convenience.In this section we demonstrate that the BO equation could be obtained as a special kind of a long-wave limit from the ILWE.
In the limit h → ∞ which corresponds to an infinitely deep lower layer we have and the equation (59) becomes the well known Benjamin-Ono (BO) equation [2,46], see also [13], Like the ILWE, the BO is an integrable equations whose solutions can be obtained by the Inverse Scattering method [27,39,41].
The limit to the one soliton solution of the BO equation (from (63) with h → ∞ but k 0 h finite, k 0 h = π − k 0 /q where q is a constant) is described in the Appendix Appendix A.2 cf.(A.12) and can be written in the form where the constants are the initial position X 0 of the soliton and its amplitude η 0 .The relation to the constant q (and hence k 0 ) is

Connection to the KdV equation
The KdV model of waves in the presence of uniformly-sheared currents has been derived previously in [12].The special situation of KdV with h 1 /h ∼ δ has been provided in Appendix Appendix A.3 for convenience.The provided analysis shows that it could be obtained as a special kind of a long-wave limit from the ILWE.
The KdV limit could be obtained assuming hk ≪ 1 or |hD| ≪ 1.The operator Then from (61) we have Noting that h 1 /h ≃ δ ≪ 1 the correction to c in the second term is of order δ 2 and should be neglected.Thus we obtain the KdV equation in the form which coincides with (A.17), Appendix Appendix A.3 since in the ILWE setup δ ≃ ε and .
The one-soliton solution of (67) can also be obtained from (63).For k 0 h ≪ 1 we have sin(k 0 h) ≈ k 0 h, cos(k 0 h) ≈ 1, and also the identity The first term does not depend on k 0 and represents the small O(δ 2 ) correction to the constant wave speed c.The second term is proportional to k 2 0 , and the approximation leads to Introducing a new constant K = k 0 /2, we finally have which coincides with the KdV one-soliton solution (A.18), Appendix Appendix A.3.

ILWE BO KdV (A.13)
Table 1: The scales of the three approximation models.

Conclusions
We have derived the integrable ILWE for the situation of solitary waves on the interface of two fluids with constant vorticities, modeling equatorial internal waves interacting with uniformly sheared currents.The surface waves which are usually of much smaller amplitude are neglected, so a "rigid lid" approximation for the upper fluid is assumed.The ILWE is an integrable model and the inverse scattering method or other methods like Darboux transforms and Hirota's method allow for the derivation of explicit multisoliton solutions.The one-soliton solution for example is (63).The two important integrable limits leading to the BO and KdV equations (for specific choice of the physical quantities) are explained in details.In addition these two limits have been applied to the ILWE one-soliton solution, and the BO and KdV one-soliton solutions have been recovered.The limits of course exist for the multisoliton ILWE solutions as well and in principle for all types of solutions.It has to be noted that the ILWE, BO and KdV correspond to different scales of the physical quantities, as it could be seen from Table 1.The BO limit corresponds to a limit to an infinitely deep lower layer.The limit from ILWE to KdV gives a KdV equation in its particular form (67). Schematically the relations between the ILWE and the KdV equations (A.13), (A.17This shows that the ILWE should be used as a "master" equation with some caution.An interesting aspect for further studies is the development of theoretical models for internal waves with currents over variable bottom, extending the results from [49,50]. with , moreover A coincides with the expression from (62) and B 1 = hB/3.The one-soliton solution is (see for example [34] for the KdV solitons) where K and X 0 are constants.

Figure 1 :
Figure 1: System setup.The function η(x, t) describes the elevation of the internal wave.