Large gyres as a shallow-water asymptotic solution of Euler’s equation in spherical coordinates

Starting from the Euler equation expressed in a rotating frame in spherical coordinates, coupled with the equation of mass conservation and the appropriate boundary conditions, a thin-layer (i.e. shallow water) asymptotic approximation is developed. The analysis is driven by a single, overarching assumption based on the smallness of one parameter: the ratio of the average depth of the oceans to the radius of the Earth. Consistent with this, the magnitude of the vertical velocity component through the layer is necessarily much smaller than the horizontal components along the layer. A choice of the size of this speed ratio is made, which corresponds, roughly, to the observational data for gyres; thus the problem is characterized by, and reduced to an analysis based on, a single small parameter. The nonlinear leading-order problem retains all the rotational contributions of the moving frame, describing motion in a thin spherical shell. There are many solutions of this system, corresponding to different vorticities, all described by a novel vorticity equation: this couples the vorticity generated by the spin of the Earth with the underlying vorticity due to the movement of the oceans. Some explicit solutions are obtained, which exhibit gyre-like flows of any size; indeed, the technique developed here allows for many different choices of the flow field and of any suitable free-surface profile. We comment briefly on the next order problem, which provides the structure through the layer. Some observations about the new vorticity equation are given, and a brief indication of how these results can be extended is offered.

Starting from the Euler equation expressed in a rotating frame in spherical coordinates, coupled with the equation of mass conservation and the appropriate boundary conditions, a thin-layer (i.e. shallow water) asymptotic approximation is developed. The analysis is driven by a single, overarching assumption based on the smallness of one parameter: the ratio of the average depth of the oceans to the radius of the Earth. Consistent with this, the magnitude of the vertical velocity component through the layer is necessarily much smaller than the horizontal components along the layer. A choice of the size of this speed ratio is made, which corresponds, roughly, to the observational data for gyres; thus the problem is characterized by, and reduced to an analysis based on, a single small parameter. The nonlinear leading-order problem retains all the rotational contributions of the moving frame, describing motion in a thin spherical shell. There are many solutions of this system, corresponding to different vorticities, all described by a novel vorticity equation: this couples the vorticity generated by the spin of the Earth with the underlying vorticity due to the movement of the oceans. Some explicit solutions are obtained, which exhibit gyre-like flows of any size; indeed, the technique developed here allows for many different choices of the flow field and of any suitable free-surface profile. We comment briefly on the next order problem, which provides the structure through the layer.  The main thrust of this work is to show how the Euler equation (coupled with the equation of mass conservation and the relevant boundary conditions) can lead directly to a solution that recovers the essential structure of a gyre of any size on the surface of the planet. Furthermore, the results allow, in principle, for any choice of velocity profile around the centre of the gyre (and so any observed variations from one gyre to another can be accommodated). We now outline the plan that we shall follow here. First, we present the governing equation of motion, and the equation of mass conservation, written in a rotating spherical coordinate system. This system is non-dimensionalized and scaled in order to produce a version that is relevant to a (relatively) thin layer of fluid on the surface of the Earth, whose dominant motion is purely rotational, so that any vertical motions are small (but still present). Then, we describe the resulting leading-order asymptotic solution which allows for any background vorticity and, at higher order, we indicate how this leads to a corresponding vertical structure in the flow field. However, it is the generality and accessibility of the leading-order problem that is the primary success of the work that we report here. In this model, we can have no surface stresses (wind) to drive the flow nor any decay mechanism-the fluid is inviscid-and so we must, perforce, allow the flow field to have been generated, and then this is to be maintained by a suitable pressure distribution at the free surface. Indeed, the pressure can be chosen in order to accommodate a (locally) non-flat surface, so that the observed domed structure in some gyres can be included. Only the assumption of a thin layer of fluid, which is naturally coupled with weak vertical motion, is needed in order to make headway; this is the problem that we describe in detail here.

Governing equations
We introduce a set of (right-handed) spherical coordinates (r , θ, ϕ): r is the distance (radius) from the centre of the sphere, θ (with 0 ≤ θ ≤ π ) is the polar angle (and then π/2 − θ is conventionally the angle of latitude); ϕ (with 0 ≤ ϕ < 2π ) is the azimuthal angle, i.e. the angle of longitude. (We use primes, throughout the formulation of the problem, to denote physical (dimensional) variables; these will be removed when we non-dimensionalize.) The unit vectors in this (r , θ, ϕ)system are (e r , e θ , e ϕ ), respectively, and the corresponding velocity components are (u , v , w ); e ϕ points from west to east, and e θ from north to south (figure 2).  . The spherical coordinate system, where θ ∈ [0, π ] is the polar angle, ϕ ∈ [0, 2π ) is the azimuthal angle and r is the distance from the origin: (π/2 − θ ) is the conventional angle of latitude, so that the North Pole corresponds to θ = 0, the South Pole to θ = π , and the Equator to θ = π/2.
The Euler equation (describing an inviscid fluid) and the equation of mass conservation are respectively, where p (r , θ , ϕ) is the pressure in the fluid, ρ is the density and (F r , F θ , F ϕ ) is the body-force vector. In view of the particular structure that we shall describe here, the density is taken as constant (although an extension of this work can readily accommodate a fluid which has piecewise constant density with depth). These equations describe the problem in a coordinate system with its origin at the centre of the sphere; we now associate (e r , e θ , e ϕ ) with a point fixed on the sphere which is rotating about its polar axis. Thus, on the left of equation (2.1), we must introduce the additional terms 2Ω × u' (Coriolis) and Ω × (Ω × r') (centripetal acceleration), where Ω = Ω (e r cos θ − e θ sin θ ), u' = u e r + v e θ + w e ϕ , r' = r e r , with Ω ≈ 7.29 × 10 −5 rad s −1 , the constant rate of rotation of the Earth; these two contributions can be written as 2Ω (−w sin θ , −w cos θ , u sin θ + v cos θ) − r Ω 2 (sin 2 θ, sin θ cos θ, 0). (2. 3) The body-force vector, At the free surface, r = R + h (θ , ϕ), where R ≈ 6378 km is the (mean) radius of the Earth, we impose a surface pressure and the kinematic boundary condition respectively. At the bottom of the ocean, r = R + d (θ, ϕ), which we take to be an impermeable, solid boundary, we have the corresponding kinematic condition 3. The problem in spherical coordinates with a thin surface layer We start the discussion of this aspect of the problem by defining a set of non-dimensional variables. To this end, we introduce the length scale H and a suitable speed scale, c ; the length scale is that associated with the average depth of the ocean, and so it is convenient to set The specific choice of speed scale is not critical to the analysis presented here (because of the particular limiting process that we invoke), so any suitable speed will suffice; for example, the average speed of currents that contribute to the motion of a gyre; the speed at the surface of the Earth, at the Equator, due to the rotation of the Earth about its axis; in the light of the hydrostatic pressure distribution (see (3.5) below), perhaps c = H g . The original (dimensional) variables are now transformed, thereby producing the non-dimensional version of all the variables (no primes) where we have used the constant density (ρ ) in the non-dimensionalization of the pressure; the scaling factor, k, associated with the vertical component of the velocity is yet to be chosen. On setting ε = H /R , the governing equations (2.1)-(2.2), for a steady flow, become and Now it is typical of thin-layer approximations-and quite obviously so for physically realistic flows-that the velocity component across the layer is small compared with the speeds along the layer. In the context of the motions observed in gyres, typical horizontal velocities are of the order of 0.01 m s −1 and the ratio of vertical speed to horizontal speed is about 10 −4 to 10 −5 ; see [4][5][6][7].
On this basis, we elect to examine the asymptotic structure of the problem on the assumption that the motion in the vertical direction is suitably weak. The most natural mathematical choice is given by k = ε 2 (and we introduce the parameter ω = Ω R /c , which is held fixed as ε → 0, and so the full effects of rotation are retained). We have now introduced a characterization of the mathematical problem that we shall discuss; of course, it is not intended that the values of ε and k should match a particular pair encountered in the field. It is sufficient to formulate the problem that incorporates the property that k is smaller than ε, which is essentially what is observed. The corresponding equations and boundary conditions, all written in non-dimensional variables, with are therefore and where we have also transformed (h , d ) = H (h, d), consistent with (3.2). We set P = Π − 1 4 ω 2 cos(2θ), (3.11) and then the leading-order problem in ε, defined here by simply setting ε = 0, is given by (3.14) and The boundary conditions appropriate at this order imply that we require the condition on the bottom is automatically satisfied with d(θ, ϕ) = constant, which is the simplification that we invoke in this first phase of the investigation. (We will make some further observations about the bottom boundary condition later.) It is now convenient to redefine the pressure at the surface so that it takes the form which is regarded as the prescribed or chosen pressure that is to be associated with the flow (and corresponding surface distortion) that we describe here. At leading order, this is the pressure that must exist throughout the depth of the fluid layer, by virtue of equation (3.12) (but note that (3.5) shows the variation with depth associated with the hydrostatic pressure distribution); any additional variation of pressure with depth arises at the next order.

Solution of the leading-order problem
The main thrust of this analysis is to construct solutions, at this order, which recover any number of large flow patterns that mirror what is observed in our oceanic gyres. To see how this comes about, we use equation (3.15) to introduce the stream function ψ where subscripts denote partial derivatives; then, from (3.13) and (3.14), the compatibility condition generated by the elimination of Π produces the vorticity equation Here, the vorticity in the flow, at leading order, expressed in spherical coordinates, is The solution of equation (4.2) is constructed by first writing where Ψ is associated with the vorticity of the underlying motion of the ocean (relative to the Earth's surface and not driven by the rotation of the Earth); equation (4.2) then becomes Throughout regions where ∇(Ψ − ω cos θ) = 0, the rank theorem (see [8]) permits us to express the solution of (4.5) in the form F(Ψ − ω cos θ) + 2ω cos θ = total vorticity, where and F is an arbitrary function; this is a new variant of the vorticity equation which, however, does not capture the solutions Ψ = ω cos θ + A of (4.5), representing stationary flows ψ = A, where A is an arbitrary constant. (Note that the classical form of the vorticity equation is the familiar .) The vorticity, which is aligned with the r coordinate, to leading order, therefore comprises two distinct components: that solely due to the rotation of the Earth (ω e = 2ω cos θ: 'spin vorticity') and that due to the motion of the ocean (ω s = F(Ψ − ω cos θ): 'oceanic vorticity'), but this, we note, is coupled to the Earth's rotation; the total vorticity ω T is the sum of these two: ω T = ω e + ω s . Although one of these contributions (ω e ) is completely prescribed, that associated with the movement of the ocean may be chosen to represent a particular gyre. We might expect ). It is obtained by first mapping P to the intersection point P (r cos ϕ, r sin ϕ, r cot θ ) of the ray from the centre of the sphere to P with the cylinder tangent to the sphere along the Equator, then moving P vertically towards the equatorial plane to the point P (r cos ϕ, r sin ϕ, r ln[tan(θ/2)]) to bring about the same north-south length distortion rate as the east-west distortion rate (1/ sin θ ) of the cylindrical projection, and subsequently unrolling the cylinder. The Mercator projection is conformal (preserves angles, transforming meridians into vertical lines and parallels into horizontal lines) but distorts areas, inflating them according to their distance from the equator, so that at the poles the scale becomes infinite (see [9]).
that the vorticity inherent in the motion of the ocean is significantly larger than that provided by the rotation of the Earth; however, a suitable model for a gyre could be to allow parts of the flow field to be irrotational (oceanic vorticity ω s = 0), and even for constant oceanic vorticity (ω s = F = constant) we see that there is no coupling.
We seek a solution of equation (4.6) in the form where A is an arbitrary constant; we observe that the logarithmic term here is that associated with the Mercator projection (figure 3) and with a Legendre function of the second kind; we define ϕ so that the centre of the gyre is always along ϕ = 0. On using this transformation of variables, equation (4.6) becomes The general procedure, at this stage, is to choose or identify those functions F that are likely to be relevant or interesting or, at the very least, enable us to produce some simple, useful, closedform solutions. We mention three cases that are readily accessible (although the third will not be pursued in this initial phase of the work). where α is an arbitrary constant; this expression corresponds precisely with the classical solution for irrotational flow in two-dimensional, planar geometry. (We will describe the properties of this, and other solutions, in the next section.) In order to examine the flow in a gyre in more detail, we should allow the (oceanic) flow to have some vorticity; the simplest choice is-quite obviously-constant vorticity (F = γ ) and then a suitable solution is which can be written as to within an additive constant, with D = γ /β; the strength of the velocity field is now proportional to γ /β, and the choice of β (for a given constant vorticity γ ) controls the type of solutions available (as we will show in the next section).
In the above two examples, the oceanic-flow component uncouples from the flow that is driven by the rotation of the Earth; we now briefly mention an example where this is not the case. This more sophisticated choice involves setting F proportional to its argument, i.e. F is linear in (Ψ − ω cos θ); we write equation (4.6)-the more natural choice, as will immediately become clear-in the form 1 where λ is an arbitrary constant. This equation can be solved by using, for example, spherical harmonic functions; this is a routine but rather tiresome procedure, and not developed here. We do, however, make one observation: a particular integral of equation (4.12) is and so there exists a solution of the homogeneous equation (for λ = −2) which is resonant with this particular rotation-driven solution. This appears to be worthy of further investigation. Finally, we observe that there is a choice of oceanic vorticity which ensures that the total vorticity (ω T ) is zero; this requires F(Ψ − ω cos θ ) = −2ω cos θ, which corresponds to the choice If we ignore the planetary (spin) vorticity (by setting ω = 0), equation

Some properties and examples of the leading-order solution
In order to examine the relevance of the solutions presented in the preceding section, we must develop some detailed properties of them. Most particularly, we must confirm that solutions do exist which describe the rotational motion evident in gyres. The first, and most illuminating property, is the streamline pattern generated by our choices of stream function; these provide confirmation that we have found solutions that correspond to large gyres.
The simplest case arises from the choice of an irrotational flow; some typical streamlines (defined by ξ = constant, for a given A; see (4.7)) for two such flows (A = 0) are shown in figure 4.   (All the figures that we present show the streamlines plotted in (θ, ϕ)-coordinates, but represented in the (flat) Cartesian plane; the 'physical' streamlines are the corresponding curves in (θ, ϕ)coordinates on the surface of the sphere; note that θ is measured downwards, north to south, from zero in these figures.) Figure 4 shows the streamlines associated with the 'oceanic' part of the flow field, i.e. based solely on the stream function given in (4.9); the streamlines as seen on the surface of the ocean, however, are generated by the stream function ψ, given in (4.4). Certainly, we must expect that the shape of these streamlines is dominated by the term in Ψ because if ω cos θ, the term due solely to the rotation of the Earth, was the dominant term, this would correspond to streamlines θ = constant. In figure 5, we show the effects of including this term at a moderate  level (ω/α = 1 2 ); the streamlines are compared with those that correspond to the flow shown in figure 4 (no spin vorticity), the distortion being evident but not overwhelming.
In the case of constant-vorticity flow, we have-not surprisingly-a more complex structure; see (4.11). For β ≤ 2, there exist solutions with closed streamlines, but these necessarily extend across the whole globe north-to-south (0 < θ < π), unless we arbitrarily restrict the region of applicability; there are also solutions with open streamlines. An altogether more satisfactory situation arises for β > 2; closed streamlines exist and, for a range of values of A and β, there are two families: one in the Northern Hemisphere and one in the Southern, the switch from one to the other requiring a change of sign in A. Some of these solutions do extend across the Equator (for the choice A = 0, for example)-so we reject these as not being relevant in this context-but for larger βs and suitable A, we find many solutions that are restricted to just one hemisphere and, furthermore, their extent is limited by a bounding streamline; two examples are shown in figure 6. These two families are necessarily limited by the Equator, across which neither exist; again, we could arbitrarily limit the region of applicability to ensure that the gyre stops some distance short of the Equator or, more satisfactorily, choose the parameters so that the bounding streamline is positioned as we wish it (as we have done here). Indeed, this is precisely the model we would choose to use when we combine together regions of different vorticity: perhaps, an outer region of constant vorticity and an inner one that is irrotational or, possibly, stationary. As the parameter β increases, so more variation is possible in the shape of the gyre; again, with a bounding streamline that ensures the gyre sits wholly within one hemisphere, we show another two solutions (figure 7). It is evident that there is considerable freedom available within the structure of these solutions, allowing the possibility of producing models for gyres that extend over considerable regions of the oceans and which correspond, roughly, with the overall shapes that are observed. The four examples shown in figures 6 and 7 indicate quite clearly how the shape and the position of the gyre can be adjusted within one hemisphere. (All these streamline patterns will be (slightly) distorted by the inclusion of the spin contribution, as shown for the irrotational case in figure 5).
On the basis of the explicit solutions that we have obtained for the streamlines (and so for the associated velocity field), we can determine the pressure distribution that is required to maintain the motions described here. In the case of the irrotational flow, represented by the stream function (4.9), the pressure, Π (θ , ϕ) (see (3.11)), is readily obtained; from this we can find the pressure to within an arbitrary additive constant. The free surface, z = h(θ, ϕ), is an arbitrary function of the stream function, ψ: this automatically satisfies the second equation in (3.16). We may therefore choose h to accommodate any observed surface profiles (such as those with a domed structure), and (5.1) then gives the surface-pressure distribution necessary to maintain the flow and the shape of the surface. The corresponding result for constant oceanic vorticity is to within an arbitrary additive constant, where h is again an arbitrary function of ψ.

Some observations concerning the solution at the next order
The structure of the governing equations, and associated boundary conditions, (3.6)-(3.10), indicates that it is reasonable to seek an asymptotic solution based on the asymptotic sequence {ε n } n≥0 . Indeed, to some extent, this guided our choice of parameter dependence that describes the velocity component in the z-direction: k = ε 2 . On this basis, we seek a solution expressed in the form where q represents each of the unknowns u, v, w, P, h (although, in the case of h, the dependence on z is omitted). Thus, the solution designated by the zero subscript is precisely that obtained in (We will provide some comments concerning the validity of this asymptotic solution later.) The system of equations and boundary conditions that describe the O(ε) problem are v 2 0 + w 2 0 + 2ωw 0 sin θ + ω 2 sin 2 θ = P 1z , (6.1) and u 0z + (v 1 sin θ) θ + w 1ϕ sin θ = 0, (6.4) with and In this formulation of the problem that describes the perturbation of the solution found earlier, we have invoked a Taylor expansion about z = h 0 (θ , ϕ) in the surface boundary conditions; this is a permitted procedure here because, as is already evident, the problem is polynomial in z (which is the typical behaviour in the shallow-water approximation). The inclusion of both h 1 and P 1s is to allow some freedom (hopefully leading to a suitable simplification) in the small correction terms to the surface profile and the surface pressure, respectively. From equation (6.1), and the associated boundary condition (6.5), we obtain the correction to the pressure which provides a contribution to the hydrostatic pressure distribution. We could, at a later stage, make choices relating to the (small) adjustments to both the surface profile and the surface pressure distribution. However, an obvious route, allowing some simplification, is to choose the relation between h 1 and P 1s so that pressure perturbation takes the form P 1 = (v 2 0 + w 2 0 + 2ωw 0 sin θ + ω 2 sin 2 θ)z. (6.8) It then follows that there is a solution of equations (6.2)-(6.3), for v 1 and w 1 , for which these two functions are proportional to z; we set v 1 = zV 1 (θ , ϕ), w 1 = zW 1 (θ, ϕ), (6.9) and then, from (6.4) and (6.7), we obtain In this development, u 0 at the surface is known (from (6.10)), and then the kinematic condition, (6.6), is satisfied by a suitable choice of h 1 (θ , ϕ); this, with the choice leading to (6.8), fixes the pressure perturbation at the surface. The problem providing the construction of solutions for V 1 and W 1 , associated with appropriate choices for the background flow (represented by Ψ ), can be formulated, although explicit, closed-form, solutions of the underlying second-order partial differential equation, for general Ψ , cannot be written down (although this is possible for the case of irrotational flow). Nevertheless, our outline demonstrates the essential character of the problem at this next order, involving a polynomial structure in z (and no obvious non-uniformities exist provided that the gyres sit strictly within a hemisphere, and not too close to a pole or the Equator). Furthermore, with this additional information, we are able to provide a representation of the way in which the direction of the horizontal component of the velocity vector varies with depth; this is controlled, at this order of approximation, by v 0 (θ , ϕ) + εzV 1 (θ, ϕ) w 0 (θ , ϕ) + εzW 1 (θ, ϕ) although the details will depend on how all these functions depend on θ and ϕ.

Discussion
The ideas presented here have demonstrated, based on Euler's equation with vorticity included, that models for gyres can be constructed. This has been possible by invoking a thin-layerequivalently a shallow-water-approximation, and nothing more. This particular asymptotic interpretation has the benefit of generating a problem, at leading order, which is both nonlinear and in a rotating, spherical coordinate system. Furthermore, this system can be solved in terms of the underlying vorticity and, in some cases, explicit, simple solutions can be obtained. It is clear, however, that the success of this approach is down to the exclusion of any variation in the vertical coordinate at leading order; this is, of course, precisely the hallmark of classical shallowwater theory. Indeed, consistent with this, the structure in the vertical (z) direction is also classical: polynomial in z. Because gyres are predominantly circulatory in motion, with very weak vertical motion (as confirmed in field data), we submit that we have generated a useful and relevant model. The difficulties encountered at higher order (i.e. being unable to construct simple, closedform solutions) are certainly irritating, but this aspect is relatively insignificant when set against the type of flows available-and there are many possibilities-in our leading-order solution. This is the main thrust of the work presented here.
The structure of gyres, as generated by our vorticity equation (for the oceanic component of the flow field), (4.6), can take many and varied forms, the extent of the possibilities is, we suggest, worthy of further investigation. We have presented some examples of the streamlines that describe the oceanic component, both for irrotational flow and constant vorticity. (One intriguing and unlooked for ingredient is the appearance-the natural choice when solving-of a transformation of coordinates that involves the essential component of the Mercator projection; see (4.7).) In the case of irrotational flow (given by (4.9)), the flow can be centred anywhere on the surface of the sphere (by adjusting the constant A), and extends in closed curves around this centre; to be physically relevant and appropriate, the region of the flow must be restricted. This ensures that the gyre can be allocated to one of the hemispheres, and outside this region we presumably impose a different type of flow field (which could be no motion at all). It is gratifying to obtain streamlines that correspond, in their general shape, to the flow configuration of gyres, but the irrotational case is far from satisfactory. Firstly, there is no natural boundary to the extent of such a flow (and we might hope and expect this to be the case for a gyre); secondly, we might reasonably argue that, in any event, the gyres are not likely to be irrotational structures. Either there is some appropriate vorticity distribution (within some regions of our oceans), or the flow is essentially stationary. So although it is reasonable to note the existence of solutions for irrotational flow, it would seem that they do not provide useful models for gyres.
On the other hand, flows with constant (oceanic) vorticity offer far better prospects. For sufficiently large β-indeed, for β > 2-and then for suitable A (see (4.11)), there exist solutions that sit in only one hemisphere and which have a bounding streamline (so there is no question of arbitrarily truncating the region). As β and A are varied, so the centre of the gyre, the shape of the streamlines and the extent of the region all change. We have not carried out a comprehensive examination of exactly how β and A affect the properties of the gyre, so this is something else for future study. We have, however, explored the nature of the problem sufficiently to confirm the general picture described above, and to produce a number of examples.
It is clear, using our description of the underlying flow (as implied by equations (3.12)-(3.16)), that we must have some oceanic vorticity, for otherwise the streamlines would be simply θ = constant, which do not represent gyres. Because we may, in principle, assign any vorticity to the flow, and we have no mechanism in our model for initiating or driving the gyre, we must choose the sign of the vorticity to match the observed clockwise (Northern Hemisphere) or counter-clockwise (Southern Hemisphere) rotations. This fixes the overall flow structure; in addition, we may impose any free-surface profile, by selecting a suitable h(ψ). Then, in order to maintain this motion (and surface), we require a surface-pressure distribution; this we may regard as a modelling of-and a replacement for-the effects of pressure and wind-action at the surface; see (5.1) and (5.2).
We believe that our approach to a better understanding of some of the underlying principles and properties of gyres is, within the confines of an inviscid theory, built on a firm foundation. The asymptotic method used here-the shallow-water approximation-that leads to the basic flow patterns that we have described, is the most natural one in this context, and it ensures that all the relevant contributions are retained. However, what we describe here can be no more than the first stage in an investigation that must, we suggest, be taken further. There is a lot of freedom afforded by our formulation of the problem, and so there are many avenues that need to be explored; we mention some of the more obvious ones.
Clearly, it would be wise to seek solutions associated with other choices of vorticity, and one obvious candidate here is to examine how expansions in spherical harmonics might play a role. With other solutions available (and this also applies to the zero and constant vorticity examples found here), we should investigate how regions of different vorticities can be combined; this is not a simple exercise. The difficulty arises because, based on what we have discovered so far, different vorticity distributions lead to different shapes for the streamlines; cf. the classical Rankine vortex, for which all streamlines are circles with a common centre. It is therefore impossible, across a streamline, to switch from one vorticity to another. We will need some mechanism for producing a transition between two sets of streamlines associated with two different vorticity distributions. The only possibility, in general, is to allow a region of no flow (stationary conditions) between the two areas of different vorticity; but perhaps other special flows can be discovered which overcome this difficulty. As our results currently stand, we can certainly have a constant vorticity flow and introduce a stationary flow both inside and outside two specific streamlines. (The former would correspond to the stagnant areas observed in the central regions of some of our gyres.) Related to this, the new type of vorticity equation, (4.6), certainly needs a detailed investigation, which might throw some light on the choices of vorticity distribution that are allowed and reasonable; a few initial observations have been provided in appendix A.
It is evident that some work needs to be done on the higher-order terms in our asymptotic solution. In a more complete discussion, we should aim to develop a fairly comprehensive description of the O(ε) term (which we have managed to do in the irrotational case, this being reduced to a quadrature; the complexity of the resulting solution, and its questionable relevance to the observed motions of gyres, persuaded us not to reproduce it here). With a clearer idea of the structure of the O(ε) solution, we can then attempt some general observations about the higherorder terms. We can report that our experience with this problem suggests that the asymptotic expansion is uniformly valid provided that the flows are not too close to the poles or the Equator; the polynomial behaviour in z precludes any difficulties in this variable.
All the comments thus far relate, more-or-less, to aspects of the model and the approach adopted here. However, we should consider whether there are any reasonably accessible avenues that might lead to an improvement in the model. One area of significance, that is relevant to gyres that we observe, is the role of the topography (both undersea and the surrounding land masses). Indeed, it is generally accepted that these are fundamental to the definition of the regions occupied by gyres. The solution presented earlier used the simplest bottom boundary condition: d = constant; but our formulation admits a more general choice. As with the surface kinematic condition (see the comments after (5.1)), we may introduce a bottom topography that is described by d = d(ψ). A viable model, within our scenario, is to introduce precisely this form of bottom condition, and to regard this as the primary driver for the geometric constraints on the gyre. We could then argue that this leads to the description of the flow field in terms of our stream function, with the associated structure of the free surface: h = h(ψ). Of course, to accommodate other behaviours of the bottom, we would need to specify appropriate functions d(θ, ϕ); see (2.7) and (3.10). However, in order to progress the general approach that we have developed here, we cannot arbitrarily assign d(θ , ϕ). Nevertheless, it might be possible to use the method of multiple scales, and introduce a depth variation (which could be large, to the extent of allowing for the appearance of shores and land masses) but which varies slowly in θ and ϕ, e.g. d = d(εθ, εϕ), or something similar. Such a possibility needs further investigation.
In conclusion, we have demonstrated that conventional ideas of fluid mechanics, with associated standard approximation procedures, can make a contribution to the study, and understanding, of gyres. As we have made clear, this can only be the start: there are many avenues that need considerably more work done on them, involving deeper and more extensive examination.
Data accessibility. All data are provided in full in the paper. Authors' contributions. The authors contributed equally to this study. Competing interests. We declare we have no competing interests. Funding. R.S.J. gratefully acknowledges the financial support of the Erwin Schrödinger Institute, Vienna, as an ESI Senior Research Fellow, during the time that the work reported here was undertaken.
Acknowledgements. The authors are grateful for helpful comments from the referees.
Appendix A. Some general considerations about the nonlinear problem for Ψ Outside the realm of explicit solutions to (4.6) one can expect to identify key qualitative features of the flow dynamics from specific properties of the vorticity F by using techniques from the theory of elliptic partial differential equations, as described below. To gain further insight, one can rely on approaches specific to Lagrangian coherent structures (see [11]), a direction of research that we plan to develop in the near future.
Let us perform the stereographic projection of the unit sphere centred at at origin, from the North Pole to the equatorial plane (figure 8): where (r, φ) are the polar coordinates in the equatorial plane. We have After several cancellations one can see that equation (4.6) takes the form For constant vorticity (F = γ ), the general solution of (A 2) is seen to be where Φ is an arbitrary harmonic function. Note that the solution (4.10) is among the family (A 3). It is convenient to use (4.4) to write (A 2) equivalently as or, with the understanding that (x, y) are the Cartesian coordinates in the complex ξ -plane, as the semi-linear elliptic equation (i) Existence of solutions for a given smooth vorticity F in a specified bounded region O with smooth boundary ∂O: Writing the lower-order part of (A 4) as F (x, y, ψ), we say that a smooth function ψ ± , defined on an open planar set whose interior comprises O ∪ ∂O, is a sub/super-solution if ±[ ψ ± + F (x, y, ψ ± )] ≤ 0 throughout O. For example, if F(0) ≥ 0, then ψ + ≡ 0 is a supersolution in the Northern Hemisphere as (A 1) yields x 2 + y 2 > 1.
On the other hand, if F(0) ≤ 0, then ψ − ≡ 0 is a sub-solution in the Southern Hemisphere because (A 1) now ensures x 2 + y 2 < 1; more sophisticated sub/super-solutions are provided by the functions ψ associated with the special solutions (4.9) and (4.11), or by seeking radial functions (in which case one has to solve a second-order ordinary differential equation). If F is smooth and if one can find sub/super-solutions ψ ± such that ψ − ≤ ψ + throughout O and such that ψ − ≤ 0 ≤ ψ + on ∂O, then a functional-analytic approach (see [12]) shows that there exists a solution ψ to (A 4) such that ψ − ≤ ψ ≤ ψ + throughout O. (ii) The comparison method and uniqueness: If the function F is non-decreasing, using the strong maximum principle (see [13]), we can reduce the above requirements that ψ − ≤ ψ + throughout O and ψ − ≤ 0 ≤ ψ + on ∂O to merely ψ − ≤ 0 ≤ ψ + on ∂O by taking advantage of the comparison principle: if ψ is a solution and ψ + a supersolution with 0 ≤ ψ + on ∂O, then the mean-value theorem yields 0 ≤ (ψ − ψ + ) − 4[F(ψ) − F(ψ + )] (1 + x 2 + y 2 ) 2 = (ψ − ψ + ) − 4F (ψ) (1 + x 2 + y 2 ) 2 (ψ − ψ + ), and the strong maximum principle ensures ψ ≤ ψ + in O, if this inequality holds on ∂O; a similar comparison can be made with a sub-solution. In particular, in this setting, we therefore have uniqueness. The above approach is also valid without the monotonicity assumption on F, but only if the domain O is sufficiently narrow or of sufficiently small area-see [14]. Moreover, if there exist constants M > 0 and q ≥ 1 with |F(s)| ≤ M(1 + |s| q ), s ∈ R, ( A 6 ) then the solution between the sub-and super-solution is a local minimum of the functional where F (t) = t 0 F(s) ds; see [12]. (iii) Numerical methods. The variational formulation (A 7) permits the implementation of the finite-element method, providing a versatile computational approach towards gaining quantitative and qualitative insight, but one can also use finite-difference methods-see the discussion in [15].