Free surface equatorial flows in spherical coordinates with discontinuous stratification depending on depth and latitude

We derive and subsequently analyze an exact solution of the geophysical fluid dynamics equations which describes equatorial flows (in spherical coordinates) and has a discontinuous fluid stratification that varies with both depth and latitude. More precisely, this solution represents a steady, purely–azimuthal equatorial two-layer flow with an associated free-surface and a discontinuous distribution of the density which gives rise to an interface separating the two fluid regions. While the velocity field and the pressure are given by means of explicit formulas, the shape of the free surface and of the interface are given in implicit form: indeed we demonstrate that there is a well-defined relationship between the imposed pressure at the free-surface and the resulting distortion of the surface’s shape. Moreover, imposing the continuity of the pressure along the interface generates an equation that describes (implicitly) the shape of the interface. We also provide a regularity result for the interface defining function under certain assumptions on the velocity field.


Introduction
Flow stratification represents a pronounced feature of geophysical ocean flows that particularly applies to large scale ocean movements [17,29]. To exemplify the previous aspect we point out the existence of a band of 150 km on each side of the Equator, stretching out longitudinally over about 16.000 km throughout the extent of the Pacific Ocean, where the stratification gives rise to a sharp interface (called pycnocline of thermocline) [16]. The thermocline separates a shallow near-surface layer of relatively warm water from a deep layer of colder and denser water.
Another important aspect that one needs to consider in the investigation of geophysical ocean flows refers to the effects generated by the presence of Coriolis forces-created by the Earth's rotation. Indeed, together with the prevailing westward wind pattern, Coriolis forces generate an underlying current field that presents flow-reversal [2,12]: in a band of about 2 • latitude around the Pacific Equator the current field shifts from a westward flow near the surface to an eastward flowing jet called the Equatorial Undercurrent (EUC) whose core resides approximately on the thermocline.
The interaction of currents and (surface and internal) waves which also takes into account the density stratification is a recent avenue of research [5-7, 11, 27] performed within the setting of nonlinear geophysical governing equations. We find important to mention that the issue of stratification was considered from a rigorous analytical perspective only (relatively) recently: a selective list of works dealing with two-dimensional gravity water flows (without Coriolis effects) consists of [13,14,18,20,[35][36][37].
Rigorous mathematical analyses of the aspects mentioned earlier were started by Constantin and Constantin & Johnson who constructed exact solutions exhibiting various features shared by geophysical ocean flows, like three-dimensionality [4,10,26], equatorially-trapped waves [3], presence of underlying currents [19], or a preferred propagation direction, cf. [8,9], where an approach that takes into account the spherical shape of the Earth was initiated. The latter approach (by means of spherical coordinates) was extended by Henry & Martin [21][22][23][24] to construct exact azimuthal equatorial flows allowing for continuous stratification (depending on depth and latitude) and by Martin & Quirchmair [30] to devise exact continuously stratified solutions concerning the Antarctic Circumpolar Current (ACC). Moreover, exact azimuthal solutions with a discontinuous density stratification, giving rise to internal waves and pertaining to EUC and ACC were presented recently by Martin [31] and Martin & Quirchmair [32], respectively.
In this paper we broaden the choice for the density function from [31] to the case of a discontinuous density that varies with both depth and latitude. The relevance of this choice stems from the existence of a zonal tilting of the thermocline, as noted in Constantin & Johnson [10]. With these aspects in mind we proceed to construct a family of exact solutions, presented in a rotating framework (by means of spherical coordinates), that represent incompressible azimuthal equatorial water flows with a free surface and a free interface. We also want to emphasize that the solutions we provide possess a depth-dependent velocity field. For studies concerning stratification, as well as the presence of an interface, in the equatorial scenario we refer the reader to [5,11,15,27,31,34]. Concerning recent progress toward the understanding of several intricate features underlying the dynamics of coupled surface and internal waves we point out the excellent work by Henry & Villari [25].
The layout of this paper is as follows: after introducing the physical problem in Sect. 2 we derive in Sect. 3.1 formulas for the velocity field and the pressure in the two domains separated by the interface arising as a result of jump in the density. An implicit formula for the interface (which plays the role of an internal wave) is obtained in Sect. 3.2 by utilizing the balance of forces along it. Moreover, the dynamic boundary condition allows us to find a relation between the enforced pressure on the surface and the deviation of the surface from a surface following the Earth's curvature. While the two implicit relations defining the free surface and the interface are quite involved they can be studied by means of the implicit function theorem. We conclude the paper with Sect. 4 where reasonable and expected physical properties of the exact solutions are presented. Indeed, it is proved that a growth in the pressure along the The variable r represents the distance from the origin, θ is the polar angle (π/2 − θ being the angle of latitude), and ϕ is the angle of longitude (azimuthal angle) surface causes a decrease in the free surface height. Moreover, a regularity result concerning the interface defining function is also presented.

Preliminaries
This section is concerned with the statement of the governing equations (together with their boundary conditions) for geophysical fluid dynamics (GFD). These are formulated in a spherical coordinate system which is fixed at a point on the Earth's surface as will be detailed shortly. As a result the fine features of the Earth's spherical geometry are reflected in the structure of our solutions. The latter are assumed to have a jet-like structure, capturing the observed azimuthal character of some equatorial flows. In the formulation of the governing equations we take into account the observation by Maslowe [33] that the Reynolds number is, in general, extremely large for the type of flows described above.
The right-handed coordinate system we will be using is denoted with (r , θ, ϕ), where r is the distance from the centre of the earth, θ (with 0 ≤ θ ≤ π) is the polar angle, and ϕ (with 0 ≤ ϕ < 2π) is the azimuthal angle, see Fig. 1. The location of the North and South poles are at θ = 0, π, respectively, while the Equator is situated at θ = π 2 . The unit vectors in this (r , θ, ϕ) system are (e r , e θ , e ϕ ), respectively, and the corresponding velocity components are (w, v, u); e ϕ points from West to East, and e θ points from North to South. More precisely, e r = (sin θ cos ϕ, sin θ sin ϕ, cos θ), e θ = (cos θ cos ϕ, cos θ sin ϕ, − sin θ), A noteworthy aspect is that the fluid domain is stratified where the stratification is brought about by the changes in the (discontinuous) density. To be more specific, denoting with R ≈ 6378 km the Earth's radius and with r 0 > 0, r 1 > 0 some constants, we assume that the flow consists of a upper layer of density ρ 1 (r , θ) that lies above a bottom layer of density ρ 2 (r , θ), with R 0 := R + r 0 and R 1 := R + r 1 . While above θ → d(θ ) is a given function, h(θ, ϕ) and k(θ, ϕ) are unknowns of the problem and denote the interface and the free surface defining functions, respectively.

Remark 2.1
Based on data collected on physical grounds, it is appropriate to assume the relation Denoting by u = ue r + ve θ + we ϕ the velocity field we have that the governing equations in the rotating (r , θ, ϕ) coordinate system are the Euler's equations, where p(r , θ, ϕ) is denotes the pressure in the fluid and (F r , F θ , F ϕ ) is the body-force vector, and the equation of mass conservation The description of the water wave problem is completed by the specification of the associated boundary conditions. These are as follows. At the free-surface r = R 0 + k(θ, ϕ) we require the dynamic condition involving the surface pressure (for some given function P 1 (θ, ϕ)) and the kinematic condition to hold. At the interface, r = R 1 + h(θ, ϕ), we ask that the normal components of the velocity fields from the upper and lower layer, respectively, are the same. The latter is equivalent with the condition Moreover, to ensure the balance of forces at the interface, we also require that the pressure from the upper layer coincides with the pressure from the bottom layer, that is, At the bottom of the ocean, which is an impermeable, solid boundary described by the equation r = d(θ, ϕ), the associated kinematic condition is

Explicit solutions
We aim at finding solutions of the Eqs. (1a), (1b) and (2) that represent purely-azimuthal steady flows with no variation in the azimuthal direction. That is to say that the velocity field satisfies θ). Moreover, the other unknows of the problem are characterized by the properties Throughout the section the range for the θ variable will be the interval π 2 , π 2 + ε where ε = 0.016 determines a strip of about 100 km about the Equator.

The velocity and the pressure
We notice that a flow with the previous features automatically satisfies the boundary and interface conditions (2) as well as the equation of mass conservation (1b). Assuming that the only body-force is due to gravity alone, that is the body-force vector is −ge r , the Euler equations are written as in the domain D i for i = 1, 2. In the first step we simplify the system above by eliminating the pressue. So we obtain that the azimuthal components of the velocity field w i satisfy for i = 1, 2. Utilizing the method of characteristics and adapting the approach from [23] we obtain that the azimuthal velocities w i (i = 1, 2) are expressed in terms of the formulas Toward deriving the formula for the pressure we first infer from (3) that the pressure gradient is given as and for i = 1, 2. An integration with respect to r yields that for all r ∈ [R + d(θ ), where and θ → C 2 (θ ) is a function such that We proceed now with the determination of the pressure in the upper layer D 1 . We recall that the interface (which is the lower boundary of D 1 ) is given as r = R 1 + h(θ ). An integration with respect to r in formula (7) yields that for all θ and all r ∈ [R 1 + h(θ ), R 0 + k(θ )] it holds where and

The interface and the free surface
Having the velocity and the pressure determined we pass now to the (implicit) determination of the two interface defining functions h and k. To this end we first exploit the balance of forces at the interface (2d), written now (in the context of azimuthal flows) as which can be detailed as We provide now a functional analytic setting for our problem by nondimensionalizing the interface defining function. That is, we set and so write (14) as where the operator G 2 acts from the Banach space C 1 π 2 , π 2 + ε into itself and is given as To recover the shape of the free surface (at least in an implicite form) we use the dynamic boundary condition (2a) and so obtain that a prescribed pressure at the surface P 1 (θ ) has to satisfy the equation called the Bernoulli relation which establishes a connection between the pressure applied on the free surface on one hand and, on the othe hand, the shape of the free surface and the shape of the interface, where the latter is determined by Eq. (15). Setting we recast the previous equation in nondimensional form as the abstract operatorial equation where G 1 is an operator from the Banach space C π 2 , π 2 + ε × C 1 π 2 , π 2 + ε × C π 2 , π 2 + ε into itself and is given through The problem of finding (k, h) that satisfy (15) and (17) can be written now as In order to prove the existence of nontrivial solutions h, k to (20) we are going to utilize the Implicit Function Theorem, cf. [1]. In order to do so, we need first to identify trivial solutions (k 0 , h 0 ) of (20). The natural candidate for a trivial solution is the flow having an undisturbed free surface and an undisturbed interface following the Earth's curvature. That is, we set k = h = 0 in (20) and find that G 1 (0, 0, P 0 1 ) = G 2 (0) = 0 if and only if and 1 P atm We evaluate To begin with we compute lim s→0 1 s θ π/2 Moreover, and lim s→0 θ π/2 Collecting now (24), (25) and (26) and recalling formula (13) we have Moreover (1+k(θ ))R 0 sin θ (1+sh(θ ))R 1 sin θ Appealing to the mean value theorem for integrals we obtain lim s→0 1 s and lim s→0 1 s Therefore, taking into account formulas (27), (29), (30) as well as the definition of the operator G 1 in (19) we obtain Similarly, we compute Performing calculations analogous to (29) we infer that If we are to consider G 2 as an operator of k, h and P then clearly, G 2k (0, 0, P 0 1 )k = 0 for all k. Therefore, as a linear operator from C [π/2, π/2 + ε] × C 1 [π/2, π/2 + ε] into itself.
Proof To begin with we note that, due to Remark 2.1, we can write On the other hand From (33), (35) and (36) we conclude that there is a negative constant a such that is a linear homeomorphism. On the other hand, since the typical values for velocity in the ocean do not exceed 4 m·s −1 we see that w(R 0 , θ) + R 0 sin θ) 2 is much smaller than g R 0 . Hence, we can conclude that there exists b < 0 such that (G 1k (0, 0, P 0 is a linear homeomorphism. In view of (34), (37) and (38) we can now conclude that (G 1 , G 2 ) k,h (0, 0, P 0 1 ) is homeomorphism from C [π/2, π/2 + ε] × C 1 [π/2, π/2 + ε] into itself. Resorting now to the implicit function theorem [1] delivers the assertion made in the statement of the theorem.

Properties of the exact solutions
This section is devoted to the analysis of the exact solutions found in Sect. 3.2. The first result indicates that our solutions in spherical coordinates present expected physical properties: while a growth in pressure on the surface leads to a decrease in the free surface height, the increase in the latter brings about a attenuation of the former.

Remark 4.1
For the next result we will assume that the given pressure P 1 on the free surface is a differentiable function of θ . Then an iterative bootstrapping procedure, cf. [1], ensures the differentiability of the surface defining function k. Theorem 4.2 Denoting, as before, with k the function defining the free surface and with P 1 the pressure on the latter, we have that and for all θ ∈ π 2 , π 2 + ε .
The last part of the proof consists in showing that A(θ ) < 0 for all θ ∈ π 2 , π 2 + . To this end, we note that assumption (42) yields The expression in the bracket above is negative as we can deduce by taking into account the size of the physical quantities involved. Hence A(θ ) is negative for all θ ∈ π 2 , π 2 + . An iterative argument shows now the asserted differentiability of h.