Stratified equatorial flows in cylindrical coordinates

We construct an exact solution modelling the geophysical dynamics of an inviscid and incompressible fluid which possesses a variable density stratification, where the fluid density may vary with both the depth and latitude. Our solution pertains to the large-scale equatorial dynamics of a fluid body with a free surface propagating steadily in a purely azimuthal direction, and is expressed in terms of cylindrical coordinates. Allowing for general fluid stratification greatly complicates the Bernoulli relation—which relates the imposed pressure to the reciprocal fluid distortion at the free-surface—thereby acting as a constraint on the existence of a solution. Employing the implicit function theorem, we establish the existence of solutions and determine that the requisite monotonicity properties hold for the flow solutions we found. Furthermore, since the fluid velocity and pressure are prescribed by explicit formulae in the framework of cylindrical coordinates, our solution is amenable to analysis by the short-wavelength stability approach, which we investigate.


Introduction
We construct an exact solution for the geophysical dynamics of an inviscid and incompressible fluid in the equatorial region which exhibits a general stratification, whereby the fluid density 3 Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. may vary with both the depth and latitude. Fluid stratification is an inherent feature of physical oceanography, with variable fluid density being particularly pertinent to the consideration of large-scale ocean motions [20,38,49,50]. In most of the ocean, and particularly in the equatorial region, temperature is the primary factor determining the fluid density [36,48]. Salinity is another appreciable factor, however low salinity contributes to density fluctuations in the ocean primarily in high-latitude regions where precipitation is high; in coastal regions where freshwater run-off mixes with surface waters; or in polar regions where sea-ice processes predominate [18,20,38].
The exact solution we present in this paper admits a fluid stratification whereby the density may vary both with fluid depth, and with latitude: this form of density distribution is sufficiently general to capture and model most geophysical fluid dynamics phenomena [20,49,50]. It is known that much of the ocean can be divided into three distinct zones, or layers, with characteristic density distributions. The surface zone (or mixed layer) consists of fluid of roughly uniform density, with the fluid temperature and density relatively constant with depth due to the mixing effects induced by wave-current interactions. This layer typically extends to a depth of approximately 150 m (however, depending on local conditions, the depth may vary from 1000 m to a complete absence). Beneath this layer is the pycnocline, a zone in which there is a rapid increase in the density with increasing depth. The pynocline isolates surface water from the denser layer below, and serves to limit vertical communication within an ocean column. As this density variation is usually linked to a similarly marked variation in water temperature, as is certainly the case in the equatorial region being considered, in many cases the pycnocline coincides with the thermocline. Beneath the pynocline lies the deep-zone (at depths below around 1000 m in mid-latitudes) where there is again relatively little change in water density with increasing depth. Although the most dramatic density variations occur in the pycnocline, and this layer serves to limit vertical communication within an ocean column, it is noteworthy that the effects of ocean stratification are not confined solely to this, or the surface layer. Indeed, ocean stratification generates the slow thermohaline circulation which is intrinsic to the ocean's deepest waters [18,20,38,49]. With regard to latitudinal density variation, we note that the exact solution we construct pertains to flows in the equatorial region where the thermocline is deeper and more pronounced, due to the solar energy available, than at higher latitudes. It is clear that attaining a detailed theoretical understanding of, and a model for, stratified flows is of the utmost practical importance.
A number of recent developments have produced an assortment of exact solutions to the GFD governing equations in various forms, see [2-6, 8-15, 21, 22, 29, 30, 39-47]; surveys of these results can be found in [23,35]. Exact solutions are extremely rare in fluid mechanics, in general, and they offer an invaluable insight into the mathematical structure of a given problem. They are also regarded as a framework from which more physically realistic and observable flows may be constructed. The works referenced above primarily involve homogeneous fluids, and accommodating variable density in the fluid vastly complicates the mathematical analysis of an already intractable problem. Nevertheless, some recent developments have been initiated which construct exact equatorial flow solutions incorporating fluid stratification. In [26,27] exact solutions are constructed which determine equatorial flows whereby the fluid density varies linearly with depth, and is independent of the latitude. In [28] an exact equatorial flow solution is presented in the framework of spherical coordinates which permits general stratification.
The solution we construct below is presented in terms of cylindrical coordinates. The cylindrical coordinate approximation offers a more transparent insight into the properties of the fluid flow compared to the framework of spherical coordinates (see [28]) while still retaining quite an amount of the mathematical structure of the full problem. Cylindrical coordinates, unlike the f− and β−tangent plane approximations, capture in detail the latitudinal variation intrinsic to the Earth's curved geometry [10,11]. The transformation from spherical to cylindrical coordinates involves appreciable modifications to the zonal structure of the solution, however we remark that the azimuthal invariance of the solution we present below may serve to lessen any deficiencies induced by this approximation. Aside from this, the form of the exact solution we construct possesses much freedom in its prescription, and it offers a model which captures the salient features of the EUC [9,10,12,48], as we describe in remark 3.3 below.
The construction of such a physically rich and complex solution generates a number of technical complications in the mathematical analysis that must be rigorously addressed. Chief among these is the Bernoulli-type relation (3.18), which provides an implicit prescription of the relationship between the imposed pressure, and the resulting distortion, at the free-surface. Relation (3.18) is greatly complicated by fluid stratification, thereby providing a constraint on the existence of a solution. We apply a functional analytic argument to establish that the Bernoulli relation ensures a well-defined relationship between the imposed pressure at the freesurface and the resulting distortion of the surface's shape for a wide-range of equatorial flows, and that this relationship exhibits the physically-expected monotonicity properties. The paper is concluded by subjecting the exact solution to a short-wavelength stability analysis. Although the mathematical implementation of the involved machinery is significantly complicated by the presence of fluid stratification, since the fluid velocity and pressure are prescribed by explicit formulae in the framework of cylindrical coordinates our solution is amenable to this technique. We find that, depending upon the initial conditions, some flows are linearly stable, but others are not.

Preliminaries
Let us begin by defining the cylindrical geometry and variables that we use to describe fluid motion in an associated rotating system. In the cylindrical coordinate system we consider, the equator is 'straightened' and replaced by a line parallel to the z-axis, while the body of the sphere is represented by a circular disc described in the corresponding polar coordinates. Thus, in a right handed system, our coordinates are (r, θ, z), where r is the distance to the centre of the disc (representing the Earth's interior), θ ∈ (−π/2, π/2) is increasing from north to south and measures the deflection from the equator, and the positive z-axis points from west to east [10,11]. The line of the equator is chosen to be associated with θ = 0. The corresponding unit vectors in the (r, θ, z) system are (e r , e θ , e z ) and the velocity components with respect to (e r , e θ , e z ) are (u, v, w).
The equations governing geophysical fluid dynamics (GFD) are given by the Euler equation, which in cylindrical coordinates with the origin at the centre of the Earth are expressed as where p(r, θ, z) denotes the pressure in the fluid and (F r , F θ , F z ) is the body-force vector, and the equation of mass conservation Subsequently, the fluid density ρ = ρ(r, θ) is permitted to vary with respect to both the fluid depth and the latitude. To examine fluid motion at fixed locations on the Earth's surface we must incorporate the effects of the Earth's rotation. We associate (e r , e θ , e z ) to a point fixed on the sphere which is rotating about its polar axis, and introduce the Coriolis, 2Ω × u, and centripetal, Ω × (Ω × r), accelerations. Here r = re r , u = ue r + ve θ + we z , and where Ω ≈ 7.29 × 10 −5 rad s −1 is the constant rate of rotation of the Earth. Therefore the total contribution of the Coriolis and centripetal acceleration, which must appear on the left-hand side of the Euler equation (2.1), is given by Additionally, we assume that the body-force is due to gravity alone. Thus, the water motion is governed by the system together with the equation of mass conservation (2.2), and is additionally subjected to the kinematic boundary conditions on the free surface r = R + h(θ, z), and on the bed r = d(θ, z), respectively. Moreover, at the free surface r = R + h(θ, z), we require the dynamic boundary condition

Construction of exact stratified solutions
We now construct an exact solution to the GFD governing equations representing a stratified, steady equatorial flow. As the flow is taken to be purely azimuthal, and invariant in this direction, we seek solutions for which the velocity field (u, v, w) satisfies while the free surface is described by for some unknown function θ → h(θ), and the bed is represented by for some given function θ → d(θ). Furthermore, due to the azimuthal character of the flow, the z-dependence is suppressed and condition (2.6) reads for some given function θ → P(θ). With the velocity field prescribed by (3.1), the system (2.3) becomes Eliminating the pressure between the two equations of the system above we obtain r sin θ(ρrΩU cos θ) r + cos θ(ρrΩU cos θ) θ = (r cos θ)gρ θ .
Denoting Z(r, θ) := ρrΩU cos θ we arrive at the equation from which we infer that c = c 1 c 2 . One can easily check that the latter relation is also sufficient forr,θ to satisfy the system (3.6). We will determine further the constants c 1 , c 2 and c. To this end we impose without loss of generality the conditioñ θ(0) = 0, which implies that c = 1. Therefore, c 1 = c 2 . Moreover, for given r, θ, we seek s 0 ∈ R such thatθ (s 0 ) = θ,r(s 0 ) = r. (3.10) The first equation of (3.10) gives while the second equation of (3.10) implies that In summary, denoting the characteristics (3.6) that satisfy (3.10) by r, θ, we have r(s) = r(cos θ)cosh(s), θ(s) = arcsin(tanh(s)) for all s.  in the neighbouring equatorial region is prescribed by (3.13). A concrete example of this is provided by the equatorial undercurrent (EUC). The EUC is a current which runs the extent of the Pacific equator whereby the surface flow is predominantly westward (due to the prevailing trade-winds) and a flow reversal occurs beneath the surface leading to an eastward-flowing jet which is confined to depths of 100-200 m, see [9,10,12,36,48]. If we suppose that the westward surface velocity has the constant magnitude W w , and the maximal eastward velocity at the core of the EUC jet has magnitude W e , then setting leads to an elementary model for the EUC. Here r = R represents the free-surface at the equator, r = R 0 is the depth at which the EUC attains its maximum speed (this typically coincides with the thermocline) and R = R 0 − (R − R 0 ) W e /(W e + W w ) represents a depth beneath which there is essentially no fluid motion.
In order to determine the pressure p, we note that (3.4) implies that  Observing that Therefore, − H(r cos θ, θ)r sin θ + H(a cos θ, θ)a sin θ.
Comparing the latter formula with (3.15) we find that and thus, the pressure p is given by the formula where A is some constant. This equation, in combination with the surface dynamic boundary condition (3.2), implicitly prescribes a Bernoulli-type relation between the imposed pressure at the surface of the ocean and the resulting deformation of that surface, namely Allowing for general fluid stratification greatly complicates this Bernoulli relation, thereby acting as a constraint on the existence of a solution: we establish the existence of solutions by implementating the implicit function theorem. This is achieved by first passing to the nondimensionalisation of the (3.18) by dividing throughout by the atmospheric pressure P atm to obtain where H, P are non-dimensional functions defined through The non-dimensional equation (3.19), relating the free surface and the pressure, can be written as the functional equation We immediately see that (3.22) stands for the (non-dimensional) surface pressure distribution required to maintain an undisturbed free surface. To prove the existence of non-trivial flow solutions we rely on the implicit function theorem. Proof. To carry out the proof we will rely on the implicit function theorem whose employment requires the computation of the derivative of f with respect to its first variable around the special solution (0, P 0 ). The latter derivative is defined as A calculation shows that (3.23) Due to the size of the involved physical quantities, clearly there is a constantc < 0 such that −gR + ΩR(cos θ)(2w + ΩR cos θ) c < 0.
We infer from the latter inequality that is an onto linear homeomorphism, a fact which enables the implementation of the implicit function theorem.
In order that the solution determined by (3.13) and (3.17) generates physically meaningful equatorial flows, it is necessary to ensure that the Bernoulli relation (3.18) predicates the expected monotonicity properties involving the surface pressure and the resulting surface geometry. The next result achieves this by establishing monotonicity properties shared by the functions P and H. Theorem 3.6. The pressure function P and the height function H satisfy P (θ) < 0 i f H (θ) 0 for some θ ∈ (0, ε), (3.24) and H (θ) < 0 i f P (θ) 0 for some θ ∈ (0, ε). (3.25) Proof. We differentiate now by θ and get The conclusions formulated in (3.24) and (3.25) follow immediately when taking into account that for geophysically plausible values of w.

Local stability of certain azimuthal flows
Since the azimuthal flow solutions (3.13) and (3.17) are represented by explicit, if complicated, formulae, they are amenable to a stability analysis. More precisely, we investigate the time growth of the amplitude of perturbations to basic flows having a velocity field which satisfies the Euler equation (2.1) and the equation of mass conservation (2.2). To reach our goal we will use the short-wavelength perturbation method for general three-dimensional flows developed by Bayly [1], Friedlander and Vishik [17] and Lifschitz and Hameiri [37], which aims to establish the uniform boundedness in time of the amplitude of the perturbation. The shortwavelength perturbation method is a particularly elegant analytical approach which has proven to be highly applicable to a variety of recently derived exact solutions of the GFD governing equations, as can be seen in [7,16,19,24,25,[31][32][33][34]. Let U = Ue r + Ve θ + We z be a perturbation of the azimuthal flow whose velocity field is u = we z with w given by (3.13). We seek U, V, W and a pressure function P such that U + u, P + p satisfy (2.1), (2.2). Ignoring quadratic terms, we see that U and P satisfy subjected to an initial disturbance U| t=0 = U 0 . We employ the (WKB) ansatz, that is, we seek U and P solutions of (4.1) and (4.2) having the specific form where A(t, r, θ, z) = A 1 (t, r, θ, z)e r + A 2 (t, r, θ, z)e θ + A 3 (t, r, θ, z)e z , and f = f(t, r, θ, z) is a scalar function and plays the role of a small parameter.
Since u = v = 0 and w z = 0 we have that and Inserting the WKB ansatz (4.3) and (4.4) in (4.1), and taking into account (4.5) and (4.6), we obtain, after suitable identifications of the coefficients, and Since the vector A is not zero, we have from (4.8) that To solve for f in (4.9) we notice that the position vector of a particle is given by Therefore, u = ue r + ve θ + we z = d dt (r(t)e r + z(t)e z ) =ṙe r + rė r +że z =ṙe r + rθe θ +że z .
Thus, the equations of a streamline for the azimuthal flow with u = v = 0, w = w(r, θ) and passing through (r 0 , θ 0 , z 0 ) are whose general solution is Moreover, the equation (4.9) satisfied by the phase f has the general solution We see now immediately from (4.10) and (4.11) that f is constant along the streamlines of the azimuthal flow. Therefore, the system (4.7) satisfied by the amplitude A of the perturbation U simplifies to The previous system can be rewritten as (4.12) We remark now that Consequently, the system (4.12) becomes ⎛ where d is a constant depending on the initial data, and (4.13) Remark 4.1. The study of the stability of the azimuthal flow with w given by formula (3.13) hinges upon the possibility to find the eigenvalues of the matrix (4.13). More precisely, the previously named flow is stable if the eigenvalues of M are purely imaginary. The computation of those eigenvalues is, for a general flow, not possible. However, for certain choices of the density function ρ and of the function F involved in the formula for w the computation of the eigenvalues becomes feasible. We make the choice ρ(r, θ) = (b − ar) cosθ and F(x) = cx, where a, b, c are some constants with a, b > 0 and b − ar > 0 for all r that describe points in the flow. with r, θ given by (3.11) is linearly stable under short-wavelength perturbations.

Remark 4.4.
On the other hand, if r 0 , θ 0 are such that ar 0 − b sin 2 θ 0 < 0, then, as it is evident from (4.18), there are values of c for which E(r 0 , θ 0 ) < 0. Thus, for such a, b, c, the flow (4.14) is not stable under short wave-length perturbations.