Characteristics of a precessing flow under the influence of a convecting temperature field in a spheroidal shell

We present results from direct numerical simulations of flows in spherical and oblate spheroidal shells, driven both by precession and thermal convection, with Ekman number $Ek=10^{-4}$, non-diffusive Rayleigh numbers from $Ra=0.1$ to $Ra=10$ and unity Prandtl number. The applied precessional forcing spans seven orders of magnitude. Our experiments show a clear transition between a convective state and a precessing flow that can be approximated by a reduced dynamical model. The change in the flow is apparent in visualizations and a decomposition of the velocity into symmetric and antisymmetric components. For the flow dominated by precession, some parameter combinations show two stable solutions that can be realized by a hysteresis or a strong thermal forcing. An increase of the Rayleigh number at a constant precession rate exhibits established scaling properties for the heat transfer, with exponents $2/7$ and $6/5$.


Introduction
It is well known that the interior of the Earth (Dziewonski & Anderson 1981) and its moon (Weber et al. 2011) consist, in large part, of a liquid layer consisting of a molten iron alloy. This scenario is also likely for other terrestrial planets. Movements in this electrically conducting fluid are assumed to be the source of planetary magnetic fields (Stevenson 2003). However, the energy source of the flow is not entirely certain. While models considering thermal and chemical convection had many successes in explaining key features of Earth's field, such as the frequent reversals and dominantly bipolar structures (Christensen & Wicht 2015), several studies from different fields of the geosciences have shed doubt on the assumption (Dwyer, Stevenson & Nimmo 2011;Olson 2013;Andrault et al. 2016;Fuller 2017). Mechanical driving mechanisms have been proposed as a viable alternative, especially the precession of the rotation axis, which has been explored in numerous studies.
A fluid in an uniformly rotating container would settle to move with the boundaries as a constant vorticity flow when no other forcings are present. Precession describes the slower rotation of the diurnal rotation axis around another axis orthogonal to the ecliptic. When precession acts on the fluid, it tries to keep its original motion due to inertia, but the forcing constantly changes its direction. Viscous and topographic torques exerted by the boundaries work to align the flow with the time-dependent rotation of the boundaries. The base precessional bulk flow with the velocity u in a sphere at position r and time t therefore consists of a fluid motion that is similar to the rotation of a rigid body with a fluid rotation axis ω f : u(r, t) = ω f (t) × r. The best known theoretical determination of ω f was made by Busse (1968), who considered the advection of a viscous boundary layer into the non-viscous interior to derive an implicit equation for the fluid rotation. This result was reobtained by Noir et al. (2003) with a torque balance approach, which was later generalized by Noir & Cébron (2013) to derive a differential equation describing the behaviour of ω f . As a general trend, the fluid rotation lags behind the rotation of the container, a behaviour that has been observed in both laboratory experiments and numerical studies, e.g. Vanyo & Dunn (2000), Tilgner & Busse (2001), Ernst-Hullermann, Harder & Hansen (2013) and Noir et al. (2003). Secondary flow structures and various viscous and inertial instabilities can develop in addition to this laminar flow, such as inertial modes, shear layers, parametric or triadic resonances and the development of turbulence. We refer to Le Bars (2016) for a recent review on the subject. Several studies have shown, with the help of direct numerical simulations, that a planetary dynamo driven by precession is possible (Tilgner 2005;Ernst-Hullermann et al. 2013;Lin et al. 2016;Cébron et al. 2019).
As mentioned in the first paragraph, many studies consider buoyancy as the main agent behind core flows. A buoyancy driven flow in a rotating shell, near the onset of convection, is characterized by two-dimensional structures that are aligned with the axis of rotation. In addition to these small scale convective columns, larger scale zonal flows develop (Dormy et al. 2004;Aurnou 2007). It has been hypothesized by Cheng et al. (2015) that these structures are not stable at parameters relevant for planetary cores. Of special interest for the study of thermal convection is the amount of heat that can be transported in addition to the basic conductive profile. Scaling laws connect the heat flux to the thermal forcing and other parameters, where different scalings characterize different flow regimes (Aurnou 2007;King et al. 2010). A great advantage of such scalings is that they can be extrapolated to regimes that are (currently) unreachable by either laboratory or numerical experiments. We refer to the recent review by Plumley & Julien (2019) for an overview on this subject with a focus on plane layer convection.
In this work, we combine a buoyancy forcing with a precessing rotation of the boundaries. A similar model was employed by Wei & Tilgner (2013), who considered a fixed, uniform background stratification in a spherical shell with a small stress-free inner core. They found that a stable stratification suppresses possible precessional instabilities, whereas an unstable stratification can lead to both stabilization and destabilization, depending on the parameters. This approach was later extended to the hydromagnetic case in Wei (2016), where it was shown that dynamos with a combined forcing are possible even when one forcing by itself leads to a failed dynamo. Here, we consider a non-uniform temperature stratification that is imposed by the temperature boundary conditions and we neglect the magnetic field. Our simulations have a larger inner core with no-slip boundaries. Most importantly, we look at flows in spheroidal shells, were theoretical models (Busse 1968; Precessing flow under the influence of a temperature field 891 A15-3 2013; Cébron 2015) have shown that the assumption of a linear, rigid-rotation-like bulk flow leads to multiple stable solutions at certain parameters, particularly low Ekman numbers or strong deformation of the boundary. This behaviour has been recreated in laboratory (Malkus 1968) and numerical experiments (Lorenzani & Tilgner 2003;Goepfert & Tilgner 2016;Vormann & Hansen 2018). We conduct a range of numerical simulations at a moderate Ekman number Ek = 10 −4 and show that this bistable behaviour persists when a temperature field is present. A jump to the second stable solution can not only be realized through a hysteresis, as for the purely mechanical forcing (Vormann & Hansen 2018), but also by the influence of thermal convection. Further, we show which of the two flow types described above dominates dependent on the parameters and point out how they can be distinguished by decomposing the flow field into symmetric and asymmetric components.
We start by introducing the mathematical model of a precessing, fluid-filled spheroid in § 2, where some space is dedicated to concisely discussing the direction of gravity in a spheroidal shell. In § 3, we summarize the torque balance model by Noir & Cébron (2013) that will be used for comparison to our numerical experiments. The results section is in two parts, where we first discuss the case of a varying precessional forcing at a constant Rayleigh number ( § 5.1) and then a smaller set of simulations where the Rayleigh number was increased at a constant precession rate ( § 5.2). The final § 6 summarizes and contextualizes the results.

Mathematical model of a fluid filled spheroid
Our model consists of a spherical or oblate spheroidal shell with outer radius r o and inner radius r i in the horizontal plane. In Cartesian coordinates, the boundaries are defined by where a and c are the horizontal major and vertical minor axes (c/a 1). The shell rotates with the frequency Ω d around the unit vectorẑ parallel to the minor axis and precesses with Ω p around another axisk p that forms an angle α withẑ. We search for the velocity u = (u, v, w), pressure p and temperature T of a fluid that is described by the Boussinesq approximation, where density variations are only retained for the buoyancy term and depend linearly on the temperature difference. In a rotating reference frame where the boundaries are at rest (mantle frame), the relevant equations are the Navier-Stokes equation the incompressibility condition ∇ · u = 0, (2.3) and the heat equation In (2.2), the fourth term on the right-hand side, the Poincaré acceleration, includes the time-dependent precession axis (2.5) The equations have been non-dimensionalized with the diurnal rotation period Ω −1 d , the distance between the outer and inner core boundary d = r o − r i and the temperature difference between inner and outer boundary T = T i − T o . With this scaling, the problem is governed by the Ekman number Ek = ν/(d 2 Ω d ), giving the ratio of viscous to rotational forces, the Poincaré number Po = Ω p /Ω d , which controls the relative strength of precession, the non-diffusive rotational Rayleigh number Ra = γ g 0 T/(Ω 2 d r o ) (with the expansivity γ and gravitation g 0 at r o ), controlling the strength of buoyancy and the Prandtl number Pr = ν/κ, which gives the ratio of viscosity ν to thermal diffusivity κ. A negative Po < 0 denotes retrograde precession, which is the case relevant for planets. The Rayleigh number can also be interpreted as the squared free-fall convective Rossby number, calculated from the ratio of a rotational to a free-fall time scale (Julien et al. 1996). We use no-slip boundary conditions (u = 0) for the velocity field and set the temperature to T i = 1 at the inner and T o = 0 at the outer boundary of the shell. Internal heat sources are not considered.
For the case of a spherical shell, the non-dimensional gravity vector in (2.2) is proportional to the radius: g = r. Due to the broken symmetry, the formulation is more involved in a spheroidal shell. For simplicity, we assume that the inner part of the shell has the same density as the fluid filled volume. We then use the derivation by Hvoždara & Kohút (2012), where the gravity field on the inside of an oblate spheroidal body is calculated using fundamental solutions of the Laplace and Poisson equations in oblate spheroidal coordinates. These coordinates (η, θ , φ) are connected to Cartesian coordinates (x, y, z) via (2.6) y = f cosh η sin θ sin φ, (2.7) z = f cosh η cos θ , (2.8) with the geometry parameter f = √ a 2 − c 2 . The solution also requires Lame's coefficient h η = h θ = f cosh 2 η − sin 2 θ, the functions and the multiplication constant where η 0 is the constant value of η on the outer boundary. With the gravitational constant G and density ρ 0 , the dimensional solution for the inside of the shell in spheroidal coordinates in a two-dimensional (2-D) plane (g φ = 0) is then g η = − 2πGρ 0 f 2 h η cosh η sinh η(sin 2 θ + E 2 P 2 (cos θ )), (2.13) Precessing flow under the influence of a temperature field 891 A15-5  which can be transferred to Cartesian coordinates in the xz plane via g x = (g η sinh η sin θ + g θ cosh η cos θ )(cosh 2 η − sin 2 θ ) −1/2 , (2.15) g z = (g η cosh η cos θ − g θ sinh η sin θ)(cosh 2 η − sin 2 θ ) −1/2 sgn (z). (2.16) For positions outside the 2-D plane (g y = 0), the solution can be obtained by vector rotation aroundẑ, since the gravity field is symmetrical. For easier comparison to the spherical case, we normalize (and non-dimensionalize) the spheroidal gravity field so that both cases have the same value at the outer boundary in the horizontal plane ( x 2 + y 2 = r o , z = 0). Since the base gravity enters into Ra, we should remember that the values are not directly comparable. Figure 1 shows an example of the resulting field in a 2-D plane for a strong flattening of c/a = 0.6. The streamlines show that gravity is not directed directly at the centre of the mass distribution (at (x, y, z) = (0, 0, 0)). Instead, the tangential streamlines exhibit a slightly curved path. When looking at the magnitude of the non-dimensional gravity field, we observe that the field is strongest at the north and south poles of the spheroid, were the distance to the centre of mass is smaller than at the boundary in the horizontal plane. Additionally, contour lines of gravity are not parallel to the inner or outer boundary of the volume. These differences from the spherical case become less significant as c/a moves closer to 1. Of course, in a real planet, the shape is influenced by pressure, rotation and gravity, which is a complex problem (Tricarico 2014). Our simplified model can be interpreted as a planetary body which shape was frozen in during an initial period of fast rotation and now has a rigid inner core and mantle. Studies of the early Earth suggest rotation rates up to ten times as high as those of today (Agnor, Canup & Levison 1999;Ćuk & Stewart 2012). For the case of strong rotation (low Ek) one may also consider the inclusion of centrifugal effects, which we do not attempt here. We mainly consider a strong deformation of c/a = 0.8 to reach a parameter region with possible bistability at acceptable computational cost -at smaller deformations, a much lower Ek would be required.

Dynamical model of the fluid rotation axis
In a purely precession driven flow (Ra = 0 in (2.2)), the main component of the flow outside of viscous boundary layers is often described as a linear rigid rotation around a fluid rotation axis ω f = (ω x , ω y , ω z ): u(r, t) = ω f (t) × r. Busse (1968) derived an implicit equation for the components of ω f based on the asymptotic advection of a viscous boundary layer flow to the non-viscous interior. Later, Noir et al. (2003) derived the same result starting from a torque balance formulation of the Navier-Stokes equation. Both approaches assumed a time-independent, linearized equation (∂u/∂t + u · ∇u = 0) and only considered the equatorial component of the viscous torque acting in the boundary layer, leading to a spin over mode. In their appendix A, Noir & Cébron (2013) present an extended model that does not restrict the time dependency of the flow and allows for a spin up or spin down of the fluid from axial differential rotation. We use this model as a point of comparison for our numerical simulations and therefore concisely repeat the result. In a precessing frame of reference, where the precession axis Ω p = (Ω x , 0, Ω z ) = Po(sin α, 0, cos α) is at rest, their dynamical model reads The generalized viscous torque is with the spin up decay factor (derived from the calculation of Greenspan & Howard (1963) for an axisymmetric container) Here, Γ (z) is the Eulerian gamma function and F(a, b; c; z) the hypergeometric series (Abramowitz & Stegun 1972). For the spin over decay factors λ r so and λ i so , we use the result from Zhang, Liao & Earnshaw (2004) and correct it for the presence of an inner core by multiplying with (3.5) The classical model by Busse (1968) as well as the dynamical model (3.1) by Noir & Cébron (2013) allow multiple solutions for some parameter combinations, especially at low Ek or smaller c/a. This behaviour was further examined by Cébron (2015), who gave approximate formulations for the range of parameters with several stable solutions. For the purpose of our study, we look for solutions to (3.1) by numerically searching for fixed points: starting with a fluid resting in the mantle frame of reference (ω f = (0, 0, 1) in the precessing frame) at Po = 0, we integrate the equations forward in time with a Runge-Kutta fourth-order time stepping scheme until a constant solution is reached. This solution is then used as a starting condition for the next larger value of |Po|, up to a maximum value. In order to find parameter combinations with multiple solutions, |Po| is then again gradually decreased while taking the final ω f as the starting condition for the next smaller |Po|. For comparison to the results from the full numerical simulations, we do not consider the three components of ω f separately but instead look at the angular kinetic energy density in the mantle frame resulting from the rotation of the fluid about ω f : The substraction ofẑ corrects for the rotation of the boundaries in the precessing frame of reference. This formulation ignores contributions from viscous boundary layers of an approximate thickness 1.4 √ Ek (Lorenzani & Tilgner 2001). In our numerical simulations, the total kinetic energy in the boundary layers is only approximately 1 % of the energy in the fluid bulk for all studied cases, leading us to neglect this inaccuracy.

Numerical method
We perform the direct numerical simulations with the spectral element code Nek5000 developed at Argonne National Laboratory (nek5000.mcs.anl.gov), in which the variables are represented in a weak formulation by high-order Lagrange polynomials defined on the Gauss-Lobatto-Legendre points in each element of a grid. The method combines high accuracy with geometric flexibility and very good parallel scaling capabilities and has been used for fluid mechanics problems both with mechanical forcings (Favier et al. 2015;Lemasquerier et al. 2017;Reddy, Favier & Bars 2018;Vormann & Hansen 2018) and convection (Scheel & Schumacher 2014. The equations are integrated forward in time with a third-order accurate backward differentiation formula. A splitting method is used, where the viscous and pressure steps are handled as sub-problems with a Jacobi preconditioned conjugate gradient solver and an additive overlapping Schwarz method. We refer to Fischer (1997), Deville, Fischer & Mund (2002), Fischer & Lottes (2005) and Karniadakis & Sherwin (2013) for details on the method.
The equations are solved on a cubed spheroidal grid as in Vormann & Hansen (2018). We project a Cartesian grid with points (x i ,ỹ i ,z i ) from each surface of a cube in the radial directionr onto a spherical shell. The corner points of the spectral elements are placed at (x ij , y ij , z ij ) = r jri , where r j is a list of positions between r i and r o placed according to the Gauss-Legendre-Lobatto distribution, refining the grid in the radial direction at the outer and inner boundaries to better resolve viscous boundary layers. The resulting grid is then compressed in the z-direction by a factor c/a to create a spheroid. The simulations reported here use a grid with 6144 elements and polynomial orders 6-12.

Results from direct numerical simulations
As mentioned in the introduction, this section on the results of our numerical experiments is split into two parts. All simulations have Ek = 10 −4 , a relative inner core size of 0.35 and a precession angle α = 23.5 • . We first discuss simulations in different geometries (c/a = 1, 0.96 and 0.8) at fixed Rayleigh numbers Ra = 0.1 and Ra = 1 ( § 5.1). For c/a = 0.96, we also used a computationally more demanding Ra = 10. In a spherical shell, Ra = 0.1 is approximately nine times overcritical, which is probably similar for a spheroidal shell. We started with the non-precessing case (Po = 0) of rotating convection and then gradually increased |Po| at constant Ra. In the second set of experiments ( § 5.2), we considered a spheroidal shell with c/a = 0.8 at a constant Po = −0.075 and, starting from a purely precessing flow (Ra = 0), increased the thermal forcing stepwise up to Ra = 5. For all simulations, we examined retrograde precession with Po < 0.
The important diagnostic parameters we will present are the kinetic energy density and the Nusselt number with the surface normal n for the outer surface S o . The value for the integral in the denominator of (5.2) for each geometry was obtained numerically in a simulation with Ra = 0 and Po = 0, since analytical solutions for spheroidal geometries are not available.
Studies of precession driven flows often decompose the velocity field into symmetric and antisymmetric components, since the basic precessional flow u = ω f × r is symmetric with respect to reflection at the origin of a Cartesian coordinate system u(r) = −u(−r). The velocity field u is therefore separated into a symmetric component u s = 1/2(u(r) − u(−r)) and an antisymmetric component u a = 1/2(u(r) + u(−r)) from which corresponding symmetric and antisymmetric kinetic energies are derived An important diagnostic parameter is then the relative antisymmetric energy An increase in E rel is often considered as a sign of instability of the base precessional flow. We ought to remember that the concept of antisymmetric energy is not as meaningful for a convection driven flow, where symmetry is not expected from the basic equations due to the buoyancy term. Still, we will see that the decomposition is useful to analyse our results. For future reference, we have compiled most of the relevant data plotted in the figures in the results section in tables 1-7 and table 8, presented in the Appendix.

Increasing the precessional forcing
The kinetic energy density (equation (5.1) . We do not show the full range of values for solutions of (3.1) to keep the figures compact. The kinetic energy density E kin is not the same for different geometries but decreases with c/a, i.e. the spherical case shows the largest absolute and relative kinetic energies (but remember that the definition of Ra depends on g 0 ).
Approximately at the value for Po were the base kinetic energy becomes smaller than the rotational energy (3.6) predicted by (3.1) (which does not take buoyancy into account), at around Po = −10 −3 to −10 −2 , the kinetic energy of the flow increases. Since this transition occurs at smaller precessional forcing for lower Ra, the kinetic energy increase happens for smaller |Po| here. When we raise |Po| further, we find the energy to be in good accordance with (3.6) for all values of Ra and c/a, with a tendency to be slightly larger for a stronger thermal forcing. In figure 5, we see three exemplary visualizations of the velocity and temperature fields. For Po = 0 and small Po, the flows are similar to typical patterns of rotating convection, with long structures in the velocity field parallel to the rotation axis. At an increased precessional forcing, a switch to a cylindrical, rigid-rotation-like flow structure occurs. For larger Ra, the situation is very similar, but smaller structures are formed. The purely precession driven flow admits two stable solutions around Po = −0.1 for c/a = 0.8, depending on the starting field, as is predicted in (3.1), which does not include buoyancy effects. Our numerical experiments with a temperature field show a similar behaviour, as is evident in figure 4, where arrows illustrate the use For the case of a larger thermal forcing, the flow directly switches to the branch of increased kinetic energy that is reached for lower Ra only when coming from a larger precessional forcing, i.e. the bistability is not realized by the starting condition but by the additional buoyancy term. When further decreasing the precessional forcing, the behaviour is analogous to Ra = 0.1, only the transition occurs earlier. Figure 6 shows two exemplary visualizations of flows in the bistable region, which are very similar in appearance to the figures discussed above.      up to the transition between the base thermal and rigid rotation solutions for the kinetic energy, as we would expect for a flow with no inherent symmetry. When reaching the Poincaré number beyond which E kin is predicted by the solution to the dynamic system (3.1), we observe a sharp drop in E rel by approximately two to three orders of magnitude. For lower Ra, the decrease in E rel is more pronounced: we can roughly estimate from the limited data that the decrease is inversely proportional to Ra. For example, in the spheroidal shell with c/a = 0.96, E rel reduces to approximately 2 × 10 −1 at Ra = 10, to 3 × 10 −2 at Ra = 1 and goes down to E rel = 3 × 10 −3 at Ra = 0.1. The plots to the right in figures 8-10 further explore the behaviour of the two flow components: there, E a and E s are shown separately as functions of Po. The antisymmetric energy is larger for larger Ra and remains approximately independent of Po, only a slight increase for larger precessional forcing is observed. The exception here is E a for the hysteresis loop at Ra = 0.1, were it increases by approximately one order of magnitude and becomes nearly as large as for Ra = 1. The symmetric part E s of the energy increases by several orders of magnitude for all values of Ra when reaching the value of Po for which the overall kinetic energy E kin also increases. For large precessional forcings, E s becomes nearly independent of the Rayleigh number, while the differences in E a remain. An increase in Ra by one order of magnitude seems to produce the same increase in E a , explaining the differences in E rel . Generally, the heat transfer represented by Nu is very similar between different geometries, as can be seen in figure 11. We observe that the heat flux is slightly lowered by the introduction of precession for Ra = 0.1 and 1, with a greater decrease for a larger Ra = 10. A clear increase of Nu above the base value is only observed for larger |Po| at c/a = 0.8 and Ra = 0.1. There, Nu stays at the higher value when the forcing is decreased again to investigate the hysteresis behaviour and then returns to its original value. 2)) for this case. While we observe oscillations of approximately ±10 % around the mean value, no long-term positive or negative trend is visible and no notable outliers occur. Note that in cases with a stronger precessional forcing, as in figure 7, the oscillations are much smaller after the initial transient.

Varying thermal forcing
In this section, we look at a set of numerical simulations were the Poincaré number was held constant at Po = −0.075 (c/a = 0.8 and Ek = 10 −4 ) and the Rayleigh number was gradually increased starting from a pure precessional flow at Ra = 0 up to Ra = 5. The value Po = −0.075 lies inside a bistable region, as was shown in the previous section. We started with a solution on the lower branch of E kin , i.e. |Po| was increased to reach the starting field. Note that, where possible, the figures also include data from simulations with Po = 0 and the respective Ra for comparison. We performed one additional experiment at Ra = 5 and Po = 0 as a comparison for the highest thermal forcing. Simulations with increasing Ra at Po = 0, where we searched for Nu > 1, indicate a critical Rayleigh number of Ra c ≈ 4.75 × 10 −3 . Compared to the non-precessing flow, our values for Ra are therefore 5 to 1000 times overcritical. For comparison, the critical Ra for c/a = 0.96 and c/a = 1 is approximately 5 × 10 −3 . Figure 13 shows results for the kinetic energy density E kin (5.1) and Nusselt number Nu (5.2). The kinetic energy shows that the base value for Ra = 0 is slightly below the value of E kin ≈ 0.024 for a solution of the dynamical model (3.1). When increasing Ra, the energy increases but stays close to this value. A fit to the data is relatively difficult due to the small range of values for E kin . Up to the transitional Rayleigh number Ra t = Ek −7/4 Ek 2 Pr −1 (1 − r i /r o ) −1 ≈ 0.15 (King et al. 2010), a quadratic function E kin = 0.0381Ra 2 + 0.0187 describes the data well, while a linear function E kin = 0.0068Ra + 0.0193 fits for larger Ra > Ra t . The energies for the same value of Ra at Po = 0 (see figure 4) are smaller by orders of magnitude for small Ra, approximately 6.5 × 10 −3 (Ra = 1) and 2.5 × 10 −4 (Ra = 0.1). For Ra = 5, the kinetic energy for the precessing flow is approximately 28 % larger than for the non-precessing experiment. The Nusselt number also increases from a base value of Nu = 1.014, starting with a steep increase that decreases as Ra becomes larger. We have Nu > 1 at Ra = 0 since the precessing flow by itself transports a small amount of heat. Here, the fit is clearer than for E kin . In the considered range of Ra, Nu can be fitted by functions Nu = 58Ra 6/5 − 1.06 for Ra < Ra t and Nu = 14.8Ra 2/7 + 1.1 for Ra > Ra t . We tried different exponents from the literature on (non-)rotating convection (Ahlers, Grossmann & Lohse 2009;Cheng et al. 2015;King et al. 2010) and found that these two exponents and the transitional Ra t provide the best explanation of the data. Here, the data for Po = 0 are very similar, as was already seen in figure 11, where we found Nu to be nearly independent of Po. Figure 14 shows how the symmetric and antisymmetric components of the flow behave as Ra is increased. The relative energy E rel changes by several orders of magnitude as Ra rises, from O(10 −8 ) up to O(10 −1 ), with the largest increase when going from Ra = 0 to Ra = 0.025. The separate plots of E a and E s show that this change is mainly due to an increase in E a . The value of E s shows a linear increase with Ra that is very similar to the overall behaviour of E kin (figure 13), while the antisymmetric component E a increases exponentially fast, although it stays below E s in the considered range. A fit of a linear function to the E s data gives E s = 0.0042Ra + 0.020, showing that the symmetric energy increases slower with Ra than the overall energy. The comparison with non-precessing experiments show that, there, both components are smaller but approximately equal,  figure 15, where we see that it remains largely unchanged as Ra increases. A minor trend is observable, where the magnitude of the vector components slightly increases at first and then decreases for the largest Ra. We further explore the behaviour of ω f with the four visualizations in figure 16. Just as in figure 5, they depict characteristic isosurfaces of the velocity magnitude and temperature. For Ra = 0, the precession flow is dominant and only leads to a small deformation of the conductive temperature profile. As the Rayleigh number increases, smaller structures emerge, but the linear rotation clearly dominates for Ra = 0.1 and 1. At the largest Ra = 5, the large roll structure is hardly visible on the small scales, though the measurements of ω f still indicate its presence.

Discussion
This study presented numerical experiments on convection in precessing spherical and spheroidal shells. We found that, at a constant Rayleigh number Ra, the flow remains largely independent of the precession over a wide range of values for Po and can be classified as rotating convection. A change occurs when the kinetic energy predicted by the dynamical system (3.1) by Noir & Cébron (2013) becomes larger than the base energy for Po = 0. After that, the flow can be described as precession dominated and is well described by the model (3.1), although a small influence of Ra remains. The transition is accompanied by a strong increase in the symmetric component of the velocity field, which becomes independent of Ra, while the antisymmetric component remains approximately constant. As was shown in Lorenzani & Tilgner (2001), the precessional forcing only drives the symmetric component of the flow, while buoyancy has no preference for either component. A convection dominated flow therefore results in E rel ≈ 0.5, while a precession dominated flow results in much smaller values. The bistability predicted by Noir & Cébron (2013) is reproduced in our simulations, again for a strong deformation of c/a = 0.8. Here, the resulting solution depended on the starting condition for a small thermal forcing Ra = 0.1. This behaviour is very similar to the purely precession driven flows reported in Vormann & Hansen (2018). For an increased thermal forcing of Ra = 1, the simulation switches to a different branch by itself, showing that an external influence can change a precessing flow in the bistable parameter range. We note that the values of E rel for the purely mechanically driven flow (Ra = 0) reported in Vormann & Hansen (2018) are several orders of magnitude smaller, even when the total kinetic energy is very similar. Here, and in the experiments discussed in the next paragraph, we only observe a minor influence of precession on the heat flux.
In the second numerical experiment with an increasing Ra, it was found that both E kin and Nu increase with Ra, with an approximately quadratic and linear dependence for E kin and a proportionality to Ra 6/5 and Ra 2/7 for Nu. Malkus (1954) first proposed an exponent of 1/3 for the non-rotating Rayleigh-Bénard system, while Shraiman & Siggia (1990) gave an exponent of 2/7 based on a study of the nesting of the thermal inside the viscous boundary layer. King et al. (2010) studied the heat transfer for dynamos in rotating shells with the help of direct numerical simulations and found scaling exponents of 6/5 (rapidly rotating regime) and 2/7 (weakly rotating regime), which validity regions are separated by a transitional Rayleigh number Ra t = Ek −7/4 Ek 2 Pr −1 (1 − r i /r o ) −1 (in our notation). A similar result was found for convection in a plane layer in King et al. (2009). They argue that the transition between the two regimes occurs when the size of the thermal and viscous boundary layers are equal. We note that the scaling in the rapidly rotating system also depends on the critical Rayleigh number, which in turn depends on the Ekman number. Since we only studied Ek = 10 −4 , we cannot make any further comments on this. King, Stellmach & Aurnou (2012) argued that the scaling may not hold at more extreme Ra and Ek. For both rotating and non-rotating convection in cylindrical containers (laboratory) and periodic Cartesian boxes (direct numerical simulations), Cheng et al. (2015) found similar relations of Nu ∝ Ra 0.322 and Nu ∝ Ra 1.29 . However, they studied a much wider range of Ek, Ra and Nu than our simulations and found other relations of the form Nu ∝ Ra β between the Rayleigh number and the heat flux, for example a steep scaling regime of Nu ∝ Ra 3.6 . At Ek = 10 −4 , they also found a larger exponent of 1.7, while the exponent 1.29, which is closest to our findings, appears at Ek = 10 −3 . For large enough Ra, the scaling converges towards the non-rotating case with exponent 0.322. We need to point out that a clear distinction between the exponents 2/7 and 1/3 for Ra > Ra t is not possible with our limited data, since the resulting fit is of similar quality for both cases. Still, it is encouraging to see that precessing convection seems to share basic features with rotating convection. Cheng et al. (2015) hypothesize a relation of the change in heat transfer regimes to the breakdown of coherent columnar structures. We have seen in our experiments that strong precessional forcing can also lead to the breakdown of such structures. At the parameter combinations shared between both sets of experiments, Po = −0.075 and Ra = 1 or Ra = 0.1, the diagnostic parameters differ by less than 1 %, which is well within the range of the fluctuations of the numerical solution. It is left to future studies to find out if a further increase in thermal forcing can trigger a switch to the branch of increased E kin for a precessional flow, as happened for Po = −0.1 and Ra = 1 when increasing the precessional frequency, although it seems unlikely since the bistability is predicted by a model that ignores buoyancy. Since the transitional Ra t as given above is dependent on Ek, studies with varying Ek and a measurement of the layer thickness are necessary to elucidate the heat transfer behaviour. Studies of rapidly rotating Rayleigh-Bénard convection at smaller Ek have shown a dependence of the heat transfer exponent on Ek. For a stronger rotation, an increased heat transfer occurs Plumley et al. 2016;Cheng et al. 2018). Since our results are similar at a moderate Ek, we might expect respective results for a precessing spheroid when approaching planetary relevant values. Cébron, Maubert & Le Bars (2010) also combined thermal and mechanical forcings in a numerical study of the tidal instability in an ellipsoidal shell. They focussed on the growth of the tidal instability and also measured the resulting heat flux. The results indicate a transitional Rayleigh number proportional to Ek −8/5 , based on a competition between the heat transfer by the tidal instability and natural convection. Similar arguments might be made for the precessional flow, though this requires information on the heat flow generated by precessional instabilities. The results of our simulations indicate that the modelled flow can be described as a competition between precession and convection, were established features of both flow types are present, such as the bistability of the precessing flow and the heat transport scaling of thermal convection. The dominance of either type of flow can be shown by studying the relative strengths of symmetric and antisymmetric components in the velocity field, were the symmetric component is well described by the model of precessing flow constructed by Noir & Cébron (2013). When this component is dominant, the precession roll component of the flow overlays the convective columns.

Implications for planetary flows
For a real planet, we can assume that the precessional forcing slowly declines due to tidal friction. If |Po| decreases below a certain value, the kinetic energy of the precessional flow becomes smaller than the value for E kin of a purely thermal flow at the core's Ra. At this point, we might expect the core flow to switch from a precessing to a convecting state and remain there as |Po| decreases further. If bistable solutions are possible, this change might occur very suddenly as the flow 'drops' to a convective state at the boundary of the bistable region. Of course, if the kinetic energy of the precessional flow is always larger than the convective flow (i.e. if Ra is too small), convective processes might not be very important at all. If, on the other hand, Ra is large enough, the convective flow can completely overlay the precessional solution. However, smaller contributions of the opposing flow forcings are of course possible. Figure 17 shows the kinetic energy density derived from the model by Busse (1968) for representative planetary parameters from the literature and c/a = 299/300 for varying retrograde precessional forcing (see the figure caption for details). We note that solutions for the dynamical system (3.1) at small Ek ≈ 10 −15 for the full range of Po are computationally expensive and have therefore been omitted. Tests for some selected values show that the results are similar enough for this approximate comparison (although we have observed that qualitative differences may occur at larger forcing). Horizontal lines show an approximate range for the kinetic energy density of the flow in Earth's outer core. For this estimate, we use approximate velocities for the core flow of 17 to 40 km a −1 , based on the values given for magnetic flux patches in Stefan et al. (2017). The kinetic energy predicted for the purely precessional flow for the Earth at Po ≈ −10 −7 is clearly below the estimates for the core flow, which might indicate that the core is not in a state dominated by precession. But, as mentioned in the introduction, other studies point in the opposite direction, so that our simple ad hoc model can certainly not give a definitive answer. Lower estimates for the core's Rossby number of 10 −7 -10 −6 , as given in Aurnou (2007) and Christensen & Wicht (2015), would actually place the prediction for the precessing flow above the flow energy for the core. We see that the Earth at Po = −10 −7 is not in a region of possible bistability, as are the other planets. Still, the possibility of multiple solutions cannot be easily refuted. The large relative precession rates at which the bistability occurs might be more likely to be realized in smaller satellites that are under the influence of a more massive planet, as for example |Po| for the Moon is slightly above the bistability range, while it is clearly below for the planets.

Outlook
Further work on this topic could include the examination of other geometries (ellipsoids, different inner core shapes, differential rotation) and an expansion to more realistic parameters (lower Ek and Pr). Apart from the bistability of the precessing flow, the spheroidal shape is also of interest for the study of convection. To our knowledge, only Evonuk (2015) considered the effect of an ellipsoidal deformation on the patterns of convection in a 2-D equatorial plane. While we did not focus on this topic, our model of a 3-D spheroid can act as a starting point for further research. Also note that we did not determine the critical Rayleigh number for the onset of convection as a function of the geometry and precession rate, but always took clearly overcritical Ra.
gratefully acknowledge the Gauss Centre for Supercomputing (GCS) for providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS share of the supercomputer JUQUEEN at Jülich Supercomputing Centre (JSC). GCS is the alliance of the three national supercomputing centres HLRS (Universität Stuttgart), JSC (Forschungszentrum Jülich) and LRZ (Bayerische Akademie der Wissenschaften), funded by the German Federal Ministry of Education and Research (BMBF) and the German State Ministries for Research of Baden-Württemberg (MWK), Bayern (StMWFK) and Nordrhein-Westfalen (MIWF).