The present numerical study aims at shedding light on the mechanism underlying the precessional instability in a sphere. Precessional instabilities in the form of parametric resonance due to topographic coupling have been reported in a spheroidal geometry both analytically and numerically. We show that such parametric resonances can also develop in spherical geometry due to the conical shear layers driven by the Ekman pumping singularities at the critical latitudes. Scaling considerations lead to a stability criterion of the form |Po| > O(E4/5), where Po represents the Poincaré number and E represents the Ekman number. The predicted threshold is consistent with our numerical simulations as well as previous experimental results. When the precessional forcing is supercriticial, our simulations show evidence of an inverse cascade, i.e., small scale flows merging into large scale cyclones with a retrograde drift. Finally, it is shown that this instability mechanism may be relevant to precessing celestial bodies such as the earth and earth’s moon.

Precession corresponds to the gyroscopic motion of a rotating object due to a torque orthogonal to its spin axis. When an internal liquid layer is present, such as the liquid core or the subsurface ocean of a planet or the fuel tank of a spinning spacecraft, energy can be dissipated in the liquid due to the induced flows. In line with this idea, it has been proposed that orbital perturbations, in particular precession and nutation, could be used to probe the interior of the planet1–3 and possibly generate magnetic fields through the so-called dynamo process.4–8 

It is well-established that the primary response in the bulk of a rapidly rotating fluid cavity subject to precession is of uniform vorticity, i.e., the fluid rotates along an axis tilted with respect to the mean rotation axis of the container.2,9–12 Although viscosity plays an important role mostly in the viscous boundary layer next to the solid wall of the cavity, it has been shown that singularities at the so-called critical latitudes excite conical shear layers in the interior that are superimposed on the uniform vorticity flow.13–16 In addition, weakly non-linear interactions in these singular regions of the boundary layer generate steady geostrophic shears,9 which were observed both numerically17–19 and experimentally.10,20–25

In a pioneering piece of experimental work, Malkus20 showed that precession driven flows can become unstable and even create space filling turbulence. These initial findings were re-confirmed by later experiments.10,21,23,26 A theoretical justification was first proposed by Kerswell,27 who argued that the precession driven uniform vorticity flow in a spheroidal cavity is inertially unstable due to the constant strain field exerted by the solid wall. Indeed, in a spheroid, the uniform vorticity flow is a superposition of a solid body rotation and a gradient flow, which is required to fulfil the non-penetration condition at the wall. This gradient flow can be decomposed into two components, one leading to elliptically deformed streamlines and the other to a shear of the centres of the streamlines. Both parts can interact with a pair of free inertial modes of the rotating cavity through parametric couplings. This mechanism has been confirmed numerically28 for oblate spheroids, while similar inertial instabilities were reported experimentally in cylindrical precessing tanks.29–31 

Both numerical simulations18,32,33 and laboratory experiments23,26 in spherical geometries have shown very rich dynamics ranging from laminar flows to fully developed turbulence. However, in the case of a spherical cavity, the uniform vorticity solution reduces to a purely solid body rotation, preventing the aforementioned inertial instability to develop. Hence, instabilities in a sphere can only originate from the viscous correction to the solid body rotation flow.34 

The present study aims at shedding light on the onset of the unstable flows in a precessing sphere. Our numerical results show that the shear in the conical structures spawned from the critical latitudes is key to the destabilization process, echoing the shear of the centres of the streamlines in the case of a spheroidal cavity. We conjecture that the conical shear layers couple with two free inertial modes leading to a parametric resonance. The known scalings of the conical shear layers allow us to derive a stability criterion for the onset of the unstable flow that is in agreement with our numerical simulations. Furthermore, we show that for large enough forcing, small scale vortices merge into large scale cyclonic vortices with a retrograde drift.

The paper is organized as follows. We first introduce the mathematical background and the numerical model in Sec. II, and then, we present our numerical results in Sec. III. Using heuristic arguments, we derive a scaling for the onset of the observed instability in Sec. IV. Finally, we discuss our findings in the context of planetary dynamics in Sec. V.

We consider a sphere of radius R filled with a homogeneous and incompressible fluid of density ρ and kinematic viscosity ν. The sphere rotates at Ω o = Ω o k ˆ and precesses at Ω p = Ω p k ˆ p , where k ˆ and k ˆ p are unit vectors along the spin and precession axes, respectively (Figure 1). Hereinafter, we refer to the mantle frame as the frame attached to the container, the precession frame as the frame rotating at Ωp and the fluid frame, the frame attached to the solid body rotation of the fluid. In this paper, all the numerical calculations are implemented in the mantle frame, while the precession frame and the fluid frame are used from time to time to discuss the results.

FIG. 1.

Sketch of the problem.

FIG. 1.

Sketch of the problem.

Close modal

Using the radius R as the length scale and Ω o 1 as the time scale, the dimensionless Navier-Stokes equations describing the fluid velocity u and pressure p in the mantle frame become35 

u t + u u + 2 ( k ˆ + P o k ˆ p ) × u = p + E 2 u P o ( k ˆ p × k ˆ ) × r ,
(1)
u = 0 ,
(2)

where the time-dependent precession vector k ˆ p is given by

k ˆ p = ( i ˆ cos t j ˆ sin t ) sin α p + k ˆ cos α p .
(3)

Here, ( i ˆ , j ˆ , k ˆ ) are unit vectors of Cartesian coordinates (x, y, z) whose z-axis is along the rotation vector k ˆ . In Eq. (1), the last term on the left hand side is the Coriolis force due to the rotation and precession, the last term on the right hand side is the so-called Poincaré acceleration, and p is the reduced pressure.

In the present study, the angle αp between the spin and precession axes is kept constant and equal to 60°. The evolution of the system is then governed by two dimensionless parameters, the Poincaré number Po and the Ekman number E,

P o = Ω p Ω o , E = ν Ω o R 2 ,
(4)

which measure the dimensionless rate of precession and the ratio of the typical viscous force and Coriolis force, respectively. Negative (positive) values of Po correspond to retrograde (prograde) precession; we consider only retrograde precession in the present study.

Equations (1) and (2), together with no slip boundary condition u = 0 on the wall, are solved using the fully spectral code developed by Marti.36 Using the so-called toroidal/poloidal decomposition in a spherical coordinate system (r, θ, ϕ)

u = × ( T r ) + × × ( P r ) ,
(5)

the velocity field is then represented by two scalar fields T and P. The incompressibility condition is automatically satisfied by such a decomposition. The scalar fields T and P are then expanded as

T ( r , θ , ϕ ) = n = 0 N l = 0 L m = l l t l , n m W n l ( r ) Y l m ( θ , ϕ ) ,
(6)
P ( r , θ , ϕ ) = n = 0 N l = 0 L m = l l p l , n m W n l ( r ) Y l m ( θ , ϕ ) ,
(7)

where Y l m ( θ , ϕ ) are the spherical harmonics of degree l and order m and W n l ( r ) are the so-called Worland polynomials. For a given harmonic degree l, the Worland polynomials are a combination of the rl factor and the Jacobi polynomials P n ( α , β ) ( x ) , i.e., W n l ( r ) = r l P n 1 / 2 , l 1 / 2 ( 2 r 2 1 ) , which exactly satisfy the parity and regularity at the origin of the sphere.37 

By using the decomposition from Eqs. (5)–(7), the radial component of the ∇ × Eq. (1) and ∇ × ∇ × Eq. (1) yield decoupled evolution equations for the spectral coefficients t l , n m and p l , n m . The time integration is implemented using a second order predictor-corrector scheme.

The numerical code has been widely benchmarked in several contexts including that of precession driven flows.32,38 We use a typical truncation of the spectral expansion up to N = 63, L = M = 127 at moderate Ekman numbers (E ⩾ 3 × 10−5). For a few calculations at lower Ekman numbers (E ⩽ 10−5), we truncated as high as N = 255, L = 511, and M = 255 to well resolve the thin boundary layer and instabilities.

In this subsection, we introduce some useful quantities derived from our numerical simulations performed in the mantle frame.

The solution to the linearized equations (1) and (2) is symmetric around the origin due to the parity of the precessional forcing32,35 (see also in Appendix  A), namely, u(−r) = − u(r), where r is the position vector. Any breaking of the symmetry must be caused by an instability. Hence, it is of interest to decompose the total velocity into its symmetric and antisymmetric parts, i.e., u = us + ua with

u s = u ( r ) u ( r ) 2 , u a = u ( r ) + u ( r ) 2 .
(8)

Such a decomposition can easily be performed in the spectral space by taking advantage of the parity of the spherical harmonics. Specifically, the symmetric part includes only odd degree l in the toroidal field T and even degree l in the poloidal field P, while the antisymmetric part is the other way round.

Following the symmetry decomposition, we define the total energy, symmetric energy, and antisymmetric energy as

E k = 1 2 V | u ( r ) | 2 d V , E k s = 1 2 V | u s ( r ) | 2 d V , E k a = 1 2 V | u a ( r ) | 2 d V .
(9)

The growth of an anti-symmetric energy component is a sufficient condition to identify the development of an instability. However, it should be noted that a symmetry-breaking is not a necessary condition, indeed instabilities with negligible anti-symmetric component have been reported in numerical studies.28 In all cases, the instabilities will lead to time variations of the energy in the system that is otherwise steady.

The main observations from previous experimental studies are reported in the precession frame;10,20,21 it is therefore natural to provide the reader with an equivalent point of view in our numerical simulations. To that end, when needed, after we perform our simulations in the mantle frame, velocities are simply transformed in the precession frame as

u p = u + z ˆ × r .
(10)

Anticipating the peculiar role of the conical shear layers which are coaxial with the rotation axis of the fluid,10,21 we extract the rotation vector of the fluid ωF from the mean vorticity in the bulk18 

2 ω F = < × u P > ,
(11)

where < ⋅ > denotes averaging in the fluid volume excluding a thin boundary layer (10E1/2).

Once ωF is determined, the spectral coefficients with respect to the system of coordinates aligned with the rotation axis of the container are converted into a new set associated with the coordinate system aligned with ωF, from which the velocities are reconstructed. This transformation is made using a Matlab subroutine39 based on the formula in Ref. 40. The angle αf between the rotation axes of the container and the fluid is defined as cos α f = k ˆ ω F . Finally, ε = | k ˆ ω F | represents the differential rotation between the container and the bulk of the fluid.

It is now well established that the base flow in a weakly precessing sphere is in the form of a solid body rotation along an axis inclined to the rotation axis of the container.9,16,35 Figure 2 represents the angle αf between the rotation axes of the container and the fluid as a function of the Poincaré number Po and Ekman number E. We observe a quantitative agreement between our simulations (symbols) and the theoretical predictions9 (solid lines) that have been thoroughly validated.10,18 We note in particular that the scaling αfE−1/2 at fixed Po expected from the resonance at |Po|cosαpE1/2 is well recovered.9 Figure 2 also includes data points for which an instability occurs (▴), showing that even in these cases the asymptotic theory of Busse9 remains valid at first order.

FIG. 2.

The angle αf (in radians) between rotation axes of the container and the fluid (a) as a function of the Poincaré number Po at fixed Ekman number E = 3.0 × 10−5 and (b) as a function of the Ekman number E at fixed Po = − 1.0 × 10−4. The solid lines represent Busse’s theory9 and symbols correspond to our numerical results. The open circles (∘) and the solid triangles (▴) represent stable and unstable flows, respectively.

FIG. 2.

The angle αf (in radians) between rotation axes of the container and the fluid (a) as a function of the Poincaré number Po at fixed Ekman number E = 3.0 × 10−5 and (b) as a function of the Ekman number E at fixed Po = − 1.0 × 10−4. The solid lines represent Busse’s theory9 and symbols correspond to our numerical results. The open circles (∘) and the solid triangles (▴) represent stable and unstable flows, respectively.

Close modal

A dominant feature of the viscous corrections to the solid body rotation flow is the conical shear layers spawned from the critical latitudes. These structures, coaxial with ωF, correspond to the so-called characteristic surfaces of the hyperbolic inertial wave equation, i.e., the unforced inviscid Navier-Stokes equations in the rotating frame. The evidence of such oblique shear layers has been reported in numerical and experimental studies,17,19,22 with a typical velocity of OE1/5) over a typical width of O(E1/5).14,16,41 A typical example at Po = − 1.0 × 10−4 and E = 1.0 × 10−6 is represented in Figure 3, showing the conical shear layers in the fluid frame. As we shall see later on in this paper, they play an important role in the destabilization mechanism of the flow.

FIG. 3.

Contours of the velocities in the fluid frame in the meridional plane across both ωF and Ωo at Po = − 1.0 × 10−4 and E = 1.0 × 10−6.

FIG. 3.

Contours of the velocities in the fluid frame in the meridional plane across both ωF and Ωo at Po = − 1.0 × 10−4 and E = 1.0 × 10−6.

Close modal

As we increase the precession rate at a given Ekman number E, the base flow becomes unstable. A typical example just above the threshold is presented in Figs. 4–6 for Po = − 7.0 × 10−3 and E = 3.0 × 10−5. Figure 4 represents the time evolution of the antisymmetric energy Eka in the mantle frame. After a transient stage, 0 < t/2π < 200, the antisymmetric energy Eka grows exponentially until saturation is reached at t/2π ∼ 1000.

FIG. 4.

Anti-symmetric energy Eka as a function of time in the mantle frame. Po = − 7.0 × 10−3, E = 3.0 × 10−5.

FIG. 4.

Anti-symmetric energy Eka as a function of time in the mantle frame. Po = − 7.0 × 10−3, E = 3.0 × 10−5.

Close modal
FIG. 5.

Contours of the antisymmetric velocities ura (a)-(c) and uθa (d)-(f) in the fluid frame at t/2π = 973, Po = − 7.0 × 10−3, and E = 3.0 × 10−5. The left column is on a spherical surface of r = 1 − 10E1/2 and the black dot denotes the rotation axis of the container. The middle column is in the meridional plane across both ωF and Ωo. The vertical dashed line in the meridional plane represents the cylinder associated with the critical latitude. The right column is in the equatorial plane perpendicular to the rotation axis of fluid. Two dashed circles are r = 0.6 and r = 0.8.

FIG. 5.

Contours of the antisymmetric velocities ura (a)-(c) and uθa (d)-(f) in the fluid frame at t/2π = 973, Po = − 7.0 × 10−3, and E = 3.0 × 10−5. The left column is on a spherical surface of r = 1 − 10E1/2 and the black dot denotes the rotation axis of the container. The middle column is in the meridional plane across both ωF and Ωo. The vertical dashed line in the meridional plane represents the cylinder associated with the critical latitude. The right column is in the equatorial plane perpendicular to the rotation axis of fluid. Two dashed circles are r = 0.6 and r = 0.8.

Close modal
FIG. 6.

Symmetric energy Eks (△) and antisymmetric energy Eka (∘) in individual mF modes for the same instant and parameters as in Figure 5.

FIG. 6.

Symmetric energy Eks (△) and antisymmetric energy Eka (∘) in individual mF modes for the same instant and parameters as in Figure 5.

Close modal

Figure 5 shows a snapshot (t/2π = 973, just before the saturation) of the antisymmetric velocities ura and uθa in the fluid frame, on a spherical shell at 10E1/2 below the surface (a) and (d), in a meridional cross section (b) and (e), and in the equatorial plane perpendicular to ωF (c) and (f). The antisymmetric flow is mostly confined between two cylinders (0.6 < s < 0.8) co-axial with ωF, which correspond to latitudes higher than the critical latitudes in the Ekman layer at ±30°. In the equatorial plane perpendicular to ωF, the antisymmetric velocities are dominated by mF = 17 for ura and mF = 18 for uθa, where mF is the azimuthal wavenumber with respect to ωF. This is consistent with the symmetry of the velocity components, indeed for an antisymmetric flow,

u r a ( r ) = u r a ( r ) , u θ a ( r ) = u θ a ( r ) , u ϕ a ( r ) = u ϕ a ( r ) .
(12)

Hence, only modes with odd mF contribute to ura and modes with even mF to uθa in the equatorial plane.

Figure 6 shows the contribution of individual mF modes to the symmetric and antisymmetric energy at the same instant and parameters as in Figure 5. The symmetric energy Eks is dominated by the mF = 1 component, whereas the antisymmetric energy is dominated by the mF = 17 and mF = 18 components. These observations in the spectrum are consistent with the observed velocities in Figure 5. The less significant peaks at higher wavenumbers may be attributed to secondary bifurcations.34 

In Figure 7(a), we extract time series of the velocities corresponding to mF = 17 and mF = 18 at a fixed position in the fluid frame during the growth phase of the antisymmetric energy. The associated Discrete Fourier Transform (DFT) is shown in Figure 7(b) where frequencies are normalized by |ωF|. The animation (Movie 1 in the supplementary material50) of the velocities in the equatorial plane indicates that the mF = 17 mode is prograde and the mF = 18 mode is retrograde. So the frequencies of the two modes are ω17/|ωF| = − 0.34 and ω18/|ωF| = 0.67, respectively, which satisfy ω18/|ωF| − ω17/|ωF| ≈ 1.0. We also see some small amplitude high frequency fluctuations in the time series which may be due to uncertainties in the determination of ωF.

FIG. 7.

(a) Time series of ur of mF = 17 mode (in blue) and uθ of mF = 18 mode (in red) at a fixed position (r = 0.7, θ = π/2, ϕ = 0) in the fluid frame and (b) the corresponding discrete Fourier transform. Po = − 7.0 × 10−3, E = 3.0 × 10−5.

FIG. 7.

(a) Time series of ur of mF = 17 mode (in blue) and uθ of mF = 18 mode (in red) at a fixed position (r = 0.7, θ = π/2, ϕ = 0) in the fluid frame and (b) the corresponding discrete Fourier transform. Po = − 7.0 × 10−3, E = 3.0 × 10−5.

Close modal

The observed azimuthal wavenumbers and frequencies are consistent with a parametric resonance mechanism similar to the shear instability described by Kerswell27 in precessing spheroids. Such an instability can arise when the background flow couples with two free inertial modes that satisfy the so-called parametric resonant conditions, ω2ω1 = 1.0, m2m1 = 1, and l2 = l1, where ω is the eigen-frequency, m is the azimuthal wavenumber, and l is the spherical harmonic degree of the inertial modes (see Appendix  C). To fully characterize the possible inertial modes interacting in our simulations, we calculate the eigen-frequencies and the velocity structures of the inviscid inertial modes in a sphere following the analytical approach in Ref. 42 (see also Appendix  C). Among all possible inertial modes, we only have to consider the antisymmetric ones (odd l) with azimuthal wavenumber m = 17 and m = 18 based on our observations. In Figure 8, we plot the eigen-frequencies ω17 + 1.0 and ω18 as a function of l for inertial modes with m = 17 (×) and m = 18 (∘). We note that there are several combinations nearly satisfying ω18ω17 = 1.0, i.e., collocated symbols in Figure 8. However, only one combination (in the dashed ellipse), with l = 27, is found to match the observed frequency (the dashed line).

FIG. 8.

Eigen-frequencies of inertial modes with m = 17 (×) and m = 18 (∘) as a function of modal degree l. The horizontal dashed line indicates the observed frequency ω18. A combination in the dashed ellipse indicates a pair of modes closely matching the observations.

FIG. 8.

Eigen-frequencies of inertial modes with m = 17 (×) and m = 18 (∘) as a function of modal degree l. The horizontal dashed line indicates the observed frequency ω18. A combination in the dashed ellipse indicates a pair of modes closely matching the observations.

Close modal

Figure 9 compares the velocity structures of the observed unstable modes (a) and (b) in the numerical simulations with the identified inviscid inertial modes (c) and (d) matching the observations. We observe that the inviscid inertial modes have the same structure in the unstable region 0.6 < s < 0.8 as in our numerical simulations. In the most outer region, the agreement is more qualitative as we do not observe significant flows in our simulations while the inviscid modes exhibit a well defined pattern with significant velocities. This discrepancy may be attributed to viscous effect that is more influential in this region.

FIG. 9.

Contours of velocities of the unstable modes in the numerical simulations (a) and (b) and possible inviscid inertial modes (c) and (d). Snapshots of ur of the mF = 17 mode (a) and uθ of the mF = 18 mode (b) are taken at the same instant as in Figure 5. In each plot, the left half is in the equatorial plane and the right half is in a meridional plane. The red (blue) contours represent positive (negative) values. Two dashed lines in the equatorial plane represent two circles with radius of 0.6 and 0.8.

FIG. 9.

Contours of velocities of the unstable modes in the numerical simulations (a) and (b) and possible inviscid inertial modes (c) and (d). Snapshots of ur of the mF = 17 mode (a) and uθ of the mF = 18 mode (b) are taken at the same instant as in Figure 5. In each plot, the left half is in the equatorial plane and the right half is in a meridional plane. The red (blue) contours represent positive (negative) values. Two dashed lines in the equatorial plane represent two circles with radius of 0.6 and 0.8.

Close modal

As we increase the precession rate well above threshold, the flow exhibits an inverse cascade of the vorticity. This is illustrated in Figure 10, which shows a sequence of snapshots of the axial vorticity in the equatorial plane and the axisymmetric azimuthal velocity in a meridional section in the fluid frame at Po = − 1.35 × 10−2 and E = 3.0 × 10−5. At the beginning (Figure 10(a)), we see small scale vortices similar to the flow observed during the growth of the parametric instability in Figure 5. The small scale structures then start to merge into large scale elongated cyclonic structures closer to the rotation axis of the fluid (Figures 10(b)10(d)). In addition, we observed a retrograde (westward) drift of the pattern, as seen in the animation (Movie 2 in the supplementary material50). These observations are consistent with an inverse cascade of energy characteristic of two dimensional turbulence. Although not yet formally established, the retrograde drift seems to result from the conservation of angular momentum.

FIG. 10.

(a)-(f) Snapshots of the axial vorticity in the equatorial plane and axisymmetric azimuthal velocity in a meridional section in the fluid frame. (g) Total kinetic energy Ek in the mantle frame as a function of time. Red dots correspond to snapshots in (a)-(f). Po = − 1.35 × 10−2, E = 3.0 × 10−5.

FIG. 10.

(a)-(f) Snapshots of the axial vorticity in the equatorial plane and axisymmetric azimuthal velocity in a meridional section in the fluid frame. (g) Total kinetic energy Ek in the mantle frame as a function of time. Red dots correspond to snapshots in (a)-(f). Po = − 1.35 × 10−2, E = 3.0 × 10−5.

Close modal

In the case of Figure 10, the large scale cyclones break down to small scales and form again (Figures 10(e) and 10(f)); in some cases, the cyclones can be sustained for more than 200 rotation periods of the container. Similar large scale vortices were also observed experimentally in a precessing cylinder43 and numerically in rotating Rayleigh Bénard convection,44 yet the underlying merging mechanism remains poorly understood.

Figure 11 shows the regime diagram in the (Po, E)-parameter space accessible in our numerical study. The flow is characterized as stable (open circles) when it is centrosymmetric and has a steady total kinetic energy after the transient stage. In contrast, symmetry-broken flows (filled triangles) or centrosymmetric flows but with time-varying total kinetic energy (open triangles) indicate the presence of an instability. We note that although the unstable flows are antisymmetric for the most part, we observe a few cases where the velocity field remains centrosymmetric, as for instance at Po = − 0.06 and E = 10−4. We also observe that the flow in a precessing sphere can be stable when the precession is sufficiently strong, as previously reported in laboratory experiments.23,26 In these cases, the fluid axis is nearly aligned with the precession vector.

FIG. 11.

Stability diagram in the plane of (Po, E). Circle and triangle symbols represent cases with steady and unsteady kinetic energy, respectively. Open symbols represent flows with centro-symmetry and filled symbols represent symmetry broken flows. The solid line represents the instability threshold from laboratory experiments in a precessing sphere.45 

FIG. 11.

Stability diagram in the plane of (Po, E). Circle and triangle symbols represent cases with steady and unsteady kinetic energy, respectively. Open symbols represent flows with centro-symmetry and filled symbols represent symmetry broken flows. The solid line represents the instability threshold from laboratory experiments in a precessing sphere.45 

Close modal

The solid line in Figure 11 represents the lower instability threshold of laboratory experiments in a precessing sphere.45 We can see that our numerical simulations are in good agreement with the experimental results particularly when E ≤ 10−4. Finally, both numerical and experimental results are consistent with a critical Poincaré number scaling as |Po|sinαp = O(E4/5) that we shall now derive on the basis of a parametric instability mechanism.

Our observations are consistent with a parametric instability mechanism similar to the shear instability described by Kerswell.27 However, in contrast with the spheroidal geometry, the shear cannot result from topographic effects. Instead, we argue that the conical shear layers spawned from the critical latitudes in the boundary layers can also induce a parametric instability. Indeed, the conical shear layers driven by the linear viscous interactions in the boundary layer are mF = 1 and ω/|ωF| = 1.0 in the fluid frame and steady in the precession frame. Thus, they satisfy the parametric resonant conditions with the observed modes in our simulations. Hereinafter, we will refer to this instability as a CSI, for Conical-Shear-Instability.

Figure 12(a) shows the streamlines in the frame of precession just before the growth of the unstable modes illustrating the distortion along the conical shear surfaces, here represented in light grey. Figure 12(b) shows the traces of the conical shear layers (left half) and velocity vectors (right half) in the meridional plane (Ωo, ω).

FIG. 12.

(a) Streamlines of the velocities in the precession frame at t/2π = 200. The gray line represents the rotation axis of the mantle and the black line represents the rotation axis of the fluid. Shadowed conical cones represent the characteristic surfaces with respect to the rotation axis of the fluid. (b) Contours of uθ (left half) and in-plane velocity vectors (right half) in the meridional plane across both the rotation axis of fluid and the rotation axis of the fluid at the same instant as in (a). Dashed lines are intersections of the characteristic surfaces and the meridional plane. Po = − 7.0 × 10−3, E = 3.0 × 10−5.

FIG. 12.

(a) Streamlines of the velocities in the precession frame at t/2π = 200. The gray line represents the rotation axis of the mantle and the black line represents the rotation axis of the fluid. Shadowed conical cones represent the characteristic surfaces with respect to the rotation axis of the fluid. (b) Contours of uθ (left half) and in-plane velocity vectors (right half) in the meridional plane across both the rotation axis of fluid and the rotation axis of the fluid at the same instant as in (a). Dashed lines are intersections of the characteristic surfaces and the meridional plane. Po = − 7.0 × 10−3, E = 3.0 × 10−5.

Close modal

Following the well established theory,46,47 we shall now derive the scaling for the onset condition. For an exact resonance, the growth rate of the parametric instability can be written as

σ = C 1 C 2 ( κ 1 + κ 2 ) 2 E 1 / 2 ,
(13)

where κ1 and κ2 are the real decay rate factors of two inertial modes. The interaction parameters between the conical shear ushear and two inertial modes u1, u2 are given by

C 1 = < u 1 , u 2 × ( × u shear ) > + < u 1 , u shear × ( × u 2 ) > < u 1 , u 1 > ,
(14)
C 2 = < u 2 , u 1 × ( × u shear ) > + < u 2 , u shear × ( × u 1 ) > < u 2 , u 2 > ,
(15)

with <A, B > = ∭VAB dV and A is the complex conjugate of A.

A complete derivation of the growth rate requires evaluation of the volume integrals in Eqs. (14) and (15). Although an explicit expression of ushear was derived recently,16 a set of partial differential equations have to be solved numerically to get ushear. Nevertheless, a scaling law for the onset of the parametric instability can be established using heuristic arguments. It has been demonstrated that the amplitude of the conical shear layers scales as |ushear| = OE1/5) over a width O(E1/5), where ε represents the differential rotation between the fluid and the surrounding solid shell.14,16,19 Therefore, we estimate |∇ × ushear| = OE1/5/E1/5) = O(ε), while the integration volume is proportional to E1/5. It follows

C 1 C 2 = O ( ε E 1 / 5 ) .
(16)

Combining Eqs. (16) and (13), we obtain the instability threshold (Re(σ) = 0)

ε E 1 / 5 ( κ 1 + κ 2 ) 2 E 1 / 2 ,
(17)

leading to

ε = O ( E 3 / 10 ) .
(18)

In a precessing sphere, the differential rotation ε is a function of Po, αp, and E. At low precession rate, i.e., |Po|sinαpE1/2, the direct resonance mechanism between the tilt-over mode and the precessional forcing leads to ε = O(|Po|sinαp/E1/2).9,16 This scaling is also confirmed by our numerical simulations in Figure 2. Thus, the lower bound for the threshold of the CSI in a sphere can be written as

| P o | sin α p = O ( E 4 / 5 ) .
(19)

This scaling is in quantitative agreement with our numerical simulations (Figure 11) and previous experimental results.45 

The conical shear layers can also be excited in precessing spheroidal cavities48 and thus the CSI should be induced as well. In this geometry, we must calculate the differential rotation ε following Busse9 which also depends on the ellipticity η of a spheroid and then apply the instability criterion of Eq. (18).

The present numerical study investigates the stability of precession-driven flows in a full sphere at moderate to low Ekman numbers. It is found that at low precession rate the flow is of uniform vorticity with a viscous correction superimposed. As the precession rate increases, the internal conical shear layers driven by the linear interactions in the boundary layer induce a parametric instability. The threshold conditions have been established using heuristic arguments leading to a scaling |Po|sinαp = O(E4/5) in quantitative agreement with both our numerical simulations and former experimental observations. At onset, the inertial modes involved in the destabilization mechanism concentrate the energy in an annular region between two cylinders coaxial with the rotation axis of the fluid (0.6 < s < 0.8). Above threshold, the flow evolves from its onset geometry with relatively large azimuthal wavenumbers to a simpler structure with few cyclonic vortices, closer to axis of rotation of the fluid and traveling in a westward direction.

Our results may allow us to shed light on the long standing ill-understood wave-like instabilities reported by Vanyo et al.21 almost 20 yr ago in a precessing spheroid with ellipticity η = 1/100. In one of their experiments, dye was injected next to the boundary layer of a precessing fluid cavity. Pictures taken over more than one hour of experimentation reveal travelling wave-like structures organized on a circular path concentric with the rotation axis of the fluid. Figure 13 compares the observations from dye injection by Vanyo et al.21 with the vorticity in our numerical simulations. We argue that the spiraling structure observed by Vanyo et al.21 can be explained as a pair of traveling inertial modes excited through a parametric resonance. In Vanyo’s experiments, which used a spheroid, three mechanisms may generate a parametric instability: the two mechanisms proposed by Kerswell27 and the conical shear layer driven one investigated in the present paper. To establish the dominant mechanism, we first calculate the differential rotation ε using Busse’s theory9 in the spheroid of Vanyo et al.,21 and then we use Eq. (18) for the CSI and the derivation of Kerswell27 for the two other mechanisms. The results summarized in Table I clearly show that the experimental setup of Vanyo et al.21 is more prone to CSI rather than the classical topography driven parametric instabilities of Kerswell.27 In the original experiment in precessing spheroid by Malkus,20 the wave-like instability was also reported (see Figure 2(b) of Ref. Ref. 20). Based on our estimates in Table I, both the topographic effect and the conical shear layers may lead to the parametric instability in the experiment by Malkus.20 

FIG. 13.

(a) Snapshot of total vorticity on the surface at r = 1 − 10E1/2 in the precession frame at the same instant as Figure 5. Three white dotted lines represent the latitudes of 0°, 30°, and 60° with respect to the rotation axis of the fluid. Po = − 7.0 × 10−3, E = 3.0 × 10−5. (b) Advection of dye injected in the boundary layer of a precessing spheroid. (Reprinted with permission from Vanyo et al., “Experiments on precessing flows in the Earth’s liquid core,” Geophys. J. Int. 121, 136–142 (1995). Copyright 1995 Oxford University Press.) η = 1/100, αp = 23.5°, Po = − 5.0 × 10−4, and E = 8.0 × 10−7.

FIG. 13.

(a) Snapshot of total vorticity on the surface at r = 1 − 10E1/2 in the precession frame at the same instant as Figure 5. Three white dotted lines represent the latitudes of 0°, 30°, and 60° with respect to the rotation axis of the fluid. Po = − 7.0 × 10−3, E = 3.0 × 10−5. (b) Advection of dye injected in the boundary layer of a precessing spheroid. (Reprinted with permission from Vanyo et al., “Experiments on precessing flows in the Earth’s liquid core,” Geophys. J. Int. 121, 136–142 (1995). Copyright 1995 Oxford University Press.) η = 1/100, αp = 23.5°, Po = − 5.0 × 10−4, and E = 8.0 × 10−7.

Close modal
TABLE I.

Comparison of different instability mechanisms in spheroids used in experiments by Vanyo et al.21 and Malkus20 and in planetary liquid cores. The growth rates of the first two mechanisms are given by Kerswell27 and the CSI discussed in this study. The prefactors of the scaling laws are assumed to be around O(1). The differential rotation δω is calculated using Busse’s theory9 for the spheroid.

Parameter Malkus Vanyo Earth Moon
Ellipticity  η  4.0 × 10−2  1.0 × 10−2  ∼2.5 × 10−3  ∼10−5 
Ekman number  E  1.0 × 10−5  8.0 × 10−7  ∼10−14  ∼10−12 
Differential rotation  ε  3.3 × 10−1  2.0 × 10−2  ∼1.7 × 10−5  ∼3.0 × 10−2 
  Inviscid growth rate       
Elliptical instability  O2η ∼4.4 × 10−3  ∼4.0 × 10−6  ∼7.2 × 10−13  ∼9.0 × 10−9 
Shear instability  Oη ∼1.3 × 10−2  ∼2.0 × 10−4  ∼4.3 × 10−8  ∼3 × 10−7 
CSI  OE1/5 ∼3.3 × 10−2  ∼1.2 × 10−3  ∼2.7 × 10−8  ∼1.2 × 10−4 
  Viscous decay rate       
Viscous damping  O(E1/2 ∼3.2 × 10−3  ∼8.9 × 10−4  ∼10−7  ∼10−6 
Parameter Malkus Vanyo Earth Moon
Ellipticity  η  4.0 × 10−2  1.0 × 10−2  ∼2.5 × 10−3  ∼10−5 
Ekman number  E  1.0 × 10−5  8.0 × 10−7  ∼10−14  ∼10−12 
Differential rotation  ε  3.3 × 10−1  2.0 × 10−2  ∼1.7 × 10−5  ∼3.0 × 10−2 
  Inviscid growth rate       
Elliptical instability  O2η ∼4.4 × 10−3  ∼4.0 × 10−6  ∼7.2 × 10−13  ∼9.0 × 10−9 
Shear instability  Oη ∼1.3 × 10−2  ∼2.0 × 10−4  ∼4.3 × 10−8  ∼3 × 10−7 
CSI  OE1/5 ∼3.3 × 10−2  ∼1.2 × 10−3  ∼2.7 × 10−8  ∼1.2 × 10−4 
  Viscous decay rate       
Viscous damping  O(E1/2 ∼3.2 × 10−3  ∼8.9 × 10−4  ∼10−7  ∼10−6 

Finally, we shall discuss our results in a planetary context. Using a hydrostatic model for the shape of the Lunar’s Core-Mantle boundary (CMB), the differential rotation ε due to the 18.6 yr precession has been estimated to be of order 3% of the mantle rotation.3,11 In line with this conclusion, Williams et al.3 argued that the resulting flow should be turbulent in order to account for the large dissipation inferred from the Lunar Laser Ranging (LLR) measurements. Yet, if one considers the two classical precessional instabilities of Kerswell27 due to the polar ellipticity, the inviscid growth rate remains smaller than the viscous decay rate leading to a stable flow. In contrast, considering the effects of the conical shear layers, we obtain an inviscid growth rate 100 times larger than the viscous decay rate assuming that the Ekman number is around 10−12 (Table I). These highly super critical conditions are likely to be associated with more complex flows than the one reported here at moderate Ekman number simulations, leading to significant dissipation as observed through inversion of the LLR time series.

In contrast with the moon, the earth’s polar flattening (η ≈ 1/400 on the CMB21) and the large period of precession (26 000 yr) result in a liquid core strongly coupled to the mantle. As a consequence, the differential rotation is small at present day, ε ≈ 1.7 × 10−5. As seen in Table I, all three mechanisms may lead to a stable precession-driven flow in the outer core assuming the Ekman number of the outer core is around 10−14. In addition, it has been shown that another set of conical shear layers emanate from the inner core boundary,17,41 with an amplitude of OE1/6) over a thickness of O(E1/3). In principle, these conical shear layers could also participate to a CSI. A stability criterion can then be derived using the same heuristic argument presented herein leading to ε = O(E1/3). Hence, CSI in the earth’s liquid core may occur if the Ekman number of the outer core is below 10−15.

In planetary condition, however, several other parameters have to be considered to draw a more reliable conclusion such as the effect of stratification and magnetic fields, as well as the interaction of precession driven flows with other sources of motions such as thermo-chemical convection.

We would like to acknowledge Andrew Jackson for very useful discussions and suggestions on this study. We thank Andreas Fichtner for computational assistance. Simulations were run on Swiss National Supercomputing Center (CSCS) under the Project No. s369. Figure 13(b) is reproduced from the Figure. 5(b) in Ref. 21, for which we thank the Oxford University Press and the original authors for permission to reprint it. Y.L. and J.N. are supported by ERC Grant No. 247303 (MFECE) at ETH Zurich. P.M. is supported by the National Science Foundation EAR CSEDI Grant No. #1067944.

The linear viscous solution of the governing Eqs. (1) and (2) is symmetric around the origin, namely, u(r) = − u(−r). Neglecting the nonlinear term in Eq. (1) and changing the sign of the position vector, we have

u ( r ) t + 2 ( k ˆ + P o k ˆ p ) × u ( r ) = p ( r ) + E 2 u ( r ) P o ( k ˆ p × k ˆ ) × r
(A1)

and

u ( r ) t + 2 ( k ˆ + P o k ˆ p ) × u ( r ) = p ( r ) + E 2 u ( r ) + P o ( k ˆ p × k ˆ ) × r .
(A2)

Adding and subtracting Eqs. (A1) and (A2), we obtain

u a t + 2 ( k ˆ + P o k ˆ p ) × u a = p a + E 2 u a
(A3)

and

u s t + 2 ( k ˆ + P o k ˆ p ) × u s = p s + E 2 u s P o ( k ˆ p × k ˆ ) × r ,
(A4)

where

u s = u ( r ) u ( r ) 2 , p s = p ( r ) + p ( r ) 2
(A5)

and

u a = u ( r ) + u ( r ) 2 , p a = p ( r ) p ( r ) 2
(A6)

are symmetric and antisymmetric parts, respectively. We can see that only the symmetric solution is forced by precession.

In order to provide different views of precession driven base flow as suggested by one referee, we show the velocities in the mantle frame and the precession frame in Fig. 14 for Po = − 1.0 × 10−4 and E = 1.0 × 10−6. In the mantle frame, the base flow is a combination of the Poincaré mode, i.e., a solid body rotation about an axis in the equatorial plane, and the secondary flow due to the viscous correction. The secondary flow is dominated by the conical shear layers spawned from the critical latitudes, which is hidden behind the solid body rotation. Note that the Poincaré mode is a travelling inertial mode in the mantle frame. In the view of the precession frame, we see both global rotation of the container and the solid body rotation around an axis in the equatorial plane (the Poincaré mode), leading to a tilted rotation axis of the fluid with respect to the rotation axis of the container. The flow is steady in the precession frame. Again, the secondary flow due to the viscous correction is hidden beneath the solid body rotation.

FIG. 14.

Contours of the velocities in the meridional plane across both Ωo and Ωp in the mantle frame (a)-(c) and the precession frame (d)-(f) at Po = − 1.0 × 10−4 and E = 1.0 × 10−6.

FIG. 14.

Contours of the velocities in the meridional plane across both Ωo and Ωp in the mantle frame (a)-(c) and the precession frame (d)-(f) at Po = − 1.0 × 10−4 and E = 1.0 × 10−6.

Close modal

Inertial modes are solutions of the so-called Poincaré equation which can be written as42 

1 s p s + 1 s 2 2 p s 2 + 1 s 2 2 p ϕ 2 ( 4 ω 2 ω 2 ) 2 p z 2 = 0
(C1)

in the cylindrical coordinates (s, ϕ, z) with the no-penetration boundary condition

s p s + 2 i ω p ϕ + ( 1 4 ω 2 z p z ) = 0 ,
(C2)

on the spherical surface s2 + z2 = 1. Solutions can be found using separation of variables and they are

p k l m ( s , ϕ , z ) = s m z γ e i m ϕ j = 1 N ( x j 2 ( x j 2 1 ) + x j 2 ( 1 ω k l m 2 4 ) s 2 + ω k l m 2 4 ( 1 x j 2 ) z 2 ) ,
(C3)

where γ = 0 if (lm) is even and γ = 1 if (lm) is odd, and xj are the N = 1 2 ( l m γ ) zeros of the Legendre polynomial P l m ( x ) of degree l and order m. The eigen-frequency ωklm is the kth solution of the transcendental equation

2 ( 1 ω 2 4 ) d d ω P l m ( ω 2 ) = m P l m ( ω 2 ) .
(C4)

For each given integer of l and m, there are lm solutions if m ≠ 0 and l − 1 solutions if m = 0. Each eigen-frequency ωklm and the associated eigenfunction are identified by three indexes (k, l, m). Here, we only consider m ⩾ 0 and the solution can be extended to m < 0 by the relationship ω(k, l, − m) = − ω(k, l, m).

Finally, the velocity field of the each eigenmode is obtained from

u k l m = i ω ( 4 ω 2 ) [ 4 ( z ˆ p k l m ) z ˆ ω 2 p k l m 2 i ω z ˆ × p k l m ] .
(C5)

The explicit expression of velocities are given by Zhang et al.49 

1.
W.
Thomson
, “
On the thermodynamic acceleration of the earth’s rotation
,”
Proc. R. Soc. Edinburgh
11
,
396
405
(
1882
).
2.
H.
Poincaré
, “
Sur la précession des corps déformables
,”
Bull. Astron.
27
,
321
356
(
1910
).
3.
J. G.
Williams
,
D. H.
Boggs
,
C. F.
Yoder
,
J. Todd
Ratcliff
, and
J. O.
Dickey
, “
Lunar rotational dissipation in solid body and molten core
,”
J. Geophys. Res.: Planets
106
, 27933, doi:10.1029/2000JE001396 (2001).
4.
E. C.
Bullard
, “
The magnetic field within the earth
,”
Proc. R. Soc. A
197
,
433
453
(
1949
).
5.
A.
Tilgner
, “
Precession driven dynamos
,”
Phys. Fluids
17
,
034104
(
2005
).
6.
A.
Tilgner
, “
Kinematic dynamos with precession driven flow in a sphere
,”
Geophys. Astrophys. Fluid Dyn.
101
,
1
9
(
2007
).
7.
C. C.
Wu
and
P. H.
Roberts
, “
On a dynamo driven by topographic precession
,”
Geophys. Astrophys. Fluid Dyn.
103
,
467
501
(
2009
).
8.
C. A.
Dwyer
,
D. J.
Stevenson
, and
F.
Nimmo
, “
A long-lived lunar dynamo driven by continuous mechanical stirring
,”
Nature
479
, 212–4 (
2011
).
9.
F. H.
Busse
, “
Steady fluid flow in a precessing spheroidal shell
,”
J. Fluid Mech.
33
,
739
751
(
1968
).
10.
J.
Noir
,
P.
Cardin
,
D.
Jault
, and
J. P.
Masson
, “
Experimental evidence of non-linear resonance effects between retrograde precession and the tilt-over mode within a spheroid
,”
Geophys. J. Int.
154
,
407
416
(
2003
).
11.
J.
Noir
and
D.
Cébron
, “
Precession-driven flows in non-axisymmetric ellipsoids
,”
J. Fluid Mech.
737
,
412
439
(
2013
).
12.
K.
Zhang
,
K. H.
Chan
, and
X.
Liao
, “
On precessing flow in an oblate spheroid of arbitrary eccentricity
,”
J. Fluid Mech.
743
,
358
384
(
2014
).
13.
H.
Bondi
and
R. A.
Lyttleton
, “
On the dynamical theory of the rotation of the earth. II. The effect of precession on the motion of the liquid core
,”
Math. Proc. Cambridge Philos. Soc.
49
,
498
(
1953
).
14.
K.
Stewartson
and
P. H.
Roberts
, “
On the motion of liquid in a spheroidal cavity of a precessing rigid body
,”
J. Fluid Mech.
17
,
1
20
(
1963
).
15.
R. R.
Kerswell
, “
On the internal shear layers spawned by the critical regions in oscillatory Ekman boundary layers
,”
J. Fluid Mech.
298
,
311
325
(
1995
).
16.
S.
Kida
, “
Steady flow in a rapidly rotating sphere with weak precession
,”
J. Fluid Mech.
680
,
150
193
(
2011
).
17.
R.
Hollerbach
and
R. R.
Kerswell
, “
Oscillatory internal shear layers in rotating and precessing flows
,”
J. Fluid Mech.
298
,
327
339
(
1995
).
18.
A.
Tilgner
and
F. H.
Busse
, “
Fluid flows in precessing spherical shells
,”
J. Fluid Mech.
426
,
387
396
(
2001
).
19.
J.
Noir
,
D.
Jault
, and
P.
Cardin
, “
Numerical study of the motions within a slowly precessing sphere at low Ekman number
,”
J. Fluid Mech.
437
,
283
299
(
2001
).
20.
W. V.
Malkus
, “
Precession of the earth as the cause of geomagnetism: Experiments lend support to the proposal that precessional torques drive the earth’s dynamo
,”
Science
160
,
259
264
(
1968
).
21.
J.
Vanyo
,
P.
Wilde
,
P.
Cardin
, and
P.
Olson
, “
Experiments on precessing flows in the Earth’s liquid core
,”
Geophys. J. Int.
121
,
136
142
(
1995
).
22.
J.
Noir
,
D.
Brito
,
K.
Aldridge
, and
P.
Cardin
, “
Experimental evidence of inertial waves in a precessing spheroidal cavity
,”
Geophys. Res. Lett.
28
, 3785–3788, doi:10.1029/2001GL012956 (2001).
23.
S.
Goto
,
N.
Ishii
,
S.
Kida
, and
M.
Nishioka
, “
Turbulence generator using a precessing sphere
,”
Phys. Fluids
19
,
061705
(
2007
).
24.
S. A.
Triana
,
D. S.
Zimmerman
, and
D. P.
Lathrop
, “
Precessional states in a laboratory model of the Earth’s core
,”
J. Geophys. Res.: Solid Earth
117
, B04103, doi:10.1029/2011jb009014 (2012).
25.
J.
Boisson
,
D.
Cébron
,
F.
Moisy
, and
P.-P.
Cortet
, “
Earth rotation prevents exact solid-body rotation of fluids in the laboratory
,”
Europhys. Lett.
98
,
59002
(
2012
).
26.
S.
Goto
,
A.
Matsunaga
,
M.
Fujiwara
,
M.
Nishioka
,
S.
Kida
,
M.
Yamato
, and
S.
Tsuda
, “
Turbulence driven by precession in spherical and slightly elongated spheroidal cavities
,”
Phys. Fluids
26
,
055107
(
2014
).
27.
R. R.
Kerswell
, “
The instability of precessing flow
,”
Geophys. Astrophys. Fluid Dyn.
72
,
107
144
(
1993
).
28.
S.
Lorenzani
and
A.
Tilgner
, “
Inertial instabilities of fluid flow in precessing spheroidal shells
,”
J. Fluid Mech.
492
,
363
379
(
2003
).
29.
R.
Lagrange
,
C.
Eloy
,
F.
Nadal
, and
P.
Meunier
, “
Instability of a fluid inside a precessing cylinder
,”
Phys. Fluids
20
,
081701
(
2008
).
30.
R.
Lagrange
,
P.
Meunier
,
F.
Nadal
, and
C.
Eloy
, “
Precessional instability of a fluid cylinder
,”
J. Fluid Mech.
666
,
104
145
(
2011
).
31.
Y.
Lin
,
J.
Noir
, and
A.
Jackson
, “
Experimental study of fluid flows in a precessing cylindrical annulus
,”
Phys. Fluids
26
,
046604
(
2014
).
32.
R.
Hollerbach
,
C.
Nore
,
P.
Marti
,
S.
Vantieghem
,
F.
Luddens
, and
J.
Léorat
, “
Parity-breaking flows in precessing spherical containers
,”
Phys. Rev. E
87
,
053020
(
2013
).
33.
X.
Wei
and
A.
Tilgner
, “
Stratified precessional flow in spherical geometry
,”
J. Fluid Mech.
718
,
R2
(
2013
).
34.
S.
Lorenzani
and
A.
Tilgner
, “
Fluid instabilities in precessing spheroidal cavities
,”
J. Fluid Mech.
447
,
111
128
(
2001
).
35.
K.
Zhang
,
K. H.
Chan
, and
X.
Liao
, “
On fluid flows in precessing spheres in the mantle frame of reference
,”
Phys. Fluids
22
,
116604
(
2010
).
36.
P.
Marti
, “
Convection and boundary driven flows in a sphere
,” Ph.D. thesis (
ETH Zurich
,
2012
).
37.
P. W.
Livermore
,
C. A.
Jones
, and
S. J.
Worland
, “
Spectral radial basis functions for full sphere computations
,”
J. Comput. Phys.
227
,
1209
1224
(
2007
).
38.
P.
Marti
,
N.
Schaeffer
,
R.
Hollerbach
,
D.
Cebron
,
C.
Nore
,
F.
Luddens
,
J.-L.
Guermond
,
J.
Aubert
,
S.
Takehiro
,
Y.
Sasaki
,
Y.-Y.
Hayashi
,
R.
Simitev
,
F.
Busse
,
S.
Vantieghem
, and
A.
Jackson
, “
Full sphere hydrodynamic and dynamo benchmarks
,”
Geophys. J. Int.
197
,
119
134
(
2014
).
39.
F. J.
Simons
and
A.
Plattner
, “
Scalar and vector slepian functions, spherical signal estimation and spectral analysis
,” in
Handbook of Geomathematics
, edited by
W.
Freeden
,
M. Z.
Nashed
, and
Thomas
Sonar
(
Springer
,
Berlin, Heidelberg
,
2014
), pp.
1
42
.
40.
F. A.
Dahlen
and
J.
Tromp
,
Theoretical Global Seismology
(
Princeton University Press
,
1998
).
41.
R. R.
Kerswell
, “
Upper bounds on the energy dissipation in turbulent precession
,”
J. Fluid Mech.
321
,
335
370
(
1996
).
42.
H. P.
Greenspan
,
The Theory of Rotating Fluids
(
Cambridge University Press
,
London
,
1968
).
43.
W.
Mouhali
,
T.
Lehner
,
J.
Léorat
, and
R.
Vitry
, “
Evidence for a cyclonic regime in a precessing cylindrical container
,”
Exp. Fluids
53
,
1693
1700
(
2012
).
44.
A. M.
Rubio
,
K.
Julien
,
E.
Knobloch
, and
J. B.
Weiss
, “
Upscale energy transfer in three-dimensional rapidly rotating turbulent convection
,”
Phys. Rev. Lett.
112
,
144501
(
2014
).
45.
S.
Kida
, “
Instability by weak precession of the flow in a rotating sphere
,”
Procedia IUTAM
7
,
183
192
(
2013
).
46.
R. R.
Kerswell
, “
Elliptical instability
,”
Annu. Rev. Fluid Mech.
34
,
83
113
(
2002
).
47.
L.
Lacaze
,
P. L.
Gal
, and
S. L.
Dizès
, “
Elliptical instability in a rotating spheroid
,”
J. Fluid. Mech.
505
,
1
22
(
2004
).
48.
A.
Tilgner
, “
Non-axisymmetric shear layers in precessing fluid ellipsoidal shells
,”
Geophys. J. Int.
136
,
629
636
(
1999
).
49.
K.
Zhang
,
P.
Earnshaw
,
X.
Liao
, and
F. H.
Busse
, “
On inertial waves in a rotating fluid sphere
,”
J. Fluid. Mech.
437
,
103
119
(
2001
).
50.
See supplementary material at http://dx.doi.org/10.1063/1.4916234 for Movie 1 and Movie 2. Movie 1 shows velocities of two unstable modes in the fluid frame in the equatorial plane for the case in Fig.7. Movie 2 shows axial vorticity in the equatorial plane in the fluid frame for the case in Fig.10.

Supplementary Material