EXACT SOLUTION AND INSTABILITY FOR GEOPHYSICAL WAVES AT ARBITRARY LATITUDE

. We present an exact solution to the nonlinear governing equations in the β -plane approximation for geophysical waves propagating at arbitrary latitude on a zonal current. Such an exact solution is explicit in the Lagrangian framework and represents three-dimensional, nonlinear oceanic wave-current interactions. Based on the short-wavelength instability approach, we prove criteria for the hydrodynamical instability of such waves.


1.
Introduction. The aim of this paper is to present an exact solution as well as criteria for its instability, to the β-plane approximation of the governing equations for geophysical fluid dynamics in a relatively narrow ocean strip at an arbitrary latitude. Geophysical fluid dynamics is the study of fluid motion where the Earth's rotation plays a significant role, the Coriolis forces are incorporated into the governing Euler equations, and applies to a wide range of oceanic and atmospheric flows [15,21,49]. Geophysical fluid dynamics contains high complexity, which leads to an inherent mathematical intractability (see the discussions in [4], in which the mathematical intractability of fluid dynamics is illustrated even in the absence of geophysical considerations). In order to mitigate the complexity, it is natural and common to derive simpler approximate equations of Euler equations. There are two approximate models which have been typically employed in the oceanographic considerations. One is the β-plane approximation, which introduces a linear variation with latitude of the Coriolis parameter to allow for the variation of the Coriolis force from point to point. This approximation applies in regions within 5 o latitude, either side of the Equator (see [15,21]). Modifications of the standard geophysical equatorial β-plane model equations which incorporate in addition centripetal forces, or a gravitational-correction term in the tangent plane approximation, were derived in [11,25,26]. Another approximate model is the f -plane approximation, which takes a constant Coriolis parameter into account and does not consider the latitudinal variations. This approximation has been applied to oceanic flows within a restricted meridional of approximately 2 o latitude, either side of the Equator [6,10,15]. An f -plane approximation at an arbitrary latitude was taken into account in [45]. In the paper at hand we will consider a β-plane approximation at an arbitrary latitude.
Exact solutions play an important role in the study of geophysical flows because many apparently intangible wave motions can be regarded as the perturbations of them, and relevant information about the dynamics of more complex flows can be extracted by controlling the perturbations. The approach pioneered by Gerstner [20] (for modern detailed descriptions see [3,22]) of finding explicit exact solutions for gravity fluid flows within the Lagrangian framework [2], was extended to geophysical flows too (see the survey [27]). Gerstner-like three-dimensional solutions were obtained in the f -plane approximation at an arbitrary latitude in [45] and very recently in [14,17], the internal wave in [38,39], in the β-plane approximation in [5,7,8,23] or in a modified β-plane approximation [25,26]; for other studies, we refer the reader to [24,30,31,32,43,44]. We mention that the above exact solutions fail to capture, for example, strong depth variations of the flows. A mixture of analytical and topological methods were applied in [46,47,48] to prove that the Lagrangian flow-map describing a number of three-dimensional geophysical exact solutions is a global diffeomorphism, thereby rigorously establishing that the flow description is globally dynamically possible. Very recently, [11,12,13,36] investigate nonlinear three-dimensional models for flows with sufficient freedom and find, in the Eulerian framework, exact solutions that capture the most relevant geophysical features (see also the survey [37]).
Once an exact solution is available, the study of its stability and instability becomes an important issue. Instability occurs when the effect of some disturbance of the forces acting on the fluid grows as time progresses, and it is important for understanding the factors that might trigger the transition from the large-scale coherent structure represented by the exact solution to a more chaotic fluid motion. Although it is very important, from the viewpoint of mathematics, the study of the hydrodynamical stability or instability of a flow is difficult, because the fully nonlinear governing equations for fluid motion are highly intractable and there are only a handful of explicit exact solutions (see the discussions in [40]). A rigorous mathematical approach to the problem of stability for general three-dimensional inviscid incompressible flows is the short-wavelength method, which was developed independently by Bayly [1], Friedlander & Vishik [18] and Lifschitz & Hameiri [42]. For Gerstner-like geophysical surface waves in various settings, the short wavelength instability analysis is suitable and elegant; these solutions have been shown to be unstable when the travelling wave profiles are steep enough, the critical steepness being very close to 1 3 -instability results were established in [9,16,19,28,29,33,34] (see also the survey [35]).
The paper is organized as follows. In Section 2, we present the β-plane governing equations for geophysical waves at an arbitrary latitude in the presence of an underlying current. In Section 3, we obtain a Gerstner-like solution to this problem; our The rotating framework with the origin at a point on the Earth's surface with latitude φ, the x-axis chosen horizontally due east, the y-axis horizontally due north (in the tangent plane) and the z-axis vertically upward. results can be regarded as an extension of the exact solution for equatorial waves in the β-plane approximation. The geophysical effects and the constant uniform zonal current influence the dispersion relation of the waves. The waves propagate both eastward and westward. We mention that we do not deal with the meridional decay of solutions, and such consideration will limit the choice of wavespeeds which are admissible, in spite of the apparently general form of the dispersion relation. We make a detailed discussion of the situations encountered in the Northern Hemisphere and the Southern Hemisphere, for admissible following as well as adverse currents. Finally, in Section 4, we prove a wave-steepness instability criterion for these exact waves: if their steepness exceeds a specific threshold, very close to 1 3 , then, certain localised small perturbations grow at an exponential rate and the flow is consequently unstable. The waves which travel from east to west are more prone to instability than those which travel from west to east. In particular, we deduce that an adverse current favours instability in the sense that the threshold on the steepness for the wave to be unstable is decreased compared to the case without current. Conversely, this threshold is increased by a following current.
2. The governing equations. We assume that the shape of the Earth is a perfect sphere of radius R =6378 km which rotates with a constant rotational speed Ω = 7.29 × 10 −5 rad/s round the polar axis towards east. A natural framework to describe our problem is a rotating one. We adopt a local "flat" coordinate system with the origin at a point on the Earth's surface with latitude φ, − π 2 ≤ φ ≤ π 2 , the x-axis chosen horizontally due east, the y-axis horizontally due north (in the tangent plane) and the z-axis vertically upward (see Figure 1). In this coordinate system, Ω = (0, Ω cos φ, Ω sin φ) is the rotation vector of the Earth along the polar axis. The Coriolis parameters are defined by  [21]). For surface water waves propagating zonally in a relatively narrow ocean strip less than a few degrees of latitude wide, it is adequate to use the f -or β-plane approximations. Within the f -plane approximation the Coriolis parameters are treated as constants, but within the β-plane approximationf is constant and for f a linear variation with the latitude is introduced, that is, f + βy, with at the fixed latitude φ. In the regions close to the line of the Equator, β = 2Ω R = 2.28 · 10 −11 m −1 s −1 .
In terms of the Cartesian coordinate system, the full governing equations for geophysical fluid dynamics in the β-plane approximation are the Euler equations together with the equation of mass conservation and with the condition of incompressibility Here t is the time, U = (U, V, W ) is the fluid velocity, P is the pressure, g = 9.8ms −2 is the standard gravitational acceleration at the Earth's surface and ρ is the water's density (which we take to be constant). Denoting the free surface by η(x, y, t) and letting P atm be the constant atmospheric pressure, the relevant boundary conditions at the free surface are the kinematic boundary condition which implies that fluid particles on the free surface remain on the surface for all time, and the dynamic boundary condition which decouples the water flow from the motion of the air above. Finally, we assume that the water is infinitely deep, with the flow converging rapidly with depth to a constant uniform zonal current of strength c 0 , that is, The direction of the current is connected with the dynamics of the exact solution that we present below.
3. Exact solution. We will use the Lagrangian framework for the exact solution.
In the Lagrangian framework, the Eulerian coordinates of fluid particles x = (x, y, z) at the time t are expressed as functions of Lagrangian labelling variables (q, s, r) which specify the fluid particle. We suppose that the position of a particle at time t is given as EXACT SOLUTION AND INSTABILITY FOR GEOPHYSICAL WAVES 4403 in which k is the wave number, c is the wave speed and c 0 is the constant underlying current such that for cc 0 > 0 the current is adverse, while for cc 0 < 0 the current is following. For later considerations, we take, on physical grounds, By remembering that g 2Ω ≈ 6.7 × 10 4 m/s, and for π 2 ≤ φ ≤ π 2 , we have g 2Ω < ĝ f , this is indeed a plausible range for c 0 .
We will prove that the system (7) defines an exact solution of the β-plane governing equations (1)- (6), where the travelling speed c and the function m depending on s are determined below. In (7), the Lagrangian labelling variables take the following values: q ∈ R, s ∈ [−s 0 , s 0 ], to show the restricted region for which the β-plane approximation is applicable, and r ∈ (−∞, r 0 ] such that to ensure that the flow has the appropriate decay properties.
For notational convenience, we set The Jacobian matrix of the transformation (7) is given by therefore, its determinant is 1 − e 2ξ , which is non-zero, and hence the transformation (7) is well defined. Furthermore, as the determinant of the Jacobian is time-independent, the flow defined by (7) is volume preserving and (3) holds.
where M (s) = f c + cβs −f c 0 m s − gm s . Now we give some suitable conditions on the pressure function P such that (13) holds and (7) is indeed an exact solution of the governing equations (1)- (6). Since the condition (5) enforces a time independence in the pressure function at the surface, it is necessary to eliminate terms containing θ in (13) by setting and From (8), we have c 0 > − ĝ f , and since − π 2 ≤ φ ≤ π 2 , we get With the constraints (14) and (15), we can solve the pressure function as The constant terms in (18) have been chosen to ensure the conditions (4) and (5) hold on the free surface, see the subsection 3.2 below. In addition, the condition (14) leads us to the dispersion relation for the wave motion prescribed by (7). By regarding (14) as a quadratic in c, if c 0 = c, we get which are well defined due to (17). The case c 0 = c gives us the classical c = ± g k , that is, the geophysical effects have no influence on the dispersion relation of the wave. If c = c + > 0, the wave travels from west to east and if c = c − < 0, the wave travels from east to west. We observe that for our solution (7), the Coriolis parameter f does not appear in the dispersion relation (14) and its solutions (19), but it appears in the expression (16) of the function m(s).
In the absence of an underlying current, for c 0 = 0, (14), (19) and (16) reduce to For equatorial waves f = 0 andf = 2Ω we recover the result obtained in [23], [26], that is, the dispersion relation (14) becomes and the function m(s) becomes 3.2. The free-surface interface. We now focus on the expression (18) of the pressure function and the fulfilment of the boundary conditions (4) and (5). In Lagrangian variables, the kinematic boundary condition (4) holds if at each fixed latitude s the free surface is given by specifying a value of r, the label q being the free parameter of the curve that represents the wave profile at this latitude. With (18) in view, this is achieved if we show that, at each fixed s, there exists a unique solution r(s) ≤ r 0 < 0 such that P (r(s), s) = P atm , which is equivalent to in which For s = 0, the choice r(0) = r 0 works in (20). For the case s = 0, by condition (9) The inequality (21) takes the form We will look for the geophysical waves satisfying Differentiating ( (9), (14), the inequalities (23) become From (23) and r(0) = r 0 , it follows that r(s) < r 0 for s = 0. Since (9) holds and m(0) = 0, we obtain that, for s = 0, m(s) > 0, or m(s) < 0 close enough to zero.
Let us see for which values of the uniform current c 0 the required restrictions (22) and (24) are satisfied. We will make a separate case discussion in the Northern Hemisphere and in the Southern Hemisphere.

I. Northern Hemisphere
For the Northern Hemisphere, we have thatf > 0, f > 0 and β > 0. We further distinguish the following cases.
Case N1. m(s) > 0 and cc 0 < 0, that is, the current is following.
In this case, (22) (16) and (17) it follows that m(s) > 0 if and only if s > 0 or s < − 2f β . Physically, the possibility s < − 2f β can be excluded, and therefore s > 0. We thus get that inequality (24) is satisfied for all s > 0.
The condition cc 0 > 0 is obtained if: (a) the solution c = c + > 0 and the uniform current 0 < c 0 < ĝ f . Just as in the case N1(a), we can exclude s < − 2f β , and we get s > 0. The inequality (24) is reduced to This condition breaks down for large enough values of s > 0. We are interested in water waves propagating in a relatively narrow ocean strip. We show that (22) We note that A(0) = 0. Thus, to ensure the validity of the inequality A(s) < 0, s > 0, the necessary condition is for s > 0 small enough, which is equivalent to for s > 0 small enough. Therefore, for a given 0 < c 0 < ce 2kr0 , (27) holds for some s ∈ (0, s 1 ] and accordingly (25) holds for s ∈ (0, s 0 ] with s 0 < s 1 . (b) the solution c = c − < 0 and the uniform current − ĝ f < c 0 < 0. In this case, − 2f β < s < 0. Since we consider a relatively narrow ocean strip, we can restrict − f β < s < 0, and therefore βs + f > 0. The inequalities (22) and (24) become now with A(s) defined in (26), and respectively. For water waves propagating in a relatively narrow ocean strip, in the same manner as in the case N2(a), we can show that (26) and (29) are satisfied for small values s depending on the size of the current c 0 . To ensure the validity of (28), the necessary condition is A (s) > 0 for − f β < s < 0 with |s| small enough, which is equivalent to for − f β < s < 0 with |s| small enough. Therefore, for a given ce 2kr0 < c 0 < 0, (30) holds for s ∈ [−s 1 , 0) and accordingly (29) holds for s ∈ [−s 0 , 0) with s 0 < s 1 .
Case N3. m(s) < 0 and cc 0 < 0, that is, the current is following.
In this case, we see easily that the inequality (22) can not be satisfied. Thus, this is a nonvalid case.

II. Southern Hemisphere
For the Southern Hemisphere, we havef > 0, f < 0 and β > 0. The four cases above can be handled in much the same way, the only difference being the range for s. Now m(s) > 0 if and only if s < 0 or s > − 2f β . The possibility s > − 2f β can be excluded for physical reasons. Therefore s < 0, which also yields βs + f < 0. For m(s) < 0, we get 0 < s < − 2f β . This range can be restricted to 0 < s < − f β for narrow ocean strips and again βs + f < 0. The conclusions in the four cases are: Case S1. m(s) > 0 and cc 0 < 0, that is, the current is following. In this case, (22) is obviously satisfied for all s = 0, and if: (a) the solution c = c + > 0 and the uniform current − ĝ f < c 0 < 0, then, the inequality (24) is satisfied for all s < 0.
Summing up, we have got the following theorem: • if the current c 0 with |c 0 | < ĝ f is following, that is, cc 0 < 0 and the function m(s) is positive; for flows with positive wave speed c = c + > 0, the range for the fixed latitude is s > 0, and for flows with negative wave speed c = c − < 0, we are restricted to latitudes in the region − f β < s < 0. • if the current c 0 with |c 0 | < ĝ f is adverse, that is, cc 0 > 0 and the function m(s) is positive, for that values of the current |c 0 | < |c|e 2kr0 , or m(s) is negative, for that values of the current |c|e 2kr0 < |c 0 | < ĝ f ; in this case, the fixed latitude s belongs to the interval [−s 0 , s 0 ], with s 0 > 0 small enough.

II. Southern Hemisphere
• if the current c 0 with |c 0 | < ĝ f is following, that is, cc 0 < 0 and the function m(s) is positive; for flows with positive wave speed c = c + > 0, the range for the fixed latitude is s < 0, and for flows with negative wave speed c = c − < 0, we are restricted to latitudes in the region 0 < s < − f β . • if the current c 0 with |c 0 | < ĝ f is adverse, that is, cc 0 > 0 and the function m(s) is positive, for that values of the current |c 0 | < |c|e 2kr0 , or m(s) is negative, for that values of the current |c|e 2kr0 < |c 0 | < ĝ f ; in this case, the fixed latitude s belongs to the interval [−s 0 , s 0 ], with s 0 > 0 small enough. The free surface z = η(x, y, t) is prescribed at s = 0 by setting r = r 0 in (7), and for any other fixed latitude s = 0 in the intervals mentioned above, there exists a unique value r(s) < r 0 which implicitly prescribes the free surface z = η(x, s, t) by setting r = r(s) in (7).
Finally, we note that the form of the surface wave, for fixed values of s and t, is an inverted trochoid (see [4,5]).
Thus, the vorticity ω = (w y − v z , u z − w x , v x − u y ) takes the form We can see that the vorticity (34) is three-dimensional away from the Equator, although the velocity field (11) is two-dimensional. Moreover, the first and third components in (34) depends explicitly on the latitude φ and the underlying current c 0 . The underlying current c 0 does not feature directly in the second component of the vorticity (34), accounting for the local spin in the y-direction, it plays an implicit role as determined by the dispersion relation (19) for the wave speed c. We also observe that the sign of this component is opposed to propagation wave speed.

4.
Instability. Since the stability issue concerns the evolution of small perturbations with time, it is reasonable to pursue its investigation within a linear framework by neglecting nonlinear terms arising from products of the perturbed quantities. Small perturbations u(t, x), p(t, x) from the geophysical flow (U(t, x), P (t, x)) (see (11), (18)) which satisfies the problem (1)- (6), are governed by the following linearized equations with L f,f ,β given by We examine the evolution in time of the localised and rapidly-varying solutions of the linearized system (35)- (36) in the wave packet form with a sharply-peaked initial condition u 0 of the form where is a small parameter, A, A are vector functions and Φ, B, B are scalar functions. Plugging (37) and (38) into the linearized equations (35)- (36) and equating the terms of the same order in , we get for the leading order terms the relations A · ∇Φ = 0 for all t ≥ 0, B = 0, and the following coupled system, consisting of the eikonal equation for the wave phase and the transport equation for the wave amplitude of the velocity, with the initial conditions Since the next-order terms A and B depend only on A and ∇Φ, and the remainder terms u rem , p rem are bounded, in an appropriate norm, at any time t by functions that can depend on t but are independent of , then, in (37), (38) these terms can be made as small as we want by multiplication with and they are dominated by the growth of the leading order terms. See, for example, [33] for details. The system of partial differential equations (39)- (40) can be written as a system of ordinary differential equations along the trajectories of the basic flow U(t, x) (11). Consider the trajectory (7) passing through a point x 0 : then, straightforward computations show that the system (39)- (40) can be rewritten along the trajectories of the basic flow in the following form where ξ := ∇Φ is the local wave vector, and ∇U is the velocity gradient matrix (33). The initial conditions for the ODE system (41)- (42) are To prove the instability of the geophysical fluid flow (7), it is not necessary to investigate (41), (42) for all initial data. We only need to chose an initial disturbance which can lead to an exponentially growing amplitude A. Let us choose the latitudinal wave vector ξ 0 = (0 1 0) T . For this initial condition, taking into account the expression (33) of the velocity gradient matrix, (41) yields Hence, from (42), it follows that For the chosen initial vector ξ 0 = (0 1 0) T and the condition of orthogonality in (43), we must have A 2 (0) = 0. Thus, the second equation in the above system yields A 2 (t) = 0 for all t ≥ 0. Thus, the system reduces to the following two-dimensional system , which can be written as with By rotating the canonical Cartesian basis with the angle α = kct 2 about the vector ξ(t) from (44), the system (45) can be transformed to an autonomous linear system. We denote by P (t) P (t) = cos(kct/2) sin(kct/2) − sin(kct/2) cos(kct/2) .
the α-rotation change-of-basis matrix, and by Ã 1 In the rotating basis, the system (45) becomes the autonomous system Ȧ where D the time-independent matrix D = .
By (46), the solution to the non-autonomous system (45) is obtained by multiplying the rotation matrix P (t) with the solution to the autonomous system (47). The matrix P (t) being time periodic, the behaviour in time of the amplitude vector A is determined by the eigenvalues of the matrix D. The amplitude grows exponentially with time if D has a positive eigenvalue. The eigenvalues of D satisfy the following equation Therefore, taking into account (9), if e ξ > 2f + kc 2f + 3kc (19) = 3f ± f 2 + 4k(f c 0 + g) then, the amplitude A increases unboundedly in time, the exponential growth rate being λ = 1 2 (2f + 3kc) 2 e 2ξ − (2f + kc) 2 1 − e 2ξ .
Due to (7), the steepness of the longitudinal wave profile, defined as the amplitude multiplied by the wave number, is e ξ . We have proved the following wave-steepness instability criterion: Theorem 4.1. At arbitrary latitude, if the geophysical waves (7) are sufficiently large -their steepness exceeds a specific threshold given by the right-hand side of (48) -then, certain localised small perturbations grow at an exponential rate and the flow is consequently unstable.
In the absence of the underlying current, for c 0 = 0, the right-hand side of (48) becomes 3f ± f 2 + 4kĝ f ± 3 f 2 + 4kg . ( Sincef f 2 + 4kg, the threshold (49)  These considerations suggest that waves which travel from east to west are more prone to instability than those which travel from west to east.
In the presence of the current c 0 = 0, we get In particular, we deduce that an adverse current with cc 0 > 0 favours instability in the sense that the threshold on the steepness for the wave to be unstable is decreased compared to the case without current. Conversely, this threshold is increased by a following current with cc 0 < 0.
For waves near the North Pole,f = 0, the right-hand side of (48) becomes 1 3 and we recover the result [41] for Gerstner's wave.