1. Introduction

Geophysical fluid dynamics (the study of fluid motion in which the effect of the Earth’s rotation plays a significant role) in the equatorial region is a highly complex subject of vast importance which has experienced a number of interesting developments in recent times from a mathematical point of view, cf. [21, 22, 25] for surveys of some of these recent developments. Of particular note has been the successful application of a variety of classical applied mathematical methodologies and techniques in deriving new approaches to the modelling of oceanographic processes [7, 8, 10, 25]. This strategy, which is based upon the assumption that the observed movement of the oceans is an intrinsic property of a fluid which can be captured by the Euler equations (rendered more tractable by way of relevant asymptotic approximations), is quite opposed to the frequently-adopted approach which consists of working with purely ad-hoc oceanographic models to capture given physical processes [9]. However, the complexity inherent in this approach is apparent since even in the setting of an inviscid, incompressible and laminar fluid, the GFD equations of motion are nonlinear and highly intractable, with Coriolis forces incorporated in the Euler equation.

In the equatorial region, whereby latitudinal variation is necessarily restricted, the governing equations are typically simplified by invoking tangent plane approximations, the classical forms being the f-plane, and the β-plane, approximations. The equatorial f-plane approximation replaces the exact spherical relations for the Coriolis parameters 2Ω(cosθ, sinθ) near the reference latitude θ = 0 by 2Ω(1, 0), while in the β-plane approximation a truncated Taylor expansion in the variable y = θ/R is used. Here R is the planetary radius, so that 2Ωsin(y/R) ≈ βy while cos(y/R) ≈ 1, where β = 2Ω/R; see [12, 15, 27]. The use of the f-plane gives insight into the equatorial ocean dynamics (see the papers [17, 18]) and the discussion in the paper [3] provides an estimate of the meridional extent of the validity of the f-plane approximation near the Equator. Nevertheless, the β-plane approximation is more accurate.

The “traditional β-plane approximation”, ubiquitous in modelling equatorial flows [12, 15, 27], neglects terms in the Coriolis force which feature the vertical velocity, and also the vertical component of the Coriolis force. The overwhelming factor behind its widespread utilisation is the convenience of this approximation, and while the traditional approximation may often be applicable in GFD, this is not always the case. Its relevance or otherwise is a highly intricate, and as yet largely not understood, matter [14]. Within the Lagrangian framework, the traditional β-plane approximation proved useful in accommodating solutions for equatorially trapped surface and internal waves (see [2,4–6,16]) which, however, do not capture the existence of strong non-uniform underlying currents. Studies have revealed that the neglected “non-traditional” components of the Coriolis force have significant effects on the equatorial ocean dynamics, something which is evident in the recent papers [8, 19, 20], and particularly in [10], whose approach motivates the subject of this paper.

The dynamics of the ocean near the equator presents some unique and complex characteristics from a modelling perspective [7, 8, 10, 13, 23, 24], among these being: a pronounced stratification with a pycnocline/thermocline interface demarcating the corresponding fluid regions; a breakdown in mid-latitude geostrophic balance; additionally the equator acts as a natural waveguide leading to equatorially-trapped zonally propagating waves. Furthermore, the Equatorial flow and its dynamics are dominated by the presence of non-uniform underlying currents, with the most prominent one being the equatorial undercurrent (EUC). Superimposed on this flow are regions of upwelling and downwelling which effectively convert a one-dimensional flow into a three-dimensional one. For the ocean in a neighbourhood of the Pacific equator, the three-dimensionality results in a flow that rises to the surface, all along the equator, and moves away from the equator close to the surface (although this upwelling may be restricted only to regions quite near the surface).

The complexity of motion exhibited by equatorial flows is remarkable, and tends to suggest that any attempt to capture the salient features of the flow by way of a tractable mathematical model would be expected to encounter significant (if not insurmountable) difficulties. Hence, it is quite striking that in [10] the authors manage to derive from the Euler equation a nonlinear and three-dimensional system of model equations which can be solved exactly. This is achieved by applying an asymptotic procedure which merely assumes slow variations in the azimuthal direction in a steady two-layer flow which is symmetric about the equator. It is shown in [10] by way of specific examples that this dynamical model is sufficiently complex to simultaneously capture a variety of oceanic phenomena such as upwelling/downwelling, zonal depth-dependent currents with flow reversal, and poleward divergence along the equator.

While the derivation of the dynamical model presented in [10] is elegant, owing to the level of complexity inherent in the GFD governing equations it is by necessity quite convoluted and involved. The three-dimensionality of the dynamical model is achieved by employing a quite ingenious, “non-traditional” form of the governing equations which retains aspects of the earth’s curvature in the β-plane approximation, to which a mixture of asymptotic techniques are applied. The resulting three-dimensionality of the solution, and it ability to capture salient features such as upwelling / downwelling, hinges on this novel formulation of the governing equations.

The aim of this note is to demonstrate that the retention of nonlinear effects is fundamental to the success of achieving such a rich description of the oceanic flow processes in any mathematical model. In particular, it is shown that any linearisation of the model will fail to capture the emergence and persistence of large coherent structures that are representative of the oceanic equatorial flows described above. This result suggests that the model presented in [10] is more-than-likely as simple a model as can possibly be described which achieves such a rich level of dynamical behaviour simulating the equatorial flow, and furthermore recommends the approach in [10] as a starting point for future studies on this issue.

2. A non-traditional β-plane formulation

In this section we give a brief overview of the derivation of the approximate model considered in [10] which motivates this study, and which is particularly novel in its retention of the effect of the curvature of the earth’s surface in a “flat” Cartesian coordinate system. The Euler equation describing the fluid motion in GFD is expressed as

where u¯=(u¯,v¯,w¯) is the fluid velocity corresponding to the coordinate variables x¯=(x¯,y¯¯,z¯¯), Ω¯=Ω¯(0,cosθ,sinθ) is the angular velocity vector of the earth’s rotation, where θ is the angle of latitude (Ω¯=73×10−6rad/s the constant rotational speed), F¯ is the external body force, ρ¯ is the water density, and p¯ is the pressure. It is assumed that the earth is a perfect sphere of radius R¯=6378km. The (x¯,y¯¯,z¯¯)-variables are taken to represent cylindrical coordinates in a coordinate system which rotates with the earth, as follows. The x¯-axis, which generates the cylinder, corresponds to the “flattened–out” great circle of the equator (the interior of the cylinder corresponding to the interior of the sphere). The positive x¯-direction denotes azimuthal flow from West to East. In terms of cyclindrical coordinates, y¯¯ is the curvilinear coordinate following the circumference of the cylinder expressed as the arc length, and z¯¯ is the radial coordinate pointing vertically upwards from the earth’s surface.

For the β-plane approximation, whereby the earth’s curved surface is approximated by a tangent plane, the Coriolis force terms in (2.1) are linearised by way of the approximations sinθ ≈ θ, cosθ ≈ 1 (since the latitude θ is small), leading to

An associated Cartesian coordinate system {x¯,y¯,z¯} is now introduced whose origin is located at a point on the earth’s surface at the equator. The x¯-axis is pointing horizontally due east (zonal direction), the y¯-axis is due north (meridional direction), and the z¯-axis is pointing vertically upwards. For the tangent plane approximation it can be shown that y¯¯=y¯, and z¯¯=z¯+y¯22R¯, representing the drop in the earth’s surface beneath the tangent place cf. [10]; note that over a distance of 10 km (commensurate with the 0.25 degree grid size for the ocean data which are typically provided by the NOAA) we have a departure of about 7.5m from local spherical geometry on the flat-space setting inherent to the f-plane or β-plane approximation (see the discussion in the paper [11]). Hence the full governing equations (2.1) for the steady flow reduce to the following “non-traditional” β-plane formulation,

and the equation of incompressibility

Here g¯=9.8ms−2 is the constant acceleration of gravity, and ρ0¯ is the constant density of the upper fluid layer: identical equations of motion hold in the layer beneath the thermocline except with ρ0¯ replaced by (1+r)ρ0¯. Typically ρ0¯=1027 km m−3 and r = 5 × 10⁻³ for the Pacific EUC. We note that it was recently proven that while nonlinear, (weakly) three-dimensional travelling wave solutions exist for the traditional β-plane approximation [2, 4, 5], in certain cases “non-traditional” components of the β-plane formulation are vital for admitting the existence of solutions with more complex dynamics, cf. [19, 20]. These solutions are Gerstner-like in the sense that they are exact and explicit in terms of Lagrangian variables [1,21]. The non-traditional element of the formulation (2.3) is the presence of the z¯¯, as opposed to the z¯, term.

Given the equations of motion (2.3), it is shown in [10] that invoking a non-dimensionalisation of the variables of the form

where (L¯,l¯,h¯) are length scales and Ū is an appropriate speed scale, and taking the limiting procedure (h¯/L¯)2→0, (l¯/L¯)2=(l¯/h¯)2(h¯/L¯)2→0, leads to the reduced system of equations

Here ω=Ω¯h¯/U¯, and P(x,y,z)=p(x,y,z)+g¯h¯U¯2z is the pressure relative to the hydrostatic pressure. There are various other conditions derived relating to interfaces at the free-surface and thermocline, but in order to present a picture of the underlying flow we merely note that the authors solve the system (2.4) explicitly as follows. Assuming that the azimuthal velocity component u has the structural form u(x, y, z) = u(x, z − y²/2), but otherwise u may be prescribed remarkably arbitrarily, the remaining velocity components can be expressed exactly as

Here ζ = z − y²/2 and ϕ_ζ (x, ζ) = u(x, z − y²/2). Taking, as examples, a number of quite elementary flow profiles for u which model the depth-dependant behaviour of the EUC the authors demonstrate that the resulting flow field (2.5) for the reduced equations (2.4) exhibits a structure which simulates upwelling, downwelling, and various other intriguing flow configurations.

3. Linearised three-dimensional flows

In this section we show that it is not possible to describe the complex flow characteristics outlined above by way of a linearised system. In particular, modelling the upwelling process — whereby in the neighbourhood of the Pacific equator the oceanic flow rises to the surface, all along the equator and, close to the surface, moves away from the equator poleward — is not achievable in the linear regime. If the fluid motion is three-dimensional and inherently linear then its velocity field has a representation of the form u = ẋ = Ax, cf. [1] for a detailed exposition of aspects of linear fluid motion. Note here that, for convenience of notation, the variables x = (x, y, z) represent any appropriate coordinate system and do not necessarily correspond to the notation of the previous section. Accordingly, the three-dimensional flow can be expressed component-wise as

and the system (3.1) has the fundamental solution given by x(t) = x₀e^At where x₀ is the initial position of the fluid parcel [26]. It is a characteristic of equatorial flows that they are typically symmetric with respect to the equator: since the y-axis is oriented in the meridional direction, with y = 0 denoting the location of the equator, we infer that and upon differentiation with respect to t, and setting t = 0, we conclude that e₀ = (0, 1, 0) is an eigenvector of A with eigenvalue 0. Conservation of mass, which is prescribed by equation (2.3d), implies that trA = 0. Consequently, the characteristic equation for A must be of the form where q must be real as the coefficients of A are real. We therefore must consider three different scenarios: q = 0, q < 0 and q > 0.

If q = 0 then the Jordan canonical form of A is either 0, or else it takes the form

in which case the fundamental matrix [26] is given by

Hence the solution can be expressed as

Consequently x″(t) = (c, 0, 0), implying a uniform zonal acceleration for each fluid particle. This evidently cannot occur physically, must less correspond with the equatorial upwelling process.

If instead we suppose that q < 0 then the eigenvalues must be

where p² = −q for p > 0. This scenario implies the existence of both a stable and unstable manifold, which cannot model the steady flow exhibited in [10] which, in the near-surface layer, moves simultaneously poleward and towards the west along the orthogonal (horizontal) direction.

Finally, supposing that q > 0 leads to the eigenvalues

where p > 0 is such that p² = q. This scenario corresponds to all paths being periodic, which is manifestly not the case for the equatorial flow we are interested in modelling.

Hence we conclude that linear fluid motion is not sufficiently rich structurally to model the large-scale coherent flow structures we are seeking, and these must depend fundamentally on nonlinear laws for their emergence and persistence in the regime which is described by geophysical fluid dynamics.

Acknowledgements

The author would like to thank the referees for their helpful and constructive comments. The author would like to thank the Faculty of Mathematics at the University of Vienna, at which part of this research was undertaken, for a very stimulating research environment. This was supported by the WWTF research grant MA16-009.