## Abstract

In rapidly rotating spheres, the whole fluid column, extending from the southern to northern spherical boundary along the rotation axis, moves like a single fluid element, which is usually referred to as geostrophic flow. A new Legendre-type polynomial is discovered in undertaking the asymptotic analysis of geostrophic flow in spherical geometry. Three essential properties characterize the new polynomial: (i) it is a function of *r* and *θ* but takes a single argument , which is restricted by 0≤*r*≤1 and 0≤*θ*≤*π*, where (*r*,*θ*,*ϕ*) denote spherical polar coordinates with *θ*=0 at the rotation axis; (ii) it is odd and vanishes at the axis of rotation *θ*=0, and (iii) it is defined within—and orthogonal over—the full sphere. As an example of its application, we employ the new polynomial in the asymptotic analysis of forced geostrophic flows in rotating fluid spheres for small Ekman and Rossby numbers. Fully numerical analysis of the same problem is also carried out, showing satisfactory agreement between the asymptotic solution and the numerical solution.

## 1. Introduction

Many planets and stars share several common features: fluidic, rotating rapidly and spherical geometry. In modelling fluid motion in the interiors of rapidly rotating planets or stars, we may take a nearly homogeneous fluid with small kinematic viscosity *ν* that is confined within a sphere of radius *r*_{0} and that rotates uniformly with a constant angular velocity , where is a unit vector and *Ω* is constant. In the understanding of various fluid dynamical problems in rapidly rotating planets or stars, we often encounter the following dimensionless equations in the co-rotating reference frame:
1.1
and
1.2
where **u** is the three-dimensional fluid velocity, *p* is a reduced pressure, the Ekman number provides the measure of relative importance between the viscous force and the Coriolis force, and the Rossby-type number *ϵ* quantifies the strength of the nonlinear effect. In equation (1.1), **f** may represent the Reynolds stress **u**·∇**u** (e.g. Busse 1983; Hollerbach & Kerswell 1995; Gillet & Jones 2006) or the magnetic Lorentz force (∇×**B**)×**B**, where **B** denotes the magnetic field (see Moffatt 1978; Hollerbach *et al.* 1992; Pais & Jault 2008; Livermore *et al.* 2009). Here, we have assumed ∂**u**/∂*t*=0 since the focus of this paper is primarily on the geostrophic flow. In the context of geophysical fluid dynamics for the Earth’s liquid core, both the Ekman number *E* and the Rossby-type number *ϵ* are small. An important mathematical problem of geophysical relevance is: for a given **f**, find the asymptotic solution to equations (1.1) and (1.2) for *E*≪1 and *ϵ*≪1 together with the non-slip boundary condition
1.3
on the bounding surface, , of the sphere, where denotes the vector normal to .

A well-known, profoundly important result in rapidly rotating systems is the Taylor–Proudman theorem that can be readily derived from equations (1.1) and (1.2) by postulating that the fluid motion is sufficiently slow *ϵ*→0 and the effect of viscosity is sufficiently weak *E*→0. To leading-order approximation, the momentum equation (1.1) reduces to
1.4
subject only to the inviscid condition
1.5
Taking the curl of equation (1.4) results in a fundamental result, the Taylor–Proudman theorem,
1.6
which states that slow steady motions in a rotating inviscid fluid cannot vary along the spin axis, i.e. are two dimensional with respect to the direction of the rotation axis. Such two-dimensional flows are said to be geostrophic, an intriguing and counterintuitive phenomena in rotating fluids. In spherical geometry, solutions satisfying both equations (1.4) and (1.5) are either given by
1.7
or by
1.8
where *U*(*x*) denotes an arbitrary function of *x*, (*r*,*θ*,*ϕ*) are spherical polar coordinates with the corresponding unit vectors and with *θ*=0 at the axis of rotation , whereas (*s*,*ϕ*,*z*) denote cylindrical polar coordinates with *s*=0 at the rotation axis. It is of importance to emphasize that the geostrophic flow given by (1.7) or (1.8), though axisymmetric with ∂/∂*ϕ*=0, is confined within the whole fluid sphere.

A deep root of mathematical difficulties in solving the problem of geostrophic flows in spherical geometry lies in the conflict between the spherical and cylindrical polar coordinates. On the one hand, the Taylor–Proudman theorem (1.6) is strongly indicative of adopting cylindrical polar coordinates. On the other hand, the boundary condition (1.3) on the spherical bounding surface clearly suggests the use of spherical polar coordinates. This spherical/cylindrical conflict in the analysis of spherical geostrophic flows suggests the necessity of a special polynomial, which will be denoted as *G*_{n}(*x*) in this paper, that should possess the following three properties. First, to represent the geostrophic flow in a sphere, the polynomial must be a function of two variables *r* and *θ*, restricted by 0≤*r*≤1 and 0≤*θ*≤*π*, but taking only a single argument . Second, since the geostrophic flow vanishes at the axis of rotation, the polynomial must be odd, i.e. with *k*=1,2,3,… and (2*k*−1) denoting its order. Third, since the geostrophic flow is defined within the sphere, the polynomial should be orthogonal over the sphere, i.e.
1.9
We shall refer to as the geostrophic polynomial. If such a polynomial meeting the three requirements exists and can be found, it resolves the conflict between the spherical and cylindrical coordinates and can be employed in the representation of a geostrophic flow within a rapidly rotating fluid sphere for the asymptotic analysis of equations (1.1) and (1.2) in the limit *E*→0.

In what follows, we shall begin in §2 by presenting a newly discovered polynomial meeting the above three requirements. An asymptotic solution to equations (1.1) and (1.2), subject to the non-slip boundary condition (1.3), is then derived by making use of the new polynomial in §3. This is followed by fully numerical analysis of the same problem in §4, which validates the asymptotic analysis. The paper closes in §5 with a summary and some remarks.

## 2. A new Legendre-type polynomial

Prior to discussing the new geostrophic polynomial, it is illuminating to look at the Legendre polynomial of odd order with *k*≥1. The odd Legendre polynomial is an eigenfunction of the Legendre differential equation
2.1
where we take as its argument in order to emphasize that in the problem of spherical geometry is usually defined on a spherical surface. The general form of the Legendre polynomial of odd degree (2*k*−1) is given by
2.2
which is *orthogonal over a spherical surface*,
2.3
In geophysical and astrophysical applications, the Legendre polynomial (2.2) is usually employed to represent an axisymmetric flow on a spherical surface.

In order to represent the geostrophic flow (1.7) in a rotating sphere, we require a new polynomial that is defined within the sphere, dependent only on and orthogonal over the sphere. In comparison with the Legendre differential equation (2.1), which is derivable from a Laplace equation via the method of separation of variables, it seems that there are no standard procedures for deriving a differential equation that gives rise to the required polynomial. After making a huge effort with many different attempts, we have successfully found a second-order ordinary differential equation
2.4
whose eigenfunction gives rise to the geostrophic polynomial, where *k* is an integer. In equation (2.4), we take as the argument of the polynomial to stress that it is defined within the sphere. The general form of the geostrophic polynomial of order (2*k*−1) satisfying equation (2.4) is
2.5
where (*r*,*θ*) are restricted by 0≤*r*≤1 and 0≤*θ*≤*π*. By a lengthy mathematical manipulation using equation (2.4), we obtain
for *k*≥1 with 0!!=1, which can be used to normalize the geostrophic polynomial *G*_{2k−1}. It should be noted, because of the term in equation (2.4), that the Legendre differential equation (2.1) cannot be transformed into the second-order differential equation (2.4) for the geostrophic polynomial.

In contrast to the Legendre polynomial (2.2), it is of primary importance to note that the geostrophic polynomial (2.5) is *orthogonal over the full sphere*. The orthogonality, similar to that of the Legendre polynomial, can be readily shown directly from the differential equation (2.4) via the standard procedure:
if *k*≠*n*, where . It may be of interest to note that Chebyshev polynomials of the second kind, *U*_{n}(*s*), defined as
have a different weight function with a different interval. It appears that this particular weight function, *w*(*s*)=*s*(1−*s*^{2})^{1/2}, for the geostrophic polynomial makes it difficult to apply the classical Rodrigues’ formula to the present problem.

Parallel to a recurrence relation for the odd Legendre polynomial (2.2), which is
2.6
there also exists a recurrence relation connecting three successive odd geostrophic polynomials:
2.7
which starts with
and
Furthermore, we can also derive a recurrence relation for the geostrophic polynomial involving its derivative:
2.8
for *k*≥2, which is also analogous to, though different from, that of the Legendre polynomial.

To illustrate the similarity, as well as the key difference, the structures of both the geostrophic polynomial *G*_{2k−1} and the Legendre polynomial *P*_{2k−1} are displayed in figure 1 for *k*=1,2,3,4,5,6. Displayed in figure 1*a* is the profile of the Legendre polynomial computed directly from the expression (2.2) for which *P*_{2k−1}(1)=1. However, the geostrophic polynomial in figure 1*b* is scaled by . In contrast to *P*_{2k−1}(1)=1, we have
for *k*≥1. Although the two polynomials, and , in figure 1 show a broadly similar profile, the most significant difference is that the geostrophic polynomial is implicitly three dimensional, defined in—and orthogonal over—the full sphere, whereas the Legendre polynomial is implicitly two dimensional, defined on—and orthogonal over—a spherical surface.

Since the existing orthogonal polynomials, such as the Legendre and Chebyshev polynomials, represent the special cases of the Jacobi polynomial family that contains two free parameters (Courant & Hilbert 1953), it would be no surprise that the new geostrophic polynomial is also connected with the Jacobi polynomial family by making a proper choice of the two free parameters together with a suitable transformation. Moreover, it should be noted that the geostrophic polynomial is fundamentally different from the Worland polynomial (Livermore *et al.* 2007), which has been used as a set of radial basis functions in full sphere computations.

Finally, since the powers 1,*x*,*x*^{2},… that are orthogonalized in a given interval form a complete system of functions (Courant & Hilbert 1953), any profile of the geostrophic flow equation (1.7) that is continuous in a fluid sphere may be approximated by a linear combination of the geostrophic polynomials. In the next section, we shall use in the representation of geostrophic flows to solve equations (1.1) and (1.2) with the non-slip condition asymptotically for *E*≪1 and *ϵ*≪1.

## 3. Asymptotic analysis

It is known that geostrophic flows (or the *z*-independent differential rotation) in spherical geometry may be generated and maintained in a number of ways. For example, the columnar convection rolls in a rotating sphere may be approximately written as (Busse 1970), where is complex and *m* is an azimuthal wavenumber of the rolls. The nonlinear interaction of the rolls would yield the axisymmetric forcing , as well as other non-axisymmetric terms, in equation (1.1) (see Busse 1983; Zhang 1992; Gillet & Jones 2006). For the purpose of illustrating an application of the new polynomial *G*_{2k−1}, we shall undertake an asymptotic analysis of equations (1.1) and (1.2) for *E*≪1 and *ϵ*≪1 by assuming with 0≤*r*≤1,0≤*θ*≤*π*,0≤*ϕ*≤2*π*. Note that in this paper is arbitrary, which is not derived from the problem of nonlinear convection or magneto-convection in a rotating fluid sphere.

Since there is indication that the effect of spatial non-uniformities of a weakly viscous flow may generally obscure the form of the asymptotic expansion even in the first viscous corrective terms in rotating fluid systems (e.g. Zhang & Liao 2008), we shall not take *E*^{1/2} as an expansion parameter in our asymptotic analysis. We assume that an asymptotic solution to equations (1.1) and (1.2), when *E*≪1 and *ϵ*≪1, can be expanded in the following manner
3.1
where **u**_{0} and *p*_{0} denote the leading-order inviscid interior solution, and are the corresponding viscous boundary correction and **u**_{1} and *p*_{1} represent the secondary interior solution with
In the asymptotic expansion (3.1), weakly viscous action on the leading-order flow **u**_{0} induces a thin viscous boundary layer on . By ejecting the normal mass flux from the thin boundary layer, the viscous effect drives the secondary interior flow **u**_{1}, communicates to the interior fluid and enables us to determine **u**_{0} via the higher order mathematical problem governing **u**_{1}.

After substitution of the expansion (3.1) into the governing equations (1.1) and (1.2), the first problem in the asymptotic sequence is concerned with the zero-order interior solution of the inviscid motion, which is governed by the equations
3.2
and
3.3
subject to
3.4
on the bounding surface of the sphere. The general solution satisfying equations (3.2)–(3.4) can be expressed in terms of the geostrophic polynomial *G*_{2k−1}
3.5
and
3.6
where the polynomial , which is obtained by solving equation (3.2), takes the form
3.7
Note that is an even polynomial and of one order higher than *G*_{2k−1} and that coefficients in equations (3.5) and (3.6) are unknown and cannot be determined in the leading-order problem.

Since the geostrophic flow **u**_{0} in equation (3.5) does not satisfy the non-slip boundary condition, a viscous boundary-layer correction, , is required to ensure that the tangential velocity vanishes at the bounding surface of the sphere. By virtue of the properties of viscous boundary layers in a rapidly rotating sphere, we can derive a fourth-order differential equation from equation (1.1) for *E*≪1 and *ϵ*≪1:
3.8
where *ξ* denotes a stretched boundary-layer variable *ξ*=*E*^{−1/2}(1−*r*). The non-slip boundary condition (1.3) implies that equation (3.8) should be solved subject to the following four boundary conditions:
The fourth-order differential equation (3.8) together with the four boundary conditions determines the solution of a spherical Ekman boundary layer. A straightforward analysis gives rise to
3.9
It is anticipated that the singularity in equation (3.9) at the equatorial plane, where the boundary-layer solution breaks down, would not significantly affect the key properties of the asymptotic solution at an asymptotically small *E* (e.g. Greenspan 1968). The value of coefficients , however, still remains undetermined.

It becomes imperative to consider the third problem in the asymptotic sequence describing the secondary interior flow **u**_{1}, which would enable us to determine the value of the coefficients . The small interior flow **u**_{1} is governed by the equations
3.10
and
3.11
The partial differential equations (3.10) and (3.11) represent an inhomogeneous boundary-value problem that requires a solvability condition. It is worth mentioning that, because *E*^{1/2} is not used as an expansion parameter in the expansion (3.1), the term *E*∇^{2}**u**_{0}, in connection with the interior viscous dissipation, must be retained in equation (3.10) and can be significant if either *E* is not sufficiently small (when comparing with the numerical results) or the scale of **f** is not sufficiently large.

After obtaining the boundary-layer flux from solution (3.9) via a standard analysis (see Greenspan 1968), we derive the solvability requirement for the inhomogeneous equation (3.10), which is
3.12
where *n*=1,2,3,…, giving rise to a system of linear equations for , *k*=1,2,3,…, which can be readily solved. Several important features should be noted in the solvability condition (3.12). First, for a given small *ϵ* and if **f** has large scale, the final term in equation (3.12) may be neglected when *E* is sufficiently small. Second, the integral of the final term in equation (3.12) has the property that
when *n*≥*k*. Third, while an analytical expression can be found for the final term in equation (3.12), the integrals in the first term, because of the factor , have to be obtained by numerical integration.

After are found through the solution of the solvability condition (3.12) for a given *ϵ* and **f**, the leading-order solution to equations (1.1) and (1.2) describing the geostrophic flow and satisfying the non-slip boundary condition can be written as
3.13
whose total kinetic energy, *E*_{kin}, is defined as
In equation (3.13), the size of *K*, representing the total number of polynomials *G*_{2k−1} required in the asymptotic expansion, is only dependent on the spatial profile of **f**. For a moderate scale of **f**, as discussed in the next section, typically *K*≤*O*(10). By a straightforward manipulation, we obtain, to leading approximation, the analytical expression for *E*_{kin}:
3.14
Here, the two asymptotic expressions, equations (3.13) and (3.14), describe the primary features of a geostrophic flow in a rapidly rotating fluid sphere.

In general, three steps are required to compute an asymptotic solution of the geostrophic problem—governed by equations (1.1) and (1.2) with the non-slip condition (1.3)—for a given **f** with *E*≪1 and *ϵ*≪1: (i) carry out the relevant integrals in equation (3.12) analytically or numerically, (ii) solve the solvability condition (3.12) to determine , and (iii) finally, substitute into the asymptotic expressions (3.13) and (3.14) for the geostrophic flow and its kinetic energy. Several examples of the asymptotic solution at various values of *E* for given **f** are shown in figure 2 and table 1, which will be validated by fuller numerics of the same problem at exactly the same parameters, and discussed in detail in §4.

## 4. Numerical analysis

The primary purpose of our numerical analysis, which is valid for both large and small *E*, is to compare with the asymptotic analysis that is valid only for *E*≪1. In the numerical analysis, we expand the velocity **u** as a sum of poloidal (*v*) and toroidal (*w*) components by writing
4.1
where the condition ∇·**u**=0 is automatically satisfied. Applying **r**·∇× and **r**·∇×∇× onto equation (1.1), we obtain the two independent scalar equations,
4.2
and
4.3
where the differential operators and are defined as
and
In terms of *v* and *w*, the non-slip boundary condition (1.3) becomes
4.4
Here, we have assumed that **f** is axisymmetric (∂/∂*ϕ*=0) and, consequently, the corresponding solution to equations (1.1) and (1.2) is also axisymmetric.

Equations (4.2) and (4.3) are then solved numerically by expanding the velocity potentials in terms of Legendre polynomials and of radial functions that satisfy the no-slip boundary condition
4.5
and
4.6
where *T*_{n}(*x*) denotes the standard Chebyshev function, and *w*_{ln} and *v*_{ln} are unknown coefficients to be determined. Note that the factor *r*^{l} in the expansions (4.5) and (4.6) is required to remove the numerical singularity at the origin *r*=0 in a fluid sphere (see Livermore *et al.* 2007). On substitution of equations (4.5) and (4.6) into equations (4.2) and (4.3) and application of the standard numerical procedure, we are able to derive a system of linear algebraic equations for *w*_{ln} and *v*_{ln} that can be found by an iterative numerical procedure. In any practical computation, we have to specify a profile for the forcing term **f** in equation (1.1). For the purpose of simplicity without loss of key physics, we shall take
4.7
where *N* may be regarded as a parameter controlling the spatial scale of **f**.

If the numerical analysis is carried out at a sufficiently small *E*, its solution should be comparable directly with the asymptotic solution (3.13) for *E*≪1. Kinetic energies, *E*_{kin}, of several asymptotic and numerical solutions obtained at exactly the same parameters are shown in table 1 for three different Ekman numbers, *E*=10^{−3},10^{−4} and 10^{−5}. In general, the asymptotic expression (3.14) yields a satisfactory agreement with the numerical simulation. In particular, the agreement between the asymptotics and numerics, as expected, becomes increasingly better with decreasing value of *E*. For example, the difference for *N*=2 between the asymptotics (*E*_{kin}=4.752×10^{−7}) and numerics (*E*_{kin}=4.121×10^{−7}) at *E*=10^{−3} is more than 10 per cent. At a smaller *E*=10^{−5}, however, the discrepancy between the asymptotics (*E*_{kin}=1.166×10^{−4}) and the numerics (*E*_{kin}=1.142×10^{−4}) reduces to about 2 per cent. Moreover, it is worth mentioning that, in calculating the kinetic energy *E*_{kin} from the asymptotic formula (3.14), the contribution from the viscous boundary flow is not included. This exclusion may also contribute to small but noticeable differences between the asymptotic and numerical values. It is also interesting to look at the spatial structure of the geostrophic flow that is depicted in figure 2*a*,*c* calculated from the asymptotic solution (3.13) and in figure 2*b*,*d* calculated from fuller numerics. Not surprisingly, the asymptotic and numerical solutions, obtained for *E*=10^{−4} and *E*=10^{−5} in figure 2, unveil nearly the same profiles without significantly noticeable differences.

## 5. Summary and remarks

This paper reports a new Legendre-type polynomial that is discovered in the asymptotic analysis of geostrophic flows in spherical geometry for *E*≪1 and *ϵ*≪1. The new polynomial is a function of two variables but takes a single argument and is orthogonal over the full sphere. This paper also presents an application of the new polynomial in deriving an asymptotic solution for forced geostrophic flows in rapidly rotating fluid spheres. Moreover, numerical analysis of the same problem is also carried out, showing a satisfactory agreement between the asymptotic solution and the numerical solution.

It should be pointed out that, although equations (1.1) and (1.2) can be solved numerically by expanding **u** in terms of the Legendre polynomial and the standard Chebyshev function *T*_{n}(*r*), there exists a profound practical difference between the numerical and asymptotic solutions. An extremely small value of *E* in the Earth’s liquid core, a consequence of the Earth’s rapid rotation and the small viscosity (Gubbins & Roberts 1987), causes severe difficulties in the numerical solution of core geostrophic flow. Even with modern powerful supercomputers, most numerical solutions are usually restricted to the range of moderately small Ekman number *E*≥10^{−5}, which is far from the geophysically realistic value *E*≤10^{−9}. Our analytical expression (3.13), which is validated by numerical analysis for moderately small *E*, is valid for an arbitrarily small Ekman number and, thus, can be employed in modelling geophysically realistic flow (e.g. Jault *et al.* 1988).

We have only discussed one possible application of the geostrophic polynomials *G*_{2k−1}. It is anticipated that the finding of these new orthogonal polynomials will open an exciting new line of research for their application to many geophysical and astrophysical fluid dynamical problems in rotating spherical systems. For example, the polynomials , together with a complete set of spherical inertial modes (Zhang *et al.* 2001), may be employed to solve the asymptotic problem of fully nonlinear convection in rapidly rotating fluid spheres. Unlike the spherical harmonics that are strongly coupled by the rotational differential operator , they would offer a more effective and efficient way to solve/understand the problem of nonlinear convection in rapidly rotating spherical systems with an asymptotically small *E* that is unreachable in the usual numerical analysis.

## Acknowledgements

X.L. is supported by NSFC/10633030, MOSTC863/2006AA01A125, STCSM 08XD14052 and CAS KJCX2-YW-T13 grants and K.Z. is supported by UK NERC and STFC grants. We would like to acknowledge constructive comments from an anonymous referee and from Prof. Andy Jackson who points out a possible connection between the new geostrophic polynomial discussed in this paper and the Jacobi polynomial family.

## Footnotes

- Received November 7, 2009.
- Accepted January 15, 2010.

- © 2010 The Royal Society