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.

Topics
Rotational dynamics, Computational fluid dynamics, Turbulence theory and modelling, Flow boundary effects

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 m2 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 m2 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ν ∝ E1/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 104, 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 21/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.

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 = ν/(ΩmR2), 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 Δ Ω = Ω m − Ω f ⁠. A boundary layer exists, the main focus of our study, with a nominal thickness L ν = ν / Ω f ⁠.

In the fluid frame, we define two Cartesian coordinate systems: one global ( X ̂ , Y ̂ , Z ̂ ) and one local ( x ̂ , y ̂ , z ̂ ) (Fig. 1). For the global coordinate system, Z ̂ is aligned with Ω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 u satisfies the following Navier–Stokes equation:
(1)
where P is the reduced pressure including the centrifugal force and r′ is the position vector referring to any point inside the cavity. The last term, r ′ × Ω p × Ω 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 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
(2)
where Re = ULν/ν is the Reynolds number and R ̂ = R / R is a unit vector. A superscript (*) denotes dimensionless variables (the variable u* is not to be confused with the friction velocity u* introduced later). In the local coordinate system, the vector k ̂ f is given by cos θ z ̂ + sin θ y ̂ ⁠, where θ is the co-latitude.
FIG. 1.
Sketch of a spherical cavity with a local box. Two coordinate systems are shown: global ( X ̂ , Y ̂ , Z ̂ ) and local ( x ̂ , y ̂ , z ̂ ) . The two vectors x ̂ and y ̂ denote the unit vectors along the directions of longitude ϕ and (negative) co-latitude θ, respectively.
View largeDownload slide

Sketch of a spherical cavity with a local box. Two coordinate systems are shown: global ( X ̂ , Y ̂ , Z ̂ ) and local ( x ̂ , y ̂ , z ̂ ) ⁠. The two vectors x ̂ and y ̂ denote the unit vectors along the directions of longitude ϕ and (negative) co-latitude θ, respectively.

FIG. 1.
Sketch of a spherical cavity with a local box. Two coordinate systems are shown: global ( X ̂ , Y ̂ , Z ̂ ) and local ( x ̂ , y ̂ , z ̂ ) . The two vectors x ̂ and y ̂ denote the unit vectors along the directions of longitude ϕ and (negative) co-latitude θ, respectively.
View largeDownload slide

Sketch of a spherical cavity with a local box. Two coordinate systems are shown: global ( X ̂ , Y ̂ , Z ̂ ) and local ( x ̂ , y ̂ , z ̂ ) ⁠. The two vectors x ̂ and y ̂ denote the unit vectors along the directions of longitude ϕ and (negative) co-latitude θ, respectively.

Close modal
Before moving on to the boundary conditions, we find it necessary to specify the orientations of Ω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 ̂ p , Y ̂ p , Z ̂ p ) ⁠, is defined such that it coincides with ( X ̂ , Y ̂ , Z ̂ ) at t = 0; otherwise, the latter is rotating with Ωf in the precessing frame. Moreover, Ωm is tilted toward negative X ̂ 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: Δ Ω = − Δ Ω X ̂ p ⁠. The coordinate is then transformed to ( X ̂ , Y ̂ , Z ̂ ) in the fluid frame and further to ( x ̂ , y ̂ , z ̂ ) for the computation of boundary velocity, which is ΔΩ × R. Explicit expressions of the transformations are given in  Appendix A. In the local coordinate system, the boundary condition at the lower end (z* → −∞) is u* = 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 z* = −30. On the fluid surface (z* = 0, corresponding to the CMB),
(3a)
(3b)
(3c)
where ϕ 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 Ωft 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 ̂ p ′ , Y ̂ p ′ , Z ̂ p ′ ) is defined such that the Z ̂ p ′ -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 ̂ p ′ -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.

We use the lunar parameters listed in Table I as an end-member case to estimate the amplitude of the Poincaré term, where k ̂ p × k ̂ f = sin γ ⁠, and find that the amplitude is 3 × 10−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 × 104. 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 Ω m / Ω f ⁠, which is equal to cos−1 θf under the no spin-up condition. Given θf = 1.543° at Re = 1.67 × 104, 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 Ω f / Ω 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),
(4)
For comparison, we numerically solve both Eqs. (2) and (4) at 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.
TABLE I.

Lunar parameters used in this study. R is the radius of the lunar core-mantle boundary (Viswanathan , 2019). α is the angle between the symmetry axis and the ecliptic normal. We use Eqs. (5)–(7) of the study by Cébron  (2019) to compute Ωf; the orientation (θf, ϕf) is given by the coordinate system ( X ̂ p ′ , Y ̂ p ′ , Z ̂ p ′ ) whose Z ̂ p ′ is aligned with Ωm (see  Appendix A). The angle between Ωf and Ωp is denoted by γ.

Parameter Definition Units Value
ν    m2 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ν  ν / Ω 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  ν/(R2Ωm)    3 × 10−12 
Re  ULν/ν    1.67 × 104 
Parameter Definition Units Value
ν    m2 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ν  ν / Ω 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  ν/(R2Ωm)    3 × 10−12 
Re  ULν/ν    1.67 × 104 
View Large
Next, we compute the dissipation in the boundary layer. In the local system, the input (output) power from the surface shear stress τ is equal to the increase (decrease) of internal and kinetic energies [see, e.g., Sec. 16 in the study by Landau and Lifschitz (2013)],
(5)
where 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 θ R2dθdϕ, we compute the total dissipation via D ν = − ∫ ( τ ⋅ u ) d S ⁠. The local viscous power per unit area is defined as P visc = τ ⋅ u 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)]
(6)
where I = (8π/15)ρR5 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.
In the turbulent regime, the surface shear stress no longer follows the analytical solution, so we prescribe the turbulent stress as kρ|u|u [e.g., in the study by Williams  (2001)], where 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 Δ Ω = Δ Ω Z ̂ 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 ⁡ cos β ⁠, where u * = | τ | / ρ 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 * = ν U / L ν 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
(7)
In principle, one can match P visc at any θ. 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.

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.

FIG. 2.
Box model with Cartesian coordinates. The typical resolution is E = 768 and N = 9. The dimension of the box is 30 × 30 × 30 in the unit of Lν. The element density is higher near the top boundary to resolve the boundary flow.
View largeDownload slide

Box model with Cartesian coordinates. The typical resolution is E = 768 and N = 9. The dimension of the box is 30 × 30 × 30 in the unit of Lν. The element density is higher near the top boundary to resolve the boundary flow.

FIG. 2.
Box model with Cartesian coordinates. The typical resolution is E = 768 and N = 9. The dimension of the box is 30 × 30 × 30 in the unit of Lν. The element density is higher near the top boundary to resolve the boundary flow.
View largeDownload slide

Box model with Cartesian coordinates. The typical resolution is E = 768 and N = 9. The dimension of the box is 30 × 30 × 30 in the unit of Lν. The element density is higher near the top boundary to resolve the boundary flow.

Close modal

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 30Lν 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.

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.

The friction velocity u* is a quantity commonly used to monitor the amplitude of surface shear stress |τ|: u * = | τ | / ρ [e.g., in the study by Tennekes  (1972)]. At θ = 0, Eqs. (B19) and (B20) give the laminar solution u * lam = R e − 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).

FIG. 3.
Friction velocity ratio u * / u * lam as a function of time with different values of Re at θ = 0. The means deviate appreciably from unity when Re ≥ 300 and almost reach 1.3 when Re = 500. The Re = 5, 100, 200 curves overlap nearly completely.
View largeDownload slide

Friction velocity ratio u * / u * lam as a function of time with different values of Re at θ = 0. The means deviate appreciably from unity when Re ≥ 300 and almost reach 1.3 when Re = 500. The Re = 5, 100, 200 curves overlap nearly completely.

FIG. 3.
Friction velocity ratio u * / u * lam as a function of time with different values of Re at θ = 0. The means deviate appreciably from unity when Re ≥ 300 and almost reach 1.3 when Re = 500. The Re = 5, 100, 200 curves overlap nearly completely.
View largeDownload slide

Friction velocity ratio u * / u * lam as a function of time with different values of Re at θ = 0. The means deviate appreciably from unity when Re ≥ 300 and almost reach 1.3 when Re = 500. The Re = 5, 100, 200 curves overlap nearly completely.

Close modal
TABLE II.

Main results of the present study. Time-averaged values are obtained from the range 1 ≤ t⊕ ≤ 5, except for the case Re = 300 where the range used is 3 ≤ t⊕ ≤ 5. The values of the friction velocity and the veering angle are obtained from the cases with θ = 0.

Re u * / u * lam u * / U β (deg) D ν / D ν 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 β (deg) D ν / D ν 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 
View Large

Figure 4 shows the total kinetic energy, Ekin, 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 ̇ 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 ̇ int ⁠, E ̇ 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 Ekin 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)].

FIG. 4.
Total kinetic energy Ekin with different values of θ at Re = 500. Most cases reach steady state within t⊕ = 5, while a long-period modulation is observed for θ = 50° and an indefinite increase is observed for θ = 60°.
View largeDownload slide

Total kinetic energy Ekin with different values of θ at Re = 500. Most cases reach steady state within t⊕ = 5, while a long-period modulation is observed for θ = 50° and an indefinite increase is observed for θ = 60°.

FIG. 4.
Total kinetic energy Ekin with different values of θ at Re = 500. Most cases reach steady state within t⊕ = 5, while a long-period modulation is observed for θ = 50° and an indefinite increase is observed for θ = 60°.
View largeDownload slide

Total kinetic energy Ekin with different values of θ at Re = 500. Most cases reach steady state within t⊕ = 5, while a long-period modulation is observed for θ = 50° and an indefinite increase is observed for θ = 60°.

Close modal

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.

FIG. 5.
Viscous power per unit area P visc (light blue curve) for the cases θ = 50° with longer simulation time. Blue dot: time-averaged value over four oscillation periods. Orange dot: value obtained in the default range 1 ≤ t⊕ ≤ 5. The default range is sufficient to obtain a steady value. The dashed line shows the laminar solution obtained in Appendix B.
View largeDownload slide

Viscous power per unit area P visc (light blue curve) for the cases θ = 50° with longer simulation time. Blue dot: time-averaged value over four oscillation periods. Orange dot: value obtained in the default range 1 ≤ t⊕ ≤ 5. The default range is sufficient to obtain a steady value. The dashed line shows the laminar solution obtained in  Appendix B.

FIG. 5.
Viscous power per unit area P visc (light blue curve) for the cases θ = 50° with longer simulation time. Blue dot: time-averaged value over four oscillation periods. Orange dot: value obtained in the default range 1 ≤ t⊕ ≤ 5. The default range is sufficient to obtain a steady value. The dashed line shows the laminar solution obtained in Appendix B.
View largeDownload slide

Viscous power per unit area P visc (light blue curve) for the cases θ = 50° with longer simulation time. Blue dot: time-averaged value over four oscillation periods. Orange dot: value obtained in the default range 1 ≤ t⊕ ≤ 5. The default range is sufficient to obtain a steady value. The dashed line shows the laminar solution obtained in  Appendix B.

Close modal
FIG. 6.
Same as Fig. 5 but for θ = 60°. The simulation runs longer than the default range required for the viscous power to reach a steady value, which seems to approach the laminar solution.
View largeDownload slide

Same as Fig. 5 but for θ = 60°. The simulation runs longer than the default range required for the viscous power to reach a steady value, which seems to approach the laminar solution.

FIG. 6.
Same as Fig. 5 but for θ = 60°. The simulation runs longer than the default range required for the viscous power to reach a steady value, which seems to approach the laminar solution.
View largeDownload slide

Same as Fig. 5 but for θ = 60°. The simulation runs longer than the default range required for the viscous power to reach a steady value, which seems to approach the laminar solution.

Close modal

Figure 7 shows the horizontal velocity Q + = u x ′ ̄ 2 + u y ′ ̄ 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 ′ = u i − 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 Ξ = ( | z + | d Q + / d | z + | ) − 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.

FIG. 7.
Mean horizontal velocity Q + = u x ′ ̄ 2 + u y ′ ̄ 2 / u * at different Reynolds numbers with θ = 0 and t⊕ = 5. Dotted line: Q+ = |z+| for |z+| ≤ 10 and Q+ = (1/0.41) log |z+| + 5 for |z+| > 10. (Inset) Kármán measure Ξ = ( | z + | d Q + / d | z + | ) − 1 for Re ≥ 300. Dashed line: Ξ = 0.41. The x-axis is |z+|, and the y-axis is Ξ.
View largeDownload slide

Mean horizontal velocity Q + = u x ′ ̄ 2 + u y ′ ̄ 2 / u * at different Reynolds numbers with θ = 0 and t⊕ = 5. Dotted line: Q+ = |z+| for |z+| ≤ 10 and Q+ = (1/0.41) log |z+| + 5 for |z+| > 10. (Inset) Kármán measure Ξ = ( | z + | d Q + / d | z + | ) − 1 for Re ≥ 300. Dashed line: Ξ = 0.41. The x-axis is |z+|, and the y-axis is Ξ.

FIG. 7.
Mean horizontal velocity Q + = u x ′ ̄ 2 + u y ′ ̄ 2 / u * at different Reynolds numbers with θ = 0 and t⊕ = 5. Dotted line: Q+ = |z+| for |z+| ≤ 10 and Q+ = (1/0.41) log |z+| + 5 for |z+| > 10. (Inset) Kármán measure Ξ = ( | z + | d Q + / d | z + | ) − 1 for Re ≥ 300. Dashed line: Ξ = 0.41. The x-axis is |z+|, and the y-axis is Ξ.
View largeDownload slide

Mean horizontal velocity Q + = u x ′ ̄ 2 + u y ′ ̄ 2 / u * at different Reynolds numbers with θ = 0 and t⊕ = 5. Dotted line: Q+ = |z+| for |z+| ≤ 10 and Q+ = (1/0.41) log |z+| + 5 for |z+| > 10. (Inset) Kármán measure Ξ = ( | z + | d Q + / d | z + | ) − 1 for Re ≥ 300. Dashed line: Ξ = 0.41. The x-axis is |z+|, and the y-axis is Ξ.

Close modal
The similarity theory (Csanady, 1967; Spalart, 1989) establishes a relationship between u* and Re of the following form:
(8)
(9)
(10)
where 0.41 is known as the von Kármán constant. The constants A, B, and C5 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 C5 = −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 ≥ 103. 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.
FIG. 8.
Ratio u*/U as a function of the Reynolds number. Blue circles: numerical solution; black line: best fit with A = 5.74, B = 1.55, and C5 = −25; dashed line: laminar solution. Results for a steady Ekman layer from previous studies are plotted with gray circles (Spalart, 1989; Coleman , 1990; Coleman, 1999; Shingai and Kawamura, 2004; Marlatt , 2012; and Braun , 2020); dotted-dashed line is the corresponding best fit with A = 5.40, B = 0.26, and C5 = −52. In the turbulent regime, u* is about 8% larger than the one for an oscillatory Ekman layer.
View largeDownload slide

Ratio u*/U as a function of the Reynolds number. Blue circles: numerical solution; black line: best fit with A = 5.74, B = 1.55, and C5 = −25; dashed line: laminar solution. Results for a steady Ekman layer from previous studies are plotted with gray circles (Spalart, 1989; Coleman , 1990; Coleman, 1999; Shingai and Kawamura, 2004; Marlatt , 2012; and Braun , 2020); dotted-dashed line is the corresponding best fit with A = 5.40, B = 0.26, and C5 = −52. In the turbulent regime, u* is about 8% larger than the one for an oscillatory Ekman layer.

FIG. 8.
Ratio u*/U as a function of the Reynolds number. Blue circles: numerical solution; black line: best fit with A = 5.74, B = 1.55, and C5 = −25; dashed line: laminar solution. Results for a steady Ekman layer from previous studies are plotted with gray circles (Spalart, 1989; Coleman , 1990; Coleman, 1999; Shingai and Kawamura, 2004; Marlatt , 2012; and Braun , 2020); dotted-dashed line is the corresponding best fit with A = 5.40, B = 0.26, and C5 = −52. In the turbulent regime, u* is about 8% larger than the one for an oscillatory Ekman layer.
View largeDownload slide

Ratio u*/U as a function of the Reynolds number. Blue circles: numerical solution; black line: best fit with A = 5.74, B = 1.55, and C5 = −25; dashed line: laminar solution. Results for a steady Ekman layer from previous studies are plotted with gray circles (Spalart, 1989; Coleman , 1990; Coleman, 1999; Shingai and Kawamura, 2004; Marlatt , 2012; and Braun , 2020); dotted-dashed line is the corresponding best fit with A = 5.40, B = 0.26, and C5 = −52. In the turbulent regime, u* is about 8% larger than the one for an oscillatory Ekman layer.

Close modal
FIG. 9.
Following Fig. 8. The veering angle β between surface shear stress and boundary velocity seems to reach an asymptotic limit when Re ≥ 103.
View largeDownload slide

Following Fig. 8. The veering angle β between surface shear stress and boundary velocity seems to reach an asymptotic limit when Re ≥ 103.

FIG. 9.
Following Fig. 8. The veering angle β between surface shear stress and boundary velocity seems to reach an asymptotic limit when Re ≥ 103.
View largeDownload slide

Following Fig. 8. The veering angle β between surface shear stress and boundary velocity seems to reach an asymptotic limit when Re ≥ 103.

Close modal

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.

FIG. 10.
Time-averaged viscous power per unit area P visc * multiplied by Re at different co-latitudes θ (°). The gray dash line is computed via P visc = τ ⋅ u , where τ and u are given by Eqs. (B6)–(B8), (B17), and (B18); the solutions obtained in the laminar regime where the non-linear term is neglected correspond to the cases of small Re.
View largeDownload slide

Time-averaged viscous power per unit area P visc * multiplied by Re at different co-latitudes θ (°). The gray dash line is computed via P visc = τ ⋅ u ⁠, where τ and u are given by Eqs. (B6)–(B8), (B17), and (B18); the solutions obtained in the laminar regime where the non-linear term is neglected correspond to the cases of small Re.

FIG. 10.
Time-averaged viscous power per unit area P visc * multiplied by Re at different co-latitudes θ (°). The gray dash line is computed via P visc = τ ⋅ u , where τ and u are given by Eqs. (B6)–(B8), (B17), and (B18); the solutions obtained in the laminar regime where the non-linear term is neglected correspond to the cases of small Re.
View largeDownload slide

Time-averaged viscous power per unit area P visc * multiplied by Re at different co-latitudes θ (°). The gray dash line is computed via P visc = τ ⋅ u ⁠, where τ and u are given by Eqs. (B6)–(B8), (B17), and (B18); the solutions obtained in the laminar regime where the non-linear term is neglected correspond to the cases of small Re.

Close modal
Figure 11 shows the total global dissipation D ν in the boundary layer at different values of 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, Ekin 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 − P visc over the spherical surface to obtain D ν ⁠. 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 ν 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
(11)
We find that χ = 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.
FIG. 11.
Boundary layer dissipation D ν as a function of Re. The dissipation is normalized by its laminar solution D ν lam . Blue circle: present study; gray square: Cébron  (2019). Dashed line: turbulence model D ν turb with χ = 1.0 and χ = 0.85.
View largeDownload slide

Boundary layer dissipation D ν as a function of Re. The dissipation is normalized by its laminar solution D ν lam ⁠. Blue circle: present study; gray square: Cébron  (2019). Dashed line: turbulence model D ν turb with χ = 1.0 and χ = 0.85.

FIG. 11.
Boundary layer dissipation D ν as a function of Re. The dissipation is normalized by its laminar solution D ν lam . Blue circle: present study; gray square: Cébron  (2019). Dashed line: turbulence model D ν turb with χ = 1.0 and χ = 0.85.
View largeDownload slide

Boundary layer dissipation D ν as a function of Re. The dissipation is normalized by its laminar solution D ν lam ⁠. Blue circle: present study; gray square: Cébron  (2019). Dashed line: turbulence model D ν turb with χ = 1.0 and χ = 0.85.

Close modal

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 ⁡ cos β ⁠, 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 × 104. We obtain the moment of inertia of the fluid core by taking the moment ratio Cf/C = 7 × 10−4 (Williams , 2014). With our best fit parameters (A, B, C5) and the reduced turbulence model [Eq. (11)], we obtain D ν 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 × 103 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, C5) shows that (1) β is sensitive to A and also sensitive to C5 when Re < 103 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).

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 θ.

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).

The authors have no conflicts to disclose.

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).

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

The solution of Ωf from the study by Cébron  (2019) is given in the precessing frame where the coordinate system ( X ̂ p ′ , Y ̂ p ′ , Z ̂ p ′ ) is defined such that Z ̂ p ′ is aligned with Ωm and Ωp is on the X ̂ p ′ – Z ̂ p ′ 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 ̂ , Y ̂ , Z ̂ ) defined in the present study has two features— Z ̂ 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 ̂ p ′ in the counterclockwise direction by ϕf and (2) rotation about the intermediate axis Y ̂ p ′ ′ in the counterclockwise direction by θf [see, e.g., the study by Thornton and Marion (2003)]. The transformation matrices are
(A1)
and
(A2)
The coordinate of a vector x in the new coordinate system ( X ̂ p , Y ̂ p , Z ̂ p ) is given by [ x ] p = λ 2 λ 1 [ x ] p ′ ⁠. The rotations are chosen so that Ωm is on the X ̂ p – Z ̂ p plane and leaning toward negative X ̂ 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)]
(A3)
Finally, in order to solve the governing equations in the local coordinate system, one further transformation is needed, which is given by
(A4)
where r ̂ = z ̂ ⁠, θ ̂ = − y ̂ ⁠, and ϕ ̂ = x ̂ (see Fig. 1). For example, the coordinates of the fluid rotation vector Ωf, which in the ( X ̂ p ′ , Y ̂ p ′ , Z ̂ p ′ ) system is given by Ωf(sin θf cos ϕf, sin θf sin ϕf, cos θf), in the local coordinate system ( x ̂ , y ̂ , z ̂ ) are (0, Ωf sin θ, Ωf cos θ). The coordinates of Ωp and R ̂ × ( k ̂ p × k ̂ f ) are computed in the same way.
FIG. 12.
Schematic diagram for the associated vectors and transformations. (Left) The angle between Ωm and Ωp is α. The co-latitude and longitude of Ωf are θf and ϕf, respectively. Note that the magnitudes and orientations of Ωm, Ωp, and Ωf are not to scale. (Middle) Counterclockwise rotation through ϕf about Z ̂ p ′ . (Right) Counterclockwise rotation through θf about Y ̂ p ′ ′ .
View largeDownload slide

Schematic diagram for the associated vectors and transformations. (Left) The angle between Ωm and Ωp is α. The co-latitude and longitude of Ωf are θf and ϕf, respectively. Note that the magnitudes and orientations of Ωm, Ωp, and Ωf are not to scale. (Middle) Counterclockwise rotation through ϕf about Z ̂ p ′ ⁠. (Right) Counterclockwise rotation through θf about Y ̂ p ′ ′ ⁠.

FIG. 12.
Schematic diagram for the associated vectors and transformations. (Left) The angle between Ωm and Ωp is α. The co-latitude and longitude of Ωf are θf and ϕf, respectively. Note that the magnitudes and orientations of Ωm, Ωp, and Ωf are not to scale. (Middle) Counterclockwise rotation through ϕf about Z ̂ p ′ . (Right) Counterclockwise rotation through θf about Y ̂ p ′ ′ .
View largeDownload slide

Schematic diagram for the associated vectors and transformations. (Left) The angle between Ωm and Ωp is α. The co-latitude and longitude of Ωf are θf and ϕf, respectively. Note that the magnitudes and orientations of Ωm, Ωp, and Ωf are not to scale. (Middle) Counterclockwise rotation through ϕf about Z ̂ p ′ ⁠. (Right) Counterclockwise rotation through θf about Y ̂ p ′ ′ ⁠.

Close modal
We analytically solve Eq. (4) in its dimensional form together with the continuity equation,
(B1)
(B2)
where Ω f = Ω f k ̂ f is the fluid rotation vector and P is the reduced pressure, including the centrifugal force. In the local coordinate system, the vector k ̂ f is given by cos ⁡ θ z ̂ + sin ⁡ θ y ̂ ⁠, where θ is the co-latitude. Following Eq. (3), the boundary conditions require u = 0 at z = −∞; at z = 0, the velocities are ux = U cos θ cos(Ωft + ϕ), uy = −U sin(Ωft + ϕ), and uz = 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.
Analytical solutions are available to validate our numerical results in the laminar regime. Equation (B1) can be solved analytically by making some assumptions. First, we assume that the system is horizontally uniform, so the solution must be independent of x and y. From the continuity equation ∇ · u = 0, we further have ∂uz/∂z = 0. Given that no vertical motion is allowed on the top boundary, we end up with uz = 0 everywhere. Equation (B1) then becomes
(B3)
(B4)
(B5)
At infinity, given ux = uy = 0, the horizontal pressure gradient should be zero: ∂P/∂x = ∂P/∂y = 0. Since ux 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.
From the boundary conditions, we assume that the solutions for the oscillatory Ekman layer are in the form
(B6)
(B7)
(B8)
where qx(z) and qy(z) are the functions to be solved. Substituting the solutions into Eqs. (B3)–(B5) and solving for qi(z), we end up with
(B9)
(B10)
where C1 = U(cos θ − 1)/2, C2 = U(cos θ + 1)/2, and σ ± = Ω f / ν ( 1 ± 2 cos θ ) 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
(B11)
(B12)
where L ν = ν / Ω f is the thickness of the steady Ekman layer. Given t0 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)],
(B13)
(B14)
The depth of penetration is usually defined as 2 ν / Ω f = 2 L ν ⁠. At critical latitudes where one of σ± becomes zero, the velocity profile shows a non-decaying component; here is the case of θ = 60°,
(B15)
(B16)
That means the flow will penetrate all the way down into the bulk of the fluid and the boundary layer has infinite thickness, which is against the very first assumption of horizontal uniformity—the vertical gradient in velocity is relatively larger than the horizontal ones in the thin boundary layer. This breakdown in the local approximation, already pointed by Bondi and Lyttleton (1953), was shown numerically to spawn internal shear layers, which has negligible contribution to the total dissipation (Hollerbach and Kerswell, 1995).
The shear stress can be computed directly from the velocity. The viscous stress tensor is given by sik = μ(∂iuk + ∂kui), 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 ̂ = z ̂ ⁠, the surface shear stress exerted on the fluid is τi = sijnj, i.e., τx = μ∂zux and τy = μ∂zuy,
(B17)
(B18)
in the local Cartesian coordinate system. Similarly, we give the solutions at the pole and the equator. At θ = 0, we have
(B19)
(B20)
Comparison to the boundary velocity shows that the velocity is leading the stress by 45°. We define the veering angle β 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 ̲ = u x + i u y ⁠, and define the veering angle β as the argument difference between τ ̲ and u ̲ ⁠, i.e., β = A r g [ τ ̲ ] − A r g [ u ̲ ] ⁠. 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
(B21)
(B22)
Comparison with the boundary velocity at θ = 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°.
From this Cartesian model, we compute the viscous torque on a spherical surface 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 × τ = − R τ ϕ θ ̂ + R τ θ ϕ ̂ ⁠. The total torque is given by integrating over the whole spherical surface. After changing to a global Cartesian coordinate system ( X ̂ , Y ̂ , Z ̂ ) ⁠, we have the torque exerted on the core,
(B23)
(B24)
(B25)
The integrals can be computed numerically. In the end, we have
(B26)
where I = (8π/15)ρR5 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 Γ = − I Δ Ω ν Ω f R 2.62 X ̂ p + 0.258 Y ̂ p ⁠. We note that the imposed boundary velocity in our local model implies that the vector of differential rotation is given by Δ Ω = − Δ Ω X ̂ 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)].
The basic theory is given by Eqs. (8) and (9) with ϑ = β [or set C5 = 0; see, e.g., the study by Coleman  (1990)]. To avoid confusion, we label the constants in the study by Sous  (2013) as A ̃ and B ̃ ⁠. Comparing to their Eqs. (12) and (13), we find that
(C1)
and
(C2)
where z0 = 0.11ν/u* for the smooth surface case. Their best fit result gives ( A ̃ = 3.3 , B ̃ = 3.0 ) ⁠, which in our convention is (A = 7.32, B = −2.67).
1.
Bondi
,
H.
 and
Lyttleton
,
R. A.
, “
On the dynamical theory of the rotation of the Earth. II. The effect of precession on the motion of the liquid core
,”
Math. Proc. Cambridge Philos. Soc.
49
,
498
–
515
(
1953
).
https://doi.org/10.1017/s030500410002867x
Google Scholar
Crossref
Search ADS
 
2.
Braun
,
L.
,
Younis
,
B. A.
, and
Weigand
,
B.
, “
A turbulence closure study of the flow and thermal fields in the Ekman layer
,”
Boundary-Layer Meteorol.
175
,
25
–
55
(
2020
).
https://doi.org/10.1007/s10546-019-00495-8
Google Scholar
Crossref
Search ADS
 
3.
Buffett
,
B. A.
, “
Constraints on magnetic energy and mantle conductivity from the forced nutations of the Earth
,”
J. Geophys. Res.: Solid Earth
97
,
19581
–
19597
, (
1992
).
https://doi.org/10.1029/92jb00977
Google Scholar
Crossref
Search ADS
 
4.
Buffett
,
B. A.
, “
Conditions for turbulent Ekman layers in precessionally driven flow
,”
Geophys. J. Int.
226
,
56
–
65
(
2021
).
https://doi.org/10.1093/gji/ggab088
Google Scholar
Crossref
Search ADS
 
5.
Buffett
,
B. A.
 and
Christensen
,
U. R.
, “
Magnetic and viscous coupling at the core-mantle boundary: Inferences from observations of the Earth’s nutations
,”
Geophys. J. Int.
171
,
145
–
152
(
2007
).
https://doi.org/10.1111/j.1365-246x.2007.03543.x
Google Scholar
Crossref
Search ADS
 
6.
Buffett
,
B.
,
Mathews
,
P.
, and
Herring
,
T.
, “
Modeling of nutation and precession: Effects of electromagnetic coupling
,”
J. Geophys. Res.: Solid Earth 107, ETG 5-1-ETG
5-14
, (
2002
).
https://doi.org/10.1029/2000jb000056
Google Scholar
 
7.
Busse
,
F. H.
, “
Steady fluid flow in a precessing spheroidal shell
,”
J. Fluid Mech.
33
,
739
–
751
(
1968
).
https://doi.org/10.1017/s0022112068001655
Google Scholar
Crossref
Search ADS
 
8.
Caldwell
,
D. R.
,
Van Atta
,
C. W.
, and
Helland
,
K. N.
, “
A laboratory study of the turbulent Ekman layer
,”
Geophys. Fluid Dyn.
3
,
125
–
160
(
1972
).
https://doi.org/10.1080/03091927208236078
Google Scholar
Crossref
Search ADS
 
9.
Cébron
,
D.
,
Laguerre
,
R.
,
Noir
,
J.
, and
Schaeffer
,
N.
, “
Precessing spherical shells: Flows, dissipation, dynamo and the lunar core
,”
Geophys. J. Int.
219
,
S34
–
S57
(
2019
).
https://doi.org/10.1093/gji/ggz037
Google Scholar
Crossref
Search ADS
 
10.
Coleman
,
G. N.
, “
Similarity statistics from a direct numerical simulation of the neutrally stratified planetary boundary layer
,”
J. Atmos. Sci.
56
,
891
–
900
(
1999
).
https://doi.org/10.1175/1520-0469(1999)056<0891:ssfadn>2.0.co;2
Google Scholar
Crossref
Search ADS
 
11.
Coleman
,
G. N.
,
Ferziger
,
J. H.
, and
Spalart
,
P. R.
, “
A numerical study of the turbulent Ekman layer
,”
J. Fluid Mech.
213
,
313
–
348
(
1990
).
https://doi.org/10.1017/s0022112090002348
Google Scholar
Crossref
Search ADS
 
12.
Csanady
,
G. T.
, “
On the ‘resistance law’ of a turbulent Ekman layer
,”
J. Atmos. Sci.
24
,
467
–
471
(
1967
).
https://doi.org/10.1175/1520-0469(1967)024<0467:otloat>2.0.co;2
Google Scholar
Crossref
Search ADS
 
13.
Deleplace
,
B.
 and
Cardin
,
P.
, “
Viscomagnetic torque at the core mantle boundary
,”
Geophys. J. Int.
167
,
557
–
566
(
2006
).
https://doi.org/10.1111/j.1365-246x.2006.03180.x
Google Scholar
Crossref
Search ADS
 
14.
Deusebio
,
E.
,
Brethouwer
,
G.
,
Schlatter
,
P.
, and
Lindborg
,
E.
, “
A numerical study of the unstratified and stratified Ekman layer
,”
J. Fluid Mech.
755
,
672
–
704
(
2014
).
https://doi.org/10.1017/jfm.2014.318
Google Scholar
Crossref
Search ADS
 
15.
Fischer
,
P. F.
,
Loth
,
F.
,
Lee
,
S. E.
,
Lee
,
S.-W.
,
Smith
,
D. S.
, and
Bassiouny
,
H. S.
, “
Simulation of high-Reynolds number vascular flows
,”
Comput. Methods Appl. Mech. Eng.
196
,
3049
–
3060
(
2007
).
https://doi.org/10.1016/j.cma.2006.10.015
Google Scholar
Crossref
Search ADS
 
16.
Goldreich
,
P.
, “
Precession of the Moon’s core
,”
J. Geophys. Res.
72
,
3135
–
3137
, (
1967
).
https://doi.org/10.1029/jz072i012p03135
Google Scholar
Crossref
Search ADS
 
17.
H. P.
Greenspan
, “
The theory of rotating fluids
,” in
Technical Report
(
Massachusetts Institute of Technology, Department of Mathematics
,
Cambridge
,
1968
).
Google Scholar
 
18.
Gwinn
,
C. R.
,
Herring
,
T. A.
, and
Shapiro
,
I. I.
, “
Geodesy by radio interferometry: Studies of the forced nutations of the Earth: 2. Interpretation
,”
J. Geophys. Res.: Solid Earth
91
,
4755
–
4765
, (
1986
).
https://doi.org/10.1029/jb091ib05p04755
Google Scholar
Crossref
Search ADS
 
19.
Herring
,
T. A.
,
Gwinn
,
C. R.
, and
Shapiro
,
I. I.
, “
Geodesy by radio interferometry: Studies of the forced nutations of the Earth: 1. Data analysis
,”
J. Geophys. Res.: Solid Earth
91
,
4745
–
4754
, (
1986
).
https://doi.org/10.1029/jb091ib05p04745
Google Scholar
Crossref
Search ADS
 
20.
Hollerbach
,
R.
 and
Kerswell
,
R. R.
, “
Oscillatory internal shear layers in rotating and precessing flows
,”
J. Fluid Mech.
298
,
327
–
339
(
1995
).
https://doi.org/10.1017/s0022112095003338
Google Scholar
Crossref
Search ADS
 
21.
Ichikawa
,
H.
 and
Tsuchiya
,
T.
, “
Atomic transport property of Fe–O liquid alloys in the Earth’s outer core P, T condition
,”
Phys. Earth Planet. Inter.
247
,
27
–
35
(
2015
).
https://doi.org/10.1016/j.pepi.2015.03.006
Google Scholar
Crossref
Search ADS
 
22.
Landau
,
L. D.
 and
Lifshitz
,
E. M.
,
Fluid Mechanics: Landau and Lifshitz: Course of Theoretical Physics
(
Elsevier
,
2013
), Vol.
6
.
Google Scholar
 
23.
Lorenzani
,
S.
 and
Tilgner
,
A.
, “
Fluid instabilities in precessing spheroidal cavities
,”
J. Fluid Mech.
447
,
111
–
128
(
2001
).
https://doi.org/10.1017/s002211200100581x
Google Scholar
Crossref
Search ADS
 
24.
Lums
,
L. I.
 and
Aldridge
,
K. D.
, “
On viscosity estimates for the Earth’s fluid outer core and core-mantle coupling
,”
J. Geomagn. Geoelectr.
43
,
93
–
110
(
1991
).
https://doi.org/10.5636/jgg.43.93
Google Scholar
Crossref
Search ADS
 
25.
Malkus
,
W. V. R.
, “
Precession of the Earth as the Cause of Geomagnetism: Experiments lend support to the proposal that precessional torques drive the Earth’s dynamo
,”
Science
160
,
259
–
264
(
1968
).
https://doi.org/10.1126/science.160.3825.259
Google Scholar
Crossref
Search ADS
PubMed
 
26.
Marlatt
,
S.
,
Waggy
,
S.
, and
Biringen
,
S.
, “
Direct numerical simulation of the turbulent Ekman layer: Evaluation of closure models
,”
J. Atmos. Sci.
69
,
1106
–
1117
(
2012
).
https://doi.org/10.1175/jas-d-11-0107.1
Google Scholar
Crossref
Search ADS
 
27.
Mathews
,
P.
 and
Guo
,
J.
, “
Viscoelectromagnetic coupling in precession-nutation theory
,”
J. Geophys. Res.: Solid Earth
110
,
B02402
, (
2005
).
https://doi.org/10.1029/2003jb002915
Google Scholar
Crossref
Search ADS
 
28.
Mathews
,
P. M.
,
Herring
,
T. A.
, and
Buffett
,
B. A.
, “
Modeling of nutation and precession: New nutation series for nonrigid Earth and insights into the Earth’s interior
,”
J. Geophys. Res.: Solid Earth 107, ETG 3-1-ETG
3-26
, (
2002
).
https://doi.org/10.1029/2001jb000390
Google Scholar
 
29.
Noir
,
J.
 and
Cébron
,
D.
, “
Precession-driven flows in non-axisymmetric ellipsoids
,”
J. Fluid Mech.
737
,
412
–
439
(
2013
).
https://doi.org/10.1017/jfm.2013.524
Google Scholar
Crossref
Search ADS
 
30.
Ozdemir
,
C. E.
,
Hsu
,
T.-J.
, and
Balachandar
,
S.
, “
Direct numerical simulations of transition and turbulence in smooth-walled Stokes boundary layer
,”
Phys. Fluids
26
,
045108
(
2014
).
https://doi.org/10.1063/1.4871020
Google Scholar
Crossref
Search ADS
 
31.
Pizzi
,
F.
,
Giesecke
,
A.
, and
Stefani
,
F.
, “
Ekman boundary layers in a fluid filled precessing cylinder
,”
AIP Adv.
11
,
035023
(
2021
).
https://doi.org/10.1063/5.0037922
Google Scholar
Crossref
Search ADS
 
32.
Pozzo
,
M.
,
Davies
,
C.
,
Gubbins
,
D.
, and
Alfe
,
D.
, “
Transport properties for liquid silicon-oxygen-iron mixtures at Earth’s core conditions
,”
Phys. Rev. B
87
,
014110
(
2013
).
https://doi.org/10.1103/physrevb.87.014110
Google Scholar
Crossref
Search ADS
 
33.
Shingai
,
K.
 and
Kawamura
,
H.
, “
A study of turbulence structure and large-scale motion in the Ekman layer through direct numerical simulations
,”
J. Turbul.
5
,
013
(
2004
).
https://doi.org/10.1088/1468-5248/5/1/013
Google Scholar
Crossref
Search ADS
 
34.
Sous
,
D.
,
Sommeria
,
J.
, and
Boyer
,
D.
, “
Friction law and turbulent properties in a laboratory Ekman boundary layer
,”
Phys. Fluids
25
,
046602
(
2013
).
https://doi.org/10.1063/1.4802045
Google Scholar
Crossref
Search ADS
 
35.
Spalart
,
P. R.
, “
Theoretical and numerical study of a three-dimensional turbulent boundary layer
,”
J. Fluid Mech.
205
,
319
–
340
(
1989
).
https://doi.org/10.1017/s0022112089002053
Google Scholar
Crossref
Search ADS
 
36.
Stys
,
C.
 and
Dumberry
,
M.
, “
A past lunar dynamo thermally driven by the precession of its inner core
,”
J. Geophys. Res.: Planets
125
,
e2020JE006396
, (
2020
).
https://doi.org/10.1029/2020je006396
Google Scholar
Crossref
Search ADS
 
37.
Tennekes
,
H.
,
Lumley
,
J. L.
,
Lumley
,
J. L.
 et al,
A First Course in Turbulence
(
MIT Press
,
1972
).
Google Scholar
Crossref
Search ADS
 
38.
Thornton
,
S. T.
 and
Marion
,
J. B.
,
Classical Dynamics of Particles and Systems
(
Cengage Learning
,
2003
).
Google Scholar
 
39.
Tilgner
,
A.
, “
Rotational dynamics of the core
,”
Treatise Geophys.
8
,
183
–
212
(
2015
).
https://doi.org/10.1016/b978-0-444-53802-4.00140-8
Google Scholar
Crossref
Search ADS
 
40.
Triana
,
S. A.
,
Trinh
,
A.
,
Rekier
,
J.
,
Zhu
,
P.
, and
Dehant
,
V.
, “
The viscous and ohmic damping of the Earth’s free core nutation
,”
J. Geophys. Res.: Solid Earth
126
,
e2020JB021042
, (
2021
).
https://doi.org/10.1029/2020JB021042
Google Scholar
Crossref
Search ADS
PubMed
 
41.
Viswanathan
,
V.
,
Rambaux
,
N.
,
Fienga
,
A.
,
Laskar
,
J.
, and
Gastineau
,
M.
, “
Observational constraint on the radius and oblateness of the lunar core-mantle boundary
,”
Geophys. Res. Lett.
46
,
7295
–
7303
, (
2019
).
https://doi.org/10.1029/2019gl082677
Google Scholar
Crossref
Search ADS
 
42.
Williams
,
J. G.
 and
Boggs
,
D. H.
, “
Tides on the Moon: Theory and determination of dissipation
,”
J. Geophys. Res.: Planets
120
,
689
–
724
, (
2015
).
https://doi.org/10.1002/2014je004755
Google Scholar
Crossref
Search ADS
 
43.
Williams
,
J. G.
,
Boggs
,
D. H.
,
Yoder
,
C. F.
,
Ratcliff
,
J. T.
, and
Dickey
,
J. O.
, “
Lunar rotational dissipation in solid body and molten core
,”
J. Geophys. Res.: Planets
106
,
27933
–
27968
, (
2001
).
https://doi.org/10.1029/2000je001396
Google Scholar
Crossref
Search ADS
 
44.
Williams
,
J. G.
,
Konopliv
,
A. S.
,
Boggs
,
D. H.
,
Park
,
R. S.
,
Yuan
,
D.-N.
,
Lemoine
,
F. G.
,
Goossens
,
S.
,
Mazarico
,
E.
,
Nimmo
,
F.
,
Weber
,
R. C.
 et al, “
Lunar interior properties from the grail mission
,”
J. Geophys. Res.: Planets
119
,
1546
–
1578
, (
2014
).
https://doi.org/10.1002/2013je004559
Google Scholar
Crossref
Search ADS
 
45.
Yoder
,
C. F.
, “
The free librations of a dissipative Moon
,”
Philos. Trans. R. Soc. London, Ser. A
303
,
327
–
338
(
1981
).
https://doi.org/10.1098/rsta.1981.0206
Google Scholar
Crossref
Search ADS
 
46.
Zhang
,
J.
 and
Dumberry
,
M.
, “
Viscous dissipation in the fluid core of the Moon
,”
J. Geophys. Res.: Planets
126
,
e2021JE006966
, (
2021
).
https://doi.org/10.1029/2021je006966
Google Scholar
Crossref
Search ADS
 
© 2023 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Supplementary Material

Supplementary materials- zip file