A Hamiltonian Approach to Wave-Current Interactions in Two-Layer Fluids

We provide a Hamiltonian formulation for the governing equations describing the two-dimensional nonlinear interaction between coupled surface waves, internal waves, and an underlying current with piecewise constant vorticity, in a two-layered fluid overlying a flat bed. This Hamiltonian structure is a starting point for the derivation of simpler models, which can be obtained systematically by expanding the Hamiltonian in dimensionless parameters. These enable an in-depth study of the coupling between the surface and internal waves, and how both these wave systems interact with the background current.


I. INTRODUCTION
Ocean flows in a band of about 2 • latitude either side of the Equator have typically a 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 1,2 .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 waves interact with the underlying currents.The vanishing of the Coriolis parameter at the Equator distinguishes the dynamics of the equatorial zone from the ocean dynamics at higher latitudes.The strong stratification confines the wind-driven currents to a shallow near-surface region, less than 200 m deep.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.The wind-generated equatorial current is modelled using the concept of vertical eddy viscosity [3][4][5] , its main features being that in the layer above the thermocline the depth-dependence is strictly monotonic and exhibits flow-reversal, while beneath the thermocline the current decays with increasing depth, being irrelevant in the abyssal region.Thus the choice of a piecewise linear current captures the primary structure of the equatorial A typical wind-induced current (curve in bold) in a two-layered fluid with a flat bed, and its inviscid approximation -a current with piecewise constant vorticity (dashed broken line).The wind drives a return flow, initiated at some sub-surface depth and extending beneath the pycnocline.
current system (see Fig. 1).While viscous theory is essential in explaining the generation of the equatorial current induced by wind forcing, inviscid theory is adequate for the study of non-turbulent wind-current interactions since the relevant Reynolds numbers are very large 6 .A further aspect of relevance is that often internal waves are highly nonlinear and occur in wave packets (of up to 30 waves) that disintegrate into trains of solitary waves, with the distances between the waves and also their amplitude decreasing from front to rear 7 -a behaviour reminiscent of soliton theory 8 .Note that the wave heights of internal waves are generally much larger than those of surface waves: while the average wave height in the equatorial Pacific is typically less than 4 m (cf. the NOAA database), internal waves more than 20-30 m high are common in this region 1 .The setting described above applies not only to equatorial flows, being typical for the evolution of large amplitude oscillations of an interface between two internal fluid layers, and its coupling with the motion of an overlying free surface, in the presence of wind-generated currents.
Concerning the assumption of flows with constant vorticity in each layer, note that since internal waves are typically long compared with the water depth, in the description of the underlying non-uniform currents it is the existence of a non-zero mean vorticity that is important rather than its specific distribution 9 .To gain insight into the nonlinear dynamics in an inviscid medium, the most advantageous approach is the Hamiltonian formalism 10 .We confine our considerations to two-dimensional flows since this setting, while simpler than the fully threedimensional case, is nevertheless suitable to describe the long-crested waves that are ubiquitous in ocean flows.We will show that for underlying currents with piecewise constant vorticity, the governing equations for twodimensional wave-current interactions, in a two-layered fluid overlying a flat bed, have a Hamiltonian formulation.Furthermore, to be able to make some quantitative predictions, we will develop structure-preserving approximations which enable us to derive model equations that are amenable to an in-depth study.It turns out that, other than preserving the original Hamiltonian structure, the weakly nonlinear shallow-water regime brings about an additional structure: by ignoring the terms that have in this setting a negligible effect, the obtained model equations have a bi-Hamiltonian, integrable structure.This opens up new possibilities for the investigation of the long-term dynamics for weakly nonlinear shallowwater waves interacting with wind-generated currents in a two-layer fluid.

II. HAMILTONIAN STRUCTURE OF THE GOVERNING EQUATIONS
Using the overbar notation for physical variables, let y = −h be the flat bed, let y = η(x, t) specify the location of the pycnocline and let y = h 1 + η 1 (x, t) be the free surface, where h 1 is the mean depth of the pycnocline and h is the mean depth of the lower layer.The interfaces delimit two fluid layers of constant density ρ > ρ 1 , in stable stratification.Let (u(x, y, t), v(x, y, t)) and (u 1 (x, y, t), v 1 (x, y, t)) be the fluid velocity fields in these two layers; see Fig. 2.
The equations of motion are the Euler equation and the equation of mass conservation, with an underlying current of piecewise constant vorticity, that is, ( For flows with piecewise constant vorticity in two-layered fluids, the separation of the flow into a current component (defined as the average velocity) and a wave perturbation thereof holds within the framework of the fully nonlinear theory, without recourse to approximations.This feature permits us to introduce time-dependent velocity poten- tials ϕ and ϕ 1 for the wave-related components: where ψ, ψ 1 are the time-dependent stream functions and κ is the (time-independent) current at the mean level of the pyconcline, y = 0. Due to mass conservation in the corresponding layer, ϕ and ϕ 1 are harmonic functions, while ∆ψ = γ and ∆ψ 1 = γ 1 .The relevant boundary conditions are the dynamic boundary condition, requiring that the pressure equals to the (constant) atmospheric pressure at the surface, together with kinematic boundary conditions which refer to the flat bed, to the interface y = η(x, t), and to the free surface (reflecting the impermeability of these three surfaces), and the continuity of the pressure across y = η(x, t) (expressing, within the inviscid setting, the balance of forces at the interface).
A reduction in dimension for the governing equations can be achieved by means of Dirichlet-Neumann operators, which deliver normal derivatives of harmonic field quantities at the boundary ('Neumann data'), given boundary measurements ('Dirichlet data').Denoting by Φ and Φ 1 the traces of ϕ and ϕ 1 on the pycnocline y = η(x, t), and by Φ 2 the trace of ϕ 1 on the surface y = h 1 + η 1 (x, t), the Dirichlet-Neumann operator G(η) associated to the lower layer is defined as where n = (−η x , 1) is the outward normal along the interface.The Dirichlet-Neumann operator G 1 (η, η 1 ) associated to the upper layer is given by where n 1 is the outward normal along the free surface.
Note that G 1 (η, η 1 ) is a matrix-operator: On the pycnocline y = η(x, t), the definition of the Dirichlet-Neumann operators and the kinematic boundary conditions yield Adding up the previous two relations, we obtain Using (7), setting and introducing the new scalar variables we can write Φ, Φ 1 , Φ 2 , in terms of ξ and ξ 1 as follows: The Hamiltonian H, which is the total energy, modulo irrelevant constants of integration, is given by where M (η, η 1 ) is the matrix The governing equations are Hamiltonian with respect to the canonical momenta p, p 1 , and coordinates q, q 1 : where The particular case when the underlying current is absent (that is, in the irrotational case γ = γ 1 = 0 with no uniform underlying current, i.e. with κ = 0), was already investigated 11 , and it represents the generalisation to a two-layered fluid of Zakharov's celebrated result for deep-water gravity waves 12 .Another particular case is obtained in the absence of stratification but with an underlying current of constant vorticity 13 .In the setting of travelling waves propagating in irrotational flow with no stratification, the Hamiltonian structure was essential in uncovering the flow pattern beneath the surface waves: theoretical studies 14,15 were confirmed numerically 16 and experimentally 17 .In this special context, the presence of an underlying current -that can only be uniform in these circumstances -increases considerably the dynamical complexity.For a two-layer fluid with a non-uniform underlying current, numerical and experimental data suggest an even richer dynamics, highlighted by the possible appearance of critical levels.The Hamiltonian structure is expected to provide useful insight into these aspects, building on available numerical and experimental results for waves in a two-layer fluid with no underlying current 18 .

III. STRUCTURE-PRESERVING APPROXIMATIONS
The process of non-dimensionalization permits us now to use the texture revealed by the Hamiltonian formulation of the governing equations to develop a structure-preserving (actually, structure-enhancing) perturbation approach in the weakly nonlinear longwavelength regime.We will derive nonlinear integrable model equations, of KdV-type.
We use as the fundamental length scale h 1 , the mean depth of the pycnocline; typically h h 1 for ocean flows (e.g.h 1 ≈ 100 m and h ≈ 4 km for the equatorial Pacific) but the intermediate long-wave regime, in which the upper and lower layer have comparable mean depths, is also a realistic possibility, even if we do not pursue it here.Associated with this length scale is the speed scale gh 1 , and let a be the average internal wave amplitude.
These scales define the set of non-dimensional variables: with the absence of the over-bar indicative of the nondimensional counterparts of the physical variables.The amplitude parameter will provide us with the linearised equations, by neglecting contributions of order O(ε 2 ).For example, from v = ηt + ūη x it follows that v = ε(η t + uη x ) and, since v is the y-derivative of the velocity potential ϕ, we have ϕ = εh 1 g h 1 ϕ, and thus ξ = ερ h 1 g h 1 ξ.The scales for u, u 1 , v and v 1 are of leading order O(1) related to the current components, with 'wave' components (arising as x-and y-derivatives of ϕ and ϕ 1 ) of order O(ε).On the other hand, the Dirichlet-Neumann operators (i.e.any of the introduced operators G, G ij ) admit expansions of the following type 11 : where For example, the leading order, O(1), expression for the operator G 1 is where D = −i∂ x (see 11 ).With this scaling, and ignoring the linear terms, whose mean vanishes, the Hamiltonian can be expanded as the equations are The quadratic part of H (the leading order of ε) produces the linearised equations: a homogeneous system of four equations admitting, for the wavelength L, solutions ũ, ũ1 , η, η 1 proportional to e i k(x−c t) , where k = 2π/ L, with a compatibility condition producing a fourth order polynomial equation for c( k) (the dispersion relation).To gain insight into weakly nonlindear wave-current interactions we study the small-amplitude long-wave regime, defined in terms of the shallowness parameter being such that δ 2 = O(ε).This setting is appropriate for the propagation of long internal waves, in which case the coupled surface waves have very small amplitude 3 .
To capture this feature, we assume that in our approximation the pair of variables ũ1 and η 1 , associated to the free surface, are of order O(ε) if compared with the variables ũ and η (at the pycnocline).In the long wave limit L → ∞, at leading order, the propagation speed c of the linear solutions, in the regime δ 2 = O(ε), satisfies a quadratic equation with solutions where In the irrotational case the above formula recovers the classical dispersion relation for internal long waves propagating on a uniform current: However, in the presence of underlying sheared currents, there is a 'threshold phenomenon': the travelling waves do not exist always.For the equatorial Pacific the values ρ = 1000 kg/m 3 , ρ 1 = 1002 kg/m 3 , h = 4 km, h 1 = 100 m, γ = 2.5 × 10 −4 s −1 , γ 1 = −0.015s −1 , κ = 1 m/s, L = 50 km are adequate (see 3 ).The predicted (realistic) speed of the eastward propagating waves is c 1 = 2.3 m/s, while for the waves propagating to the West the speed value c 2 = −63 m/s emerges.This latter speed, of the order of hundreds of km/h, is only encountered in special circumstances, e.g. in connection with the generation of tsunami waves (for example, the 22 May 1960 Chilean tsunami 8 ).The full dispersion equation for linear waves being a polynomial equation of degree four, other than producing another large speed c 3 , it provides us with an additional realistic speed value c 4 = −0.46m/s.Since the surface wind-drift current flows westward at a speed of 0.5 m/s and flow-reversal occurs beneath the surface, the speed c 4 is a critical speed, triggering the appearance of critical levels (where the wave speed equals the meanflow speed); there is evidence 19 of thermoclinic eddies associated with internal waves in the western equatorial Pacific Ocean.The natural mechanism for describing the structure of the flow near critical levels, within an inviscid system, is to invoke nonlinear theory.
The approximate Hamiltonian of the weakly nonlinear system retains terms of orders ε 2 and ε 3 .From the linearised equations it follows that, at leading order, with c 0 = c − κ, for some yet unknown function r(x, t), while The most general form for r is for some constants b i , i = 1, 2. Now we are in a position to express the evolution equations only by means of the variable η.This way we obtain two evolutionary equations for η which should coincide up to O(ε).The equality of their coefficients allows to recover the constants: Then η satisfies the KdV-type equation which, in the irrotational limit when all vorticities are zero, becomes (with c 0 = c − κ): From the solution η of (15), one can recover ũ and as well as η 1 and ũ1 : In the setting of the equatorial Pacific, formula (18)  shows that for the slow propagation, whether eastward or westward, the coupled internal and surface waves are in phase since in both cases c 0 /c > 0. However, for other choices of the mean depth of the two layers, coupled outof-phase wave propagation is possible (see 20 ).One appealing feature of the KdV model for shallow water waves is that it can be solved exactly using inverse scattering theory: starting with arbitrary initial data that are smooth and sufficiently localised in space, the solution evolves into an ordered set of localised solitary waves (solitons), with the tallest in front, followed by an oscillatory tail.Moreover, the details of this general picture can be predicted fairly easily from detailed knowledge of the initial data 8 .This feature opens up a wide range of new possibilities for the qualitative investigation of wave-current interactions in stratified fluids.

a)
FIG.1.A typical wind-induced current (curve in bold) in a two-layered fluid with a flat bed, and its inviscid approximation -a current with piecewise constant vorticity (dashed broken line).The wind drives a return flow, initiated at some sub-surface depth and extending beneath the pycnocline.

FIG. 2 .
FIG.2.Depiction of a cross-section of the fluid domain, showing the coupled waves at the surface and on the pycnocline.