On an exact solution of a nonlinear three-dimensional model in ocean flows with equatorial undercurrent and linear variation in density

The aim of the paper is to develop an exact solution relating to a system of model equations representing ocean flows with Equatorial Undercurrent and thermocline in the presence of linear variation of density with depth. The system of equations is generated from the Euler equations represented in a suitable rotating frame by following a careful asymptotic approach.The study in this paper is motivated by the recently developed Constantin-Johnson model [ 13 ] for Pacific flows with undercurrent and the exact results provided therein. The model formulated is two-layered, three-dimensional and nonlinear with a symmetric structure about the equator. The equations contain Coriolis effect and is consistent with \begin{document}$ \beta $\end{document} - plane approximation. Exact results of the asymptotic system of equations have been derived in a region close to the equator.

1. Introduction. Equatorial Undercurrent plays an important role in ocean dynamics and has a controlling influence on the temperature, salinity, and the nutrient enrichment of equatorial thermocline. It also impacts on biological productivity, atmospheric carbon exchange, and the El-Niño Southern Oscillation.
Since the observation of these undercurrents in 1950s by Cromwell, several investigations have been carried out [29]. Although significant progress has been made to explain features, primarily based on linear theory, which has been successful in modelling a number of observations, there are still some open questions yet to be resolved. In this context, suitable nonlinear theories are expected to provide more insight as the ocean dynamics is essentially three-dimensional and nonlinear [15].
The following three properties of the equatorial dynamics are: (1) A shallow westward flow of surface current (typically within the top 50m depth). This is confined mainly to a few degrees on either side of the equator.
(2) A strong coherent eastward undercurrent flowing at a speed much faster than the surface current. This is primarily responsible for the vertically-integrated transport. Beneath this undercurrent, a relatively weak flow exists.
(3) Flow in the westward direction on either side of the EUC, with eastward countercurrents poleward of this. This current is strongest in the Northern hemisphere and meets the surface.
The flow structure is nonlinear and complex with regions of upwelling and downwelling superimposed on the EUC. Hence, the flow is essentially three dimensional, though upwelling may be restricted to certain specific zones. While flow at the surface moves predominantly to the west driven by wind, and some flows moving to the north or to the south, further down little is known about flow structures [24,25,26]. A schematic shown in Fig.1  Equatorial ocean dynamics has some striking features. Close to the Equator the meridional component of the Coriolis force is negligible (vanishing at the Equator) leading to the breakdown of geostrophy. Hence, the flow in this region is essentially driven by the easterly trade winds in the Pacific and Atlantic, and by the seasonally reversing monsoonal winds in the Indian Ocean [35]. Equatorial ocean regions also present a pronounced stratification, greater than anywhere else in the ocean [16].
The simplest type of model of the EUC consists of a homogeneous fluid subject to uniform westward stress at the surface. This local model, in a rudimentary form, is unstratified. The model can be further developed to account for the effect of horizontal viscosity and more realistic models can also consider effects of stratification leading to a layered model of undercurrents. Other non-local and physical models involve a pressure head acting on the undercurrent which has an extra-equatorial origin. Observations indicate that some of the water in the EUC flows from subtropical gyres as the temperature in the core of the undercurrent ranges between 16 − 22 0 C. This range of temperature is much lower than the surface equatorial temperatures except at the eastern end of the ocean basins where the undercurrent ends. A two-layer inertial model of undercurrents can account for equatorial and extra-equatorial dynamics.
Local theories of undercurrent were proposed by [14]. These theories were developed further by [34,36,33]. Other linear non-dissipative models can be found in [19,20,31]. Significant extensions to linear models including effects of continuous stratification are due to [28]. In addition, some nonlinear models have also been developed [5,3,4,31]. Inertial theories began with the work of [17] and were further advanced by [32]. In the investigations by [30], the EUC was considered as a part of a larger and more complex subtropical current system, with both local and inertial effects. However, a complete model of the integrated dynamics requires numerical approaches.
The basic variation with depth of the main azimuthal flow of the EUC can be modelled in the spherical geometry of a rotating Earth [11], to capture the fine structure. Approximations are possible with certain limitations; for an overview of ocean dynamics in the equatorial Pacific in the f -plane setting see [10]. Explicit solutions to the nonlinear governing equations in the equatorial β-plane, were recently obtained in the Lagrangian framework [6,8,9,21,22]. These solutions represent realistic flow only for the region near the surface or in the neighbourhood of the thermocline.
Constantin and Johnson [13] recently proposed a model based on the classical fluid-dynamics approach. The nonlinear, three-dimensional model, herein called the Constantin-Johnson model, was developed following a systematic approach which is mathematically consistent, capturing the essential properties of the flow observed in the Pacific equatorial undercurrent (EUC), and avoids any oversimplification. The model, for example, accounts for the Coriolis effects due to the rotation of the Earth and retains non-traditional terms arising in the formulation. The model is based on the assumption of slow evolution of a two-layer flow in the equatorial direction and successfully provides asymptotic solutions. Some numerical studies on this model to investigate the properties of velocity field and flow paths have been carried out in [1,2].
Motivated by the Constantin-Johnson model and the exact results obtained by them, this paper considers the case of linear variation of density along the depth from the surface down to the bottom of the undercurrent which is close to but little below the thermocline, and a constant density below this, which represents a realistic scenario. Exact results have been derived in this paper for the case described with density variation, retaining the nonlinear, three-dimensional structure of the model. This also indicates that the Constantin-Johnson model is indeed versatile enough to cater for several realistic situations.

2.
Governing equations for an EUC. We choose a coordinate system that rotates with the Earth, with x pointing due east, y due north, and z vertically upwards (see Fig.2). This chosen coordinate system is the one associated with the tangent plane (and therefore consistent with the β-plane approximation). The β-plane approximation adopts a local flat Cartesian coordinate system, while still capturing the effects of the curvature of the Earths surface. These coordinates are obtained from the original spherical coordinate system by using a two-stage manoeuvre. First, the great circle (the equator) is straightened to become the generator of a cylinder, its circular cross section now representing the interior of the sphere. The coordinate perpendicular to this, on the surface, is then be expressed as the arc length along the circular rim of the cylinder, following the circumference at a mean radius of the Earth. The coordinate y in Fig.2 is in the direction of the tangent at a point on the equator. This is, however, not the curvilinear coordinate which we denote by y. In typical Cartesian coordinates, the z-axis points vertically upwards from the centre of the Earth, for all points on the surface; we denote the corresponding  Figure 2. The rotating frame of reference based on tangent plane, with the x axis chosen horizontally due east, the y axis horizontally due north and the z axis vertically upward. cylindrical (radial) coordinate by z. The two governing equations (Euler and mass conservation) with spatially varying density may therefore be written as where u = (u, v, w) is the velocity of the fluid corresponding to x = (x, y, z), t is the time (and D/Dt the associated material derivative), while p and ρ are the pressure and the spatially varying density, respectively. The body force is represented by F (which incorporates the effect of constant gravitational acceleration only) and Ω = Ω(0, cos θ, sin θ) is the angular-velocity vector describing the rotation of the Earth (with |Ω| = Ω ≈ 7.29 × 10 −5 rad s −1 and θ the angle of latitude). We now invoke a suitable approximation of the Coriolis term close to the equator. For small values of y/R, where R is an average radius of the Earth, we write the term 2Ω × u as 2Ω(w − yv/R, yu/R, −u), since sin θ ≈ θ and cos θ ≈ 1 for |θ| 1. This is in effect the equatorial β-plane approximation, which is regarded as acceptable within about 2 0 of the equator [18]. In this context, more recent discussions on equatorial f -plane are available in [7,27,12,23]. We therefore have x measured, equivalently, along the curvature of the equator (eastwards) and y measured along a line of longitude (northwards). It may be noted that, with this same approximation, the curved surface of the Earth drops below the (x, y)-tangent-plane by an amount of y 2 /2R, approximately. This is indeed, consistent with, and a consequence of, the β-plane approximation. We thus have z = z + 1 2 y 2 /R and y = y (see Fig. 2). With this interpretation of the chosen coordinate system, the distance from the centre of the Earth is R + z, so that the gravitational potential is ρzg and thus the gravitational field (its negative gradient) becomes (0, 0, −ρg). The gravity term appears only in the z-component equation. The governing equations for a steady flow therefore become (writing y for y) with g as the constant acceleration due to gravity. Let us consider the case of a linear variation of density from the surface down to the bottom of the equatorial undercurrent (EUC) and then constant density below this. If so desired the problem can be reformulated by considering linear variation in density from the surface down to a characteristic depth corresponding to the position of the theromocline. The depth from the surface to the bottom of the EUC is denoted by h. Let the density at the surface and at the point at the bottom of the EUC be denoted by ρ s and ρ 0 , respectively. Since the warmer water with lower density is at the top of the dense, colder layer, ρ s ≤ ρ 0 .
Assume the density variation of the form is the depth of the mean ocean surface and we denote α = ρ s /ρ 0 .
Eqs. (2)-(5) have been expressed within the limitations of the β-plane approximation and, precisely as written, applied to the region above the bottom of the undercurrent. For the region below, the density is assumed to be constant, the equations similar to those given in [13] then follow directly.
The boundary conditions that we require involve the pressure at the surface (z = η(x, y)) and a kinematic condition that describes the motion of the free surface. These are, respectively, given by the following where P s (x, y) is the pressure prescribed at the surface. Further discussions on the pressure will follow later. There is a corresponding kinematic condition at the thermocline (z = T (x, y)), which is given by However, since there is no pycnocline, the pressure is automatically continuous across the thermocline. Further, the kinematic condition at the impermeable ocean bed z = −d(x, y) is given by 3. Non-dimensionalization of the equations. The next stage in our analysis is to recast the governing equations in a suitable form by introducing nondimensionalization of the variables so that we can carry out an asymptotic analysis.
To this end, we use the following non-dimensionalization of the variables x = (x, y, z) = (Lx, ly, hz), where the length scales are (L, l, h) and U is an appropriate speed scale. The velocity components have been non-dimensionalised in a form that gaurantees the existence of a stream function, via a balance of terms of equal size, irrespective of whether the flow is in the (x, z) or (y, z) frame. Normalizing with respect to ρ 0 , (2) -(5) become Rḡ yu = −p y ; and with the boundary conditions p = P s (x, y)/(ρ 0 U 2 ) = P 0 (x, y); w = uη x + vη y on z = η(x, y), and w = uT x + vT y on z = T (x, y), where η = hη and T = hT . Here, we do not use the bottom boundary condition since the flow below a certain depth is assumed stagnant and this condition is automatically satisfied by the flow field. It may be noted that if α = 1, we arrive at exactly the same equations obtained for piece-wise constant density as used in [13]. We further redefine the pressure so that the pressure is expressed relative to the ambient hydrostatic pressure. This gives Further, we setl = √Rh ,Ωh/ū = ω. With this choice for l, we invoke the limiting process described by (h/L) 2 → 0 with (l/L) 2 = (l/h) 2 (h/L) 2 → 0 and all other parameters held fixed. For some typical values of length and velocity scales, e.g. for the case of Pacific ocean, see [13]. The governing equations for momentum and mass conservation therefore take the form with the boundary conditions P (x, y, z) = P 0 (x, y) + κη(x, y) w = uη x + vη y on z = η(x, y), where κ(η) = [α−(1−α) η 2 ]gh/U 2 . The governing equations that correspond to (16)- (19) and (20)- (21), but valid below the region with the undercurrent corresponding to constant density, are written as (see [13]) The interpretation of our choice of scales is such that, in the flow direction, we have 0 ≤ x ≤ 1. In the context of the EUC, we could associate the western end, away from any shorelines, by x = 0 and correspondingly x = 1 at the eastern end. The width must be restricted (by virtue of the β-plane approximation), which we represent by −y 0 < y < y 0 , for suitable (finite) y 0 . 4. Pressure function. We now examine the form taken by the pressure in this asymptotic formulation.
Let us introduce leading to P y = −2ωφ ζ y ; It may be noted that ζ = z h is the non-dimensional coordinate measured from the tangent plane. We then have where A(x) is an arbitrary function. Once the velocity field u is prescribed, one can then determine φ and henceφ. For convenience, we define φ(x, 0, 1 2 y 2 ) = 0 which impliesφ(x, 0, 1 2 y 2 ) = 0 if u = 0 at z = 1 2 y 2 . We now investigate the pressure boundary condition in more detail. To do so, we first describe the free surface z = η(x, y) relative to the tangent plane by setting η(x, y) = 1 2 leading to the pressure at the free surface given by h(x, y), η(x, y)) (26) which is the relation between the pressure prescribed at the surface and the surface distortion.

5.
Simplifications of the model based on physical considerations. Since the free surface incorporates the curvature in the y direction as measured with respect to the tangent plane, no further dependence on y is assumed in this model, i.e. h = h(x) only. Further, we also prescribe that the free surface is flat in the x direction which is not an unreasonable simplification as the ocean rises from east to west (relative to the curvature of the Earth) by a small amount (e.g. by 0.5 m for the case of the Pacific ocean). Hence, we have h(x) = 0. The free surface is therefore the (local) undisturbed surface of the sphere within the β-plane approximation. Pressure at the surface is thus given by We now invoke the condition of constant pressure along y = 0 leading to A(x) = constant. The thermocline is identified with a plane parallel to the tangent plane (and so follows the curvature in the y-direction), but we allow it to evolve in the x-direction. Thus, we choose to set Thus the special problem that we discuss here is the one which follows the geometry of the curved surface (in the y-direction, locally) and has a (local) flat free surface associated with a constant atmospheric pressure at the surface along the equator and a thermocline at a level which may evolve in the equatorial direction. This is described by the set of equations with w = yv on z = 1 2 y 2 and w = −ut + yv on z = −t(x) + 1 2 y 2 (31) where the prime denotes the derivative with respect to x, and u is assumed given as a function of z on y = 0, and therefore it is known across the (z, y) plane (with the x-evolution superimposed on this).

Method of solution.
We now proceed to solve the problem using some similar concepts used in [13]. At the non-dimensional longitude x, let us denote the depth at which u first vanishes below the surface by ζ 0 (x) < 0. Note that field observations indicate that u(x, 0) < 0 all along the equator at the surface, while u vanishes in the abyssal layer where there is no motion.
Since φ(x, 0, 1 2 y 2 ) = 0, we conclude that now using (19) in (16) we get which can be re-expressed as where andφ Eq. (33) gives an equation of similar form as in [13], but with a different expression on the right hand side. To solve (33) we follow a similar approach to the one proposed in [13] by introducing a quasi-stream-function and setting v(x, y, z) which on integration with respect to the z variable leads to wherê Expressions for v and w are obtained as We can write (33) as and then from (39) we obtain ) and B(x, ζ) = −(1 − α)u 2 /E(x, ζ) , then from (40) we have where with It may be noted that when α = 1, we recover precisely the same equation as obtained in [13] for constant density.
We now proceed to find the characteristic curve for the solutions of (41). From (41), we have the equations for the characteristics and the compatibility condition as dz .
The trajectories of the characteristic curve are given by the differential equation Let us now consider a function H(x, ζ, z) ∈ C 1 (D) such that By appropriately choosing the vertical length scale and restricting the region of flow of interest to a certain width on either side of the equator, we can limit y to a suitable value y L , such that |y| ≤ y L < 1; and D = {(x, ζ) : 0 < x < 1 , ζ 0 < ζ = z − 1 2 y 2 < 0 with |y| < y L }. Thus, y 2 /2 1 and can be assumed to be close to zero within this region for the purpose of evaluating the variation of H(x, ζ, z) with respect to y in this particular region. Hence, in this domain we may assume leading to H y (x, ζ, z) = 0 .
In fact, numerical investigations [1] have shown that a number of interesting features are concentrated within a range of |y| ≤ 0.5, near the equatorial region, which corresponds to a total width of about 35 km (i.e., ∼ 20 km on either side of the equator). If the vertical length scale to be used for normalization is chosen to be, e.g. 2h, y would be restricted to values less than 0.35, if we keep the region of interest unchanged (i.e., to a total width of about 35 km); or if we so desire, we can cover a wider region, e.g. with a total width of about 60 km across the equator by allowing |y| to be within a value of about 0.6.
In this region close to the equator, the characteristics curves take the form where C is an arbitrary constant.
Notice that in [13], B(x, ζ) = 0, which when used in (45) leads to the characteristic curves with K an arbitrary constant, alongwhich the solutions propagate in [13].
Noting the form of the characteristics of the solution of (41) we define a function Let us seek a solution for Ψ in (41) of the form Ψ(x, y, z) = yψ(x, ξ(x, y, z)) .
This is consistent with the fact that the velocity v along y = 0 should be zero (for symmetry). On using the expressions for H z and H y from (46) and (48) respectively, we then have The expressions in (53) when substituted in (41) gives leading to Using the expression for Ψ; v and w follows. We note that there is a singular term in this solution, (55), arising from the zero of (F − B), a problem that we address in the next section. The expression for v and w are obtained as and respectively, where the subscripts in (56) and (57) denote partial derivatives. In the region below the bottom of the undercurrent, the solution given in [13] is valid. With the hypotheses that the singular term F − B = 0 is avoided, a solution of the velocity field in the fluid region {D = (x, ζ) : 0 < x < 1 , ζ 0 < ζ = z − 1 2 y 2 < 0 with |y| < y L } , is obtained. The solution is also valid for the fluid region ζ < ζ 0 as it satisfies (22) and (23). 7. Singular term in the solution. For the solution in the previous section to exist, we need F − B = 0. Consider the following relation which leads to since E = 0. Since the model is developed based on asymptotic formulation, and the choice of suitable scales and velocity (see [13]) leads to a non-dimensional value of ω ≈ 0.03, the consideration of ω → 0 is relevant and appropriate. When ω → 0, we need to avoid which on integration yields where g(x) is an arbitrary function. Hence, we need to avoid the profile given by (61). However, this is not very restrictive as we can choose various different profiles not of the form in (61), e.g., see [13] for some typical profiles. Notice that when α → 1, the profile in (61) has a singularity which results in an unbounded and unrealistic profile.
8. Remarks. We notice that for the proposed solution if u ∈ C 1 (D), then (v, w, P ) ∈ C 1 (D) where D is the closure of the fluid domain. If convergence in the sense of distributions is considered then it is convenient to work with the Sobolev space.
In particular, if we consider the Hilbert space with u ∈ H 1 (D), we can relax the regularity requirement and (v, w, P ) will belong to the Hilbert space H 1 (D) and will have to satisfy the kinematic and dynamic boundary conditions as previously discussed.
9. Discussions. This paper presents an extension of the original work by Constantin and Johnson, which was for piece-wise constant density, to a linear variation (in depth) of the density. The work carried out in this paper has shown that the Constantin-Johnson model is versatile enough and can be readily developed to include some more realistic flow properties. It may be noticed from (45) that the main contribution to the change in the flow structure due to linear variation (in depth) of density as compared to the one when the density is assumed to be constant, comes from the non-zero B(x, ζ) term. Since the variation in density is small (less than 1%) over the warm shallow layer, the term (1−α) → 0 and in turn B(x, ζ) is close to zero, in general. However, in the presence of strong undercurrents, i.e. large values of u, the term B(x, ζ) may have non-zero contribution as B(x, ζ) is proportional to u 2 ; and the flow structure may be influenced by the linear variation of density in such a situation. It would be interesting to carry out some numerical investigations, both for the case of solutions close to the equator (for which exact solutions have been developed in this paper) and solutions over a wider range on either side of the equator. In addition, it would also be interesting to observe how the streamlines and flow patterns are for the case with linear variation (in depth) of the density, and compare those with the cellular structures and their evolution along the equator as has been obeserved in previous studies [13,1]. Simulations of flow paths to investigate the presence of features such as upwelling, downwelling and extra-equatorial dynamics as observed previously [2] would also be of interest in future studies.