Equatorial Wave–Current Interactions

We study the nonlinear equations of motion for equatorial wave–current interactions in the physically realistic setting of azimuthal two-dimensional inviscid flows with piecewise constant vorticity in a two-layer fluid with a flat bed and a free surface. We derive a Hamiltonian formulation for the nonlinear governing equations that is adequate for structure-preserving perturbations, at the linear and at the nonlinear level. Linear theory reveals some important features of the dynamics, highlighting differences between the short- and long-wave regimes. The fact that ocean energy is concentrated in the long-wave propagation modes motivates the pursuit of in-depth nonlinear analysis in the long-wave regime. In particular, specific weakly nonlinear long-wave regimes capture the wave-breaking phenomenon while others are structure-enhancing since therein the dynamics is described by an integrable Hamiltonian system whose solitary-wave solutions are solitons.


Introduction
The Earth's rotation profoundly affects the dynamics of the atmosphere and of the ocean. However, an in-depth qualitative study of the governing equations for rotating fluids is analytically untractable and prohibitively expensive computationally. For this reason, it is imperative to derive simpler, manageable models which can be studied in detail. A significant feature of equatorial ocean dynamics is that the change of sign of the Coriolis force across the Equator produces an effective waveguide, with the Equator acting as a (fictitious) natural boundary that facilitates azimuthal flow propagation which is symmetric with respect to the Equator: wave-current interactions which propagate in the longitudinal direction and are symmetric with respect to the Equator, featuring so-called Kelvin-waves (large-scale wave motions affected by the Earth's rotation and trapped at the Equator) and a depth-dependent underlying current field. Moreover, due to the pronounced density stratification of the equatorial regions, greater than anywhere else in the ocean, a rather sharp interface-the thermocline-separates a shallow nearsurface layer of relatively warm water from a deep layer of colder and denser water, each layer being practically of constant density. Consequently, there are two types of equatorial Kelvin waves: surface and internal waves. These two types of waves have quite different spatial and temporal scales, the internal waves being significantly larger and slower, with much longer wavelengths. Let us also point out that, in the ocean, energy tends to be concentrated in the lower frequencies (see [16]). Internal equatorial Kelvin waves are slow (taking more than two months to cross the equatorial Pacific) and have very long wavelengths (measured in tens/hundreds of kilometres), being typically more strongly excited than other ocean waves. They play an essential role in the "El Niño"phenomenon (see [22]): a spectacular, naturally occurring anomalous appearance, every few years, of unusually warm water in the central Pacific, transported eastwards by the internal Kelvin waves until a large expanse of the equatorial Pacific becomes much warmer than the average, thus altering the usual heat exchange process with the air above it and resulting in dramatic shifts of weather patterns in a chain reaction around the world. Moreover, the fact that a Kelvin wave is uni-directional facilitates the growth and accumulation of nonlinear perturbations as the wave propagates, so that nonlinear effects can become very large even if the initial amplitude is relatively small-it is not uncommon to record internal wave heights in excess of 40 m. Nonlinearity causes the often observed distortion of some internal equatorial Kelvin waves: the leading edge steepens while the trailing edge becomes flatter. These physical motivations prompted an intensive research activity, with the available studies of equatorial ocean waves falling into one of the following categories: (I) observational data and/or numerical simulations (see the discussion in [20]); (II) linear wave perturbations of a depth-dependent underlying current (see the discussion in [7]); (III) nonlinear studies assuming a passive equatorial current field (see the discussion in [29]).
Category (I) studies convey valuable information but, given the complexity of the encountered flows, by itself it does not suffice to identify clearly the important processes that are at work. Concerning the studies in category (II), the common occurence of nonlinear equatorial phenomena provides the impetus to go beyond the linear theory. As for (III), the importance of the interactions between waves and currents is highlighted by the fact that the underlying current field in the equatorial Pacific, generated by the prevailing westerly ambient wind pattern and the forces created by the Earth's rotation (see the discussion in [3,8]), presents flow-reversal: an eastward jet whose core resides in the thermocline (the Equatorial Undercurrent) is sandwiched between a westward surface flow and an abyssal layer of practically still water. The fact that a flow with piecewise constant vorticity (negative above the thermocline and positive below it) captures this salient feature within a Hamiltonian framework (see Sect. 3) opens up the possibility to develop a rigorous in-depth study of nonlinear wave-current interactions. A detailed description of the salient physical features of the ocean flow in the equatorial Pacific and of the precise aims of our mathematical analysis is provided in Sect. 2. The analysis of the Hamiltonian structure of the governing equations is pursued in Sect. 3. In Sect. 4 we develop a Hamiltonian perturbation theory. The study of the linearised equations, performed in Sect. 4.1, reveals acute differences between the two most important physical regimes for wave propagation: • Internal linear short waves propagate slowly eastwards, while the short waves at the surface are fast and can propagate eastwards as well as westwards, with each propagation mode having no noticeable effect at the other interface. • The effects of linear long-wave propagation are practically confined to the motion of the thermocline, whose oscillations propagate slowly eastwards or westwards. The fact that the westward internal waves are slower is an outcome of the dynamic response of the ocean to the presence of strong underlying currents with flow reversal.
However, there are two major limitations of linear theory. Firstly, the ubiquitous wavebreaking phenomenon is outside the realm of linear waves. Moreover, in the context of equatorial wave-current interactions, the common occurrence of internal solitary waves-localised waves that maintain their coherence, propagating with constant speed and unaltered shape-is not captured by linear theory. In light of this, and given that most ocean energy is concentrated in the long waves, in Sect. 4.2 we lay out the foundations for carrying out a systematic weakly nonlinear analysis in the principal long-wave scaling regimes of geophysical relevance (see Table 1). Of particular interest is the derivation of model equations that capture the most striking geophysical manifestations of nonlinear phenomena: • Internal solitary waves of permanent form which owe their existence to a balance between wave-steepening effects and wave dispersion-these wave patterns are important because they are often large, very energetic events, and they have therefore a significant role in mass and momentum transport across the ocean. Actually, the corresponding model equation is not merely Hamiltonian, it is completely integrable and therefore the solitary waves present the enhanced structure of a soliton-solitary waves that have an elastic scattering property (after colliding with each other, they eventually emerge unscathed, retaining their shape and speed) and present remarkable stability properties. This opens up the possibility of a detailed analysis of the nonlinear interaction of several such wave patterns by means of an inverse scattering approach, based on an appropriate Riemann-Hilbert problem formulation. Table 1. Overview of the considered long-wave scaling regimes, expressed in terms of the waves (a similar scaling is performed on the tangential velocities at the two interfaces), and of the main qualitative features that can be accommodated within a specific regime Long-wave regime Spatio-temporal scaling of the waves at the two interfaces (internal/surface) Main qualitative feature of the wavecurrent interaction Linear η(x, t) = ε η (x, t), • Large-amplitude internal waves that break-other than their widespread occurrence, they are perceived as important for the turbulent mixing their death throes produce. This type of solutions correspond to a model similar to the classical nonlinear shallow water equations, modified to account for the presence of underlying currents in a stratified flow. We use the method of characteristics to provide insight into the fascinating process of wave breaking.
In Sect. 5 we overview the obtained results and we offer a perspective on promising future directions for related research. We conclude this introduction by pointing out that, while we made a conscious effort to avoid being pedantic, there are some aspects of our investigation where a cavalier disregard for proper mathematical rigour is counter-productive. These are best illustrated by some conceptual and technical challenges that contrast to the situation encountered in the periodic setting, being specific for localised wave-current interactions. A central issue is the fact that the Hamiltonian cannot be the total energy of the flow since the kinetic energy contribution of the current is already infinite. Instead, we have to filter out the contribution from the current field and prove that this is achievable at the nonlinear level. Also, defining canonical variables of Schwartz class in the presence of depthdependent underlying currents is not a foregone result; to the best of our knowledge, this problem has not been addressed in the research literature.

Preliminaries
The aim of this section is to briefly discuss the observational record as a physical background for the systematic mathematical approach developed in this paper.

Key features of equatorial ocean dynamics.
Within a band of about 100-150 km of the Equator and extending longitudinally over about 16,000 km, the Pacific Ocean possesses some remarkable features: a significant density stratification (that is greater than anywhere else in the ocean, see [14]), an underlying current field with flow-reversal, and a wide variety of observed wave propagation phenomena. The hallmark of the pronounced density stratification is the presence of a rather sharp thermocline that separates a shallow near-surface layer of relatively warm water from a deep layer of colder and denser water; the assumption that each layer is of constant density is reasonable (see [31]). Within a band of width about 300 km, centred on the Equator, there is, confined to a depth of about 100 m, a westward current that is driven by the prevailing trade winds, below which lies the Equatorial Undercurrent (EUC)-an eastward jet whose core resides on the thermocline. Below the EUC the flow dies out rapidly so that, at depths in excess of about 500 m there is an abyssal region of almost still water (see [32]). This equatorial background state interacts with various oceanic wave propagation modes, including long waves (with wavelengths exceeding 50 km; see [16]) and short waves (with wavelengths of a few hundreds of metres; see [34]). A comprehensive model of this wave-current interaction must accommodate the density stratification as well as the coupling between the waves at the surface and on the thermocline. Observations provide evidence for highly nonlinear regimes of internal wave motion. While explicit nonlinear solutions can be obtained, they do not cope satisfactorily with the complexity of the equatorial flow due to the limitations on the permissible underlying currents (see the discussion in [7]). For this reason, it is necessary to perform approximations that capture the relevant dynamics. The available approaches rely on linear theory, whether they resort to a numerical treatment by means of finite differences (see [32,39]) or to the method of multiple scales (see [7]). We aim to develop an approach that captures nonlinear effects which, in particular, can explain the propagation of observed equatorial solitary-like waves (see [38]). These are of special interest because they are completely missed by linear theory (see Sect. 4.1.3) and, as pointed out in [2], they are much more easily observable than nonlinear perturbations of an oscillatory wavetrain. Note also that geophysical fluid-dynamical considerations show that the Reynolds number is extremely large (see [30]), so that nonlinear effects typically dominate over viscosity. Field evidence for critical levels (locations where the wave speed equals the mean-flow speed) in the near-surface layer above the thermocline, where the flow-reversal of the underlying current occurs is available (see [34]). It is therefore advisable to rely on the nonlinear alternative to the conventional linear viscous boundary-layer approach for the description of the flow in the neighbourhood of critical levels, where Kelvin 'cat's eye' flow patterns appear (see the discussion in [9] for the somewhat simpler case of gravity water flows with constant non-zero vorticity).

Basic modelling assumptions.
A few realistic simplifying assumptions can be made. Firstly, since the Reynolds numbers are typically very large in geophysical ocean flows, it is reasonable to employ inviscid theory (see [30]). Secondly, because we are considering flows in the neighbourhood of the Equator, it is adequate to use the f -plane approximation in the governing equations (see [26]). Thirdly, since field data shows that the meridional velocities are much smaller than the zonal velocities at the Equator (see [20]), we study flow configurations which are latitude-independent and with a vanishing meridional velocity component. Consequently, we investigate two-dimensional inviscid flows which present no variation in the meridional direction, regarding them as wavecurrent interactions due to localised wave perturbations of a pure current background state. The presence of underlying depth-dependent currents places us within the framework of flows with non-zero vorticity. Since for wave-current interactions in which the waves are long compared with the mean depth of the effective flow region the importance of a non-zero mean vorticity preponderates that of its specific distribution (see [13]), the simplest realistic setting is that of flows with constant vorticity above and below the thermocline: negative above, to permit a reversal from the surface westward wind-drift to the eastward-flowing subsurface EUC, and positive below to model a flow that withers with increasing depth. This choice is propitious not just because in two-dimensional flows the vorticity of a particle remains constant as the particle moves about, as one can easily check. Remarkably, for two-dimensional stratified flows with constant vorticity in each layer, a separation of the flow into a pure current and an irrotational wave perturbation can be performed within the framework of the fully nonlinear theory, without recourse to approximations (see Sect. 3.3). This feature permits us to view the nonlinear wavecurrent interaction as an irrotational wave perturbation of a mean flow representing a pure current.

The Hamiltonian perspective.
The ability to pursue in-depth nonlinear studies is contingent upon a rich structure. Since dissipation is not important, this motivates the quest for a Hamiltonian formulation. The Hamiltonian perspective results in a significant simplification, serving as a guide to the choice of new dependent and independent variables in which the equations take their simplest form. Note that the Hamiltonian formulation of the governing equations for two-dimensional gravity water flows was pioneered in [46] for irrotational flows, and was extended to rotational flows with constant vorticity in [43,44]. A Hamiltonian approach to two-dimensional two-layer gravity water flows with a free surface was developed in [10], and recently the rotational counterpart, with constant vorticity in each layer, was obtained in [6] for periodic flows. This paper establishes the validity of a Hamiltonian perspective in the presence of Coriolis effects in the equatorial f -plane approximation, for zonal flows with no variation in the meridional direction which represent localised perturbations of an underlying pure current background state. A suitable nondimensionalisation reveals that in certain oceanographically relevant regimes the geophysical effects are not merely a small perturbation of the governing equations for gravity water waves: scaling ascertains the relative importance of the different components of the flow, and in this geophysical regime the (linear) contribution due to the Earth's rotation to the longitudinal momentum equation balances the linear part of the material derivative, both being of the same order. We will also show that the Hamiltonian framework is adequate for structure-preserving perturbations. In particular, a specific weakly nonlinear long-wave regime turns out to be structure-enhancing since in this setting the dynamics is described by an integrable Hamiltonian system (with infinitely many degrees of freedom). In other regimes one can derive nonlinear models that capture wave breaking.

The Governing Equations
The objective of this section is to present the nonlinear governing equations and to elucidate the basic structure of the associated mean flow. Also, by specifying the associated geophysical scales we can introduce a set of non-dimensional variables which are useful for ascertaining the relative importance of the different components of the flow.
The fundamental model is that of an inviscid two-layer fluid which admits non-zero vorticity. Using the over-bar to represent the physical variables, we choose a coordinate system with its origin at a point on the Earth's surface, with the x-axis horizontally due East, the y-axis horizontally due North and the z-axis upwards (see Fig. 1). Because z N S x y Fig. 1. The rotating frame of reference, with thex-axis chosen horizontally due east, theȳ-axis horizontally due north and thez-axis upward we are considering flows in the neighbourhood of the Equator, we use the f -plane approximation (see [26]). Since observations show that the meridional velocities are much smaller than the zonal velocities at the Equator (see [20]) then, by assumption, our flow configuration is latitude-independent with vanishing meridional velocity component, and so we restrict the domain to be at a fixed latitude (see Fig. 2). Let z = h 1 + η 1 (x, t) be the free surface, z = η(x, t) be the thermocline and z = −h be the flat bed, where t stands for time. Here h 1 and h are the mean depths of the near-surface layer and of the abyssal layer, with typical values about 120-250 m and 4 km, respectively (see [14]). The density of the fluid above the thermocline is constant, ρ 1 ≈ 10 3 kg m −3 , and is replaced by ρ = ρ 1 (1 + r ) for the fluid below the thermocline, where r is a small positive constant; typically r is about 10 −3 (see [22]).
Denoting by (u 1 (x, z, t), v 1 (x, z, t)) and (u(x, z, t), v(x, z, t)) the velocity fields in the layers above and below the thermocline, respectively, the equations of motion are the suitably adjusted Euler equations in D 1 (t), (3.1) where P = P(x, z, t) is the pressure, g ≈ 9.8 m s −2 is the (constant) acceleration of gravity and ≈ 7.29×10 −5 rad s −1 is the (constant) rotational speed of the Earth about the polar axis (towards the East), supplemented by the equations of mass conservation, The vorticity distribution is specified by for suitable (physical) constants γ 1 and γ . The appropriate boundary conditions are, at the free surface, the dynamic and kinematic boundary conditions 7) w 1 = η 1,t + u 1 η 1,x on z = h 1 + η 1 (x, t), (3.8) respectively, where P atm is the constant pressure of the atmosphere at the surface of the ocean. At the thermocline we have the kinematic boundary conditions (3.10) together with the requirement that the pressure is continuous across the thermocline, where P ± refer to the limits of the pressure P from above and below the interface. Finally, at the fixed, horizontal, impermeable bottom, we have the kinematic boundary condition Note that (3.9)-(3.10) ensure the continuity of the normal component of the fluid velocity across the thermocline z = η(x, t). Since within the inviscid setting the stress at an interface has no tangential component, due to the absence of any interaction derived from friction at the interface, in principle it is permissible for the tangential component of the fluid velocity to present discontinuities at the thermocline. Nevertheless, the available field data for equatorial flows suggests a continuous transition between the two layers (velocity discontinuities across the thermocline would correspond to a delta-sheet of vorticity-see also the discussion at the end of Sect. 3.3), so that we also require a tangential velocity balance: Consequently, the velocity field is continuous across the thermocline. However, since the constant vorticities in the upper and lower layer have opposite signs, this continuity property is not valid in the continuously differentiable sense. 14) the associated pressure distribution P(z) being given by (3.16) The current profile (3.14)-(3.15) captures the salient features: a westward surface drift (since |γ 1 |h 1 > γ h) below which resides an eastward jet (the EUC) which overlies Depiction of the flow model in the absence of waves (pure current): constant negative vorticity γ 1 in the near-surface homogeneous layer accommodates a flow-reversal from the westward wind-drift at the free surface z = h 1 to an eastward jet whose core resides on the thermocline z = 0, while a constant positive vorticity in the abyssal homogeneous region below the thermocline permits the adjustment to practically still water at great depths (with no motion on the flat bed z = −h). The corresponding velocity field is given by u 1 (z) = γ 1 z + γ h and w 1 = 0 in the near-surface layer 0 < z < h 1 , with u(z) = γ (z + h) and w = 0 in the abyssal layer −h < z < 0; here γ > 0 > γ 1 and |γ 1 |h 1 > γ h an abyssal layer where a gradual transition to no motion on the flat bed occurs (see Fig. 3). Note that represents the mass transport per unit width (in m 2 s −1 ), so that for there is a net eastward mass transport. Since typically (see below) h 1 /h and γ /γ 1 are both of order O(10 −2 ), an eastward net flow occurs within 1 • from the Equator; this is compensated by a return flow at higher latitudes (see [36]). We seek nonlinear wave-current interactions arising as localised wave perturbations of the background flow (3.14-3.15).

3.2.
Nondimensionalisation. The first stage in expressing the governing equations in a useful form is to introduce a suitable non-dimensionalisation of the variables; to this end we write (3.17) where the length scale is L = 500 m and U 0 = 0.5 m s −1 is an appropriate speed scale. (The omission of the overbar now indicates that we are using the non-dimensional version of the corresponding variable.) Set Above the thermocline, in the region Eqs. (3.1), (3.3) and (3.5) therefore become while below the thermocline, in the region where for a flow reversal from a westward surface wind drift of 0.5 m s −1 to a maximal eastward EUC speed of 1 m s −1 in a layer of mean depth h 1 = 120 m, with an abyssal layer of mean depth h = 4 km; these values, corresponding to h 1 ≈ 0.24 and h ≈ 8, are appropriate for the 6000 km stretch between about 140 • E and 150 • W (see [14]). For a flow that gradually dies out with depth, being motionless at the flat bed, we compute γ 1 ≈ −1.25 × 10 −2 s −1 and γ ≈ 2.5 × 10 −4 s −1 . Written in nondimensional form, the boundary conditions (3.7)-(3.12) become (3.34) where p ± in (3.32) refer to the limits of the function p from above and below the common boundary z = η(x, t) of the regions D 1 (t) and D(t).
The non-dimensional counterpart of the pure current background state (3.14)-(3.15) is with p(z), corresponding to the background pressure distribution (3.16), given by (3.37)

Nonlinear flow separation.
For wave-current interactions the "current" component of the flow is defined as the average velocity, and the localised perturbations that fluctuate around this average are ascribed to the wave motion, the direction of wave propagation being along the x-axis. We view the nonlinear equatorial wave-current interaction as a localised wave perturbation of the pure current solution (3.35)-(3.36) in a domain that is infinite in the x-direction. The localised nature of the wave disturbances is captured by assuming the smooth functions The corresponding velocity fields (u 1 , w 1 ) and (u, w) are to be smooth in the domains D 1 (t) and D(t), with both vector functions admitting continuous extensions to the closure of the domains and across their adjacent boundary. It is remarkable that the interpretation of (3.35)-(3.36) as the underlying current can be achieved within the framework of fully nonlinear theory, without performing approximations.
Combining this with the absence of a flow on the flat bed, u = 0 on z = −h, we get U(z 0 , t) = γ (z 0 + h), which coincides with (3.36). Similarly, the underlying current at a depth z 0 above the crest of the internal wave and below the trough of the surface wave is Let us point out that the velocity field is actually Hölder continuous across the smooth thermocline. Indeed, the velocity field of the pure current background state (3.35)-(3.36) is Lipschitz continuous in the horizontal strip −h ≤ z ≤ h 1 . Since the wave-current interaction is a localized perturbation of this background flow, in an open horizontal strip O containing the thermocline, the velocity field V F of the wave motion is divergence-free, has a piecewise constant curl, decays as it approaches the lateral boundary components of O at infinity, and is smooth on the upper and lower boundary components ∂O ± of O (since in the fluid domain away from the thermocline, the Laplacian of each velocity component is constant, so that the regularity for Poisson's equation [15] applies). For s > 3, the a priori elliptic estimate for the div-curl system, interpreted for V F ∈ L ∞ (O) in the sense of distributions (see [4,42]), ensures, by means of a Sobolev imbedding (see [15]), that V F is Hölder continuous with exponent α < 1/3 throughout O. Consequently the velocity field of the wavecurrent interaction is Hölder continuous across the thermocline. As a first step towards the choice of new dependent and independent variables that reveal the Hamiltonian structure, we introduce the stream function and the perturbed velocity potential. The equations of mass conservation, (3.21) and (3.25), ensure the existence of a stream function in each layer, ψ 1 (x, z, t) in D 1 (t) and ψ(x, z, t) in D(t), determined, up to an additive term that depends only on time, by Relying on (3.21), (3.25), (3.33) and (3.30), we may set The discussion after (3.13) shows that is continuously differentiable in K(t) without being twice continuously differentiable (due to discontinuities in the second order partial derivatives across the thermocline). For this reason, rather than using we prefer to use ψ and ψ 1 . With the vorticity distribution defined by (3.22) and (3.26), (3.39) yields We now introduce harmonic perturbed velocity potentials, ϕ in D(t) and ϕ 1 in D 1 (t), by requiring (3.43) More precisely, since the velocity field is a smooth localised perturbation of the background state (3.35)-(3.36), with u = 0 on the flat bed z = −h, we may set (3. 45) to see that the first integral in (3.45) is well-defined, use (3.34) and subsequently the divergence theorem applied to the vector field Both perturbed velocity potentials admit continuously differentiable extensions to the closures of their domains, ϕ ∈ C 1 (D(t)) and ϕ 1 ∈ C 1 (D 1 (t)). However, in contrast with the stream functions, any oscillation of the thermocline impedes a continuous extension of the perturbed velocity potentials to the bulk of the fluid K(t): while smooth additive functions that depend solely on time may be added to the right sides of (3.44)-(3.45), this process does not lead to an equality of ϕ and ϕ 1 on the interface. Indeed, for η ≡ 0, (3.43) and the continuity of the velocity field across the thermocline show that differentiation of a presumed equality The kinematic boundary conditions (3.29)-(3.31) can now be written as respectively, so that the expressions on which the gradient operates are just timedependent. Recalling (3.40)-(3.41), (3.44)-(3.45), and taking into account that the flow is a smooth localised perturbation of (3.35)-(3.36) with the background pressure distribution (3.37), we evaluate these expressions as x → −∞ on z = η(x, t) and on z = −h, respectively. We get Consequently, we can write the dynamic boundary condition (3.28) as while the continuity of the pressure across the thermocline, (3.32), becomes as the Hamiltonian, where are the excess kinetic and potential energies, respectively. Other than the technical spur to avoid infinite energy, the rationale for the above choice is that not all potential energy can be converted into kinetic energy, the portion that is available for this conversion being the difference in the potential energy between the perturbed wave-current state and the pure current background state. In view of (3.17) and (3.27), the nondimensional expression in (3.54) corresponds to the total excess energy per unit width, is a suitable scale, note that, over the 16,000 km length of the EUC, it corresponds to a mean total energy of about 4 kg s −2 , which is of the order of the mean kinetic energy at the surface of the equatorial Pacific (about 2 kg s −2 ; see [19]). As for the estimate of potential energy, field data gathered at 45 m depth indicates an average ratio between potential and kinetic energy of about 0.76, varying between 0.18 and 3.63 with means of the order of 10 −2 kg s −2 (see [28]). Moreover, the energy is dominated by the zonal component, which is about one order of magnitude greater than the meridional component, thus corroborating the realistic nature of an approach that invokes the f -plane approximation and the neglect of meridional variations. Using (3.43) and (3.38), we get . 4. Sketch of the geometric features that are relevant in the definition of the Dirichlet-Neumann operators associated with the two layers: at the instant t, the outward unit normals for the abyssal layer D(t) and for the near-surface layer D 1 (t) at a point X on the thermocline are N − (X ) and N + (X ), respectively, while the outward unit normal at a point X 1 on the free surface is N + (X 1 ) we define on D 1 (t) the Dirichlet-Neumann matrix operator that associates to the boundary values of the harmonic function ϕ 1 on the lower and upper boundaries, 1 and 2 , respectively, the outward normal derivatives N · ∇ϕ 1 at these boundaries. Then Similarly, the divergence theorem applied to 1 2 From the definition of the Dirichlet-Neumann operators, using (3.43) and (3.29)-(3.31), we see that Adding up the relations (4.9) and (3.65), we obtain Let us now define the operator Introducing the variables relation (3.67) enables us to express and 1 in terms of ξ and ξ 1 : Regarding h 1 , h, γ 1 , γ, ω as fixed parameters and gathering (3.54)-(3.57), (3.59), (3.61)-(3.63) and (3.70)-(3.71), we express H as a functional depending solely on the variables η, η 1 , ξ , ξ 1 : perturbations. This property is passed on to ξ and ξ 1 , while for η and η 1 it is part of our setup. It is of interest to provide a concise explicit form of the functional H in (3.72). Since 1 = (1 + r ) − ξ by (3.69), B and G are unbounded self-adjoint operators on L 2 (R) while (G 12 ) * = G 21 , and using the relation (1 + r ) We now compute variations of the functional H, with respect to the inner product associated to square-integrable real functions defined on the real line. We consider variations of the wave field, regarding the underlying current field (3.35)-(3.36) as fixed; in particular, the flat bed, the vorticity and the mean depths of the layers do not change. Note that On the other hand, respectively. For harmonic variations of ϕ, using the divergence theorem and the identity Similarly, for harmonic perturbations of ϕ 1 we get using the identities (3.85) Since differentiation of (3.38) with respect to the t-variable yields R η t dx = R η 1,t dx = 0, using integration by parts, (3.85) and (3.81), we get Cancellations occur in (3.84) due to the above two relations in (3.84), so that Therefore we proved the following result. u = (η, η 1 , ζ, ζ 1 ) . With H given by the explicit expression (3.73), the system The above variational derivative of H :

Theorem 2. Let
, with respect to the inner product ·, · in the Hilbert space L 2 (R)×

Hamiltonian Perturbations
We regard the equatorial waves as localised perturbations of the background state, typical examples being solitary waves. A reasonable measure of the "wavelength" of a solitary wave is the spatial extent of the region where the deviation of the profile from its mean level is at least 1% of the maximum height. In the geophysical regime in which the wavelength is long compared to the average depth, both gravity and the rotation of the Earth play a role in the dynamics. In our nondimensional framework x = 1 corresponds to 500 m, so that the relevant spatial scaling is with ε 10 −1 . This defines a physical regime in which the dependent and independent variables are O(1), and the relation to the nondimensional variable x is by means of (4.1); in particular, a perturbation that mainly occurs for values x ∈ (−1, 1) corresponds to x ∈ (−1/ε, 1/ε) and thus represents a physical region spreading over more than 10 km. Our approach to systematic derivations of model equations is from the point of view of Hamiltonian perturbation theory. The parameter ε from (4.1) will also be introduced through choices of scaling of the dependent variables (η, η 1 , ζ, ζ 1 ), corresponding to scaling regimes of interest, in which dispersive and nonlinear effects are brought into play. Since our main interest is the propagation of internal waves, the representative scales are determined by the motion of the thermocline, and we use these scales also for the surface waves. This results in a Hamiltonian that is a function of the small parameter ε. The approximating Hamiltonian systems are obtained by retaining a finite number of terms in the Taylor expansion in ε of the Hamiltonian.
Given a constant α > 0 and a smooth local diffeomorphism f : R 4 → R 4 , the change of variables [αt = t, u = f (u)] transforms the Hamiltonian system ∂ t u = J δ u H to the Hamiltonian system ∂ t u = J ∂ u H with Hamiltonian H (u ) = H(u) if the operator J = α(δ u f )J (δ u f ) is symplectic, that is, (J ) J J = J (see [35]); moreover, if the change of variables preserves the Hamiltonian form of all Hamiltonian equations, then J must be symplectic (see [33]). Due to the physical interpretation in which the dynamics is governed by the perturbations of the two interfaces, the relevant coordinate transformation in our setting is the spatial scaling η = aη, ζ = aζ, 2) in terms of the magnitude of the deviation of the surface/interface from the quiescent flat state and of the relative sizes of the wave/current components of the fluid velocity field at these boundaries. We will provide details of this interpretative aspect in each specific case that we discuss. The linear Dirichlet-Neumann operators G and G + are analytic in their dependence on η and (η, η 1 ), respectively, having convergent Taylor series expansions and given a smooth function m : R → C whose derivatives of any order have polynomial growth, the Fourier multiplication operator maps S(R) into S(R) and extends to a self-adjoint operator in L 2 (R) if and only if m is real-valued; see [40].

The operator (4.7) is bounded if and only if m ∈ L ∞ (R), with m(D) f real-valued or imaginary-valued, for real-valued f ∈ S(R), if m(−x) = m(x) or if m(−x) = −m(x)
for all x ∈ R, respectively. In terms of Fourier multipliers the leading order terms of the Dirichlet-Neumann operators G and G + are given by while G (00) where F(x) = f (x/ε). On the other hand, the effect of the scaling (4.2) is expressed by with more intricate but workable scaling formulas expressing G  (h D) and To obtain the linearised form of (3.86) we have to exhibit in (4.15) a functional dependence on u = (η, η 1 , ζ, ζ 1 ) . This can be achieved as follows. Note that (3.82)-(3.83) yield Using the calculus of pseudodifferential operators with smooth symbols, by means of (4.17) we can handle the last three terms in (4.15) by relying on the fact that although the functions and the operators of interest are all real-valued, the natural set-up for Fourier multipliers is the complex inner product on L 2 (R), which is also adequate for taking advantage of the self-adjointness of operators of the type (4.7). We get The linearised equations of motion are Given an initial data u 0 , representing a localised wave perturbation of the pure current background state described in Sect. 3.1, we can solve the linear system (4.18) using the Fourier transform. Indeed, settingf (k) = R f (x)e −ikx dx for f ∈ S(R) in each component of u, the system (4.18) is transformed for any fixed k ∈ R into the linear autonomous system of ordinary differential equations where (4.20) Note that in the absence of vorticity (γ = γ 1 = 0) the matrix M(k) is real but the context of equatorial wave-current interactions (for which γ = 0) brings about purely imaginary entries on the main diagonal. The complex system (4.19) can be cast in the form of a real Hamiltonian system with double degree of freedom (of dimension eight) by separating the real and complex part of each of the four components of the vectorû ∈ C 4 . However, in our context it is advisable to work with the complex linear Hamiltonian system to avoid going beyond the treshold of Galois theory for roots of polynomials. In particular, finding the spectrum of M(k) is computationally within reach for a fixed value of ξ since explicit formulas are available for its quartic characteristic polynomial. On the other hand, since the available theoretical studies of Hamiltonian systems deal with systems in R 2n , there is a need for theoretical adjustment but with guaranteed structural features. For example, the matrix M(k) being Hamiltonian means that J M(k) is self-adjoint, so M(k) = J −1 (−M(k) )J . Therefore the matrices M(k) and M(k) are similar. Thus if (k) is an eigenvalue of M(k), then so is − (k) and the whole Jordan block structure will be the same for (k) and − (k).
The unique solution to (4.19) with initial dataû 0 (k) = R u 0 (x) e ikx dx, given bŷ u(k, t) = e M(k)tû 0 (k) for t ≥ 0, corresponds, by means of the inverse Fourier transform, to the solution of (4.18) with initial data u 0 ∈ S(R). To a purely imaginary eigenvalue (k) = −ikc with c ∈ R\{0} of M(k), with corresponding eigenvector v(k) = 0, we can associate the oscillatory mode solution e M(k)t v(k) = e −ikct v(k) of (4. 19), and for the initial data  (4.19) for |k| < 64 is physically relevant. Note that an eigenvalue (k) with nonzero real part is the symptom of an oscillation mode whose amplitude grows exponentially in time since (k) or − (k) leads to a coefficient of the type e ω(k)t with ω(k) > 0. This Kelvin-Helmholtz instability corresponds to a situation when a wave is capable of extracting energy from the background pure current state, by drawing either kinetic energy from the pure current motion or potential energy from the stratification. 1 Since these instability phenomena are a prelude to mixing and turbulence, it is important to understand whether or not they are inherent to the equatorial wave-current interactions. For considerations of spectral type the explicit formulas for the roots of the quartic characteristic polynomial are too unwieldy. We advocate an approach that takes advantage of the fact that each entry of the matrix M(k) depends smoothly (actually, the dependence is real-analytic) on the parameter k. This ensures that the eigenvalues depend continuously on k, even if splitting and coalescence of eigenvalues may take place (see the discussion in [21], Chapter VII) and the eigenvalues need not be differentiable, which makes it difficult to translate this feature quantitatively. Also, note that the eigenspaces can behave quite singularly-they are not necessarily continuous. However, for simple eigenvalues the eigenvalue itself as well as the (unique) normalised eigenvector will exhibit a real-analytic dependence on k (see [25]). Note that the multiparameter version of this last result fails spectacularly: for simple eigenvalues, a locally Lipschitz-dependence holds but differentiability might fail and it may not even be possible to choose eigenvectors in a continuous way (see the discussion in [24]).
Since our main interest is the propagation of localised perturbations, it is advantageous to use the Paley-Wiener theorem to translate the feature that x → u(x, t) has compact support at every instant into the requirement that k →û(k, t) is an analytic function of exponential type. The general solution of (4.19) is a linear combination of solutions of the form t n e (k)t U, where n is a nonnegative integer, U is a constant vector and (k) is an eigenvalue of M(k). Since the fact that no entry of the matrix M(k) exhibits a superlinear growth in k towards infinity ensures, by inspection of Ferrari's classical algebraic solution to the quartic, the existence of a constant m > 0 such that | (k)| ≤ m(1 + |k|) for all k ∈ R, we deduce that k →û(k, t) remains forever of exponential type if it is initially so. Consequently the linearised system captures the propagation of localised wave perturbations of the pure current background state.
A simple calculation shows that M(0) has (0) = 0 as an eigenvalue of algebraic multiplicity four and geometric multiplicity two, with the eigenspace spanned by the eigenvectors (1 0 0 0) and (0 1 0 0) . This behaviour is quite different from that encountered for waves in the physically relevant regime 0 < |k| < 64.

Lemma. The matrix M(k) has four distinct purely imaginary eigenvalues for
To investigate the roots of the characteristic polynomial p(λ) = det(R − λI 4 ), we first perform two sets of operations to simplify its structure. Add the first row multiplied by μ to the third row, and the second row multiplied by μ 1 to the fourth row, and in the outcome add the third column multiplied by μ to the first column and the fourth column multiplied by μ 1 to the second column to obtain a determinant edxpressed in terms of X = λ − γ h, corresponding to the wave speed relative to the maximum speed of the EUC, asp (4.23) whereĝ = g+γ 2 1 h 1 −ωγ h and s = s(k) = sech(h 1 k). Expanding the above determinant by minors along the first column yieldŝ since > 0 > and the inequality s < (1 + s) tanh(s) for s > 0, in combination with the relation since the inequality tanh(s) < s < (1 + s) tanh(s) for s > 0 yields , (4.27) while μ > 0 > μ 1 , < 0, and (4.28) An expansion of (4.23) by minors along the first column yieldŝ We conclude from −μ B < 0 < − /(2μ) and (4.24), (4.26), (4.29), that the quartic (4.23), with leading term X 4 , has four distinct real roots, two positive and two negative.

Long waves.
We now describe the main features of the linear wave propagation under the over-arching assumption of long waves, characterised by wavelengths in excess of 16 km and corresponding roughly to |k| < 0.2.
The lemma in the preamble of this section validates the applicability of perturbation theory. Using the approximations for k → 0, in combination with we can approximate the determinant (4.23) by An expansion of the above determinant by minors along the first column yields the polynomial While the intricacy of the coefficients prevents us from using the explicit quartic formulas to gain insight into the location of the roots of p 0 , we can nevertheless take advantage of the approach used to reduce the solution of the quartic to solving the resolvent cubic. More precisely, we seek real numbers A, α and β such that a process factorising the quartic into a product of two quadratics whose roots are readily available. Expanding the brackets on the above right side and comparing the coefficients with (4.32), we get 36) and the two available expressions for 4α 2 β 2 yield a cubic equation for A, It is a daunting task to use Cardano's formulas to identify a real root A of the cubic (4.37) so that α 2 and β 2 , defined in (4.34)-(4.35), are positive. It turns out that a computational approach can be avoided by using the available structure. Indeed, in terms of a = A + g(h + h 1 )/2, we can write the cubic A(a) on the left side of (4.37) in the form Since . Thus A has a negative root a 0 > −(γ 1 h 1 + γ h) 2 /8 and the cubic (4.37) has a root A 0 such that Since γ 1 h 1 + γ h < 0, corresponding to A 0 we get from (4.35) a value α 0 < 0 with while (4.34) yields a value β 0 > 0 with Writing now (4.33) in the form we see that the four real roots of the quartic p 0 are given by By (3.17) and (3.27), the numerical range of the relevant physical constants is (4.43) We will first obtain approximations of the roots (4.41)-(4.42) taking only into account the relative sizes of the above parameters, more precisely, the fact that The specification (4.43) singles out the details for the ocean dynamics in the equatorial Pacific but other values, subject to (4.44), are relevant for flows in the equatorial regions of the Atlantic and Indian ocean, respectively. The estimates (4.38), (4.39) and (4.40) yield so that X 2 3 + X 2 4 ≈ γ 2 1 h 2 1 + γ 2 h 2 . Using now (4.32) to write p 0 (X ) = 0 for X = X 3 and X = X 4 as we can approximate X 3,4 by the roots of the quadratic 1+r , that is, . (4.46) Since c = λ = γ h + X and c = c U 0 , we obtain from (3.17)-(3.18) that the dispersion relations corresponding to (4.45) and (4.46) are, respectively, and .
Note that a root X of the quarticp(X ) is related by means of λ = γ h + X to an eigenvalue λ of the matrix R(k), and the corresponding eigenvector . If E is an eigenvector of R(k), then, multiplying the first equation of the system R(k)E = λE by μ and adding it to the second equation yields while the second equation multiplied by μ 1 and added to the fourth equation yields and the first equation reads Using the first two relations in the third yields Taking (4.30) and (4.31) into account, (4.49) becomes For X = X 1,2 , (4.45) and (4.44) yield E 1 ≈ gh X X 3 E 2 ≈ h h+h 1 E 2 , that is, while for X = X 3,4 the quadratic polynomial providing (4.46) yields [(γ − γ 1 )X − rg]hh 1 ≈ −(h + h 1 )X 2 , so that (4.50) simplifies to

and in this case
for example, for the values specified in (4.43), the ratio η/η 1 is about 500 (see Fig. 5).
The above considerations prove the validity of the following result.

Theorem 3.
Within the framework of linear theory, long waves in slow-mode propagate at speeds given by the formula (4.48), with the effects confined to the motion of the thermocline, due to (4.52), while the propagation speed of the fast long-waves is given by (4.47).
We showed that the Paley-Wiener theorem applied to the linear system (4.19) ensures that a localised disturbance occurring in the equatorial mid-Pacific will generate a unique solution in the form of a localised wave perturbation of the underlying pure current background state, that we can regard as a superposition of periodic modes by means of (4.22). Since in the ocean, energy tends to be concentrated in the lower frequencies, the physically most relevant regime is that of long waves, in which the asymptotic procedure implemented above supresses high frequencies by acting like a low-pass filter. From (4.47) and (4.48), using the values in (4.43), we obtain four possible propagation speeds in this regime: The speed values (4.53) correspond to the extremely rare event 2 of tsunami waves: if triggered by a submarine earthquake or by a meteorite impact, these would consist, in view of (4.51), of surface waves coupled with internal waves of roughly the same amplitude, propagating along the Equator eastwards at 727 km h −1 and westwards at 720 km h −1 . While these fast waves are exceptional, note that on Sunday, 22 May 1960, several powerful earthquakes in rapid succession, occurred along 1000 km of fault parallel to the Chilean coastline, with the epicentre at about 200 km off the coast of Central Chile (see the discussion in [5]). Tsunami waves were generated which propagated across the Pacific Ocean and records of the travel time of the tsunami indicate a speed of about 720 km h −1 , which matches well our predictions. On the other hand, the two values (4.54) represent ubiquituous slow internal waves (known as eastward Kelvin waves and westward Rossby waves): (4.52) yields that the corresponding surface waves are insignificant, given that the typical amplitude of an internal wave is about 20 m. Note that the predicted speed values (4.54) for the slow internal waves fit reasonably well with reported field data (see the discussion in [7]). They highlight an important feature of the equatorial ocean dynamics in the Pacific: the eastward internal wave speed exceeds the 2 While four out of every five tsunamis happen in the Pacific Ocean, initiated by earthquakes along the boundary of one of the Earth's main geologically active tectonic plates, those in the tropical Pacific tend to be modest in size. However, a few times every century, tsunami waves generated by great earthquakes in the North Pacific or along the Pacific coast of South America sweep across the entire Pacific. As for large asteriod impacts into the deep ocean, these are very rare-the only currently known deep-ocean impact event on Earth being the Eltanin asteroid impact and tsunami, ocurring about 2.5 million years ago when a asteriod exceeding 1 km in diameter did strike the ocean about 1500 km SSW of Chile (see [17]). thermocline free surface Fig. 5. Linear long-wave theory: coupled localised wave perturbations over a large spatial scale exist but the amplitude of the internal wave exceeds more than 500 times that of the surface wave. Typically detected oscillations of the thermocline from its mean level are about 10 m maximum speed of the EUC, while the westward internal wave is slower than the maximum speed of the surface wind-driven current. Consequently, the dynamic response of the equatorial Pacific presents an east-west asymmetry due to the interaction between waves and the depth-dependent underlying currents.

Short waves.
High-frequency ocean waves, with wavelengths of 50-100 m, correspond to the range 31 < |k| < 64, in which the approximations tanh(h 1 k) ≈ 1 and tanh(kh) ≈ 1 are reasonable, so that (k) ≈ 1 (2+r )|k| and 1 (k) ≈ 1 |k| . On the other hand, the fact that s → s/ sinh(s) is decreasing for s > 0 yields 0 < max{|μμ 1 |, |μk|, |μ 1 k|, k 2 } sech(h 1 k) (k) in particular, sinh(kh 1 ) > 817 for k > 31 yields the upper bound 0.0189 above. Retaining only the leading order of each entry, we approximate the matrix M(k) in (4.19) by with , 1 < 0. The characteristic polynomial has the four disjoint purely imaginary eigenvalues without any effect on the (flat) free surface, and two modes in which periodic harmonic oscillations of the free surface propagate at the speeds with an undisturbed (flat) thermocline (see Fig. 6). The appearance of dispersive effects is in marked contrast with the long waves: the initial disturbance can be thought of as being composed of the sum of a great many modes and, as time passes by, due to (4.57)-(4.58), the longer components (that is, those for which |k| is relatively smaller) will propagate faster than the shorter modes. Let us point out succinctly the practical meaning of the nondimensional formulas (4.57)-(4.58). In view of (3.17), the nondimensional travelling mode e ik(x−ct) , of wavelength 2π |k| , corresponds in physical variables to a harmonic oscillation of wavelength 2π |k| L, propagating at the speed c = c U 0 . More precisely, due to (3.27), (3.18), (3.17), (3.82), (3.83) and (4.20), the dimensional counterpart of (4.57) for the wavelength , (4.59) while that of (4.58) is (4.60) The previous considerations prove the following result. Note that formula (4.60) predicts that sinusoidal surface waves with a wavelength L = 100 m (corresponding to k = 10π ) propagate at speeds These values of the wavelength and of the wave speeds are realistic for wind-generated surface waves in the Pacific (see [23]). On the other hand, by Such values, close to the maximal speed of the Equatorial Undercurrent, are reported in field data (see the discussion in [7]).

Solitary waves.
For the linear system (4.18) the existence of a smooth localised wave solution u(x − c 0 t) which propagates without change of shape at the constant speed c 0 is equivalent, upon applying the Fourier transform, to (k) = −ikc 0 being an eigenvalue of M(k) with corresponding eigenvectorû(k) for every k ∈ R for whicĥ u(k) = 0. In physical terms this means that the internal solitary wave profile η ∈ S(R) is expressed in Fourier integral form η(x) = 1 2π Rη (k) e ikx dk as consisting of a superposition of harmonic wave components 1 2πη (ξ ) e ikx of wavelength 2π/|k| and amplitude 1 2π |η(k)|, each travelling out at the same velocity so that the shape of the initial localised disturbance is not altered as time passes by, the only effect being a mere translation in the direction of wave propagation. We now prove that linear theory does not capture the solitary wave phenomenon.  ∈ (a, b). Consequently it must be identically zero, so that p(c 0 ) = 0 for all k ∈ R. Since lim k→∞ sech( Consequently we should have c 0 = γ h or c 0 = 1 . Letting k → 0 in the identity (4.24) yields which rules out the possibility c 0 = γ h. On the other hand, from (4.23) we can express lim k→0 p( 1 ) expanding by minors along the last row yields to obtain the value so that c 0 = 1 . The obtained contradiction completes the proof.

Weakly nonlinear models.
Even the most powerful computers cannot handle the full range of space and time scales that must be resolved to successfully simulate numerically the real-world equatorial flows. Model equations and conceptual simplifications can be brought to bear on the problem, thereby providing us with a considerable level of understanding of various physically relevant regimes for equatorial wave-current interactions. The results of the previous subsection show that linear models are not sufficient to capture the solitary wave phenomenon. Recall that linear theory corresponds to retaining only quadratic terms in the expansion of the Hamiltonian, with the expansion parameter describing the amplitude of the wave perturbation of a background pure current state. Weakly nonlinear models can be derived by expanding the Hamiltonian up to cubic terms in the wave amplitude, thereby permitting nonlinear interactions to become relevant.
It is more convenient to write the governing Hamiltonian system (3.86) in terms of the variables (q, q 1 , η, η 1 ), where so that [q − (μ + γ − γ 1 )η|/r and q 1 − μ 1 η 1 are the tangential velocities of the wave perturbations at the interface and free surface, respectively. Therefore the variables (q, q 1 , η, η 1 ) basically encode the profiles and the propagation speeds of the interface and of the free surface, whilst taking into account the vorticity values above and below the thermocline as well as the effect of the Earth's rotation. Note that (3.38) in combination with (3.82)-(3.83) ensure q, q 1 ∈ S(R). Since for q = ζ x andq =ζ x , using integration by parts, we see that Similarly, since (3.82)-(3.83) yield ξ x = q − μη and ξ 1,x = q 1 − μ 1 η 1 , we get 1 to be small parameters, we consider the shallow-water (long wave) regime defined by the spatial scale x = εx, (4.64) and the wave perturbation scales (4.65) in which η and q, as well as η 1 and q 1 , are considered of the same order of magnitude, with the surface wave having a smaller amplitude than the internal wave-an impulse that generates wave perturbations of the background pure current state will typically have a pronounced oscillatory effect at the thermocline 3 whereas the sea surface is hardly affected at all. 4 The interpretation of the scaling (4.65), in which the typical internal wave amplitude ε 2 and wavelength ε −1 are in quadratic balance, is that x , t, η , q , with Hamiltonian H = Hε −3 . Due to (4.8), (4.9), (4.13), (4.14), under the scaling (4.65) the Taylor expansion for the Dirichlet-Neumann operator G(η) is 3 The impulse causes the denser water to rise above the flat equilibrium thermocline level, but gravity pulls it back down because the water beneath the thermocline is heavier than that above it, and the inertia acquired during the falling movement causes the water to penetrate below the level of equilibrium of the thermocline, thus setting up an oscillatory motion. This wave disturbance propagates horizontally since the local increase in pressure that occurs where a parcel of water rises above the equilibrium level pushes the fluid below it away from that place, generating in this process an acceleration (and thus also a force) which also comprises a horizontal component. 4 Since the internal restoring buoyancy force is much smaller than that at the surface, the density variation across the thermocline being about one-thousandth of that across the air-water interface. since tanh(s) = s − 1 since (4.64) and (3.82)-(3.83) yield ξ = εξ and ξ 1 = ε 2 ξ 1 , and the only pseudodifferential operators for which we have to take into account one term beyond leading order are G 22 and On the other hand, in analogy to (4.62), the relations where Therefore the wave speed c of a linear travelling wave, in which the (x , t)-dependence is solely in terms of (x − ct), satisfies the quadratic equation which, in view of (4.76), has the solutions Note that g γ hω, so that the solutions are real. Recalling from (3.36) that γ h is the speed of the underlying current at the mean level z = 0 of the thermocline, we see that in (4.88) the plus sign corresponds to the speed of waves outrunning the current (downstream linear waves), while the minus sign corresponds to the speed of waves running counter to the current (upstream linear waves). At this (leading) order c is a constant and dispersion effects are not observable. To deal with nonlinear effects, setting with c given by (4.88), we observe that (4.85)-(4.86) show that q = ((r + 1)c 1 η )/ h at the leading order. We therefore expect that for some constants b 1,2 that are yet to be determined. With the Ansatz (4.90), we substitute q in (4.83)-(4.84) and write both two equations in terms of η only up to the order O(ε 2 ), thus obtaining two evolutionary equations for η , which should coincide up to O(ε). The equality of their coefficients allows us to find the constants b 1 and b 2 , as follows: The resulting evolution equation for η is while the corresponding q can be recovered from (4.90). Although we can write (4.93) in the standard form of the KdV-equation in a moving frame of reference by means of a scaling transformation (see [35]), we prefer to work with the form (4.93) to facilitate the physical interpretation of the mathematical results. One appealing feature of the KdV equation (4.93) is its bi-Hamiltonian structure: the equation has two expressions as a Hamiltonian evolutionary equation, with the two Hamiltonian operators compatible, that is, any linear combination is again a Hamiltonian operator (see [35]). The bi-Hamiltonian structure ensures the existence of infinitely many integrals of motion that are functionally independent. Moreover, the equation can be solved exactly using inverse scattering theory: starting with arbitrary initial data in the Schwartz class, the solution that evolves from this data develops into a finite number of localised solitary waves (solitons) with speeds proportional to the amplitude, plus an oscillatory tail (see the discussion in [5]). Each solitary wave recovers its localised identity after interacting with other solitary waves (re-emerging unscated after such interaction, the only hallmark of the interaction being a phase shift), while the oscillatory tail disperses and spreads out in space. Therefore the solution evolves in to an ordered set of solitons, with the tallest in front, followed by an oscillatory tail. The details of this general picture can be predicted fairly explicitly from detailed knowledge of the initial data. For example, the internal solitary wave solution is a wave of depression, given in the original variables by with A > 0 and A = 1 > 0. We see that the propagation speed of the solitary wave, c + A|A| 3 , exceeds the speed c of the leading-order linear wave by an amount that is proportional to the wave amplitude A. Therefore larger solitary waves travel faster. With regard to the soliton interaction of two solitary waves that are initially separated, with the larger (and faster) one overtaking the smaller (and slower) one, they emerge unscated from the interaction, but each is shifted, this being the hallmark of the nonlinear nature of the interaction-in contrast to mere superposition. More precisely, the faster/slower solitary wave is shifted forward/backward , respectively, where a 1 > a 2 > 0 are the amplitudes of the solitary waves. We point out that case studies for such wave phenomena in equatorial regions are available in the literature (see the interpretation of space shuttle photographs in [47], which show the relevance of soliton theory to the propagation of equatorial internal waves, and the discussion in [2], where it is pointed out that internal solitons are likely to be generated by the relaxation of the trade winds at the western boundary of the equatorial Pacific that triggers the "El Niño" event), although in general the presence of underlying currents is neglected. We see that the effect of the currents is encoded in the coefficients of the equation, so that these currents alter the shape and speed of these solitons (in comparison to the classical irrotational setting).

Inviscid Burgers regime.
Equatorial ocean dynamics present peculiar features when compared to off-equatorial regions. A remarkable anomaly is that, despite high levels of internal wave energy and shear due to the presence of strong, zonal, basinscale currents, study observations show that internal wave breaking is less frequent than one might anticipate, with rates near the Equator less than 10% of those typical at midlatitudes for internal waves of comparable energy (see [18]). Nevertheless, internal wave breaking does occur in the equatorial Pacific: observations indicate that, in a process in which the interactions with the large shear caused by the EUC play a significant rôle (with the EUC enhancing the internal wave activity-see the discussion in [27]), a strong wind stress propagates downward and generates quite significant internal waves (see [45]) which often break. The breaking is preceded by a steepening of the internal wave profile which results in finite-time gradient blow-up. We will now indicate how the Hamiltonian perturbation approach yields in a specified physical regime model equations that capture this type of behaviour.
Taking ε 1 and δ 1 to be small parameters, we consider the shallow-water regime defined by the spatial scale x = εx, (4.94) and the wave perturbation scales (4.95) in which x , t, η , q , η 1 , and q 1 are of the order of magnitude O(1), with the surface perturbations faster but less significant than those at the thermocline (in view of the fact that δ 1). The regime (4.94)-(4.95) captures internal waves of larger amplitude than the long-wave regime (4.64)-(4.65), transforming (4.63) into the Hamiltonian system r + 1 q η + A 1 η 2 +δ 2 (g − ω 1 )(η 1 ) 2 dx + εδ h r + 1 R q q 1 − μη q 1 − μ 1 η 1 q + μμ 1 η η 1 dx while q can be recovered from (4.106). It is well-known (see, for example, the discussion in [12]) that any initial data in the Schwartz class will lead to finite-time blow-up for the solution to (4.108) in the form of wave breaking: the solution remains bounded but its slope becomes unbounded in finite time.

Concluding Discussion
Within the framework of two-dimensional equatorial flows in the f -plane, with no meridional variations, we presented a Hamiltonian formulation of the nonlinear governing equations for wave-current interactions in a two-layer inviscid fluid with a flat bed and a free surface, for flows with constant vorticity in each layer. A key step was to prove, without recourse to approximations, that the flow can be viewed as an irrotational perturbation-possibly of the same order-of a mean flow that represents the underlying current field. We have also derived several simplified models for the equatorial wave-current interaction across the Pacific through systematic structure-preserving perturbation theory. These have genuine potential as simplified diagnostic and prognostic models since they accommodate all the salient features of equatorial ocean dynamics in the Pacific (strong stratification, as well as a current field with flow-reversal) and are able to capture nonlinear effects. In particular, at some specific geophysical scales the weakly nonlinear long-wave regime turns out to be structure-enhancing since in this setting the nonlinear dynamics is described by a KdV equation-an integrable infinite-dimensional Hamiltonian system. The derived models can discriminate between different physical effects in observations and simulations, dictated by different scales. For example, the variety of oceanographically relevant scales permits us to single out various types of interplay between dispersion and nonlinearity, ranging from regimes in which dispersion and nonlinearity balance each other and allow wave solutions that propagate without change of form to dispersionless regimes that favour wave breaking. The general systematic approach developed in this paper is also useful for developing accurate numerical procedures. Note that, even within the framework of irrotational two-layer inviscid flow with a free surface, the problem of treating computationally the nonlinear interactions between the surface and internal modes of oscillation is challenging (see the discussion in [37]). For example, for fixed densities ρ > ρ 1 of the layers but allowing variations of their mean depths h and h 1 , in the rigid-lid approximation the polarity of the internal solitary waves (that is, whether they are waves of elevation or depression) changes once in the range h 1 /h ∈ (0, 1), while no such change occurs in the coupled configuration (see [1,10]). Furthermore, the presence of non-zero vorticity in each layer complicates considerably the dynamics. Showing that within the nonlinear framework it is possible to split the velocity field into an underlying steady current component and a harmonic wave velocity field opens up the possibility to find numerically solitary wave solutions to the fully-nonlinear problem using a boundary integral method based on Cauchy integral formula (see the discussion in [41]).
The Hamiltonian approach developed in this paper and the above mentioned topics open up exciting perspectives for future interdisciplinary research on equatorial wavecurrent interactions.