The energy dissipation in a uniform vorticity flow, such as the flow in a precessing spheroid or the one associated with the earth’s free core nutation, is mainly confined to the boundary layers. However, the thinness of the boundary layer renders it difficult to study the energy dissipation in the turbulent regime, either in laboratory experiments or through direct numerical simulations. Here, we use a local Cartesian model to study the energy dissipation in the boundary layer of a precessing sphere when the flow becomes turbulent, contrasting it with the laminar case. We compute the evolution of the boundary layer over time at individual co-latitudes based on direct numerical simulations using the computational fluid dynamics solver Nek5000. We then estimate the total global dissipation by summing up individual contributions. A comparison with known analytical results in the laminar case validates this approach. We briefly discuss the applications to the lunar and the earth’s core cases.

## I. INTRODUCTION

Free Core Nutation (FCN) is a rotational normal mode of the earth. Since its frequency is close to forced nutation, the retrograde annual nutation, some properties of the FCN can be determined by measuring the resonant amplitude and the phase lag of the forced nutation [e.g., Herring (1986)]. Inversions from the observed phase lag placed an upper bound for the fluid core kinematic viscosity of 0.54 m^{2} s^{−1} (Gwinn , 1986), about six orders larger than the value inferred from studies of liquid metal [see the review of Lums and Aldridge (1991)]; a low viscosity of *ν* = 10^{−6} m^{2} s^{−1} is still favored by recent studies based on *ab initio* calculations [e.g., Pozzo (2013) and Ichikawa and Tsuchiya (2015)]. Turbulence in the core’s boundary layers can lead, in principle, to an enhanced eddy viscosity. However, Ohmic dissipation in the core and mantle is a more natural way to explain the observations [e.g., the studies by Buffett (1992), Mathews and Guo (2005), Deleplace and Cardin (2006), and Buffett and Christensen (2007)]. However, accounting for the observed phase lag requires either a strong magnetic field at the core-mantle boundary (CMB) or a highly conductive layer at the base of the mantle (Mathews , 2002; Buffett , 2002). Although the flow near the CMB associated with the FCN motion is too weak to be turbulent, the flow associated with precession might be strong enough to be turbulent (Buffett, 2021; Triana , 2021). Thus, it is important to properly quantify the turbulent contribution to the energy dissipation from precessional motion to obtain tighter bounds on the Ohmic dissipation.

Instabilities in the boundary layer of a precession-driven flow have been previously identified in numerical simulations (Lorenzani and Tilgner, 2001); however, the thinness of the boundary layer makes it difficult to study experimentally [see the review of Tilgner (2015)]. That same reason, combined with current computational constraints, precludes numerical studies from reaching regimes of parameters appropriate for planetary interiors where the Coriolis force plays a central role. This limitation is reflected in the smallness of the *Ekman number* *E*, representing the ratio of viscous to Coriolis forces and controlling the thickness *L*_{ν} of the boundary layer (also known as the Ekman layer), which scales as *L*_{ν} ∝ *E*^{1/2}. We can define a local Reynolds number *Re* to quantify the condition for turbulence in the boundary layer as *Re* = *UL*_{ν}/*ν*, where *ν* is the viscosity and *U* is the velocity. For the nominal value of the earth’s precession-driven flow, the Reynolds number is about *Re* = 500, higher than the critical value of *Re* ∼ 150 for the onset of turbulence proposed by Sous (2013). With a global spherical model, Cébron (2019) showed that the total viscous dissipation increases by a factor of 1.17 compared to the laminar solution when *Re* is around 230. Buffett (2021) restricted the modeling domain to the boundary layer and broke it into pieces according to the co-latitude *θ*. Using a local Cartesian model, similar to ours, he studied two representative cases at *θ* = 0 and 60° and showed that the *friction velocity* *u*_{*}, a quantity proportional to the square root of the surface shear stress (see the definition below), increases by a factor of 1.2 at *Re* = 500 compared to the laminar value. Buffett (2021) demonstrated a promising way to study the turbulent effects in the boundary layer from a local point of view. However, the connection to the global picture is still missing. Dissipation depends not only on the surface shear stress but also on the *veering angle* *β*, the angle between the stress and the fluid velocity on the boundary. In the present work, we numerically investigate the local flow at different co-latitudes and integrate the results to obtain global estimates of the possible turbulent energy dissipation.

Energy dissipation in the lunar core due to precession is a key parameter in understanding the possible dynamo action of the past (Stys and Dumberry, 2020; Zhang and Dumberry, 2021). The *Re* for the boundary layer flow is about 10^{4}, much larger than earth’s. Therefore, it has been believed that the boundary layer is turbulent [e.g., Yoder (1981)]. The model by Buffett (2021) is, in principle, not applicable to the moon not only because of the much higher Reynolds number but also because the flattening in their model must be much larger than the ratio of precession to rotation (i.e., the Poincaré number *Po*). However, in the lunar case, the latter condition is not met as the CMB flattening of 10^{−4} (Viswanathan , 2019) is smaller than the Poincaré number, *Po* = 4 × 10^{−3}. In order to extend the result achieved by Buffett (2021), we investigate the problem by considering the case where the flattening is zero, i.e., a precessing fluid in a spherical cavity.

The scaling of quantities in the boundary layer can be used to estimate the dissipation at large *Re*. The theory for the scaling of *u*_{*} and *β* in the Ekman layer was proposed by Csanady (1967). It was later improved by Spalart (1989) by including higher order terms. Scaling has been studied via laboratory experiments (Caldwell , 1972; Sous , 2013) and numerical simulations (Coleman , 1990; Coleman, 1999; Shingai and Kawamura, 2004; Deusebio , 2014; and Braun , 2020). However, the Ekman layer in those studies is steady, i.e., it remains constant all the time, contrary to the one from precession-driven flow, which oscillates with diurnal frequency. The study of a Stokes boundary layer, an oscillating boundary layer but without Coriolis force, reveals complex behavior during the transition to turbulence: an instability growing in one oscillation cycle may eventually decay at intermediate *Re* (Ozdemir , 2014). Besides, analytical solutions show that the *u*_{*} associated with an oscillatory Ekman layer is smaller than that of a steady Ekman layer by a factor of 2^{1/4} (Buffett, 2021). Therefore, we want to know if scaling is applicable to the oscillatory Ekman layer and establish the difference with those previous studies.

In Sec. II, we present the Navier–Stokes equation in the fluid frame using local Cartesian coordinates and compare it to the one used by Buffett (2021). Next, we outline the numerical implementation in Sec. III. Numerical solutions are presented in Sec. IV. Finally, the discussion and conclusions appear in Sec. V.

## II. PROBLEM FORMULATION

We consider a spherical cavity with radius *R* filled with an incompressible and viscous fluid, rotating with angular velocity **Ω**_{m} and precessing with angular velocity **Ω**_{p}. The angle between **Ω**_{m} and **Ω**_{p} is denoted as *α*, the precession angle. In the case of weak precession (Ω_{p} ≪ Ω_{m}), the leading order solution for the flow in the cavity may be expressed as a solid body rotation **Ω**_{f}, which represents the mean rotation vector of the fluid. Since the dominant forcing in the system is solely related to precession and the fluid viscosity *ν*, the orientation of **Ω**_{f} depends on their relative amplitudes. According to Busse’s theory (Busse, 1968), **Ω**_{f} can be determined, given the Poincare number *Po* = Ω_{p}/Ω_{m}, the Ekman number *E* = *ν*/(Ω_{m}*R*^{2}), and the precession angle *α*. We choose our reference frame as the one rotating with the fluid’s mean angular velocity, and we refer to it as the *fluid frame*. In such a frame, the bulk of the fluid remains at rest, and the cavity rotates with $\Delta \Omega = \Omega m \u2212 \Omega f $. A boundary layer exists, the main focus of our study, with a nominal thickness $ L \nu = \nu / \Omega f $.

**Ω**

_{f}, and the coordinate system co-rotates with the fluid. For the local coordinate system, the

*x*-axis is eastward, the

*y*-axis is northward, and the

*z*-axis is in the radial outward direction. The fluid domain extends from

*z*= 0 toward the negative

*z*-axis. In the fluid frame, the flow velocity

**satisfies the following Navier–Stokes equation:**

*u**P*is the reduced pressure including the centrifugal force and

**′ is the position vector referring to any point inside the cavity. The last term, $ r \u2032 \xd7 \Omega p \xd7 \Omega f $, is known as the Poincaré term [coined by Malkus (1968)] or the Euler term, which is associated with the fictitious force appearing in such a non-uniformly rotating reference frame. Since the flow is mainly confined to the viscous boundary layer, we take**

*r**L*

_{ν}as the length scale to analyze Eq. (1). We further assume that

*L*

_{ν}is much smaller than the cavity radius

*R*, so the cavity curvature is neglected in the local model. Together with the differential velocity

*U*=

*R*|Δ

**Ω**| as the velocity scale, Eq. (1) can be written in dimensionless form as

*Re*=

*UL*

_{ν}/

*ν*is the Reynolds number and $ R \u0302 =R/R$ is a unit vector. A superscript (*) denotes dimensionless variables (the variable

*** is not to be confused with the friction velocity**

*u**u*

_{*}introduced later). In the local coordinate system, the vector $ k \u0302 f $ is given by $cos\theta z \u0302 +sin\theta y \u0302 $, where

*θ*is the co-latitude.

**Ω**

_{m}and

**Ω**

_{f}in a coordinate system defined in the precessing frame, which rotates with

**Ω**

_{p}. In the precessing frame, both vectors are stationary based on the theory proposed by Busse (1968). The coordinate system, denoted as $ ( X \u0302 p , Y \u0302 p , Z \u0302 p ) $, is defined such that it coincides with $ ( X \u0302 , Y \u0302 , Z \u0302 ) $ at

*t*= 0; otherwise, the latter is rotating with

**Ω**

_{f}in the precessing frame. Moreover,

**Ω**

_{m}is tilted toward negative $ X \u0302 p $, forming an angle

*θ*

_{f}relative to

**Ω**

_{f}(see below for the definition of

*θ*

_{f}). Under the no spin-up condition (i.e., Ω

_{f}= Ω

_{m}sin

*θ*

_{f}), the differential rotation vector will be lying on the equatorial plane: $\Delta \Omega =\u2212\Delta \Omega X \u0302 p $. The coordinate is then transformed to $ ( X \u0302 , Y \u0302 , Z \u0302 ) $ in the fluid frame and further to $ ( x \u0302 , y \u0302 , z \u0302 ) $ for the computation of boundary velocity, which is Δ

**Ω**×

**. Explicit expressions of the transformations are given in Appendix A. In the local coordinate system, the boundary condition at the lower end (**

*R**z** → −

*∞*) is

*** = 0 (i.e., the bulk of the fluid sufficiently away from the boundary is at rest). Note that, in the numerical implementation discussed below, the lower end is set at**

*u**z** = −30. On the fluid surface (

*z** = 0, corresponding to the CMB),

*ϕ*is the longitude. Thanks to the azimuthal

*m*= 1 symmetry of the flow, it is anticipated that the local flow field at a given

*ϕ*

_{0}is equivalent to the flow at the prime meridian (

*ϕ*= 0) with a lag time

*ϕ*

_{0}/Ω

_{f}. For the numerical implementation, we set

*ϕ*= 0, and then the azimuthal dependence is retrieved directly from Ω

_{f}

*t*for post-processing analysis, that is, a time-average can be treated as an azimuthal average.

To determine **Ω**_{f}, we use the reduced model developed by Cébron (2019) [i.e., their Eqs. (5)–(7)]. Their solutions are obtained in the precessing frame where a Cartesian coordinate system $ ( X \u0302 p \u2032 , Y \u0302 p \u2032 , Z \u0302 p \u2032 ) $ is defined such that the $ Z \u0302 p \u2032 $-axis and **Ω**_{m} are aligned. The spherical coordinates of **Ω**_{f} are denoted as the co-latitude *θ*_{f} and the longitude *ϕ*_{f}. The coordinates are transformed to the local system by three steps: (1) rotation of the coordinate system to align the $ Z \u0302 p \u2032 $-axis with **Ω**_{f}, (2) transformation from the precessing frame to the fluid frame, and (3) transformation from the global system to the local system. Details about the transformation are given in Appendix A. So far, we are ready to solve Eq. (2) numerically.

^{−10}. The ratio of the amplitudes of the Poincaré term to the viscous term is 4 × 10

^{−6}, which is the value obtained in the turbulent regime where

*Re*= 1.67 × 10

^{4}. In the laminar regime, say

*Re*= 5, we find the ratio is 7 × 10

^{−3}, thus allowing us to drop the Poincaré term in the

*Re*range discussed in the present study. Second, we consider the magnitude of $ \Omega m / \Omega f $, which is equal to cos

^{−1}

*θ*

_{f}under the no spin-up condition. Given

*θ*

_{f}= 1.543° at

*Re*= 1.67 × 10

^{4}, the magnitude is very close to unity with a difference of 0.04%. With lower

*Re*(high

*ν*), the angle

*θ*

_{f}will be even smaller, thus making $ \Omega f / \Omega m $ closer to unity. Third, the condition of weak precession (

*Po*≪ 1) allows us to drop the precession vector in the Coriolis term. With these approximations, we arrive to the reduced form of Eq. (2),

*Re*= 500 and find no significant difference between the two solutions (see Fig. S1 in the supplementary material). In other words, neglecting the precessional terms makes no difference in our case. In the following, all the discussions will be based on Eq. (4). We note that Eq. (4) is generally identical to the one derived for a precessing spheroidal cavity [see the Appendix of the study by Buffett (2021)], even though the reasons to drop the Poincaré term are different. In a precessing spheroid, the Poincaré term is balanced by the terms associated with the flattening. Nevertheless, both formulations make the governing equation resemble a simple rotating fluid. The appearance of Eq. (4) may seem counter-intuitive indeed, but it is important to remember the initial assumption that the container is precessing and that the flow is mostly a uniform vorticity flow. The local model aims to examine the flow in a small rectangular box near the boundary. Precession comes into play when we determine

**Ω**

_{f}. The only difference to the study by Buffett (2021) lies in the selection of the rotation vector, where we employ

**Ω**

_{f}in contrast to

**Ω**

_{m}used by Buffet. Considering the aforementioned assumption, that is, Ω

_{f}/Ω

_{m}≈ 1, the discrepancy between the two vectors becomes negligible. Nevertheless, we think that adopting

**Ω**

_{f}is more appropriate and self-consistent as the local model requires the bulk fluid to be at rest, a condition that may not hold when employing

**Ω**

_{m}in the mantle frame. In Appendix B, we derive the analytical solution in the laminar regime where Re is small and the non-linear term can be neglected, which matches the numerical result, as expected.

Parameter . | Definition . | Units . | Value . |
---|---|---|---|

ν | m^{2} s^{−1} | 10^{−6} | |

Ω_{m} | rad s^{−1} | 2.66 × 10^{−6} | |

Ω_{f} | Ω_{m} cos θ_{f} | rad s^{−1} | 2.659 × 10^{−6} |

Ω_{p} | rad s^{−1} | 1.07 × 10^{−8} | |

R | km | 380 | |

L_{ν} | $ \nu / \Omega f $ | m | 0.61 |

α | deg | 1.543 | |

θ_{f} | deg | 1.543 | |

ϕ_{f} | deg | 0.060 | |

γ | deg | 0.002 | |

U | RΩ_{m} sin θ_{f} | m s^{−1} | 2.72 × 10^{−2} |

Po | Ω_{p}/Ω_{m} | 4 × 10^{−3} | |

E | ν/(R^{2}Ω_{m}) | 3 × 10^{−12} | |

Re | UL_{ν}/ν | 1.67 × 10^{4} |

Parameter . | Definition . | Units . | Value . |
---|---|---|---|

ν | m^{2} s^{−1} | 10^{−6} | |

Ω_{m} | rad s^{−1} | 2.66 × 10^{−6} | |

Ω_{f} | Ω_{m} cos θ_{f} | rad s^{−1} | 2.659 × 10^{−6} |

Ω_{p} | rad s^{−1} | 1.07 × 10^{−8} | |

R | km | 380 | |

L_{ν} | $ \nu / \Omega f $ | m | 0.61 |

α | deg | 1.543 | |

θ_{f} | deg | 1.543 | |

ϕ_{f} | deg | 0.060 | |

γ | deg | 0.002 | |

U | RΩ_{m} sin θ_{f} | m s^{−1} | 2.72 × 10^{−2} |

Po | Ω_{p}/Ω_{m} | 4 × 10^{−3} | |

E | ν/(R^{2}Ω_{m}) | 3 × 10^{−12} | |

Re | UL_{ν}/ν | 1.67 × 10^{4} |

**is equal to the increase (decrease) of internal and kinetic energies [see, e.g., Sec. 16 in the study by Landau and Lifschitz (2013)],**

*τ**S*is the surface area. The increase in internal energy is further linked to the energy dissipated in the system. Therefore, for an equilibrium system where the change in kinetic energy reaches zero, the viscous power can be used to infer the dissipation. Suppose that the area

*S*corresponds to a small patch on the cavity surface

*dS*= sin

*θ R*

^{2}

*dθdϕ*, we compute the total dissipation via $ D \nu =\u2212\u222b ( \tau \u22c5 u ) dS$. The local viscous power per unit area is defined as $ P visc =\tau \u22c5u$ for later use. We can integrate the laminar analytical solutions to obtain the global laminar dissipation as [e.g., in the study by Cébron (2019)]

*I*= (8

*π*/15)

*ρR*

^{5}is the moment of inertia of the fluid. This result is consistent with the torque approach, which is given in Appendix B. A script for the computation is included in the supplementary material.

*kρ*|

**|**

*u***[e.g., in the study by Williams (2001)], where**

*u**k*is a dimensionless parameter that depends on the viscosity. The total turbulent dissipation −(45

*π*/32)

*kI*ΔΩ

^{3}is obtained by integrating the negative viscous power over the spherical surface. Since dissipation is independent of the orientation of Δ

**Ω**, one can take $\Delta \Omega =\Delta \Omega Z \u0302 $ to simplify the computation, as performed by Yoder (1981). For an arbitrary orientation, we obtain the same coefficient, −(45

*π*/32), by approximating the integral as a Riemann sum. Note that the stress prescription assumes that the stress is aligned with the velocity, which is not true for most cases (see Fig. S10 for the veering angle at different values of

*θ*), so the turbulent torque derived from this definition should thus be treated with extra care. In order to determine the parameter

*k*, we adopt the approach by Cébron (2019) and require $ P visc $, computed from turbulent and laminar stresses (B17) and (B18), to match

*θ*= 0; it gives $k= ( u * / U ) 2 \u2061cos\beta $, where $ u * = | \tau | / \rho $ is the friction velocity and

*β*is the veering angle, the angle between the surface shear stress and the fluid velocity on the boundary. We notice that, when we substitute the laminar solutions $ u * = \nu U / L \nu $ and

*β*= 45° into

*k*and compute the total turbulent dissipation, the obtained dissipation disagrees with the laminar expectation given in Eq. (6). To reconcile the inconsistency, we define the turbulent dissipation as

*θ*. We choose

*θ*= 0 because the asymptotic similarity theory is used to estimate the friction velocity and veering angle at this special case, where the presence of a logarithmic layer supports the use of the theory. The purpose of this modification is to provide a way to compute

*k*and relate it to the total dissipation, which can be compared to our numerical results.

## III. NUMERICAL IMPLEMENTATION

We use the spectral element code Nek5000 (Fischer , 2007) to solve Eq. (4) for an incompressible fluid in a Cartesian box. The domain is divided into $E$ elements, where solutions are represented by *N*-th order Lagrange interpolation polynomials (see Fig. 2). The typical resolution in the present study is $E=768$ and *N* = 9. The temporal discretization is based on a third-order backward-difference scheme for implicit terms and on an extrapolation scheme for explicit terms. We use dealiasing following the 3/2 rule.

The domain is extended downward from *z* = 0. The non-dimensional size is 30 × 30 × 30. The top boundary is a wall (representing the CMB), where the velocity is specified to mimic the precessing container (mantle), as seen from the fluid frame [see Eq. (3)]. The four lateral boundaries are periodic. The bottom boundary is stress-free (i.e., no tangential stress), and the normal velocity is zero. We have performed a height test to make sure that the domain height of 30*L*_{ν} is large enough to accommodate the boundary layer flow (see Fig. S2 in the supplementary material). To resolve the boundary layer flow, we use a higher element density near the top boundary. We will explore the parameter space: *Re* = [5, 500] and *θ* = [0, 90°]. Simulations for *Re* = 700 are added to extend the *Re* range. We also simulate the case with *Re* = 600 at *θ* = 0. The range of *Re* covers both the laminar and turbulent regimes of the oscillatory Ekman layer.

The simulation always starts from the fluid at rest, where a random perturbation in velocity with root-mean-square amplitude of about 10^{−2} is imposed. In most cases, the flow will go through a transient state for a certain period and reach a *statistically steady state* (we shall refer to it as a steady state for the sake of brevity in the following sections), which means time-averaged quantities over the oscillation period(s) remain constant. The mean values presented below are obtained by taking time averages over oscillation periods. Due to the homogeneity on *x*- and *y*-directions, all the values are also averaged on the horizontal plane. The unit of time *t*_{⊕} presented below is 2*π*/Ω_{f}, indicating the period of the oscillating boundary. In what follows, the superscript (+) denotes the quantities non-dimensionalized by the viscous wall scales, velocity *u*_{*}, and length *ν*/*u*_{*}. The notation is chosen to describe the flow near the boundary.

## IV. RESULTS

We present now our numerical results, comparing them with analytical solutions, highlighting the difference between laminar and turbulent regimes. Then we use the numerical results in the turbulent regime together with the asymptotic theory to construct a turbulence model for the energy dissipation. Finally, we compute the total dissipation based on the local model and calibrate the turbulence model accordingly.

### A. Friction velocity and kinetic energy

The friction velocity *u*_{*} is a quantity commonly used to monitor the amplitude of surface shear stress |** τ**|: $ u * = | \tau | / \rho $ [e.g., in the study by Tennekes (1972)]. At

*θ*= 0, Eqs. (B19) and (B20) give the laminar solution $ u * lam = R e \u2212 1 / 2 U$, which serves as the benchmark for low Reynolds number cases. Figure 3 shows the ratio $ u * / u * lam $ as a function of time at

*θ*= 0 for different values of

*Re*. We see that the friction velocity deviates from its laminar solution when

*Re*≥ 300 and the deviation is larger for higher

*Re*. Concerning the case of

*Re*= 300, it appears that the deviation from the laminar solution by about 10% can be attributed to turbulence, as supported by the presence of a logarithmic region presented later in this study. We observe that the transient period happens before

*t*

_{⊕}= 3 (before

*t*

_{⊕}= 1 for

*Re*= 400 and higher) and the friction velocity reaches a steady value. The result is consistent with the one from the study by Buffett (2021), where the solution at

*Re*= 500 is about 20% larger than the laminar value after the onset of turbulence. To obtain the mean value, we take a time average over the interval 1 ≤

*t*

_{⊕}≤ 5, except for

*Re*= 300 where we use the 3 <

*t*

_{⊕}< 5 interval (see Table II).

Re
. | $ u * / u * lam $ . | $ u * /U$ . | $ \beta $ (deg) . | $ D \nu / D \nu lam $ . |
---|---|---|---|---|

5 | 1.00 ± 0.00 | 0.447 192 ± 0.000 419 | 45.002 058 ± 0.095 827 | 0.96 |

100 | 1.00 ± 0.00 | 0.099 973 ± 0.000 224 | 45.010 85 ± 0.229 41 | 0.96 |

200 | 1.00 ± 0.00 | 0.070 693 ± 0.000 158 | 45.010 364 ± 0.229 03 | 0.97 |

300 | 1.13 ± 0.02 | 0.064 989 ± 0.001 234 | 29.350 199 ± 1.461 621 | 1.34 |

400 | 1.21 ± 0.02 | 0.060 584 ± 0.001 197 | 25.635 712 ± 1.273 489 | 1.61 |

500 | 1.29 ± 0.03 | 0.057 866 ± 0.001 131 | 22.978 039 ± 1.135 976 | 1.86 |

600 | 1.37 ± 0.03 | 0.055 831 ± 0.001 122 | 21.037 396 ± 1.317 849 | ⋯ |

700 | 1.43 ± 0.03 | 0.054 069 ± 0.001 123 | 20.081 535 ± 1.411 917 | 2.42 |

Re
. | $ u * / u * lam $ . | $ u * /U$ . | $ \beta $ (deg) . | $ D \nu / D \nu lam $ . |
---|---|---|---|---|

5 | 1.00 ± 0.00 | 0.447 192 ± 0.000 419 | 45.002 058 ± 0.095 827 | 0.96 |

100 | 1.00 ± 0.00 | 0.099 973 ± 0.000 224 | 45.010 85 ± 0.229 41 | 0.96 |

200 | 1.00 ± 0.00 | 0.070 693 ± 0.000 158 | 45.010 364 ± 0.229 03 | 0.97 |

300 | 1.13 ± 0.02 | 0.064 989 ± 0.001 234 | 29.350 199 ± 1.461 621 | 1.34 |

400 | 1.21 ± 0.02 | 0.060 584 ± 0.001 197 | 25.635 712 ± 1.273 489 | 1.61 |

500 | 1.29 ± 0.03 | 0.057 866 ± 0.001 131 | 22.978 039 ± 1.135 976 | 1.86 |

600 | 1.37 ± 0.03 | 0.055 831 ± 0.001 122 | 21.037 396 ± 1.317 849 | ⋯ |

700 | 1.43 ± 0.03 | 0.054 069 ± 0.001 123 | 20.081 535 ± 1.411 917 | 2.42 |

Figure 4 shows the total kinetic energy, *E*_{kin}, for *Re* = 500 at different values of *θ*. We observe that the total kinetic energies reach a steady state for *most* cases. It implies that the time-averaged change in the kinetic energy is zero, $ E \u0307 kin =0$, and the system is in equilibrium state. According to Eq. (5), we can use the viscous power to approximate the dissipation. Although we can compute the dissipation directly from our model, the reason to use $ P visc $ is two-fold: (1) it is consistent with the conventional torque approach to compute the dissipation, and (2) we want to apply the turbulence model that is associated with the surface shear stress (i.e., *u*_{*}). A figure showing the time evolution of the three quantities ($ E \u0307 int $, $ E \u0307 kin $, and $ P visc $) for the case at *Re* = 500 and *θ* = 0 is given in the supplementary material (see Fig. S3), verifying the energy balance equation [Eq. (5)]. In Fig. 4, a longer transient period is observed at *θ* = 50°, but the flow eventually reaches a steady state. At *θ* = 60°, we see that *E*_{kin} increases without reaching saturation within the default simulation time *t*_{⊕} = 5. This phenomenon is directly related to the well-known anomalous behavior of the Ekman boundary layer thickness at the so-called *critical latitude*, which in our case is precisely at *θ* = 60° [see, e.g., Sec. 8.07.1.4 in the study by Tilgner (2015)].

Figures 5 and 6 show the viscous power per unit area $ P visc $ for the cases *θ* = 50° and 60° with longer simulation time. The time-averaged value over four oscillation periods is plotted in dots. The orange dot indicates the value obtained in the default range, 1 ≤ *t*_{⊕} ≤ 5. As can be seen, for *θ* = 50°, the power reaches a steady value within *t*_{⊕} ≤ 5, suggesting that *t*_{⊕} = 5 is sufficient for most cases. On the other hand, *θ* = 60° requires a longer transient before reaching a steady value, which approximates the laminar solution. This finding suggests that we should replace the value obtained in the default range by its laminar solution when computing the total dissipation.

### B. Mean velocity profiles

Figure 7 shows the horizontal velocity $ Q + = u x \u2032 \u0304 2 + u y \u2032 \u0304 2 / u * $ at different Reynolds numbers with *θ* = 0 and *t*_{⊕} = 5. The overbar indicates the spatial average on the horizontal plane, and the prime symbol denotes the difference with respect to the boundary velocity, i.e., $ u i \u2032 = u i \u2212 u i ( z = 0 ) $. We adopt this notation in order to compare our results to the law of wall formalism [see the study by Tennekes (1972)]. According to the formalism, the velocity profile is divided into three parts: a viscous sublayer *Q*^{+} = |*z*^{+}|, a buffer layer, and an inertial sublayer *Q*^{+} = (1/0.41) log|*z*^{+}| + 5. The inertial sublayer is also known as the logarithmic region. As shown in Fig. 7, the velocity profiles approach the one predicted by the law of the wall when *Re* ≥ 300. A constant slope reveals the presence of a logarithmic region. The inset of Fig. 7 shows the Kármán measure $\Xi = ( | z + | d Q + / d | z + | ) \u2212 1 $ (i.e., the reciprocal of the slope) for *Re* ≥ 300. A clear logarithmic region develops after *Re* ≈ 300. In the supplementary material, we show the velocity profiles at different time instants for the case *Re* = 300 and find that the appearance of a logarithmic layer coincides with an appreciable deviation of the dissipation from its laminar value, as expected from the presence of turbulence. Besides, the velocity profiles show a maximum above the logarithmic region, which is usually referred to as a *low-level jet* and is also observed for a steady Ekman layer (Deusebio , 2014). The fact that the maximum is not well represented in the case of *Re* = 700 implies that the domain height should be increased accordingly if we want to extend the *Re* range.

### C. Similarity theory

*u*

_{*}and

*Re*of the following form:

*A*,

*B*, and

*C*

_{5}are to be determined by fitting with either experimental or numerical results. We use our numerical solution from

*Re*= 300, 400, 500, 600, and 700 to find the best fit parameters

*A*= 5.74 and

*B*= 1.55 with

*C*

_{5}= −25 (solid lines in Figs. 8 and 9). The parameters are obtained by performing a least-squares fit. Dashed lines show analytical solutions derived in the laminar regime for comparison. The fitting is good for the friction velocity

*u*

_{*}but rather inaccurate at

*Re*≈ 200 for the veering angle

*β*. For comparison, results for an steady Ekman layer are plotted in gray:

*u*

_{*}is about 8% larger while

*β*seems to reach an asymptotic limit at

*Re*≥ 10

^{3}. Comparison to the model proposed by Sous (2013) is shown in the supplementary material (Figs. S5 and S6). Based on our best fit parameters, we compute

*u*

_{*}and

*β*for a given

*Re*in the turbulent regime and then use Eq. (7) to obtain the dissipation, comparing it with the one obtained from our local model.

### D. Viscous dissipation

Figure 10 shows the viscous power per unit area $ P visc * $ at each co-latitude *θ*, which we compute by taking a time-average in the default range, 1 ≤ *t*_{⊕} ≤ 5, except for the case *Re* = 300 (see Sec. IV A). Note that we multiply the power by *Re* to remove the viscosity dependence and to see the effect of turbulence. As it can be seen, when *Re* ≥ 300, the power increases at almost every *θ* compared to the laminar solution, while the models with *Re* ≤ 200 remain close to the laminar solution. Given *Re* ≥ 300, the increments are not the same across *θ*. No increment is observed at the equator and is barely present near the critical latitude. As shown in Sec. IV A, the small increment at the critical latitude is the consequence of insufficient simulation time; the power will eventually reach the laminar solution, given long enough simulation time. This latitude variation can be attributed to the orientation of the rotation vector of the fluid, which goes from perpendicular to the flow at the pole to parallel at the equator. When parallel, there is no Coriolis force, and the Ekman boundary layer flow will reduce to the well-studied Stokes type, whose critical Reynolds number for the fully turbulent regime *Re* = 2500 [see, e.g., the study by Ozdemir (2014)] is higher than the present setting.

*Re*. We obtain a global value by integrating the contributions from each numerical simulation at different values of

*θ*, using the trapezoidal rule. That means the dissipation in the interval, say, for example, 10° <

*θ*< 20°, is obtained by linear interpolation. Since the flow is symmetrical along the azimuthal direction and with respect to the equator, we only need to integrate over the range 0 ≤

*θ*≤ 90°. As we have seen in Sec. IV A,

*E*

_{kin}remains steady for most cases, and the viscous power from surface shear stress will be equal but with an opposite sign to the internal dissipation. Concerning the diverging case at the critical latitude, we replace the value by its laminar solution to avoid overestimation. Therefore, we integrate $\u2212 P visc $ over the spherical surface to obtain $ D \nu $. To see the effect of turbulence, we normalize the numerical solution with the laminar solution [Eq. (6)], so it should be very close to unity in the laminar regime if our numerical approach is valid. As can be seen at

*Re*= 5, the difference in the laminar solution is about −4%, which is associated with the limited number of co-latitudes used. The dissipation deviates from the laminar solution when

*Re*≥ 300. A similar behavior was also observed in the global model by Cébron (2019), where dissipation is found to be mainly confined to the boundary layer. Next, we plot the dissipation predicted by the turbulence model $ D \nu turb $ [Eq. (7)] using a dashed line (Fig. 11). It can be seen that the prediction overestimates the dissipation when

*Re*≥ 300. Since the turbulence model assumes that the parameter

*k*at different values of

*θ*is the same as the one at the pole, this overestimate is possibly due to the latitudinal dependence shown in Fig. 10. Even from the limited

*Re*range given in this study, latitudinal variation in the turbulence effect is evident, thus raising the concern of the reliability of the turbulence model based on

*k*, which is typically treated as a constant in the literature. To accommodate the discrepancy, we introduce a new parameter

*χ*and rewrite the turbulent dissipation as

*χ*= 0.85 marks the lower bound for the results obtained in this study. The dissipation obtained at

*Re*= 300, 400, and 500 falls on the curve of

*χ*= 0.85, which means the total dissipation is reduced by about 15% due to the possible laminar flow near the equator. The excess dissipation at

*Re*= 700 shows a hint of increasing dissipation when

*Re*is even larger.

## V. DISCUSSION AND CONCLUSIONS

We have derived the governing equation for the boundary layer generated by the precession-driven flow in a spherical cavity using a local Cartesian coordinate system. Both the order of magnitude estimation and the numerical solution suggest that the Poincaré term is negligible in the boundary layer if precession is weak. In order to solve the governing equation numerically, we constructed a local Cartesian box model to simulate the boundary layer flow. In the laminar regime, we found that the total energy dissipation, which is computed from the local contributions obtained from the individual box models, is consistent with the laminar value from previous studies based on a global model [e.g., Cébron (2019)]. In the turbulent regime, we have found that (1) the flow at certain latitudes (i.e., *θ* = 60° and 90°) may remain laminar and (2) the total energy dissipation deviates appreciably from the laminar value. For the flow at the pole, *θ* = 0, the energy dissipation begins to deviate from the analytical laminar expectation when *Re* ≥ 300, and the flow exhibits turbulent characteristics, such as a logarithmic layer region, allowing us to use an asymptotic similarity formalism. Using this, we have computed the asymptotic curves for both *u*_{*} and *β* as a function of *Re*, which can be used to compute the total energy dissipation in the turbulent regime via turbulence models.

In the conventional turbulence model, the total turbulent dissipation is given as −(45*π*/32)*kI*ΔΩ^{3}, where parameter *k* is supposed to quantify the asymptotic behavior of the local dissipation in the turbulent regime at high *Re*. Note that the derivation of the model implicitly assumes that *k* is independent of latitude. Following the same approach as that followed by Cébron (2019) to estimate *k* via $k= ( u * / U ) 2 \u2061cos\beta $, we have found that both (*u*_{*}/*U*) and *β* at high *Re* can be computed from the best fit asymptotic curve obtained at *θ* = 0. In order to recover the laminar solution [Eq. (6)] when *Re* goes into the laminar regime, we propose a turbulence model by writing the total dissipation as Eq. (7). The overestimation of the turbulence model, compared to the value obtained from the numerical solution, suggests that the assumption on *k* may be inappropriate. At least for the case of *Re* = 500, our numerical solution shows that the flow is not even turbulent at certain latitudes. A thorough investigation on the latitudinal variation on *k* combined with an investigation of the flow asymptotic behavior at each latitude would be helpful to establish a more reliable turbulence model. That aside, the uncertainty of the turbulence model given by Eq. (7) is about 15% in the *Re* ≤ 700 range.

The previous study by Deusebio (2014) shows that the asymptotic theory applied to the intermediate *Re* range (400 < *Re* < 775), where a logarithmic layer is observed, is sufficient to get a good prediction of *u*_{*} and *β* at high *Re* (=1600). We adopt the theory for an oscillatory Ekman layer at *θ* = 0 not only because we see the logarithmic layer in the flow (Sec. IV B) but also because the theory should be generally unaffected by imposing a time varying boundary velocity. The theory developed by Csanady (1967) and Spalart (1989) is based on the assumption that there exists an overlapping region where the velocity profiles in both the inner layer (viscous sublayer) and the outer layer (free-stream region) match, a fact unaffected by imposing time dependence.

Application to the earth’s core is possible since the numerical results presented here are based on the reduced form of the governing equation [Eq. (4)], which is identical to the one derived for a spheroidal cavity such as Earth by Buffett (2021). We are able to reach the earth’s parameter regime with the local Cartesian box model, so we do not need to rely on any turbulence model, which will be discussed in the following. Our numerical result shows that the total energy dissipation in the turbulent regime at *Re* = 500 is 1.86 times of the laminar value. If the damping of FCN shares the same increasing factor, then the effective viscosity will increase by a factor of 3.5, which is still far below the required amount (6 orders of magnitude higher) inferred from the observation. It implies that the electromagnetic coupling may play a significant role in closing the gap.

An accurate calculation of the energy dissipation at the CMB in the moon is unfeasible since we not only lack a proper turbulence model but also because interior models obtained from seismic and gravity measurements are poorly constrained. There are three relevant parameters: (1) the angle between the rotation axes of the core and the mantle, (2) the lunar core radius, and (3) the moment of inertia of the fluid core. There is no direct measurement for the orientation of the core rotation axis. According to the dynamic argument that the lunar CMB flattening is smaller than its Poincaré number, the core is unlikely to share the precession of the mantle (Goldreich, 1967); the core rotation axis should thus be closely aligned with the ecliptic normal. We use the reduced model by Cébron (2019) to compute the orientation of **Ω**_{f} (see Table I). The plausible range of the core radius obtained from gravity measurements is 200–380 km (Williams , 2014; Viswanathan , 2019). If we take the core radius as 380 km, then the relative velocity at CMB will be around 2.72 × 10^{−2} m s^{−1}, and the Reynolds number will be 1.67 × 10^{4}. We obtain the moment of inertia of the fluid core by taking the moment ratio *C*_{f}/*C* = 7 × 10^{−4} (Williams , 2014). With our best fit parameters (*A*, *B*, *C*_{5}) and the reduced turbulence model [Eq. (11)], we obtain $ D \nu turb =72$ or 84 MW for *χ* = 0.85 or 1.0, respectively. For a smaller core, say *Re* = 200 km, with the same moment ratio, the Reynolds number reduces to 8.79 × 10^{3} while the dissipation increases by about 10%. For comparison, we also compute the dissipation with the parameters (*A* = 7.32, *B* = −2.67) (Sous , 2013, the conversion to the present notation is presented in Appendix C) and the conventional turbulence model, −(45*π*/32)*kI*ΔΩ^{3} [see Sec. II; e.g., Cébron (2019)]; the obtained value is 126 MW, about 60% larger than our result. According to the analysis of lunar laser ranging (LLR) data, the dissipation at CMB is about 58 MW (Williams , 2001), 84 MW (Williams , 2014), and 73 MW (Williams and Boggs, 2015), respectively. Concerning the uncertainty of the interior model, it is still insufficient to say if our model matches the observation or not, but at least we provide a set of parameters with numerical support, which can be used for the turbulence model. A sensitivity test on the parameters (*A*, *B*, *C*_{5}) shows that (1) *β* is sensitive to *A* and also sensitive to *C*_{5} when *Re* < 10^{3} and (2) the sensitivity of *u*_{*}/*U* to both *A* and *B* decreases as *Re* increases (see Figs. S7–S9 in the supplementary material).

The implication for planetary cores with turbulent boundary layers should be interpreted with care. In this study, the quoted *Re* for Earth and the Moon is obtained based on the laminar solution [e.g., Busse (1968) and Cébron (2019)], where the boundary layer is treated in the laminar regime, even though the obtained *Re* suggests the other way around. For Earth, the conclusion that the turbulence effect is unlikely to account for the observed FCN damping remains as Earth is at the margin of turbulence. For the Moon, our purpose is to demonstrate a way to forward-compute the dissipation, and a more comprehensive model is needed for future exploration.

We aim to demonstrate that through appropriate formulation, it is feasible to complement the results from global models. By narrowing our focus to the boundary layer exclusively, we can expand the parameter range and explore the turbulent regime. Potential applications of this approach include spin-up, libration, and precessing cylinders (Pizzi , 2021).

## SUPPLEMENTARY MATERIAL

See the supplementary material for the figures of friction velocity showing the comparison between the solutions obtained from Eqs. (2) and (4), velocity profiles for the domain height test, energy balance for the case of Re = 500 and θ = 0, velocity profiles of Re = 300 and θ = 0 at different time instants, u* and β as a function of Re with the comparison to the result of Sous et al. (2013), and sensitivity test on (A, B, C5), β as a function of θ.

## ACKNOWLEDGMENTS

We are grateful for the comments and suggestions by all the anonymous reviewers. S.-A.S. would like to thank Benjamin Favier and Philipp Schlatter for the helpful discussion on the numerical methods. The numerical solutions are obtained using NEK5000 Version 20.0-dev at Argonne National Laboratory, Illinois. (available at https://nek5000.mcs.anl.gov). S.-A.S. is a Research Fellow of the F.R.S.-FNRS. V.D., J.R., and S.A.T. acknowledge the funding from the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Program (GRACEFUL Synergy Grant Agreement No. 855677).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Sheng-An Shih (施勝安)**: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Methodology (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). **Santiago Andrés Triana**: Conceptualization (equal); Methodology (supporting); Supervision (equal); Writing – review & editing (equal). **Jérémy Rekier**: Conceptualization (equal); Writing – review & editing (equal). **Véronique Dehant**: Conceptualization (equal); Funding acquisition (equal); Supervision (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX A: TRANSFORMATION BETWEEN REFERENCE FRAMES AND COORDINATE SYSTEMS

**Ω**

_{f}from the study by Cébron (2019) is given in the precessing frame where the coordinate system $ ( X \u0302 p \u2032 , Y \u0302 p \u2032 , Z \u0302 p \u2032 ) $ is defined such that $ Z \u0302 p \u2032 $ is aligned with

**Ω**

_{m}and

**Ω**

_{p}is on the $ X \u0302 p \u2032 $–$ Z \u0302 p \u2032 $ plane with the co-latitude

*α*(Fig. 12). We use (

*θ*

_{f},

*ϕ*

_{f}) to denote the orientation of the fluid rotation vector. Since the global coordinate system $ ( X \u0302 , Y \u0302 , Z \u0302 ) $ defined in the present study has two features—$ Z \u0302 $ is aligned with

**Ω**

_{f}and it is defined in the fluid frame—we need a series of transformations to adopt the solution in the present study. First, we implement a two-step rotation: (1) rotation about $ Z \u0302 p \u2032 $ in the counterclockwise direction by

*ϕ*

_{f}and (2) rotation about the intermediate axis $ Y \u0302 p \u2032 \u2032 $ in the counterclockwise direction by

*θ*

_{f}[see, e.g., the study by Thornton and Marion (2003)]. The transformation matrices are

**in the new coordinate system $ ( X \u0302 p , Y \u0302 p , Z \u0302 p ) $ is given by $ [ x ] p = \lambda 2 \lambda 1 [ x ] p \u2032 $. The rotations are chosen so that**

*x***Ω**

_{m}is on the $ X \u0302 p $–$ Z \u0302 p $ plane and leaning toward negative $ X \u0302 p $ with a co-latitude

*θ*

_{f}, which gives the boundary conditions consistent with the ones mentioned in the main text. Second, a transformation from the precessing frame to the fluid frame is given by [see, e.g., the study by Noir and Cébron (2013)]

**Ω**

_{f}, which in the $ ( X \u0302 p \u2032 , Y \u0302 p \u2032 , Z \u0302 p \u2032 ) $ system is given by Ω

_{f}(sin

*θ*

_{f}cos

*ϕ*

_{f}, sin

*θ*

_{f}sin

*ϕ*

_{f}, cos

*θ*

_{f}), in the local coordinate system $ ( x \u0302 , y \u0302 , z \u0302 ) $ are (0, Ω

_{f}sin

*θ*, Ω

_{f}cos

*θ*). The coordinates of

**Ω**

_{p}and $ R \u0302 \xd7 ( k \u0302 p \xd7 k \u0302 f ) $ are computed in the same way.

### APPENDIX B: ANALYTICAL SOLUTION

*P*is the reduced pressure, including the centrifugal force. In the local coordinate system, the vector $ k \u0302 f $ is given by $cos\u2061\theta z \u0302 +sin\u2061\theta y \u0302 $, where

*θ*is the co-latitude. Following Eq. (3), the boundary conditions require

**= 0 at**

*u**z*= −

*∞*; at

*z*= 0, the velocities are

*u*

_{x}=

*U*cos

*θ*cos(Ω

_{f}

*t*+

*ϕ*),

*u*

_{y}= −

*U*sin(Ω

_{f}

*t*+

*ϕ*), and

*u*

_{z}= 0, where

*U*is the relative velocity of the mantle and

*ϕ*is the longitude. We have to note that the governing equations generally imply a simple rotating fluid with angular velocity

**Ω**

_{f}, and this due to the fact that the Poincare term is neglected in the boundary layer of the precession-driven flow. However, precession enters the formulation through the boundary conditions. Nevertheless, the following formulation, borrowed from the classic Ekman layer analysis, serves not only to demonstrate the flow but also to validate the result against known solutions in the literature.

*x*and

*y*. From the continuity equation ∇ ·

**= 0, we further have**

*u**∂u*

_{z}/

*∂z*= 0. Given that no vertical motion is allowed on the top boundary, we end up with

*u*

_{z}= 0 everywhere. Equation (B1) then becomes

*u*

_{x}=

*u*

_{y}= 0, the horizontal pressure gradient should be zero:

*∂P*/

*∂x*=

*∂P*/

*∂y*= 0. Since

*u*

_{x}is independent of

*x*and

*y*, the derivative of Eq. (B5) shows that the horizontal pressure gradient is an invariant of height. That means the horizontal pressure gradients are not involved in the horizontal force balance.

*q*

_{x}(

*z*) and

*q*

_{y}(

*z*) are the functions to be solved. Substituting the solutions into Eqs. (B3)–(B5) and solving for

*q*

_{i}(

*z*), we end up with

*C*

_{1}=

*U*(cos

*θ*− 1)/2,

*C*

_{2}=

*U*(cos

*θ*+ 1)/2, and $ \sigma \xb1 = \Omega f / \nu ( 1 \xb1 2 cos \theta ) i $. Notice that the sign of the term (1 ± 2 cos

*θ*) changes at

*θ*= 60° and 120°, so we cannot simply expand

*σ*

_{±}at this stage. The above-mentioned expression is easy to implement for numerical computation. To get the physical picture, we expand the solution at the three representative co-latitudes:

*θ*= 0°, 90°, and 60°. At the north pole (

*θ*= 0), the solution reads

*t*

_{0}and

*ϕ*

_{0}, the solution shows an Ekman spiral, which decays exponentially to the depth. On top of that, the spiral is rotating in the clockwise direction. At the equator (

*θ*= 90°), the boundary drives the fluid in the direction parallel to the rotation vector: there is no Coriolis force acting on the flow. The type of the boundary flow then reduces to the Stokes boundary layer [see, e.g., the study by Landau and Lifshitz (2013)],

*σ*

_{±}becomes zero, the velocity profile shows a non-decaying component; here is the case of

*θ*= 60°,

*s*

_{ik}=

*μ*(

*∂*

_{i}

*u*

_{k}+

*∂*

_{k}

*u*

_{i}), where

*μ*is the dynamic viscosity (note

*ν*=

*μ*/

*ρ*is the kinematic viscosity) [see, e.g., Eqs. (15.2) and (15.8) in the study by Landau and Lifshitz (2013)]. On the top boundary

*z*= 0, given the normal vector $ n \u0302 = z \u0302 $, the surface shear stress exerted on the fluid is

*τ*

_{i}=

*s*

_{ij}

*n*

_{j}, i.e.,

*τ*

_{x}=

*μ∂*

_{z}

*u*

_{x}and

*τ*

_{y}=

*μ∂*

_{z}

*u*

_{y},

*θ*= 0, we have

*β*as the phase difference between the boundary velocity and the surface shear stress (One can also consider the velocity and the stress vectors in complex form, e.g., $ u \u0332 = u x +i u y $, and define the veering angle

*β*as the argument difference between $ \tau \u0332 $ and $ u \u0332 $, i.e., $\beta =Arg [ \tau \u0332 ] \u2212Arg [ u \u0332 ] $. However, the definition is only valid at the poles where the vector components have the same magnitude.). It determines the magnitude of the viscous power from

**: smaller**

*τ**β*means larger power. At

*θ*= 90°, we have

*θ*= 90° shows that the stress is leading the velocity by 45°. We note that the phases derived from the

*x*- and

*y*-components evolve in opposite directions when

*θ*increases from 0° (see Fig. S10 in the supplementary material). Derived from the

*y*-component,

*β*increases from −45° at the north pole to +45° at the critical latitude and then remains until the equator. The flow shows distinct features in the polar region and equatorial region, with the border being

*θ*= 60°.

*S*of radius

*R*by integrating the local contributions:

**Γ**=

*∫*

**×**

*R*

*τ**dS*. In a spherical coordinate system, we have

*τ*

_{θ}= −

*τ*

_{y}and

*τ*

_{ϕ}=

*τ*

_{x}, where

*τ*

_{x}and

*τ*

_{y}are local stresses given by Eqs. (B17) and (B18). The local torque is $R\xd7\tau =\u2212R \tau \varphi \theta \u0302 +R \tau \theta \varphi \u0302 $. The total torque is given by integrating over the whole spherical surface. After changing to a global Cartesian coordinate system $ ( X \u0302 , Y \u0302 , Z \u0302 ) $, we have the torque exerted on the core,

*I*= (8

*π*/15)

*ρR*

^{5}is the moment of inertia of the sphere. In the fluid frame, the torque circulates in the clockwise direction. Transformed back to the precessing frame [see Eq. (A3)], the torque becomes $\Gamma =\u2212I\Delta \Omega \nu \Omega f R 2.62 X \u0302 p + 0.258 Y \u0302 p $. We note that the imposed boundary velocity in our local model implies that the vector of differential rotation is given by $\Delta \Omega =\u2212\Delta \Omega X \u0302 p $. In this case, Eq. (B26) agrees with the ones derived from previous studies [e.g., the studies by Yoder (1981), Williams (2001), and Cébron (2019)], where the results are based on the solution for the spin-over mode (Greenspan, 1968). In principle, the boundary velocity can be modified to represent other types of flow (e.g., spin-up), and the above-mentioned formulation can be applied to derive the corresponding torque [e.g., the study by Noir and Cébron (2013)].

### APPENDIX C: SIMILARITY THEORY

*ϑ*=

*β*[or set

*C*

_{5}= 0; see, e.g., the study by Coleman (1990)]. To avoid confusion, we label the constants in the study by Sous (2013) as $ A \u0303 $ and $ B \u0303 $. Comparing to their Eqs. (12) and (13), we find that

*z*

_{0}= 0.11

*ν*/

*u*

_{*}for the smooth surface case. Their best fit result gives $ ( A \u0303 = 3.3 , B \u0303 = 3.0 ) $, which in our convention is (

*A*= 7.32,

*B*= −2.67).

## REFERENCES

*P*,

*T*condition

*Fluid Mechanics: Landau and Lifshitz: Course of Theoretical Physics*

*A First Course in Turbulence*

*Classical Dynamics of Particles and Systems*