Hamiltonian Formulation for Wave-Current Interactions in Stratified Rotational Flows

We show that the Hamiltonian framework permits an elegant formulation of the nonlinear governing equations for the coupling between internal and surface waves in stratified water flows with piecewise constant vorticity.


Introduction
Tropical ocean dynamics are quite intricate due to significant density stratification and to the interaction of waves with depth-dependent current fields [4,8]. The density stratification is particularly pronounced because of the presence of a sharp pycnocline that separates a shallow near-surface layer of relatively warm water from a deep layer of colder and denser water; since the decline in temperature with depth is responsible for the increase in density, the pycnocline is also a thermocline. The stirring of surface waters by the wind, typically in the ocean's uppermost 100-200 m, produces a well-mixed layer of practically uniform density; the assumption that the dark, cold deep layer below the pycnocline layer is of constant density is also physically reasonable-see [23]. While the pycnocline allows some kinetic energy to penetrate into the abyssal layer, it mainly acts as a barrier to vertical motion within the ocean. The internal waves that arise as fluctuations of this sharp interface are coupled with surface waves. Their primary sources are winds and tides [18] that alter either the free surface or the pycnocline so that, due to the coupling effect, waves propagate along the thermocline and at the ocean's surface. The strength of the coupling varies, and often very large internal waves are hardly noticeable at the surface. An observable manifestation of this feature is the dead-water phenomenon. This can occur when a ship moves in still water and its displacement sets up small surface waves, with typical amplitudes of the order of some tens of centimetres. The coupling effect triggers concealed internal waves that can reach amplitudes of up to tens of metres, until their presence disturbes the sub-surface flow to a degree that the ship becomes hard to manoeuver. On the other hand, the importance of the underlying currents is highlighted by a spectacular wave phenomenon, recently discovered in the South China Sea: field data provides evidence for tide-generated internal waves with amplitudes in excess of 100 m, extending over 100 km; see the discussion in [28]. The internal energy produced by the tides is essential in generating these waves. This is a striking illustration of the need to amend the common modelling assumptions (see [11,12,16,24] and references herein) that do not accommodate the possibility of underlying currents. While this issue is addressed in the recent research literature by means of numerical simulations, this type of approach presents the disadvantage that it makes it difficult to identify clearly the important processes at work. Therefore it is desirable to go beyond the limitations inherent to a case-by-case numerical computation. We will show that the Hamiltonian point of view can alleviate this unfortunate situation in the physically relevant case of flows with piecewise constant vorticity. The fact that, cf. [29], tidal currents are approximately uniform with depth and winddriven currents are sheared (thus being flows with non-zero vorticity), shows the oceanographic relevance of our theoretical considerations.
In this paper we will present the Hamiltonian formulation of the nonlinear governing equations, refraining from pursuing approximations, since the discussion of the wide range of physical regimes that are relevant to the derivation and interpretation of approximate models requires quantitative geophysical considerations specific to oceanography; we plan to pursue this direction in a future investigation. After presenting the governing equations in Section 2, we show in Section 3 that they can be expressed as an infinite-dimensional Hamiltonian system. In Section 4 we illustrate the benefit of the Hamiltonian formulation by a group-theoretical discussion of known conservation laws.
We would like to point out that, in addition to representing a technically more challenging problem, if compared with the classical irrotational setting of a pure wave motion with no underlying currents, wave-current interactions give rise to new dynamical phenomena associated with the possible appearance of critical layers. These are surfaces where the phase speed of wave propagation equals the mean flow speed-see [9,10,20,31] for a discussion in the simpler context of travelling waves in a homogeneous fluid. The extraction of momentum from the mean flow feeds the growth of critical layers, triggering instability mechanisms (see [22]). The fact that the typical Reynolds numbers in geophysical fluid dynamics are very large, cf. [22], precludes the standard fluid dynamics approach towards the resolution of the difficulty associated to shear flows with critical layers by means of a viscous structure across the layer. Instead, nonlinear effects dominate over viscosity. The possibility to perform in-depth nonlinear studies is therefore contingent upon finding a succinct formulation that highlights structural properties, and this is an inherent feature of Hamiltonian systems.

The Governing Equations
We consider a two-dimensional water flow with a free surface, moving under the influence of gravity and propagating in the positive x-direction, while the y axis points vertically upwards. The water flow occupies the domain bounded below by the rigid flat bed y = −h (with h > 0), and above by the free surface y = η 1 (x, t) + h 1 , the latter being a perturbation of the flat free surface y = h 1 . Here, h 1 > 0 is a constant, t stands for time and x → η 1 (x, t) is a periodic function of principal period L with mean zero, that is, L 0 η 1 (x, t) dx = 0 for all t. The density of the fluid is assumed to be distributed as follows: adjacent to the flat bed, there is a layer * := {(x, y, t) : x ∈ R, t ∈ R, −h < y < η(x, t)} of constant density ρ, separated by the interface y = η(x, t) from the free-surface adjacent layer * is, at any fixed time t, an L-periodic function with zero mean. Thus, the pycnocline and the free surface fluctuate about their mean levels y = 0 and y = h 1 , respectively; see Fig. 1. Throughout this paper some physical variables (for example, the density of the fluid and the vorticity of the flow, the perturbed velocity potentials or the components of the velocity field, but not the stream function or the pressure) may present discontinuities across the interface y = η(x, t) that separates the two layers. To draw attention to this, we use the index 1 as a label for the upper layer. Whenever we refer to the overall physical variable without specification of the layer, we shall use the hat or caret-shaped symbol placed on top of variable (for example, the densitŷ ρ takes on the value ρ 1 in the upper layer and the value ρ > ρ 1 in the lower layer).

Formulation in the Physical Variables
The assumption of inviscid flow is appropriate for ocean flows (see [22]). Therefore, denoting with (û(x, y, t),v(x, y, t)) the velocity field, the equations of motion are Euler's equations where P = P(x, y, t) is the pressure and g is the (constant) acceleration of gravity, supplemented by the equation of mass conservation Complementing the equations of motion are appropriate boundary conditions. At the free surface, the dynamic boundary condition (with P atm being the constant atmospheric pressure at the surface of the ocean) decouples the motion of the water from that of the air above. The kinematic boundary conditions, which refer to the flat bed y = −h, the pycnocline y = η(x, t) and the free surface y = h 1 + η 1 (x, t), ensure that no fluid particles cross such a boundary (and thus a particle on it will remain confined to it); they read as where (u 1 , v 1 ) and (u, v) are the velocity fields in * 1 and * , respectively, each admitting continuous extensions to the closure of its domain, Finally, another condition is required at the interface y = η(x, t), to ensure a balance of forces at this internal boundary. Since the only force exerted by an inviscid fluid on a boundary acts in the normal direction and is due to pressure, this balance is expressed by the condition that the pressure is continuous across the boundary: where the subscripts a, b in P a , P b refer to the limits of the pressure P from 'above' and 'below' the interface. A distinctive feature of the governing equations (2.1)-(2.7) is that the free surface and the interface are both unknown and must be determined as part of the solution, due to the strong interdependency of the waves and the flow in this nonlinear free-boundary problem. Let us also point out that, due to (2.5), the normal component of the fluid velocity is continuous at the interface y = η(x, t), so that the pycnocline acts as a barrier, effectively preventing cross-currents. In contrast to this, tangential velocity discontinuities may occur at the interface. To gain insight into this, let us consider the case of pure currents, in the absence of waves. A steady pure current in a two-dimensional flow over a flat bed is represented in the interior of the fluid domain by a vanishing vertical fluid velocity component and a depth-dependent horizontal fluid velocity component, the gradient of which is the vorticity: the velocity field is of the form (u 1 Note that, by (2.1), (2.3) and (2.7), in this case the pressure is hydrostatic throughout the fluid, with P = P atm + ρ 1 g(h 1 − y) for 0 y h 1 and P = P atm + ρ 1 gh 1 − ρgy for −h y < 0. In tropical ocean dynamics, ignoring geophysical Coriolis effects, the main sources for the perturbation of a static ocean are winds and tides. Nearsurface ocean currents are primarily caused by the wind. When the wind blows over the surface, there is transfer of momentum from air to water, part of which causes a net forward motion of the surface water in the direction that the wind is blowing. The surface water molecules moved by the force of the wind drag deeper layers of water molecules below them, with each deeper layer moving more slowly than the layer above it, until the movement ceases at a depth of about 100 m. While friction is essential in the generation of wind-driven currents, the fact that seasonal tropical wind pattern is characteristically regular and uniform (with occasional storms of great violence) permits us to model such currents within the inviscid setting as a permanent feature of the underlying flow, thus bypassing the mechanisms by which they come about. The simplest realistic model for a wind-drift current in the mixed layer * 1 with a flat free surface and interface, y = h 1 and y = 0, respectively, is a flow with u 1 = γ 1 y and v 1 = 0, for some constant positive vorticity γ 1 ; the abyssal layer * is still, as wind effects are confined to the mixed layer * 1 . As for the tidal currents, these are the alternating horizontal movements of water associated with the rise and fall of the tide. They are practically constant from sea surface to bottom, unless modified by stratification or topographic irregularities. Satellite altimeter data show that sub-pycnocline tidal currents decouple from those in the overlying near-surface layer, being weaker (Figs. 2, 3). Note that tidal currents flow alternately in approximately opposite directions, with a short period of little or no current, called slack water, at each reversal of the current. While during the uni-directional flow, the (uniform) speed varies from zero at the time of slack water to a maximum, about midway between the slacks, this variation in time is slow, so that for periods of several hours the unsteadiness of the current needs not be considered ( [27]). As they move over the ocean bottom, tidal currents sometimes encounters submarine seamounts, ridges, and other rugged features on the ocean floor. Deflected by such topographic obstacles, they can easily generate motions of the pycnocline, which propagate away from the source in waves. For example, the South China Sea is the site of intense internal wave motion. Internal waves are generated in the Luzon Strait, located between the islands of Luzon (Philippines) and Taiwan and connecting the South China Sea to the Pacific Ocean, where semidiurnal and diurnal tidal currents with magnitudes reaching 1 m/s flow nearly eastwest and are deflected by two north-south submarine ridges, generating internal waves that propagate westwards along the latitudinal zone of 21 • N, away from the Luzon Strait into the South China Sea. They subsequently evolve into large depicting an adverse tidal current (opposing the wind-drift). In the mixed layer above the pycnocline, the wind effect is dominant, with u 1 (h 1 ) > 0. At the flat interface y = 0, the jump discontinuity 0 > u(−h) = u(0) > u 1 (0) relates the limits u 1 (0) and u(0) from above and below the pycnocline and fast waves, attaining amplitudes in excess of 100 m and propagation speeds close to 3 m/s; see [28]. While from October to February the zonal wind pattern is oriented westward [32], the tidal currents present daily variations, alternating between the eastward ebb current and the westward flood current [14]. Since over a longitudinal zone in excess of 100 km, between 118.5 • E and 119.5 • E, the South China Sea bed is almost flat, with depth about 2.6 km (see [14]), the ocean flow in this region enters our general theoretical framework. The fact that the uniform tidal current in the abyssal layer * beneath the thermocline is weaker than the uniform tidal current in the mixed layer * 1 gives rise to a discontinuity of the fluid velocity component at the interface (of the tangential component, since the normal component vanishes throughout the flow), discontinuity that can not be mitigated by the wind-driven near-surface drift since this current dies out at the pycnocline. The possibility of accommodating discontinuities stems from the fact that the pycnocline can sustain stress, since the inviscid setting allows the upper and lower fluids to slide relative to each other without any viscous friction at the interface. Note that such discontinuities trigger Kelvin-Helmholtz instabilities (see [13]). Let us also point out that even in the absence of underlying currents, laboratory experiments in two-layer fluids consisting of fresh water above salt water show that the horizontal component of the velocity field induced by an internal wave appears to present jump-like discontinuities at the pycnocline; see [17]. Consequently, since a corresponding boundary-layer structure is not available in our inviscid settingjustifiable since in conventional geophysical-fluid-dynamical considerations, the Reynolds number is extremely large (see [22])-discontinuities of the tangential velocity field at the interface are intrinsic.
The previous considerations indicate the proper function space setting for our problem: the regular solutions of the governing equations (2.1)-(2.7) correspond to smooth functions η and η 1 , with corresponding smooth velocity fields in the two open domains * and * 1 , both admitting continuous extensions to the upper and lower boundaries, that is, The coupling between the two velocity fields is only reflected in the boundary condition (2.7), for a pressure that is thus uniquely determined, due to (2.1) and (2.3). While each velocity field proper to * and * 1 may admit a unique continuation beyond the boundary (for example, by analytic continuation in the case of travelling waves in flows of constant vorticity without stagnation points-see [6]), the continuous extension of the velocity field from the doublyconnected domain * ∪ * 1 to its closure might fail due to discontinuities across the interface.
For wave-current interactions the vorticity of the flow, defined bŷ characterizes the underlying currents (see [6]). Throughout this paper we assume that the vorticity is constant throughout each layer (but might be discontinuous across the interface), that is, there are constants γ, γ 1 ∈ R such thatγ = γ in * andγ = γ 1 in * 1 . Note that if the flow presents this feature initially, at time t = 0, the fact that the vorticity of a particle is preserved as the particle moves in a two-dimensional flow (see [6]) ensures that this feature persists at all later times.

Re-formulation in Terms of the Stream Function and the Perturbed Velocity Potential
In a situation where both waves and underlying currents are present in a flow, the "current" component at a point that is always lying within the fluid is defined as the average velocity, and the periodic components that fluctuate around this average are ascribed to the wave motion. While this definition excludes points lying at some instants out of the fluid, for example, above the trough level of a travelling surface wave, dividing the motion between waves and currents turns out to be especially useful in the context of flows with constant vorticity.

Let us denote by and 1 the time-dependent domains
representing, at a fixed time, the two fluid layers of constant density. The equation of mass conservation (2.2) ensures the existence of a stream function in each layer, ψ in and ψ 1 in 1 , determined, up to an additive term that depends only on time, by u = ψ y and v = −ψ x in , u 1 = ψ 1,y and v 1 = −ψ 1,x in 1 . (2.9) From (2.9) and (2.6) we see that the stream functions can be written as for some smooth functions ψ 0 , ψ + . Thus each stream function has a C 1 -extension to the closure of its respective domain, ψ ∈ C 1 ( ) and ψ 1 ∈ C 1 ( 1 ). From (2.5) and (2.9) we see that the stream functions ψ and ψ 1 differ only by a function of time at the interface y = η(x, t). The degree of freedom in (2.10) allows us to choose them equal on y = η(x, t), so that there exists ∈ C( ∪ 1 ) with = ψ on and = ψ 1 on 1 . The discussion in Section 2.1 shows that, in general, ∈ C 1 ( ∪ 1 ). Also, note that is periodic of period L in the x-variable: in this follows from the definition (2.10) of ψ, and due to (2.5) and the fact that the interface η has mean zero. With the vorticity defined by (2.8), equation (2.9) yields We now introduce in each layer a perturbed velocity potential, ϕ in and ϕ 1 in 1 , by requiring and v = ϕ y in , u 1 = ϕ 1,x + γ 1 y and v 1 = ϕ 1,y in 1 . (2.12) More precisely, let (2.14) Both perturbed velocity potentials admit extensions to the closures of their domains, ϕ ∈ C 1 ( ) and ϕ 1 ∈ C 1 ( 1 ). The discussion in Section 2.1 shows that we cannot expect the existence of a C 1 ( ∪ 1 )-function that coincides with ϕ in and with ϕ 1 in 1 . Actually, in contrast to the case of the stream functions, we cannot even expect the existence of a continuous function that extends both perturbed velocity potentials.
[In general, it is not possible to accommodate the equality of ϕ and ϕ 1 on the interface by taking advantage of the fact that on the right side of (2.13) and (2.14) we can add an arbitrary smooth time-dependent function. Indeed, using For γ = γ 1 , this would imply u = u 1 on the interface and consequently (2.5) forces v = v 1 on the interface, so that the velocity field would be continuous across the pycnocline, which is not granted (see Section 2.1).] Moreover, the periodicity issue for the perturbed velocity potentials is not straightforward. To see this, let us denote by U (y, t) the current underlying the interface, defined at the level y below the trough of the internal wave η by Then the periodicity of v and the definition of γ yield U y = γ , so that While it would seem that κ is time-dependent, note that the first equation in (2.1) in combination with (2.6) and the periodicity of u and P ensure κ (t) = 0. Returning to the issue of the periodicity of the perturbed velocity potential ϕ, from (2.13) and (2.15)-(2.16) is easy to see that ϕ(x + L , y, t) = ϕ(x, y, t)+κ L for all (x, y) ∈ . Therefore, the function (x, y) → ϕ(x, y, t) − κ x is periodic in the x-variable, with period L. To deal with ϕ 1 , we can proceed similarly in the upper layer 1 , where the underlying current at a level y, above the crest of the internal wave η and below the trough of the surface wave h 1 + η 1 , is defined by Using (2.1), (2.5) and the relation obtained by differentiating (2.5) with respect to the x-variable, the periodicity of P and of the velocity field (u 1 , v 1 ) yields where the subscript s stands for the evaluation of the function at (x, η(x, t), t). Therefore where κ 1 , given by (2.19), is time-independent. From (2.14), (2.18) and (2.20) we get In view of (2.12), we may write u =φ x + γ y + κ and v =φ y in , u 1 =φ 1,x + γ 1 y + κ 1 and v 1 =φ 1,y in 1 , Sinceφ andφ 1 are harmonic functions, (2.21) represents a splitting of the velocity field into an underlying steady current component and a periodic harmonic wave-velocity field. The kinematic boundary conditions (2.4) and (2.5) can now be written as and (2.23) respectively; here the subscript s 1 indicates the traces of the involved functions on the free surface y = η(x, t). We can also recast Euler's equations as for some functions f and f 1 . Changingφ 1 by a suitable additive time-dependent term-a procedure that does not alter its essential properties-and making use of (2.3), we obtaiñ Introducing the new variables and we can express (2.24) alternatively as while (2.7) can be formulated as the following constraint: (2.28) The advantage of introducing the stream functions ψ, ψ 1 and the perturbed velocity potentialsφ,φ 1 is twofold. On one hand, it permits us to eliminate the pressure from the system: once η, η 1 and these variables are known, to recover the pressure in the fluid we simply set This re-formulation deals automatically with the equations (2.1) and (2.2), and the boundary conditions (2.3) and (2.6), since these relations are innate to these auxiliary variables. Also, note that the kinematic boundary conditions (2.4) and (2.5) can be easily expressed in terms of the wave profiles η and η 1 as Consequently, the total number of unknown functions is reduced to six: η, η 1 , ψ, ψ 1 ,φ,φ 1 ; moreover, there is an interdependency between the pairs (φ, ψ) and (φ 1 , ψ 1 ) since (2.9) and (2.21) show that the complex-valued functions (x +iy) → φ +i ψ − γ 2 y 2 − κ y and (x +iy) → φ 1 +i ψ 1 − γ 1 2 y 2 − κ 1 y are analytic in and 1 , respectively. The second important aspect is that this process brings to light some physically relevant structural features; in particular, the importance of the underlying current field becomes apparent.

Main Result
We will show that the governing equations in two space variables, described in Section 2, may be reduced to a Hamiltonian system by collapsing the dynamical variables into suitable one-dimensional representations. Recall from [26] that a Hamiltonian system of partial differential equations has the form where t → ω(t) describes a path in a Hilbert space H equipped with an inner product (·, ·), the associated Hamiltonian functional H : D ⊂ H → R is defined on a dense subspace D of H, while the structure map J is a skew-adjoint (pseudo-) differential operator which defines a Poisson bracket for functionals by {F 1 , If the map J is invertible the Poisson structure is termed a symplectic structure, J is called the cosymplectic operator and its inverse K is the symplectic operator. In our context the variable ω is vector-valued, having four scalar components, and the Hilbert space is an appropriate product of spaces of square-integrable functions, with D the subspace of smooth periodic functions with an L-periodic dependence on the horizontal space variable x. The fact that, in the governing equations the partial differential equations in the interior of the domains and 1 , expressing the harmonicity of the respective perturbed velocity potentials, are supplemented by constraints on a free boundary and an interface is the hallmark of permissible variations of the fluid domain (see [15]). While in the preliminary considerations in Section 3.1 we will retain the space-time dependence of the functional involved in the nearly-Hamiltonian formulation, the subsequent analysis will show that one can change variables in such a way that the system becomes Hamiltonian, with the corresponding functional defined on a spatial domain (and not in space-time). This feature gives a compelling advantage with respect to other possible variational formulations (for example, of Lagrangian type) since the dynamical evolution will proceed smoothly in time, the Hamiltonian representation itself remaining intact, over a time-interval determined by the initial data, until the possible development of singularities (presumably with the onset of wave breaking).
There are several advantages associated with a Hamiltonian formulation. First, it represents an aesthetically pleasing procedure to reduce the number of variables, that can result in significant simplifications. Second, once a system is known to be Hamiltonian, an immediate benefit is the availability of a mechanism that provides conservation laws associated with symmetry groups of the system (like changes to moving coordinate frames, spatial rotations, time translations, scaling transformations), thus gaining insight into the dynamics. Third, the Hamiltonian perspective is very useful in perturbation analysis since it guarantees that if appropriate symmetries are preserved in the approximation process, then conservation properties of the exact system are maintained in its approximations. In our context this modelling aspect necessitates the detailed analysis of field data in the context of various physical regimes and for this reason it will be explored in a future publication. We aim to provide a Hamiltonian formulation that applies to waves of small and large amplitude, and for this reason we investigate the governing equations; the vast research literature on integrable and nearly-integrable models for water waves (see the survey [19]) is not relevant here, as the validity of these models is restricted to waves of small amplitude (and, mostly, within the setting of shallow-water theory). In the direction of Hamiltonian formulations for the governing equations for water waves, let us point out that some simpler scenarios were investigated earlier; they can be grouped in basically three categories: (i) the discovery by Zakharov [33] that the governing equations for two-dimensional irrotational deep-water waves have a Hamiltonian structure; (ii) the possibility of extending Zakharov's approach to the governing equations for two-dimensional gravity water-flows with constant vorticity [30]; (iii) the Hamiltonian formulation of the governing equations for interfacial and surface waves in the absence of underlying currents, which describe the evo-lution of an interface between two internal fluid layers, and its coupling with the motion of an overlying free surface [11].
Our setting blends the challenges presented by the amplifications (ii) and (iii) to the seminal result (i), due to the combined effect of stratification and underlying linearly sheared currents. While physical considerations suggest the total energy as an apt Hamiltonian functional, the more subtle mathematical aspect is the choice of suitable dynamical variables. In this context, we would like to point out the key contribution [2], which provides the correct dynamical variables for the Hamiltonian formulation of interfacial irrotational waves with a rigid lid.

The Nearly-Hamiltonian Formulation
Starting from the natural candidate for the Hamiltonian functional, given by the total energy of the flow, we now present a nearly-Hamiltonian formulation of the governing equations: a system that is Hamiltonian in the absence of sheared underlying currents (that is, for γ = γ 1 = 0). In Section 3.2 we will show that a suitable change of variables can be used to obtain a genuine Hamiltonian formulation.
In our considerations the following notation turns out to be convenient: (3.1) At a fixed time, consider the total energy of the flow, given by where the first term represents the kinetic energy (energy of motion) andρgy is the gravitational potential energy (energy of position). Taking into account the stratification of the fluid, we can rewrite H as which, by (2.21), equals Before computing variations of the functional H , let us collect a few useful formulas. For harmonic functions θ 1 and θ 2 the identity ∇ ·(θ 1 ∇θ 2 ) = (∇θ 1 )·(∇θ 2 ) holds, yielding, for harmonic variations δφ ofφ, that We will also employ the following rule for computing variations (see [15]): For variations that do not alter the fixed bed y = −h, and at a fixed wavelength L, by (3.4) we have Since harmonic variations of the perturbed velocity potentialφ are of interest, to explicitate the first term above we employ (3.3) and the divergence theorem in the relying in the last step on (2.6) and on the periodicity ofφ. Since so that, due to (3.6), (3.8) In the same manner we obtain We proceed now to compute the variations of the other terms in (3.2). Making use of (3.5), (3.7), (2.6) and of the periodicity of the functions η, η 1 ,φ andφ 1 , we get and similarly, Moreover, (3.12) and Finally, for m = 1, 2, 3, we have the following standard variational formulas: Note that h 1 does not vary, as y = −h is the rigid bed and h + h 1 represents the mean depth of the water-we are considering flows without sources or sinks. After these preliminary considerations, let us introduce the variables Throughout this subsection we take for granted the non-trivial fact (proven in Section 3.3) that the total energy functional H , defined in (3.2), has the form for some Hamiltonian density function h that depends solely on the four scalar variables ξ, ξ 1 , η, η 1 , and their spatial derivatives. This permits us to present the following partial result. (3.18) where (χ , χ 1 ) are defined in (2.25)-(2.26) and (ξ, ξ 1 ) in (3.16).

Theorem 1. The governing equations admit the nearly-Hamiltonian formulation
In the absence of stratification Theorem 1 particularizes to the nearly-Hamiltonian formulation for water waves with constant vorticity obtained in [7]. Indeed, setting ρ = ρ 1 and γ = γ 1 , the fact that = 1 allows us to ignore η and ξ in (3.18), thus obtaining the reduced system For irrotational flows (γ 1 = 0), the system (3.19) is Hamiltonian with phase space of smooth and L-periodic scalar functions. Therefore in this setting the water-wave problem has a symplectic structure; we recover the finite-depth analogue of Zakharov's Hamiltonian formulation for irrotational deep-water waves [33].
On the other hand, for stratified irrotational flows (γ = γ 1 = 0) the system (3.18) is again Hamiltonian; we recover the Hamiltonian formulation from [11], with phase space Using (2.22), we identify the factor of δ 2 in the above integrand with ρ 1 η 1,t , while from (2.23) we see that the factor of δ is ρη t and that of δ 1 is (−ρ 1 η t ). The simplification of the factors of δη and δη 1 is more involved. To deal with the factor of δη, let us first rely on (2.23) to write it as (3.20) We now take advantage of the following identities obtained at once by combining (2.9) and (2.21), to further transform the expression After algebraic cancellations, the above expression simplifies to Using (2.28), this becomes due to (3.1) and (3.16). This is a satisfactory expression, so now we turn to the factor of δη 1 , proceeding analogously. First, we use (2.22) to write it in the form Since the combination of (2.9) and (2.21) yields the identity we obtain the equivalent expression for the factor of δη 1 . Algebraic cancellations simplify this to and, using (2.27), we can write the above as relying in the last steps on (3.1) and (3.16).
We can summarize the previous computations as follows: (3.21) in view of (3.16). The system (3.18) now emerges from the definition of the variational derivative with respect to the inner product in the space L 2 [0, L] of squareintegrable functions.

The Hamiltonian Formulation
We now show that if we perform the change of variables then the nearly-Hamiltonian system (3.18) becomes Hamiltonian even for nonvanishing (piecewise constant) vorticity.

Theorem 2. The Hamiltonian system
is a re-formulation of the governing equations.

The Hamiltonian Functional in Terms of the Dirichlet-Neumann Operator
The purpose of this section is to write the Hamiltonian functional defined by (3.2) in a more concise form by means of the Dirichlet-Neumann operator. This process will also establish the validity of the formula (3.17) that was relied upon in Theorem 1.
Given smooth, L-periodic scalar functions and η such that η(x) > −h for all x ∈ [0, L], letφ be the (unique) solution, with an L-periodic dependence on the x-variable, of the boundary value problem where * (η) = {(x, y) : x ∈ R, −h < y < η(x)}. Denoting by n the outward unit normal along the upper boundary y = η(x) of the domain (see Fig. 4), the Dirichlet-Neumann operator G = G(η) associated to * = * (η) is defined by Since we consider spatially periodic motions, the lateral boundary conditions on a periodicity cell of * express the condition that both the data and the solution exhibit a periodic dependence on the horizontal spatial variable x. Similarly, given smooth, L-periodic functions η, η 1 , 1 , 2 , such that η(x) < h 1 + η 1 (x) for all x ∈ [0, L], we denote byφ 1 the (unique) solution, with an L-periodic dependence on the x-variable, of the Dirichlet boundary value problem The Dirichlet-Neumann operator G 1 = G 1 (η, η 1 ) associated to * 1 = * 1 (η, η 1 ) is defined by the normal derivatives of the solution on the lower and upper boundaries of the domain * 1 , that is, Here n 1 is the outward unit normal vector along the upper boundary y = h 1 +η 1 (x) of * 1 , while, for consistency with the notation used in defining G, n is the inward unit normal at the lower boundary of * 1 ; see Fig. 5. We denote the entries of the matrix operator G 1 (η, η 1 ) by There are two approaches for investigating the structural properties of Dirichlet-Neumann operators: a perturbative approach-possible due to the analytic dependence of the operator on the curves representing the upper and lower boundaries of the domain [5] and an alternative based on a boundary integral representation for the solution using single and double layer potentials (see the discussion in [1]). For background information on Dirichlet-Neumann operators we refer to [10,11,21].
Let us now show that the Hamiltonian functional H , given by the total energy (3.2), depends only on the variables ξ, ξ 1 , η, η 1 and their spatial derivatives.
From the definition of the Dirichlet-Neumann operators, using (2.23), we see that and while (2.22) yields (3.30) Adding up the relations (3.28) and (3.29), we obtain (3.31) Let us now define the operator Recalling (3.16), (3.31) enables us to express , 1 , 2 in terms of ξ and ξ 1 , as follows With the help of the Dirichlet-Neumann operators, using Green's second identity and (3.28)-(3.30), we get taking advantage of (3.31) in the last chain of equalities. Similarly, using Green's theorem for the vector field (−yφ, 0) in the simple regions and 1 , we obtain that while Green's theorem for the vector field (−φ, 0) yields Summing up the above three relations, we get (3.36) Introducing the notation and taking into account (3.33)-(3.35), in view of the equalities ργ − ρ 1 γ 1 1 = ρ(γ − γ 1 ) + γ 1 ξ , ρκ − ρ 1 κ 1 1 = ρ(κ − κ 1 ) + κ 1 ξ, we obtain that where, in the second equality, we have used (3.32) to infer that and the last equality is obtained by exploiting the relation ρG 11 = B − ρ 1 G, due to (3.32), and taking advantage of the fact that both operators B −1 and G 11 are self-adjoint, while G * 12 = G 21 , cf. [11]. Finally, due to (3.2) and the fact that η as well as η 1 have mean zero, we can sum up the previous computations by the formula This representation proves our claim about the functional dependence of the total energy H .

Symmetries and Conservation Laws
Apart from the intrinsic aesthetic appeal of the Hamiltonian formulation, there are many potential benefits. A prime example of practical importance is that it permits the application of Noether's theorem (actually, of the generalization of Noether's theorem for infinite-dimensional evolution equations, discussed in [26]), which relates conservation laws to continuous symmetries of the system of governing equations.
The two-dimensional free boundary problem for incompressible irrotational gravity waves propagating at the surface of a homogeneous inviscid fluid (one layer) is known to have eight nontrivial conservation laws [3,25]. The conservation laws are a consequence of the symmetry group of the problem. In the set-up of the present study the depth-dependent currents introduce an additional y-dependence of the velocity field, which breaks some of the symmetries which are valid in the setting of the irrotational flow of an inviscid homogeneous fluid. More specifically, from the symmetries listed in [26] only the x− and t-translations and adding constant to the potentials survive (the horizontal Galilean boost is not a symmetry anymore, because of the γ 1 ψ 1 , γ ψ terms in the Bernoulli-type equations, for example, (2.24) with stream functions ψ 1 and ψ depending on the shear flow). According to Noether's theorem (see [26]), the t-translation invariance leads to the conservation of the energy and the x-translation invariance leads to the conservation of the horizontal momentum, L 0 J 0 dx, with density where L(η, η 1 , η t , η 1,t ) is the Lagrangian density. Here δL δη t = z and δL δη 1,t = z 1 are the canonical momenta, conjugate to the canonical coordinates η and η 1 (see Section 3.2). Thus L 0 (zη x + z 1 η 1,x ) dx is time-independent. In addition to these conservation laws (energy and horizontal momentum), we have the mass conservation laws L 0 η dx and L 0 η 1 dx, as explained in [25]. For convenience we did set both equal to zero in the process of defining the mean depths h 1 and h of the two layers, as discussed in the beginning of Section 2.