Geophysical internal equatorial waves of extreme form

The existence of internal geophysical waves of extreme form is confirmed and an explicit solution presented. The flow is confined to a layer lying above an eastward current while the mean horizontal flow of the solutions is westward, thus incorporating flow reversal in the fluid.


Introduction
Hydrodynamic waves of extreme form have played a prominent role in the mathematical study of water waves since the 19th century, when Stokes [39] first predicted the existence of the wave of greatest height. Extreme waves are inherently nonlinear, for instance extreme Stokes waves are solutions of full Euler equations for inviscid, incompressible fluids, and their existence (cf. [1]) relies on the nonlinear character of these governing equations. Moreover, extreme waves incorporate numerous mathematical intricacies [3,40], especially in regions close to the wave crest, which presents itself in the form of a cusp. Nevertheless, numerous physical features of such flows have been analysed mathematically, for instance the particle trajectories and pressure in extreme Stokes waves in deep-water and water of finite depth may be found in [9,28,29,30,31]. In this work we will investigate internal geophysical flows of extreme form near the equator, and it will be found that the presence of a cusp at the wave crest of these flows once again introduces several mathematical difficulties. In particular, the necessary criteria to apply the implicit function theorem used to prove the existence of smooth geophysical flows are absent in the extreme case, while the cross-section of the wave about the equator also requires a more intricate analysis, both of which are addressed in this work.
The mathematical analysis of geophysical waves in the equatorial region has recently witnessed an explosion of research activity. A significant motivation for this interest has been the fact that such waves have an exact description in a variety of settings in the form of equatorially trapped Gerstner like solutions, see [23] for a review of three-dimensional Gerstner like solutions in the equatorial region. Equatorially trapped surface waves have been investigated in [7,13,14] while surface waves in the presence of underlying currents have been treated in [22,25]. Internal geophysical flows are investigated in [10,4,33,32,36,27], while instabilities in such flows are addressed in [12,26,18,24].
The waves we shall consider are inherently nonlinear and while this nonlinearity may preclude an exact analysis of their behaviour in many instances, this is not strictly the case, with the Gerstner solution [19] being one of the few known exact solutions of the governing equations for incompressible, inviscid fluids in deep-water. The Gerstner wave was subsequently rediscovered by Rankine [34] and of particular interest it was extended to heterogeneous fluids by Dubreil-Jacotin [17]. A more modern perspective of Gerstner's wave is found in [5,6,21] where the evolution of the fluid domain under the action of this solution is analysed using a combination of analytical and topological methods.
Of particular significance to the current work, Gerstner's solution may also be extended to hydrodynamic systems described in rotating coordinate systems. Gerstner-like solutions were initially used in [8,10] to model three-dimensional equatorially trapped geophysical flows and in [11] a Gerstner-like solution was constructed as an exact solution for internal geophysical water-waves near the equator in the β-plane formulation. This explicit solution is constructed in terms of Lagrangian variables (see [2] for a comprehensive review of the Lagrangian approach to fluid dynamics), and it is found that these waves may propagate in an eastward direction and are also found to decay significantly in the meridional direction away from the equator. Moreover, this Gerstner-like solution exists in a layer of finite depth, lying directly above a uniform current layer, which in turn lies above a transitional layer where the uniform current decays to zero when the motionless layer begins. The uniform current layer and the transitional layer are both subject to eastward currents, while it was shown in [11] that the Gersner-like solution induces a mean current towards the west, and so the model accounts for flow reversal, as observed in the Pacific equatorial undercurrent.
Notably, the results presented in [8,10,11] concern nonlinear flows which are smooth throughout. In the current work we aim to extend the results presented in [11] to incorporate the case of nonlinear extreme waves in the Gerstner-like solution, which feature cusps at the wave-troughs along the equator in the thermocline, the interface between the Gerstner-like solution and the uniform current layer. The presence of these cusps leads to a number of mathematical complications, in particular it means the criteria necessary to confirm the existence of the extreme wave profile in the thermocline by means of the implicit function theorem are not satisfied along the equator. This issue is overcome in this work thus confirming the existence of extreme waves in the thermocline. Based on this, we proceed to outline certain features the meridional cross-section of the thermocline must display. Finally, we also confirm that current flow reversal across the thermocline is preserved, even when it is disturbed to the form of an extreme Gerstner-like solution.

Preliminaries
2.1. The governing equations in a rotating frame. The rotation of the Earth about the polar axis induces a fictitious force, known as the Coriolis force, on all bodies away from the equator. In the following, the Earth is approximated to be an ideal sphere with radius R = 6731 km and rotating with angular velocity Ω = 7.29 × 10 −5 rad s −1 about the polar axis. The acceleration due to gravity is assumed uniform over Earth's surface with numerical value g = 9.8 m s −2 .
x y z Ω Figure 1. The rotating (x, y, z)-coordinate system fixed to the surface of the Earth. The x-axis points due-east, the y-axis points due-north and the z-axis points vertically upwards from the surface.
In Figure 1 we sketch the surface of the Earth with a Cartesian coordinate system, (x, y, z) constructed at some point about the equator. The x−axis of this coordinate system points eastwards parallel to the equator, while the y-axis points due-north towards the pole. The z-axis is the line perpendicular to the surface of the Earth at this point. The latitude of the equator is set to be 0 • , while the lines of latitude approximately ±5 • either side of the equator denote the boundary of the equatorial region.
In the following, we shall consider a fluid in motion relative to the rotating (x, y, z)-coordinate axes with associated orthonormal basis vectors {e x , e y , e z }. This rotation is defined in relation to a reference frame which is stationary with respect to the centre of the Earth and whose orthonormal basis vectors we denote by {i, j, k}. The coordinate transformation from the inertial to the rotating frame of reference is given by   e x = − sin φi + cos φj e y = − sin θ cos φi − sin θ sin φj + cos θk e z = cos θ cos φi + cos θ sin φj + sin θk with φ ∈ [0, 2π] and θ ∈ − π 2 , π 2 respectively denoting the zonal and meridional coordinates.
The inertial coordinate system may be chosen such that the rotation vector of the Earth is given by Ω = Ωk, where k is the unit vector projecting from the north pole along the polar axis. Inverting the transformation given by equation (1), the rotation vector Ω is given by (2) Ω = Ω cos θe y + Ω sin θe z with respect to the rotating reference frame. Moreover, given any rotating orthonormal basis {e x , e y , e z } with rotation vector Ω it is always the case that the basis vectors satisfy and the reader is referred to [20] for a comprehensive derivation of this formula. Thus a fluid particle following some trajectory will have a corresponding velocity and acceleration given by The acceleration of a fluid element, a, has two contributions due to fictitious forces arising from the rotation of the Earth. The terms proportional to Ω 2 correspond to the centripetal acceleration, which is the acceleration of a body away from the surface of the Earth due to the rotation of the surface. In practice, this term is negligible since it is effectively cancelled by the gravitational acceleration of the object, although, the effects of this force are observable in the deformation of the Earth away from a perfect sphere (see [16] for further discussion). The term proportional to Ω is the result of the Coriolis force, and in the following we shall only consider the effect of this fictitious force.
2.1.1. The Euler equation in rotating coordinates. Given equation (4) and the effective cancellation of centripetal and gravitational forces on the surface of the Earth, it is found that the familiar material derivative is modified according to thus modifying the Euler equation according to The hydrodynamic pressure in the fluid is represented by P with the fluid density denoted by ρ, while u, v and w denote the velocity field components.
In the following, the meridional coordinate is restricted to the interval −5 • < θ < 5 • , and in this equatorial region the approximation cos θ 1 sin θ y R may be used in equation (5), cf. [7,10,11,13]. Introducing the parameter β = 2Ω R this approximation yields the so-called β-plane formulation, given by where the velocity field is restricted by v = 0 and u y = w y = 0, meaning the equatorial flows under consideration are effectively restricted to zonal motion. While there is no motion in the meridional direction, this does not exclude variation in that direction, and in particular the depth of each interface between the layers of the flow will depend on latitude, thus making the flow three-dimensional. Incompressibility of the flow is ensured when while mass-conservation along the flow follows from (8) ρ t + ∇ · (uρ) = 0, and the system (6)-(8) constitute the governing equations of the model. The boundary conditions coupled to these governing equations are of the form where κ > 0 is constant, and whose relevance will be clarified in section 4. The first of these boundary conditions ensures the pressure on the free surface z = η(x, t) is given by the constant atmospheric pressure P 0 . The second boundary condition on this free surface implies that a fluid particle always remains on this surface. The final kinematic boundary condition ensures the fluid motion ceases below a certain depth, which demarcates the motionless layer.

The flow about the equator
The model resembles that presented in [11], in that the waves occur in a fluid with two regions of different densities. The interface where the density changes is the thermocline, and it is this change in density which allows internal waves to arise. These internal waves propagate towards the east and move with constant wavespeed c which will be prescribed by the dispersion relation in equation (30). The waves, as already pointed out, have zero meridional velocity, and occur beneath a surface layer L(t) whose behaviour is heavily influenced by wind effects. The behaviour of the fluid in this layer is not under consideration in this work. The lower boundary of the layer is described by the surface z = η + (x − ct, y), and beneath this boundary there are four regions wherein the fluid motion exhibits distinctive behaviour.
The upper most layer we investigate in this work is given by Wind effects on the internal flow in this work will be neglected for simplicity, however, it is known that atmospheric motion can play a role in the dynamics of internal flows in stratified layers, cf. [13,32]. The fluid density within this region is of constant value ρ 0 , while the lower boundary of this region corresponds to the thermocline. Beneath the thermocline, the uniform flow layer is bounded from below by the surface where d is constant. In this region the fluid density increases slightly to ρ + > ρ 0 . The typical relative change in fluid density observed in the equatorial Pacific is of  the order (see [37]). Beneath this uniform flow layer, the transitional layer is bounded below by the surface wherein the fluid is again of constant uniform density ρ + . The final layer, is unbounded below and corresponds to the motionless layer, wherein the fluid motion vanishes completely, and the fluid density is again ρ + throughout. The layers of the fluid domain and the boundaries between those layers are shown in Figure 2. In particular, the thermocline exhibits a cusp at the wave trough, which is one of the distinctive features of extreme waves. In the case of extreme waves, the curve described by the thermocline is a cycloid, which is the path traced out by a particle at a fixed point on the circumference of a circle which rolls without slipping. In the case of smooth internal waves investigated in [11], the thermocline wave profile is a trochoid, which is the path traced out by some point along the radius of a circle which rolls without slipping, cf. [6].
3.1. The layers beneath the thermocline. The behaviour of the flow in the layers beneath the thermocline is identical to that in [11], and we outline here the relevant features of these layers for the sake of completeness. The significant difference between the wave of extreme form and the internal waves previously investigated in that work occur in the layer M(t), in particular along the thermocline, due to the presence of cusps at the wave-trough as illustrated in Figure 2. 3.1.1. The motionless layer. At any meridional location y the upper boundary of the motionless layer is given by and beneath this boundary the fluid is still with u = v = w = 0. In this case, the hydro-static pressure is given by The constant parameter D is the depth of the upper boundary of this layer, as measured along the equator y = 0.
3.1.2. The transitional layer. The upper boundary of this region is given by the surface where the constant parameter d is the depth of the boundary as measured along the equator. Within this region the fluid is motionless in the vertical and meridional directions, however, the fluid does move in the zonal direction, with a velocity which depends on the depth and the meridional location. Specifically, the horizontal velocity in this layer is given by The vertical component of the Euler equation within this transitional layer now becomes (14) 2Ωcρ and upon integrating we have The pressure on the lower boundary of this region, P | z=−D+ κβ 4Ω y 2 , is deduced from the continuity of the pressure across the interface of two layers in the fluid domain, along with the expression (12) for the pressure in the motionless layer.
the flow is uniform, given by Replacing this velocity field in the Euler equation (6), we deduce that having integrated each component and imposed equation (15) on the lower boundary z = −d + κβ 4Ω y 2 .
3.2. The layer above the thermocline. The layer immediately above the thermocline, given by η 0 (x − ct, y) ≤ z ≤ η + (x − ct, y) is fully dynamic, in that the velocity field does not reduce to a simplified form as in all the layers beneath the thermocline. In this layer, it is convenient to analyse the flow in terms of Lagrangian coordinates, (q, s, r), which serve to label the individual fluid particles within this layer. The flow in this layer is described in parametric form by k is the wavelength of the flow, while the wave speed c is deduced from the continuity of the pressure across the thermocline. The function f := f (s)  controls the decay of the fluid disturbance away from the equator which corresponds to the line s = 0. The parameter space variables (q, s, r) ∈ R × (−s 0 , s 0 ) × (r 0 (s), r + (s)) may be interpreted as the initial locations of a motionless fluid body prior to being disturbed by a remotely generated wave propagating in the form (18). The functions r 0 (s) and r + (s) characterise the surfaces z = η 0 (x, y, t) and z = η + (x, y, t) respectively, while s 0 ≤ 250 km is the typical meridional extent of such disturbances (cf. [11]). Under the action of this flow, the particle labelled (q, s, r) describes a circular path in the plane of fixed s, with centre (q, r − d 0 ) and radius 1 k e −k(r+f (s)) which decreases with depth, as illustrated in Figure 3.
The dynamic possibility of the flow (18) has been confirmed in the context of surface waves in the work [38], and in the context of internal geophysical waves in [35]. In either regime, as time advances, the mapping (18) remains a diffeomorphism from the parameter space (q, s, r) to the fluid domain when r ≥ 0. Meanwhile, the wave profile described by (18), at any fixed (s, r, t) is characterised to a large extent by the radius e −k(r+f (s)) . Crucially, for this profile to correspond to the graph of a function it is necessary that r + f (s) ≥ 0, with the wave of extreme form emerging when r + f (s) = 0. When r + f (s) > 0 the wave-form is trochoidal, while the case r + f (s) < 0 yields a curve which intersects with itself and so does not correspond to the graph of a function, cf. [6,11].
Introducing the parameters the Jacobian matrix of the coordinate transform (18) is of the form while corresponding Jacobian is given by We note that this Jacobian is time-independent indicating that the flow is volume preserving. In the case of smooth waves studied in [11], the Jacobian satisfies J > 0 at all points in the closure of the layer M(t), and this allows the use of the implicit function theorem to infer the existence of an associated profile r 0 (s). In contrast, for extreme waves this strict inequality now becomes J ≥ 0, with the Jacobian vanishing on the thermocline and specifically along the equator, which gives rise to a number of mathematical complications. In particular it precludes the use of the implicit function theorem to confirm the existence of waves of extreme form propagating in the thermocline, an issue which is addressed in the Section 3.3.
A particularly appealing feature of the Lagrangian formalism is that the fluid velocity at any point along the flow given in equation (18)  In terms of the labelling variables (q, s, r), the gradient of the hydrodynamic pressure P is obtained via (23) ∇ (q,s,r) P = ∂(X, Y, Z) ∂(q, s, r) ∇ (x,y,z) P, and so applying equations (19)- (23) to the β-plane approximation (6), we find Integrating the expression for P q and comparing with the expression for P s , we deduce that .
Moreover the pressure distribution, written in terms of the labelling variables, is given by whereP 0 is a constant.
3.2.1. The wave phase speed. Since the thermocline z = η 0 (x, y, t) is identified with the parameterised surface r = r 0 (s), the pressure in the uniform current layer given by equation (17) combined with the third component of equation (18) evaluated on r 0 (s), give the following expression for the pressure on the thermocline where ξ 0 = k(r 0 (s) + f (s)). Alternatively, the pressure in the layer M(t) given by (26) and restricted to the thermocline r = r 0 (s) also yields Since the pressure is continuous throughout the flow at all times, it follows from expressions (27)-(28) that from which we conclude having retained only the positive root since the waves travel eastward. Crucially this dispersion relation highlights the necessity of a heterogeneous fluid body for the existence of these internal waves, since it is immediately clear that if ρ + = ρ 0 , then c = 0 and the internal waves no longer exist. Decay of the internal wave with meridional distance from the equator is assured when kc − 2Ω > 0, and so it follows from (30) that which combined with inequality (11), ensures that a maximal wavelength for internal waves is of order L = 2π k ≈ 1.2 × 10 7 m. 3.3. Existence of extreme waves in the thermocline. Equating the expressions given in equations (27)-(28) ensures the continuity of the pressure across the interface of the the layer M(t) and the uniform current layer. The equivalence of these two expression allows us to deduce an implicit relation for the function r 0 (s) of the form (32) ρ + κβc 2 s 2 + ρ 0 (kc 2 − 2Ωc) 1 2k e −2k(r0(s)+f (s)) + r 0 (s) Let us denote and differentiating with respect to r we find (33) G r (s, r) = ρ 0 (kc 2 − 2Ωc)(1 − e −2k(r+f (s)) ) > 0, when r > 0.
In particular, smooth nonlinear flows always satisfy r 0 (s) ≥ r * 0 > 0 for all s ∈ (−s 0 , s 0 ), and so under these circumstances it is clear that r → G(s, r) is a strictly increasing diffeomorphism from the open interval (0, ∞) onto the open interval e −2kf (s) , ∞ . This in turn ensures the implicit function theorem may be used to infer the existence of a solution r = r 0 (s) of equation (32).
In addition we observer that G(s, r) is effectively a function of s 2 only, and so the implicit function theorem ensures r 0 (s) is effectively a function of s 2 also, in which case it is clear that r 0 (s) is an even function of s, see [11] for further discussion. Extreme waves emerge from (18) only when r 0 (s) ≥ r * 0 = 0 with r 0 (s) = 0 along the equator s = 0 (in contrast r * 0 > 0 for smooth waves). In this scenario equation (33) now becomes G r (s, r) ≥ 0 for r ≥ r 0 and so significantly G(s, r) is no longer a strictly increasing function of r for all s ∈ (−s 0 , s 0 ) (see also the discussion after equation (20)). Equivalently, the fact that G r (s, r) is not strictly positive everywhere essentially means G(s, r) is not automatically an injective function of r for each fixed s ∈ (−s 0 , s 0 ), and so the implicit function theorem cannot be automatically applied to infer a solution of the relation (32). Nevertheless, even in the case of extreme waves it can still be shown for any fixed s ∈ (−s 0 , s 0 ) that the function r → G(s, r) is an injective map as follows: Proof. Away from the equator we see that When s = 0 we note that G(0, r) is continuous for any r ∈ [0, ∞) and differentiable for any r ∈ (0, ∞). Hence, the mean value theorem ensures that for any pair of values 0 ≤ r 1 < r 2 < ∞, there exists some γ ∈ (r 1 , r 2 ) such that (34) G(0, r 2 ) − G(0, r 1 ) since 0 ≤ r 1 < γ and kc 2 − 2Ωc > 0. Thus G(0, r) is injective for r ∈ [0, ∞). It follows that the function G(s, r) is an injective function of r for every fixed s.
Using the fact that G(s, r) is an injective function on the interval r ∈ [0, ∞) we may now show that there exists a solution of the implicit relation (32).
However, we have already shown that an interval (−s 0 , s 0 ) may always be chosen which ensures the existence of some function r(s) for each s ∈ (−s 0 , s 0 ). Consequently we may always choose an appropriate sub-interval (s 1 −δ, s 1 +δ) ⊂ (−s 0 , s 0 ) for each s 1 ∈ (−s 0 , s 0 ), so that the continuity of r(s) is ensured (see [15] for further discussion). 4. Physical features of the flow 4.1. The profile near the equator. Given any point (s, r 0 (s)) where G r (s, r 0 (s)) > 0, the derivative of r 0 (s) with respect to s is given by .
We recall that Lemma 3.1 and Theorem 3.2 require r ∈ [0, ∞), and so it follows that r 0 (s) must not decrease, at least initially, as we move away from the equator. Hence, it follows from equation (38) that r 0 (s) ≤ 0 for s < 0 and r 0 (s) ≥ 0 for s > 0, at least about the equator. This in turn requires ρ + κ − ρ 0 e −2ξ0 ≤ 0, for all s = 0 from which it follows that We note from equation (38) if κ = ρ0 ρ+ then r 0 (s) ≤ 0 when s ≥ 0 and vice-versa, meaning r 0 (s) decreases in moving away from the equator, hence the necessity for a strict inequality in (39).
When κ < ρ0 ρ+ , equation (32) has an exact solution of the form where W (·) is the Lambert W function, defined by The function W (x) is real valued for any x ∈ [− 1 e , ∞). The function is double valued on the interval [− 1 e , 0) with an increasing and decreasing branch in this interval, while the function is single valued on the interval [0, ∞). An interesting application of the Lambert W function for a convective Rayleigh-Benard model may be found in [41]. In Figure 5 in the left panel, the graph of (s, r 0 (s)) is shown for a wave length L = 200 m and as expected the function is initially increasing when moving away from the equator. In the right panel, the cross section (Y (q, s, r ( s), t), Z(q, s, r 0 (s), t)) is illustrated and again we see the shape of r 0 (s) reflected in this graph, except at the wave crest where the profile becomes smooth. In Figure 6 the thermocline is illustrated from beneath to emphasise the cusp at the wave trough, one of the hallmarks of extreme waves. We also note from the figures that the wave profile does not appear to be differentiable at the wave trough, again a common feature of extreme waves.
When the wavelength L is sufficiently small, there are values s ∈ (−s 0 , s 0 ) such that the solution r 0 (s) < 0, for instance when L = 100 m we find that r 0 (s) < 0 when |s| 234 km. However, we also note that for |s| 234 km then f (s)+r(s) > 0, meaning the details of Lemma 3.1 and Theorem 3.2 remain valid, and crucially the graph of (X(q, r 0 (s), s, t), Z(q, r 0 (s), s, t)) does not intersect itself. Thus, it is only near the equator that we must insist r 0 (s) > 0, in line with our considerations above.  4.2. Current flow reversal. In [11] it was shown that the mean Eulerian velocity of the solution (18) induces a westward flowing current in the near surface region M (t), and owing to the fact that current in the layer beneath the thermocline is eastward, it was found that such behaviour reflects the features of the equatorial undercurrent (with relevant data in [42]). We shall now demonstrate that this property of the dynamic solution is also a feature of the flow in the extreme case, thus maintaining one of the most attractive features of the flow from a geophysical perspective.
The mean horizontal flow at this latitude s and depth z 0 is then given by