1. Introduction
Buoyancy-driven convection in a rotating spherical shell is a classical framework for studying the dynamics and magnetic field generation in the cores of planets and stars. However, despite significant progress (Gastine et al. Reference Gastine, Wicht and Aubert2016; Schaeffer et al. Reference Schaeffer, Jault, Nataf and Fournier2017; Long et al. Reference Long, Mound, Davies and Tobias2020; Schwaiger et al. Reference Schwaiger, Gastine and Aubert2020; Gastine & Aurnou Reference Gastine and Aurnou2023), state-of-the-art direct numerical simulations cannot reach the parameter values that are representative of the astrophysical bodies, and are unlikely to do so in the near future (Davies et al. Reference Davies, Gubbins and Jimack2011; Roberts & King Reference Roberts and King2013). Therefore, extensive work has focused on developing scaling relations between the governing input parameters and global output diagnostic quantities of numerical simulations, which can be used to extrapolate to the conditions of planetary cores, and facilitate comparisons with available observations. These scaling relations rely on the balance of forces that determine the system dynamics (e.g. Christensen et al. Reference Christensen, Aubert and Hulot2010; King & Buffett Reference King and Buffett2013; Aubert et al. Reference Aubert, Gastine and Fournier2017), so it is crucial to quantify dynamical balances in numerical simulations accurately.
The dynamics in rotating non-magnetic convection are governed by the Ekman number
$E$
, a measure of the ratio of viscous to Coriolis forces, the Prandtl number
$Pr$
, the ratio of viscous to thermal diffusivities, and the Rayleigh number
$Ra_{T}$
, measuring the buoyancy force driving convection. Theoretical considerations suggest that the convecting system exhibits at least three distinct dynamical regimes (e.g. Aubert et al. Reference Aubert, Brito, Nataf, Cardin and Masson2001; King & Buffett Reference King and Buffett2013; King et al. Reference King, Stellmach and Buffett2013; Gastine et al. Reference Gastine, Wicht and Aubert2016; Aurnou et al. Reference Aurnou, Horn and Julien2020; Long et al. Reference Long, Mound, Davies and Tobias2020; Kunnen Reference Kunnen2021). At fixed
$E$
and
$Pr$
, raising the thermal forcing to just above the onset of convection (i.e.
$Ra_T\geqslant Ra_T^{c}$
, where
$Ra_T^{c}\propto E^{4/3}$
is the critical Rayleigh number at the onset of convection) leads to weakly nonlinear (WN) convection. In this regime, the primary dynamics is expected to be a quasi-geostrophic (QG) balance between Coriolis and pressure forces, while the residual (i.e. the ageostrophic Coriolis force) is balanced by buoyancy and viscous forces forming a secondary Viscous–Archimedean (buoyancy)–Coriolis (ageostrophic) or VAC balance. Increasing the thermal forcing increases the role of inertia, which leads to a turbulent rapidly rotating (RR) flow regime where strong nonlinearity and rotational constraint coexist. In this regime, both viscosity and inertial forces are expected to be small relative to the Coriolis force in the primary QG balance, with an Inertia–Archimedian–Coriolis (ageostrophic) or IAC balance generally assumed among the secondary forces. Further increase in thermal forcing leads to the weakly rotating (WR) regime where the primary QG balance is gradually broken, eventually resulting in non-rotating behaviour at sufficiently high
$Ra_T$
. Among these three regimes, the RR regime is thought to be more relevant for investigating planetary core convection, as compared to the other regimes (Kunnen Reference Kunnen2021).
Depending on the dynamical balance, global and local flow diagnostic quantities (e.g. average velocity, length scale, heat transport, boundary layer thickness) may exhibit distinct scaling with the governing input parameters. Scaling regimes in rotating spherical shell convection have been studied extensively by Gastine et al. (Reference Gastine, Wicht and Aubert2016) and Long et al. (Reference Long, Mound, Davies and Tobias2020) (henceforth referred to as L20). The conformity of the simulation diagnostics with theoretical scaling laws was utilized to demarcate boundaries between the dynamical regimes (i.e. WN, RR, WR). For example, Gastine et al. (Reference Gastine, Wicht and Aubert2016) and L20 defined the WN regime based on the theoretical expectation that at low supercriticality, the Nusselt number
$Nu$
, representing the ratio of the total average heat flux from the shell to the conductive flux, scales as
$Nu-1\propto Ra_{T}/Ra_T^{c} -1$
(Gillet & Jones Reference Gillet and Jones2006). They found that this heat transfer behaviour is consistent with the dimensionless flow length scale scaling
$\ell \sim E^{1/3}$
and Reynolds number scaling
$Re \sim B^{1/2} E^{1/3}$
(where
$B$
is the convective power), as expected from a VAC balance. For higher supercriticality, a length scale scaling
$\ell \sim Ro^{1/2}$
(where
$Ro$
, the Rossby number, is a measure of the ratio of inertia and Coriolis forces) and a Reynolds number scaling
$Re \sim B^{2/5} E^{1/5}$
are approached, indicating an IAC balance. The asymptotically reduced models of Guervilly et al. (Reference Guervilly, Cardin and Schaeffer2019) also found a ‘diffusion-free’
$\ell \sim Ro^{1/2}$
scaling. An asymptotically reduced plane layer model of convection (Oliver et al. Reference Oliver, Jacobi, Julien and Calkins2023) found the scaling of
$Nu$
and
$Re$
consistent with IAC predictions, while the length scale remained viscously controlled (i.e.
$\ell \sim E^{1/3}$
). A rotating cylinder experiment (Abbate & Aurnou Reference Abbate and Aurnou2023), however, was unable to distinguish between a viscously controlled regime and a ‘diffusion-free’ regime from the scaling of the length scale. These comparisons relied on assumed force balances; however, the quantitative behaviour of forces in numerical simulations of spherical shell rotating convection (RC) remains largely unexplored. A few studies have explored force balances in plane layer geometries (Guzmán et al. Reference Guzmán, Madonia, Cheng, Ostilla-Mónico, Clercx and Kunnen2020; Naskar & Pal Reference Naskar and Pal2022a
,Reference Naskar and Pal
b
), which may represent the dynamics in the tangent cylinder region of a spherical shell (Gastine & Aurnou Reference Gastine and Aurnou2023), while Schwaiger et al. (Reference Schwaiger, Gastine and Aubert2020), Teed & Dormy (Reference Teed and Dormy2023) and Nicoski et al. (Reference Nicoski, O'Connor and Calkins2024) have computed force balances in a limited number of spherical shell RC runs. Hence the primary objective of our study is to investigate the dynamical balances that emerge in numerical simulations of spherical shell RC.
Quantifying dynamical balances in RC is an intricate issue. The force representation in incompressible flows is complicated by the existence of the gradient portions of forces, which are balanced by the pressure gradient term, but are not directly relevant to the dynamics (Hughes & Cattaneo Reference Hughes and Cattaneo2019). A simple way to remove the gradient contributions is to curl the force balance (Dormy Reference Dormy2016), though the derivative operation can enhance small-scale contributions to the individual terms (Teed & Dormy Reference Teed and Dormy2023). Differences between the dynamics predicted by force and curled balances are important since some theoretical predictions consider the asymptotic behaviour of forces (Nicoski et al. Reference Nicoski, O'Connor and Calkins2024), whereas others are based on curled forces (e.g. VAC and IAC in L20). Teed & Dormy (Reference Teed and Dormy2023) found that forces and their curls predicted different dynamical balances in a single simulation of rotating spherical shell convection. Here, we investigate the consistency between force and curled depictions of RC dynamics across a broad range of parameters (
$Ra_{T}/Ra_{T}^{c}=1.2{-}297$
,
$E=10^{-3}{-}10^{-6}$
and
$Pr=1$
).
When calculating forces or curls, it is crucial to distinguish dynamics in the boundary layers and the convective bulk. Integrating forces over the entire spherical shell may overestimate the role of viscosity in the bulk dynamics (Soderlund et al. Reference Soderlund, King and Aurnou2012; Yadav et al. Reference Yadav, Gastine, Christensen, Wolk and Poppenhaeger2016). Most studies of force calculations remove a region corresponding to one velocity boundary layer (VBL) thickness at the top and bottom of the domain (e.g. Schwaiger et al. Reference Schwaiger, Gastine and Aubert2020; Teed & Dormy Reference Teed and Dormy2023; Nicoski et al. Reference Nicoski, O'Connor and Calkins2024) in order to isolate the bulk dynamics. However, we are unaware of previous systematic studies that have established the thickness of the layer near the boundaries that should be excluded from the volume-averaged forces/curls to obtain a robust estimate of the bulk dynamics.
In rotating spherical shell convection, the balance in the axisymmetric part of forces (corresponding to spherical harmonic order
$m=0$
) can be different from the balance in their non-axisymmetric counterparts. Nicoski et al. (Reference Nicoski, O'Connor and Calkins2024) partitioned the forces into azimuthally averaged and corresponding fluctuating parts, and analysed the scaling behaviour of the fluctuating radial forces as a function of
$E$
and
$Ra_{T}$
in RC with strong zonal flows, finding a primary QG balance across a broad range of parameters. In contrast to the QG force balance in the small-scale (i.e.
$m\neq 0$
) convective motions (Nicoski et al. Reference Nicoski, O'Connor and Calkins2024), Aubert (Reference Aubert2005) focused on the large-scale dynamics and used curls to identify a horizontal thermal wind (TW) balance in the large-scale azimuthal (
$m=0$
) motions for RC simulations operating in the RR regime. Here, we systematically compare mean and fluctuating components of forces and their curls across different dynamical regimes of spherical shell RC.
In a spherical shell geometry, an additional difficulty in representing the system dynamics arises due to the dependence of each force and curl component (
$\hat {r},\hat {\theta }, \hat {\phi }$
) on the spatio-temporal coordinates (
$r,\theta ,\phi ,t$
) and scale (spherical harmonic degree
$l$
and order
$m$
), as well as the governing parameters (
$Ra_{T}$
,
$E$
,
$Pr$
). Since the buoyancy force is radial while the system rotates about a vertical axis, different dynamical balances may be expected in the three orthogonal directions. Directional anisotropy of the force balance has been considered rarely in spherical shell models (for an exception, see Calkins et al. Reference Calkins, Orvedahl and Featherstone2021), and forces are usually represented by their total magnitudes (Soderlund et al. Reference Soderlund, King and Aurnou2012; Guzmán et al. Reference Guzmán, Madonia, Cheng, Ostilla-Mónico, Clercx and Kunnen2021; Orvedahl et al. Reference Orvedahl, Featherstone and Calkins2021). Studies of azimuthally averaged forces (Calkins et al. Reference Calkins, Orvedahl and Featherstone2021) and curls (Aubert Reference Aubert2005) in dynamo simulations have indeed revealed a TW balance in the meridional plane with a Coriolis–Lorentz balance in the azimuthal direction. We complement these works by considering the balance of individual (
$\hat {r},\hat {\theta }, \hat {\phi }$
) force components in rotating spherical shell convection.
The scale-dependence of the force balance in RC spherical shell simulations has been investigated by considering the force contributions from each spherical harmonic degree
$l$
. Schwaiger et al. (Reference Schwaiger, Gastine and Aubert2020) found a secondary balance between ageostrophic Coriolis and buoyancy forces at low
$l$
, and a cross-over to Coriolis–Inertia balance at higher
$l$
. They argued that the length scale at which inertia and buoyancy forces cross over is related to energetic scales of the flow, as estimated from the peak of the poloidal kinetic energy spectra. However, the scale-dependent representation of curled forces may not exhibit such a dynamically relevant cross-over scale (Teed & Dormy Reference Teed and Dormy2023). Furthermore, some simulations show that the small-scale curled balance differs markedly compared to the force balance (Teed & Dormy Reference Teed and Dormy2023). It is therefore useful to systematically compare the scale-dependent nature of force and curled balances in different dynamical regimes of RC.
In summary, previous studies have investigated various aspects of dynamical balances (e.g. forces versus curls, mean versus fluctuating balances, components versus combined terms, scale-dependence) in spherical shell RC. However, to our knowledge, no study systematically compares the different depictions of dynamical balances in different regimes of spherical shell RC. A recent related study by Nicoski et al. (Reference Nicoski, O'Connor and Calkins2024) investigated the behaviour of fluctuating radial forces in spherical shell RC with fixed temperature and free-slip boundaries for
$E=5\times 10^{-4}{-}5\times 10^{-7}$
,
$Pr=1$
and the convective supercriticality range
$Ra_{T}/Ra_{T}^{c}=3{-}161$
. They found a primary QG balance in the radial fluctuating forces at all parameters considered. At low convective supercriticalities, they obtained a secondary balance between buoyancy and ageostrophic Coriolis forces, while inertia and viscosity enter the balance at higher thermal forcing. In this paper, we systematically compare integrated and scale-dependent representations of mean and fluctuating forces and curled balances in the large suite of rotating spherical shell convection simulations from Mound & Davies (Reference Mound and Davies2017) and L20, supplemented by three new simulations. In contrast to Nicoski et al. (Reference Nicoski, O'Connor and Calkins2024), our simulations use no-slip, fixed heat flux boundaries.
This paper is organized as follows. We present the mathematical model in § 2.1, the numerical details in § 2.2, and the force calculation method in § 2.3. Results are presented in § 3, beginning in § 3.1 with an analysis of the boundary region that must be excluded to ensure a robust representation of the bulk dynamics. In §§ 3.2 and 3.3, we investigate mean and fluctuating bulk-integrated force and curl balances, respectively. The total magnitude of the forces and their scale dependence are discussed in § 3.4. In § 4, we compare transitions in the computed force and curled balanced to the regime diagram obtained by L20. We discuss our observations and summarize the findings in § 5.
2. Method
2.1. Mathematical model
We employ a numerical model of convection of a Boussinesq fluid in a rotating spherical shell. The relevant physical properties of the fluid are the kinematic viscosity
$\nu$
, thermal expansivity
$\alpha$
, and thermal diffusivity
$\kappa$
, defined as
$\kappa =k/\rho _0 c_p$
, where
$k$
is the thermal conductivity,
$\rho _0$
is the reference density, and
$c_p$
is the specific heat capacity. A spherical coordinate system
$(r,\theta ,\phi )$
is used to represent the domain bounded by the inner and outer boundaries,
$r_i$
and
$r_o$
, respectively. The system rotates with a constant angular velocity
$\boldsymbol {\Omega }=\Omega \hat {\boldsymbol {z}}$
, about the vertical, and gravity
$g$
varies linearly with radius, with
$g=g_o$
at the outer radius. The governing equations are cast in non-dimensional form using the shell gap
$h$
as the length scale, the viscous diffusion time
$h^2/\nu$
as the time scale, and
$\beta /h$
as the temperature scale, giving
\begin{align} \boldsymbol {\nabla \cdot u}=0, \qquad\qquad\qquad\qquad\qquad\end{align}
\begin{align} \frac {\partial \boldsymbol {u}}{\partial t}+ (\boldsymbol {u}\boldsymbol {\cdot} \boldsymbol {\nabla })\boldsymbol {u}+\frac {1}{E}\left (\hat {\boldsymbol {z}}\times \boldsymbol {u}\right ) =-\boldsymbol {\nabla }\widetilde {P}+\left (\frac {Ra}{Pr}\right )T\boldsymbol {r}+ {\nabla ^{2}}\boldsymbol {u}, \end{align}
\begin{align} \frac {\partial T}{\partial t}+ (\boldsymbol {u}\boldsymbol {\cdot} \boldsymbol {\nabla })(T+T_c) = \frac {1}{Pr}\,{\nabla ^{2}}(T+T_c), \quad\qquad\end{align}
where
$T$
is the temperature fluctuation relative to the conductive state
$T_c$
, given by
$\partial T_c/\partial r=-\beta /r^2$
. The parameter
$\beta$
is related to the fixed heat flow through the boundaries as
$Q=4{\unicode[Arial]{x03C0}} \beta k$
. No-slip velocity conditions and fixed heat flux temperature conditions have been used at both boundaries.
The non-dimensional numbers appearing in these equations are the Ekman number (
$E$
), Rayleigh number (
$Ra$
), and Prandtl number (
$Pr$
), which are defined as
\begin{equation} E = \frac {\nu }{2\Omega h^2}, \quad Ra=\frac {g_o\alpha \beta h^3}{\nu \kappa r^{*}_{o}}, \quad Pr=\frac {\nu }{\kappa } . \end{equation}
Here,
$r^{*}_{o}=1/(1-\eta )$
is the non-dimensional outer radius, where
$\eta =r_i/r_o$
is the radius ratio fixed at
$\eta =0.35$
.
2.2. Numerical details
The velocity field is represented by toroidal and poloidal scalar fields, which are expressed as radially varying Schmidt-normalized spherical harmonics. Radial variations are expressed using second-order finite differences on the zeros of Chebyshev polynomials. A predictor–corrector scheme is used for time stepping in spectral space that treats the diffusion terms implicitly. Further numerical details can be found in previous studies that use the same solver (Willis et al. Reference Willis, Sreenivasan and Gubbins2007; Davies et al. Reference Davies, Gubbins and Jimack2011; Matsui et al. Reference Matsui2016).
We consider all the simulations reported in Mound & Davies (Reference Mound and Davies2017) for
$E=10^{-4}$
,
$10^{-5}$
and
$10^{-6}$
, that impose homogeneous heat flux at the outer boundary. The simulations performed by L20 at
$E=10^{-3}$
,
$E=3\times 10^{-4}$
and
$E=3\times 10^{-5}$
are also included. It should be noted here that the flux Rayleigh number in this study relates to the modified flux Rayleigh number defined by Mound & Davies (Reference Mound and Davies2017) and L20 as
$\widetilde {Ra}=Ra\,E$
. To complement the database, we have run three more simulations at
$E=10^{-6}$
, for
$\widetilde {Ra}=350, 550, 30\,000$
. The details of these simulations are given in Appendix A. The Prandtl number is fixed at
$Pr=1$
for all runs. The fixed temperature Rayleigh number (
$Ra_T$
) is related to the fixed flux Rayleigh number (
$Ra$
) as
\begin{align} Ra_T=\frac {Ra}{Nu}\,\frac {(1-\eta )^2}{\eta }. \end{align}
The reduced fixed-temperature Rayleigh number used in § 3 is defined as
$\widetilde {Ra}_T=Ra_{T}\,E^{4/3}$
. The critical values for the three definitions of Rayleigh numbers (
$\widetilde {Ra}^{c}$
,
$Ra_{T}^{c}$
and
$\widetilde {Ra}_{T}^{c}$
) at the onset of thermal convection are provided in Appendix B, based on the analysis of L20.
2.3. Force calculation
We refer to the terms from left to right in (2.2) as time derivative (TD), inertia (I), Coriolis (C), pressure gradient (P), thermal buoyancy (or Archimedean, A) and viscous (V) forces. Additionally, the ageostrophic Coriolis force
$({1}/{E}) (\hat {\boldsymbol {z}}\times \boldsymbol {u} )+\boldsymbol {\nabla }\widetilde {P}$
is denoted as
$\text{C}_{ag}$
. Henceforth, we use these abbreviations to refer to the balance in the simulations. For example, a balance between inertia, Archimedean (i.e. thermal buoyancy) and Coriolis forces will be referred to as an IAC balance. Also, a TW balance refers to an ACP force balance or an AC balance of curled forces. Similarly, a QG balance refers to a CP balance of forces.
We partition dependent variables into their azimuthally averaged mean and corresponding fluctuating parts as
\begin{align} f(r,\theta ,\phi ,t)&=\overline {f}(r,\theta ,t)+f ^{\prime}(r,\theta ,\phi ,t), \notag \\ \overline {f}(r,\theta ,t)&=\frac {1}{2{\unicode[Arial]{x03C0}} }\int _{0}^{2{\unicode[Arial]{x03C0}} }f(r,\theta ,\phi ,t)\,{\rm d}\phi. \end{align}
In spectral representation, this is equivalent to partitioning into the azimuthally symmetric harmonics of order
$m=0$
(mean) and the corresponding part with order
$m\neq 0$
. Azimuthally averaging (2.2) leads to
\begin{align} \underbrace {\frac {\partial \overline {\boldsymbol {u}}}{\partial t}}_{\overline {TD}}+ \underbrace {(\overline {\boldsymbol {u}}\cdot {\boldsymbol {\nabla }})\overline {\boldsymbol {u}}}_{I^{mm}}+\underbrace {(\overline {\boldsymbol {u}^{\prime}\cdot \boldsymbol {\nabla })\boldsymbol {u}^{\prime}}}_{\overline {I^{ff}}}+\underbrace {\frac {1}{E}\left (\hat {\boldsymbol {z}}\times \overline {\boldsymbol {u}}\right )}_{\overline {C}} =\underbrace {-\boldsymbol {\nabla }\overline {P}}_{\overline {P}}+\underbrace {\frac {Ra}{Pr}\overline {T}}_{\overline {A}}\boldsymbol {r}+\underbrace {{\nabla ^{2}}\overline {\boldsymbol {u}}}_{\overline {V}}, \end{align}
where partitioning the mean inertial force term
$\overline {I}= {\overline {(\boldsymbol {u}\cdot \boldsymbol {\nabla} ) u}}$
as
$\overline {I}=I^{mm}+\overline {I^{ff}}$
leads to the mean–mean inertia term
$I^{mm}=(\overline {\boldsymbol {u}}\cdot \boldsymbol {\nabla })\overline {\boldsymbol {u}}$
and the Reynolds stress
$\overline {I^{ff}}=(\overline {\boldsymbol {u}^{\prime}\cdot \boldsymbol {\nabla })\boldsymbol {u}^{\prime}}$
. Subtracting the mean momentum equation (2.7) from (2.2) leads to the corresponding fluctuating part of the momentum equation:
\begin{align} \underbrace {\frac {\partial \boldsymbol {u}^{\prime}}{\partial t}}_{TD^{\prime}}+& \underbrace {(\overline {\boldsymbol {u}}\cdot \boldsymbol {\nabla })\boldsymbol {u}^{\prime}}_{I^{mf}}+ \underbrace {(\boldsymbol {u}^{\prime}\cdot \boldsymbol {\nabla })\overline {\boldsymbol {u}}}_{I^{fm}}+ \underbrace {(\boldsymbol {u}^{\prime}\cdot \boldsymbol {\nabla })\boldsymbol {u}^{\prime}}_{I^{ff}}-\underbrace {(\overline {\boldsymbol {u}^{\prime}\cdot \boldsymbol {\nabla })\boldsymbol {u}^{\prime}}}_{\overline {I^{ff}}}+\underbrace {\frac {1}{E}\left (\hat {\boldsymbol {z}}\times \boldsymbol {u}^{\prime}\right )}_{C^{\prime}}\notag \\ &=\underbrace {-\boldsymbol {\nabla }P^{\prime}}_{P^{\prime}}+\underbrace {\frac {Ra}{Pr}T^{\prime}}_{A^{\prime}}\boldsymbol {r}+\underbrace {{\nabla ^{2}}\boldsymbol {u}^{\prime}}_{V^{\prime}}, \end{align}
where the fluctuating inertial term
$I^{\prime}$
has four parts,
$I^{\prime}=I^{mf}+I^{fm}+I^{ff}-\overline {I^{ff}}$
, signifying the mean–fluctuating, fluctuating–mean, fluctuating–fluctuating inertial terms and the Reynolds stress. Among the fluctuating inertial terms in (2.8), the fluctuating–fluctuating term (
$I^{ff}$
) always dominates in our simulations. Therefore, we have discussed only the sum of the four terms
$I^{\prime}$
, for presentational convenience. The spherical harmonic degree
$l=0$
of the mean radial forces and
$I^{ff}$
are balanced by the hydrostatic part of the pressure gradient, which is irrelevant to the convective dynamics (Calkins et al. Reference Calkins, Orvedahl and Featherstone2021; Nicoski et al. Reference Nicoski, O'Connor and Calkins2024). Therefore, we remove the degree
$l=0$
from the radial component of
$I^{ff}$
from (2.8), and from all the radial mean forces in (2.7). In our calculations, the temporally averaged values of TD terms in (2.7) and (2.8) remain small compared to all the other terms, though they may be significant for simulations with free-slip boundaries (Nicoski et al. Reference Nicoski, O'Connor and Calkins2024). Therefore, this term is not discussed in the next section, and the nonlinear advection terms are referred to as inertia throughout the text. Also, the abbreviations of mean (fluctuating) forces are mentioned explicitly rather than using an overbar (or prime). For example, a VAC balance in the fluctuating curled forces will be referred to as a curled fluctuating VAC balance, rather than a
$\text{V}^{\prime}\text{A}^{\prime}\text{C}^{\prime}$
balance of curled forces.
The total force magnitude can be calculated from the vector components as
\begin{equation} f_{tot}=|\boldsymbol {f}|= \sqrt {f_r^{2}+f_\theta ^{2}+f_\phi ^{2}}, \end{equation}
where
$f_r$
,
$f_\theta$
and
$f_\phi$
are the components in the
$\hat {r}$
,
$\hat {\theta }$
and
$\hat {\phi }$
directions, respectively. Further, the force components (or the total force magnitude) are represented by their root mean square (r.m.s.) values, where the ‘mean’ refers to a volume average performed by first averaging over a spherical surface and then averaging in the radial direction excluding regions near the boundary:
\begin{align} \langle f_j \rangle _{\mathcal {S}}=\frac {1}{4{\unicode[Arial]{x03C0}} r^{2}}\int _{0}^{{\unicode[Arial]{x03C0}} }\int _{0}^{2{\unicode[Arial]{x03C0}} }f_j(r,\theta ,\phi ,t)\, r^2 \sin\theta \, {\rm d}\phi \, {\rm d}\theta , \end{align}
\begin{align} \langle f_j \rangle _{\mathcal {V}}=\frac {1}{h-r_{ex}\delta ^{i}_v-r_{ex}\delta ^{o}_v}\int _{r_i+r_{ex}\delta ^{i}_v}^{r_o-r_{ex}\delta ^{o}_v}\langle f_j \rangle _{\mathcal {S}}\ {\rm d}r, \end{align}
where
$j=r,\theta ,\phi \ \text {or}\ \textit {tot}$
. In (2.10b
),
$r_{ex}$
represents the multiple of the VBL thicknesses at the inner and the outer boundaries,
$\delta ^{i}_v$
and
$\delta ^{o}_v$
, respectively, that we exclude from the volume average to ensure the representation of bulk dynamics. There are two common methods of estimating the thickness of the VBL from the radial profile of horizontal velocity (L20). One is the ‘local maxima method’ (L20), where the distance of the nearest maximum in this profile from the respective boundaries is defined as the VBL thickness, whereas the distance of the intersection of the tangent to the profile at the wall and the tangent at the nearest maximum near the wall is used as VBL thickness in the ‘linear intersection method’ (Gastine et al. Reference Gastine, Wicht and Aubert2016). In our study, VBL thickness is estimated using the linear intersection method, which has been reported as a better estimate than the local maxima method (Gastine et al. Reference Gastine, Wicht and Aurnou2015). The value of
$r_{ex}$
is determined in § 3.1.
We evaluate the scale dependence of the forces as a function of spherical harmonic degree following Aubert et al. (Reference Aubert, Gastine and Fournier2017). The volume-averaged r.m.s. forces are averaged in time over at least a hundred advective time units after the simulations reach a statistically stationary state (Mound & Davies Reference Mound and Davies2017; Long et al. Reference Long, Mound, Davies and Tobias2020). The calculation of forces follows the methodology reported by Calkins et al. (Reference Calkins, Orvedahl and Featherstone2021) and Nicoski et al. (Reference Nicoski, O'Connor and Calkins2024), and a couple of cases from Nicoski et al. (Reference Nicoski, O'Connor and Calkins2024) have been reproduced to validate our force calculation.
3. Results
In this section, we investigate the force balances in the thermally driven RC simulations at
$E=10^{-5}$
reported in Mound & Davies (Reference Mound and Davies2017) and L20. We start, in § 3.1, with a systematic study to assess the thickness of the region near the boundaries that should be excluded to ensure the representation of bulk dynamics. We explore the mean and fluctuating parts of each component of forces in § 3.2, and their curl in § 3.3. The scale dependence of the forces is investigated in § 3.4, whereas the role of inertia in the force balance in various regimes of RC is considered in § 4. The thermal forcings at which transitions occur from WN to RR to WR behaviours have been predicted previously based on theoretical scaling laws (L20) and are marked by vertical lines in the volume-averaged balances presented in the next section. We use representative cases from each of these regimes at
$\widetilde {Ra}=90$
(WN),
$\widetilde {Ra}=1200$
(RR) and
$\widetilde {Ra}=13000$
(WR) to establish the required boundary layer exclusion in § 3.1 and to demonstrate the scale dependence of the balances in § 3.4. The force balances at
$E=10^{-4}$
and
$E=10^{-6}$
are presented in supplementary figures S1 and S2, respectively, for comparison.
3.1. Boundary layer exclusion
Variation of volume-averaged r.m.s. fluctuating forces (left-hand column) and their curls (right-hand column) at
$E=10^{-5}$
with thermal forcing
$\widetilde {Ra}_T$
where only single boundary layer thicknesses are excluded (
$r_{ex}=1$
) in (a,b), and with boundary layer exclusion factor (
$r_{ex}$
) in (c–h) for the annotated cases in (a,b). The representative cases from WN, RR and WR regimes correspond to (c,d)
$\widetilde {Ra}=150$
, (e,f)
$\widetilde {Ra}=1200$
and (g,h)
$\widetilde {Ra}=13\,000$
.
We begin by assessing the variation of the fluctuating forces (figures 1
c,e,g) and their curls (figures 1
d,f,h) as a function of
$r_{ex}$
, which represents the number of VBLs excluded near each boundary from the volume average. Mean quantities exhibit similar behaviour, so we focus on fluctuating quantities in this subsection. A common practice (e.g. Yadav et al. Reference Yadav, Gastine, Christensen, Wolk and Poppenhaeger2016; Aubert et al. Reference Aubert, Gastine and Fournier2017) is to exclude a single VBL before calculating the volume average (figures 1
a,b), which corresponds to
$r_{ex}=1$
. We demonstrate the impact of changing
$r_{ex}$
in our representative cases from the WN (figures 1
c,d), RR (figures 1
e,f) and WR (figures 1
g,h) regimes. Since viscous forces are strongest near the no-slip boundaries in our simulations, that force term changes most significantly as boundary layers are excluded (figures 1
c,e,g). The viscous force term decreases with increasing
$r_{ex}$
and converges to its bulk values for
$r_{ex}\geqslant 5$
for all cases. The curled viscous force requires more exclusion than the uncurled force to converge to its bulk values (e.g. compare figures 1
c,d), probably due to the presence of an extra spatial gradient that inflates sharp changes near the boundaries.
Crucially, excluding only a single VBL to calculate the forces is generally insufficient to properly capture bulk dynamics in our simulations. In the WN regime, the viscous force is larger than inertia for
$r_{ex}=1$
(figure 1
c), while it falls below inertia for
$r_{ex}\gt 2$
. Overestimation of the viscous term in the bulk with
$r_{ex}=1$
is particularly problematic when considering curled forces, where it can result in a more than an order of magnitude overestimate (e.g. figure 1
d). In the curled forces, a VC balance is obtained for the WN cases using
$r_{ex}=1$
, whereas an AC (or TW) balance can be observed for
$r_{ex}=10$
(figure 1
d). In the RR regime (figure 1
f), the curled force balance changes from VC to IVAC with increasing exclusion. For the WR case (figure 1
h), the IVC balance obtained for
$r_{ex}=1$
becomes an IVAC balance with increasing
$r_{ex}$
. Figure 1 also shows that the ordering of terms in the force and curled balances can differ significantly at low
$r_{ex}$
, but show better agreement for large
$r_{ex}$
.
Variation of fluctuating viscous to Coriolis (a) force and (b) curled force ratios, with the thickness of the excluded layers as a multiple (
$r_{ex}$
) of VBL thickness. The representative cases from WN (circles), RR (diamonds) and WR (stars) regimes for three Ekman numbers correspond to
$\widetilde {Ra}=30,900,13\,000$
for
$E=10^{-4}$
(blue),
$\widetilde {Ra}=90,1200,13\,000$
(red) for
$E=10^{-5}$
and
$\widetilde {Ra}=150,2000,18\,000$
for
$E=10^{-6}$
(green).
To test the appropriate value of
$r_{ex}$
across the entire suite of runs, we plot the ratio of volume-averaged fluctuating viscous to Coriolis terms as a function of
$r_{ex}$
in figure 2. This ratio is generally largest for low
$r_{ex}$
, and decreases with increasing
$r_{ex}$
. The value of
$r_{ex}$
required to obtain converged results (i.e. that do not change upon a further increase in
$r_{ex}$
) varies somewhat with the parameters, but is always larger for the curled terms than the forces. Across our suite of simulations, a converged bulk force balance is always obtained with
$r_{ex}\geqslant 5$
. The curled force ratio (figure 2
b) exhibits slower convergence to the bulk values than the uncurled forces, and can require
$r_{ex}\geqslant 10$
. Following the analysis in figure 2, we exclude ten VBLs (i.e.
$r_{ex}=10$
) in all volume averages. This corresponds to
$2{-}4$
VBLs if the local maxima method is used instead of a linear intersection method, as the former leads to thicker estimates of the VBL (Gastine et al. Reference Gastine, Wicht and Aurnou2015). For the thickest boundary layers at
$\widetilde {Ra}=30$
and
$E=10^{-4}$
, approximately
$34{-}38\,\%$
of the gap is excluded from the radial extent of the domain at
$r_{ex}=10$
. For the majority of simulations, the total excluded region amounts to less than 20 % of the gap width.
3.2. Force balance
Volume-averaged r.m.s. mean (left-hand column) and fluctuating (right-hand column) force components in (a,b)
$\hat {r}$
, (c,d)
$\hat {\theta }$
and (e,f)
$\hat {\phi }$
for
$E=10^{-5}$
.
The volume-averaged mean and fluctuating forces as functions of thermal forcing
$\widetilde {Ra}_{T}$
are shown in figures 3(a,c,e) and 3(b,d,f), respectively. Because buoyancy acts only radially, and the Coriolis force is horizontal, there will be a directional anisotropy in the force balances, therefore we consider the radial (
$\hat {r}$
, figures 3
a,b), co-latitudinal (
$\hat {\theta }$
, figures 3
c,d), and azimuthal (
$\hat {\phi }$
, figures 3
e,f) force components separately.
In the radial direction (figure 3
a), we find a primary TW balance between the mean Coriolis, pressure and buoyancy forces, with the buoyancy force gradually becoming subdominant at the highest values of
$\widetilde {Ra}_{T}$
considered. In the
$\theta$
direction (figure 3
c), the primary balance in the mean forces is QG (i.e. between Coriolis and pressure forces) for all
$\widetilde {Ra}_{T}$
. These mean balances in
$\hat {r}$
and
$\hat {\theta }$
are the same as the balances reported in recent dynamo simulations (Calkins et al. Reference Calkins, Orvedahl and Featherstone2021). The residuals of the primary TW balance in
$\hat {r}$
and the primary QG balance in
$\hat {\theta }$
are balanced by the Reynolds stress (
$\overline {I^{ff}}$
), which dominates the two components of the total mean advection term (
$\overline {I}=\overline {I^{mm}}+\overline {I^{ff}}$
) until the largest forcings considered. In the mean
$\hat {r}$
and
$\hat {\theta }$
balances, viscosity is always subdominant, as found by Calkins et al. (Reference Calkins, Orvedahl and Featherstone2021). The primary mean forces in the azimuthal direction (figure 3
e) have similar magnitude with the secondary forces in the
$\hat {r}$
and
$\hat {\theta }$
directions. Owing to our choice of averaging, there is no mean pressure gradient or buoyancy force in the azimuthal direction. Therefore, the Coriolis force is balanced by inertia and viscous forces at low
$\widetilde {Ra}_{T}$
(IVC balance), whereas an IC balance dominates at high
$\widetilde {Ra}_{T}$
.
In the fluctuating forces in figures 3(b,d,f), the QG balance between Coriolis and pressure forces dominates in all directions. At low
$\widetilde {Ra}_T$
in the radial direction, there is a secondary
$\text{AC}_{ag}$
balance between ageostrophic Coriolis force and buoyancy, with viscosity approximately an order of magnitude weaker. At high
$\widetilde {Ra}_T$
in the radial direction, there is a secondary
$\text{IAC}_{ag}$
balance, with viscosity still subdominant but less than an order of magnitude weaker. In the
$\theta$
and
$\phi$
directions, the ageostrophic Coriolis force is balanced by both viscous and inertial forces for low
$\widetilde {Ra}_T$
(
$\text{IVC}_{ag}$
balance), with the relative importance of viscosity weakening as
$\widetilde {Ra}_T$
increases.
The behaviour of the fluctuating radial forces can be compared with the free-slip and fixed temperature simulations of Nicoski et al. (Reference Nicoski, O'Connor and Calkins2024) (see figure 16(b) in their paper). For both studies, the primary radial balance is QG, while the secondary balance is
$\text{AC}_{ag}$
at low thermal forcings, and
$\text{IAC}_{ag}$
with a weakly subdominant viscous contribution at higher thermal forcings. Though the viscous force can be higher than inertia for the lowest
$\widetilde {Ra}_T$
values considered here (figure 3
b), it is the smallest force in both studies when the same range of
$\widetilde {Ra}_T$
is considered. However, the two studies find different dominant contributions to the inertia force. Figure 3(b) shows that
$I^{\prime}$
is dominated by
$I^{ff}$
, which balances
$\text{C}_{ag}$
and
$\text{A}^{\prime}_r$
for large
$\widetilde {Ra}_T$
. In contrast, Nicoski et al. (Reference Nicoski, O'Connor and Calkins2024) found that the mean-fluctuating inertia (
$I^{mf}$
) dominates the nonlinear advection terms, which, combined with the TD, balances
$\text{C}_{ag}$
and
$\text{A}^{\prime}_r$
for large
$\widetilde {Ra}_T$
. This difference probably arises from the strong mean zonal flows in the simulations of Nicoski et al. (Reference Nicoski, O'Connor and Calkins2024), owing to the use of free-slip conditions.
The partitioning into mean (
$m=0$
) and fluctuating (
$m\neq 0$
) forces brings out distinct balances (comparing the left-hand and right-hand plots of figure 3). For example, the radial force balance exhibits a primary mean TW balance and a primary fluctuating QG balance. The unpartitioned (i.e. mean + fluctuating) radial force balance has a primary QG balance that exceeds a subdominant
$\text{AC}_{ag}$
balance by a factor
$2{-}3$
, essentially averaging the large- and small-scale dynamics.
The azimuthal component of azimuthally averaged forces has been considered in previous studies (Sheyko et al. Reference Sheyko, Finlay, Favre and Jackson2018; Menu et al. Reference Menu, Petitdemange and Galtier2020) since it removes the pressure gradient, although there is no buoyancy force in this direction. Our analysis indicates that this representation (see figure 3(e) for an example) does not reflect the balance of mean forces in the
$\hat {r}$
and
$\hat {\theta }$
directions, and also does not correspond to the balances of the fluctuating part of the forces.
3.3. Curled force balance
Volume-averaged r.m.s. mean (left-hand column) and fluctuating (right-hand column) curled force components in (a,b)
$\hat {r}$
, (c,d)
$\hat {\theta }$
and (e,f)
$\hat {\phi }$
for
$E=10^{-5}$
.
Though the force balance provides useful insights into the dynamics, all forces in (2.2) are non-solenoidal and therefore will have gradient portions. These gradient portions of the forces are balanced by the pressure gradient term, which plays no role in the dynamics (Hughes & Cattaneo Reference Hughes and Cattaneo2019; Teed & Dormy Reference Teed and Dormy2023). Our approach to removing these gradients is to take the curl of the momentum equation. We partition the curled forces into mean (figures 4
a,c,e) and fluctuating (figures 4
b,d,f) parts, and separately consider the radial (figures 4
a,b), co-latitudinal (figures 4
c,d), and azimuthal (figures 4
e,f) components. Since the buoyancy force is radial, its curl acts only in the angular directions (
$\theta$
and
$\phi$
), making the curled balance inherently anisotropic.
In the
$\hat {r}$
and
$\hat {\theta }$
directions, we find a primary mean balance between the Coriolis, inertial and viscous terms (IVC balance), with the viscous contribution weakening at the highest values of
$\widetilde {Ra}_T$
considered. In the
$\hat {\phi }$
direction, the mean balance is a TW, except at the highest values of
$\widetilde {Ra}_T$
considered, where the mean curled inertia also enters this primary balance. A TW arises only in the mean azimuthal component because our choice of averaging causes the
$\theta$
-component of the mean curled buoyancy force to vanish (since
$\partial \overline {T}/\partial \phi = 0$
). Aubert (Reference Aubert2005) also found a TW balance in the azimuthally averaged curled force balance for non-magnetic simulations at
$E = 10^{-4}{-}10^{-5}$
.
The curled fluctuating forces in
$\hat {\theta }$
and
$\hat {\phi }$
exhibit a primary balance between Coriolis, buoyancy and viscous terms (VAC balance) at low
$\widetilde {Ra}_T$
, while the inertial force gradually enters this balance with increasing
$\widetilde {Ra}_T$
. A similar trend is observed in the radial balance, except for the omission of the buoyancy term. We note here that the radial curled force balance, as reported by Dormy (Reference Dormy2016), may not represent the curled force balance in the other directions.
As with the forces, partitioning curled terms into mean and fluctuating components brings out distinct balances that would be obscured if only the unpartitioned curled terms were considered. In particular, because the fluctuating curled quantities have much higher amplitude than the mean curled quantities, an unpartitioned (i.e. mean + fluctuation) curled force representation would simply show the fluctuating small-scale balance. The balance obtained from the total magnitude of the three mean curled components follows the balance in the
$\hat {\phi }$
direction shown in figure 4(e), while the balance obtained from the total magnitude of the three fluctuating curled components reflects the balances in the
$\hat {\theta }$
and
$\hat {\phi }$
directions shown in figures 4(d) and 4(f).
3.4. Scale-dependent force balance
Volume-averaged r.m.s. of the total magnitude of fluctuating forces (left-hand column) and their curls (right-hand column) at
$E=10^{-5}$
as a function of thermal forcing
$\widetilde {Ra}_T$
in (a,b) and with SH degree (
$l$
) in (c–h) for the exemplar cases highlighted in (a,b). The representative cases from WN, RR and WR regimes correspond to (c,d)
$\widetilde {Ra}=90$
, (e,f)
$\widetilde {Ra}=1200$
and (g,h)
$\widetilde {Ra}=13\,000$
.
In figure 5, we compare the scale dependence of the total magnitude of the fluctuating forces (left-hand column) and their curls (right-hand column) with the corresponding scale-integrated representation (figures 5
a,b). For brevity, we again select one case from each of the WN (figures 5
c,d), RR (figures 5
e,f), and WR (figures 5
g,h) regimes, corresponding to
$\widetilde {Ra}=90, 1200, 13\,000$
, respectively (these cases are highlighted in figures 5
a,b). The regime boundaries are shown using dashed and solid vertical lines, following the analysis of L20, as described in detail in the next subsection. Although the mean forces are also dependent on the spherical harmonic (SH) degree, we focus on the fluctuating forces because of their more direct correspondence with the convective motions and the heat transfer behaviour.
Based on the scale-integrated fluctuating forces, all simulations have a QG primary balance (figure 5
a). This primary force balance holds across all scales in the WN and RR regimes (figures 5
c,e), while inertia enters the primary balance at the small scales (i.e. large
$l$
) in the WR regime (figure 5
g). The secondary force balance in the WN and RR regimes is characterized by an
$\text{AC}_{ag}$
balance between buoyancy and ageostrophic Coriolis at large scales, with inertia entering at small scales (figure 5
e). In the WR regime, there is a secondary
$\text{IAC}_{ag}$
balance with a significant inertial contribution at all scales. The scale dependence of forces in the RR regime (figure 5
e) is similar to the force spectra reported in Schwaiger et al. (Reference Schwaiger, Gastine and Aubert2020) (see figure 7(a) in their paper) and is often referred to as a QG–
$\text{IAC}_{ag}$
balance.
Previous studies have attempted to relate crossings of scale-dependent forces (such as the crossing of buoyancy and inertia forces at
$l\sim 22$
in figure 5
e) to dynamically relevant length scales (e.g. Schwaiger et al. Reference Schwaiger, Gastine and Aubert2020). However, for high
$\widetilde {Ra}_T$
in the WR regime (figure 5
g), the observed crossings occur at
$l\sim 3$
, which is much smaller than the dominant wavenumber of the flow (
$l\sim 15$
, based on the peak of the kinetic energy spectra). Therefore, it may not always be possible to relate crossings of scale-dependent forces to dynamically significant length scales.
Compared to the scale-dependent forces, the scale-dependent curled balances do not produce a clear separation of balances (i.e. primary/secondary/tertiary). In the WN regime (figure 5 d), an AC balance is evident at large scales, while viscous and inertia forces enter the balance at small scales. The curled force magnitudes in the RR regime (figure 5 f) are comparable at all scales, leading to an IVAC balance in the scale-integrated representation (figure 5 b). In the WR regime (figure 5 h), an IC balance is observed at large scales, while the inertia dominates the balance at small scales.
Figure 5 compares scale-integrated and scale-dependent representations of dynamical balances in our simulations. Integrated forces can not capture cross-over scales where they exist, or the general decrease of the buoyancy force and increase of the inertial force with increasing
$l$
(e.g. figure 5
e). Furthermore, integrated quantities generally do not reflect the force balance at the smallest scales of the solution. Nevertheless, the integrated representation quantitatively captures the overall ordering of forces in a single measure that can be easily compared across a large suite of simulations. In comparison, the curled balances are comparatively less scale-dependent, exhibiting neither a clear ordering of the forces nor any distinct cross-over scales (see also Teed & Dormy Reference Teed and Dormy2023). Therefore, a scale-integrated analysis (figure 5
b) is sufficient to represent the curled forces.
3.5. Ekman number dependence
Volume-averaged r.m.s. fluctuating forces (left-hand column) and their curls (right-hand column) at (a,b)
$E=10^{-4}$
and (c,d)
$E=10^{-6}$
. The dashed vertical line represents the thermal forcing where a transition from WN to WR regimes happens according to the scaling predictions of Long et al. (Reference Long, Mound, Davies and Tobias2020). Their analysis does not predict the existence of the RR regime for
$E\geqslant 10^{-4}$
.
We can now consider the general trend in the balance of forces with increasing
$\widetilde {Ra}_T$
and decreasing
$E$
. The total magnitude of the fluctuating forces and their curls for
$E=10^{-4}$
and
$E=10^{-6}$
increases with
$\widetilde {Ra}_T$
(figure 6), exhibiting a similar dependence to their counterparts at
$E=10^{-5}$
(figures 3 and 4). The Coriolis and pressure forces in the primary QG balance become increasingly separated from the secondary forces as
$E$
is reduced (compare figures 6
a,c). Although the separation between viscous and Coriolis forces increases with decreasing
$E$
, the viscous force magnitude remains similar to the other secondary forces. According to asymptotic predictions (Nicoski et al. Reference Nicoski, O'Connor and Calkins2024), although the viscous force should decrease compared to Coriolis with decreasing
$E$
, the ratio of the viscous force to buoyancy and inertia should remain invariant, which is consistent with our results. The viscous force in the curled balance (figures 6
b,d) also has a magnitude comparable to other curled forces over the investigated range of
$\widetilde {Ra}_T$
and
$E$
. In summary, we see broadly similar primary and secondary balances in the forces and curled forces at different
$E$
, including transitions between dynamic regimes as
$\widetilde {Ra}_T$
varies.
4. Summary of dynamical balances and comparison to regimes of L20
Table 1 presents a qualitative summary of the different balances that we found within our suite of simulations depending on whether we considered forces or curled forces, total magnitudes or individual vector components, or partitioning into azimuthally mean and fluctuating contributions. We reiterate here that our analysis reflects bulk dynamics, with the volume averages obtained after removing ten VBLs from each boundary of the domain. Balances similar to those described in detail above for
$E=10^{-5}$
are also found at
$E=10^{-4}$
and
$E=10^{-6}$
; the main difference is that the separation between the primary and secondary balances, denoted by a dash (
$-$
) in the table, increases with decreasing Ekman number (e.g. compare figure 5 in this text with figure 3 in supplementary materials S1 and S2). We have tried to denote ‘balances’ that are groupings of two or more terms that are separated by an order of magnitude in amplitude from other terms; however, such a large separation is not always present. Changes in the balances with
$\widetilde {Ra}_T$
are indicated by a right arrow (
$\rightarrow$
), which indicates the general trends with increased thermal forcing but not the specific values at which the balances change. Therefore, table 1 is only a general description of the complex variations in force balances amongst our suite of simulations. In general, increased thermal forcing results in an increase in the relative importance of inertial terms in the balances.
Summary of force and curled force balances in our simulations. The ACP balance of forces (or the AC balance of curled forces) is referred to as a TW balance, while the residual of these forces is designated here as
$(\text{TW})_{res}$
. Similarly, the primary balance between Coriolis and pressure gradient forces is denoted as a QG balance. Primary and secondary force balances are separated by a dash (
$-$
), while the changes in the balance with increasing thermal forcing (
$\widetilde {Ra}_T$
) are designated with a right arrow (
$\rightarrow$
).
We now compare transitions in the force and curled force balances to previous predictions of regime transitions based on scaling laws (L20). We put emphasis on the
$E\leqslant 10^{-5}$
cases as they are more appropriate for comparing with the asymptotic scaling theories used by L20. They defined the WN–RR regime transition (dashed vertical lines in figures 3–5) as
$Ra_{T} = 8Ra_T^{c}$
based on the observed gradual departure from the linear
$Nu-1\propto Ra_{T}/Ra_T^{c} -1$
scaling expected just above onset. In the WN regime, they found that the simulated flow length scale
$\ell$
and convective Reynolds number
$Re_c$
follow the predictions of VAC theory. L20 defined the RR–WR transition (solid vertical lines in figures 3–5) based on the condition
$Ra_{T}E^{8/5}\sim O(1)$
of Julien et al. (Reference Julien, Knobloch, Rubio and Vasil2012a
) above which the thermal boundary layers lose geostrophic balance. They found scalings for
$\ell$
and
$Re_c$
close to but statistically different from the predictions of IAC theory.
Figures 3 and 4 show that mean forces exhibit no changes in primary or secondary balances over the range of
$\widetilde {Ra}_T$
considered, hence do not conform to the regime transitions found by L20. This is expected, since L20 defined transitions based on quantities that depend strongly on convective fluctuations such as
$Nu$
,
$\ell$
and
$Re_c$
. In the fluctuating forces, the WN–RR transition correlates with viscous and inertial terms coming into approximate balance with the ageostrophic and buoyancy terms that comprise the secondary balance (figure 5
a), while in the fluctuating curled forces, this transition arises when the magnitude of the nonlinear advection term becomes comparable to the Coriolis, buoyancy and viscous terms (figure 5
b). In the fluctuating forces, the RR–WR transition broadly correlates with the amplitude of the viscous term falling below that of the secondary balance, while in the fluctuating curls, this transition appears to correlate with the amplitude of the inertial term rising above the remaining terms. However, the total fluctuating forces and curled forces do not suggest an exact value of
$\widetilde {Ra}_T$
where regime transitions occur. Indeed, these quantities exhibit gradual changes with
$\widetilde {Ra}_T$
and
$E$
(e.g. figures 5(a,b) and figures 3(a,b) in both S1 and S2), hence any transition inferred from them is necessarily broad rather than abrupt.
Table 2 summarizes the general character of the balances in the regimes defined by L20, calculated using the total magnitude of the mean (supplementary material S3) and fluctuating forces and curled forces. L20 inferred a VAC balance in the WN regime, and an IAC balance in the RR regime, using scaling theory based on the curled force balance. In the WN regime, the calculated fluctuating curled force balance is VAC at low
$\widetilde {Ra}_T$
, transitioning to an IVAC balance as
$\widetilde {Ra}_T$
increases (figure 5
b). This behaviour is broadly consistent with the assumptions of L20. In the section of the RR regime accessed by our simulations, the calculated fluctuating curled force balance is IVAC rather than the IAC balance assumed by L20. The viscous force is also significant in the force balance (figure 5
a), though it remains smaller than the other forces in the RR regime. Similar behaviour of the viscous force can also be observed at
$E=10^{-4}$
and
$E=10^{-6}$
(see figures 3(a,b) in both S1 and S2). Reconciling the calculated dynamical balances in the RR regime with other flow diagnostics such as
$\ell$
,
$Nu$
and
$Re$
must await a scaling theory for the IVAC regime.
Summary of force and curled force balances in the regimes of RC simulations as predicted by L20. The balances in total force magnitudes (2.9) have been used here. The abbreviations used here are the same as described in table 1. Primary and secondary force balances are separated by a dash (
$-$
), while the changes in the balance with increasing thermal forcing (
$\widetilde {Ra}_T$
) are designated with a right arrow (
$\rightarrow$
).
In summary, the fluctuating force and curled force balances exhibit smooth variations over the range of
$\widetilde {Ra}_T$
and
$E$
considered, reflecting gradual rather than abrupt changes in the dynamics. Broadly speaking, it appears that the RR regime as defined by L20 corresponds to a range of
$\widetilde {Ra}_T$
where inertial effects enter the primary fluctuating curled force balance (or the secondary fluctuating force balance). This observation motivated us to seek a single parameter to characterize the changing dynamics with
$\widetilde {Ra}_T$
. However, it is difficult to find a single quantity that adequately represents the transitions in heat transport and flow behaviour identified by L20. This is perhaps unsurprising given that even the simple VAC scaling laws used by L20 are defined by at least two parameters, while the dominant dynamical balances identified in table 2 involve at least three terms. We therefore classify these balances by introducing two parameters based on our calculated dynamical balances.
Regime diagram with predicted regime boundaries from the analysis of L20. The primary force balance is classified as QG when
$\mathcal {F}_{I/C}\leqslant 0.1$
(filled symbols), and non-QG when
$\mathcal {F}_{I/C} \gt 0.1$
(open symbols). The markers are coloured by
$\log_{10}(\mathcal {CF}_{I/C})$
. L20 found that the WN–RR regime transition corresponded to the condition
$Ra_{T}/Ra_{T}^{c}=8$
, and the RR–WR transition corresponded to the condition
$Ra_{T}=0.6E^{-8/5}$
of Julien et al. (Reference Julien, Rubio, Grooms and Knobloch2012b
).
Figure 7 shows a quantitative comparison of our dynamical balances with the regime diagram of L20 (their figure 14). To classify the balance in the
$E{-}Ra_{T}/Ra_{T}^{c}$
parameter space, we are motivated by the observation that inertia varies most strongly with
$\widetilde {Ra}_T$
in our simulations (figure 5). We therefore introduce two new measures based on the total magnitude of the fluctuating forces and the curled forces (figure 5
a,b) that assess the role of inertia in the force balance.
To measure the degree of geostrophy in the balance we define a force ratio
\begin{equation} \mathcal {F}_{I/C} = \frac {I^{\prime}}{C^{\prime}}. \end{equation}
We find that the value
$\mathcal {F}_{I/C}=0.1$
can be used to demarcate the simulations belonging to the WR regime as demonstrated in figure 7. Here,
$\mathcal {F}_{I/C}\gt 0.1$
(open symbols) can describe almost all simulations that fall in the WR regime, as compared to the simulations in the range
$\mathcal {F}_{I/C}\leqslant 0.1$
(filled symbols) that mostly fall inside the WN and RR regimes. A comparison can be made between the force ratio
$\mathcal {F}_{I/C}$
and the Rossby number
$Ro=2\,Re\,E$
(where
$Re$
is the Reynolds number defined in Appendix A), which is often used to represent the ratio of the two forces. We find that the agreement between the magnitudes of
$\mathcal {F}_{I/C}$
and
$Ro$
worsens as
$E$
decreases, and at our lowest
$E$
,
$Ro$
underestimates
$\mathcal {F}_{I/C}$
by approximately a factor of ten.
To measure the role of inertia in the curled balance, we further define a curled force ratio
\begin{equation} \mathcal {CF}_{I/C} = \frac {\nabla \times I^{\prime}}{\nabla \times C^{\prime}}. \end{equation}
The symbols in figure 7 are coloured by
$\log_{10} (\mathcal {CF}_{I/C} )$
. We emphasize again that the roles of all terms, including inertia, in the dynamical balances change gradually, hence values of
$\mathcal {F}_{I/C}$
and
$\mathcal {CF}_{I/C}$
that reflect regime transitions must be chosen arbitrarily. Nevertheless, the RR regime can be characterized by the combination of
$\mathcal {F}_{I/C}\lesssim 0.1$
and
$\log_{10} (\mathcal {CF}_{I/C} )\approx 0$
(i.e.
$\mathcal {CF}_{I/C}\approx 1$
), which reflects the primary QG balance in the fluctuating forces and the IVAC balance in the fluctuating curled forces, respectively (table 2).
5. Discussion and conclusions
We have analysed different representations of dynamical balances in simulations of spherical shell RC. The radial, co-latitudinal and azimuthal components of the forces have been considered separately to demonstrate the anisotropic nature of dynamical balances. We also partition the forces into azimuthally averaged mean and corresponding fluctuating parts that exhibit distinct balances. Furthermore, the curled force components are also analysed to investigate the solenoidal force balance. The utility of a scale-dependent representation of curled and uncurled forces has been addressed. Our main findings are presented in figures 3–5 and table 1, and can be summarized as follows.
-
(i) The bulk curled force balance depends critically on the number of VBLs that are removed near the upper and lower boundaries. We find that removing ten VBLs from each boundary provides a robust estimate of the curled force balance that is broadly consistent with the balance obtained from calculating forces.
-
(ii) Mean and fluctuating forces and curled forces exhibit distinct balances, consistent with the results of Calkins et al. (Reference Calkins, Orvedahl and Featherstone2021) and Nicoski et al. (Reference Nicoski, O'Connor and Calkins2024). In particular, the primary mean force and curled force balances are thermal wind (TW), while the primary fluctuating force and curled force balances are quasi-geostrophic (QG) and Inertia–Viscous–Archimedian–Coriolis (IVAC), respectively.
-
(iii) Radial, co-latitudinal and azimuthal forces exhibit distinct balances as found by Calkins et al. (Reference Calkins, Orvedahl and Featherstone2021) and Aubert (Reference Aubert2005) for dynamo simulations. For example, mean forces exhibit a primary TW balance in the radial direction, QG balance in the latitudinal direction, and Inertial–Viscous–Coriolis (IVC) balance in the longitudinal direction. A total force magnitude representation underestimates the role of buoyancy compared to the radial balance.
-
(iv) In the scale-dependent balances, the separation of magnitude between the forces decreases when a curl operation is performed. Cross-over scales are observed in some but not all force balances, and are not observed in curled force balances, consistent with the results of Teed & Dormy (Reference Teed and Dormy2023). The curled forces are only weakly scale-dependent and therefore suitably represented by scale-integrated quantities.
-
(v) Transitions in fluctuating force and curled force balances are broadly consistent with the three regimes of RC obtained by Long et al. (Reference Long, Mound, Davies and Tobias2020). However, the relative importance of forces (and their curls) varies gradually with thermal forcing rather than exhibiting any abrupt changes, and therefore does not define precise values of transition parameters.
-
(vi) The rapidly rotating (RR) regime broadly corresponds to a range of thermal forcing where inertia is of comparable magnitude to the other terms in the primary curled force balance. We find an IVAC rather than an IAC balance in the fluctuating curled forces in the RR regime for the investigated parameter regime (see table 1). Also, the viscous force increases with increasing thermal forcing, and remains significant in the secondary force balance even if
$E$
is lowered. These observations are consistent with the predictions of Nicoski et al. (Reference Nicoski, O'Connor and Calkins2024).
The dynamical balances in table 1 can be compared to results obtained in previous studies. The mean force balance in the dynamo simulations of Calkins et al. (Reference Calkins, Orvedahl and Featherstone2021) with no-slip boundary conditions is TW in the
$\hat {r}$
direction and QG in the
$\hat {\theta }$
direction, consistent with our results. The mean curled forces in a dynamo simulation with no-slip conditions (Aubert Reference Aubert2005) is also in a TW balance in the azimuthal direction as in our non-magnetic simulations (figure 4
a). This indicates that the primary mean balance in
$\hat {r}$
and
$\hat {\theta }$
is consistent between non-magnetic and dynamo simulations, hence is relatively unaffected by the presence of a magnetic field.
The radial fluctuating force behaviour in our no-slip simulations is consistent with the results of Nicoski et al. (Reference Nicoski, O'Connor and Calkins2024), who observed a primary QG balance in simulations with stress-free boundary conditions. This indicates that the fluctuating balance in non-magnetic RC is not sensitive to the velocity boundary conditions. Notably, the fluctuating viscous force remains non-negligible within our suite of simulations, similar to the findings of Nicoski et al. (Reference Nicoski, O'Connor and Calkins2024). Indeed, in the fluctuating curled forces (figure 5
b), the viscous term is always part of the dominant balance, even at
$E=10^{-6}$
(see figure 3(b) in supplementary material S2).
Our scale-dependent force balance in the RR regime is consistent with the balance reported in a previous non-magnetic simulation (Schwaiger et al. Reference Schwaiger, Gastine and Aubert2020). Although the scale-dependent force balance does not always exhibit a clear cross-over between forces in all regimes in RC, we find a cross-over between buoyancy and inertia forces in the secondary balance in the RR regime (figures 5(e) and 3(e) in S1 and S2). However, as observed previously by Teed & Dormy (Reference Teed and Dormy2023), such cross-overs are not found in the scale-dependent curled balance (figures 5(f) and 3(f) in S1 and S2). This can be attributed to the separations among the terms, which reduce owing to the curl operation that removes the dynamically irrelevant gradient part. Also, the viscous force, which is important only in the small-scale force balance (figure 5
e), is significant at all scales in the curled balance (figure 5
f). This behaviour of the viscous force remains the same for a lower Ekman number (
$E=10^{-6}$
; see figure 3(f) in S2). Whether the asymptotic separation among various forces, as demonstrated in the numerical simulations of Nicoski et al. (Reference Nicoski, O'Connor and Calkins2024), can also be demonstrated for curled forces, requires future simulations at lower Ekman numbers (
$E \lt 10^{-6}$
).
We choose to look into the force components in a spherical coordinate system (
$\hat {r}$
,
$\hat {\theta }$
,
$\hat {\phi }$
), which matches the spherical shell geometry that we considered and takes advantage of the spherical harmonic representation of our numerical implementation. In rotationally constrained convection (i.e. WN and RR regimes), it can also be useful to study the forces in a cylindrical coordinate system (
$\hat {s}$
,
$\hat {\phi }$
,
$\hat {z}$
), where
$\hat {s}=\sin {\theta }\,\hat {r}+\cos {\theta }\,\hat {\theta }$
is the cylindrical radius, and
$\hat {z}=\cos {\theta }\,\hat {r}-\sin {\theta }\,\hat {\theta }$
is the vertical axis of rotation. In particular, strong rotation causes the small-scale leading-order motion to be independent of the
$\hat {z}$
direction, known as the Taylor–Proudman constraint (Proudman Reference Proudman1916; Taylor Reference Taylor1923). Asymptotically reduced models enforcing this constraint predict the horizontal (i.e.
$\hat {\phi }$
) viscous force to be large compared to the vertical component (
$\hat {z}$
). In the supplementary material, we plot the
$\hat {z}$
-component of mean and fluctuating forces (S4, figure 1) and the curled forces (S4, figure 2) for the three Ekman numbers
$E=10^{-4}$
(a,b),
$E=10^{-5}$
(c,d) and
$E=10^{-6}$
(e,f). The mean vertical (i.e.
$\hat {z}$
) viscous force (S4, figure 1
c) in the WN regime can be seen to be approximately one order of magnitude smaller than the azimuthal (
$\hat {\phi }$
) viscous force (figure 3
e), which corresponds well with asymptotic theory (e.g. Sprague et al. Reference Sprague, Julian, Knobloch and Werne2006). However, for higher
$\widetilde {Ra}_T$
in the RR and WR regimes, we find that the viscous forces have comparable magnitude in all directions. Also, the fluctuating viscous forces do not exhibit any separation between vertical and horizontal directions. It is difficult to draw any direct comparison with the multiscale asymptotic theory of Sprague et al. (Reference Sprague, Julian, Knobloch and Werne2006) since the definitions for mean and fluctuating quantities in the theory are different from our partitioning of forces as defined in (2.6), and the Ekman number range in our study may not be low enough to ensure clear separation at all scales.
Another issue related to the incompressibility condition is the vanishing of azimuthal Coriolis force (
$C_{\phi }$
) when averaged over geostrophic cylinders (i.e. over
$\phi$
and
$z$
). This leads to the zonal flow balance between Reynolds stress (
${\overline {{I}_{ff}}}^{z}$
) and viscous forces (
${\overline {V}}^{z}$
). We recover this balance when integrating over geostrophic cylinders that extend to the solid boundaries. Since the incompressibility condition that enforces the vanishing of Coriolis force requires a closed surface, this balance is not recovered for integration over cylinders constrained to the bulk (i.e. excluding the boundary regions). Therefore, this balance is not recovered for the azimuthal forces in figure 3(e), which have been volume-integrated over the bulk by excluding regions near the boundary.
Our analysis does not take into account the temporal variation of dynamical balances (see e.g. Schaeffer et al. Reference Schaeffer, Jault, Nataf and Fournier2017), though previous studies suggest that these variations are small (Aubert et al. Reference Aubert, Gastine and Fournier2017). We also do not consider dynamical balances in the boundary layers since we aim to characterize the bulk dynamics that would ultimately be used to extrapolate to the conditions of planetary interiors and can be compared to available observations. Nevertheless, the calculated bulk dynamical balances form a basis for comparison with different theoretical analyses of RC (e.g. Aubert et al. Reference Aubert, Brito, Nataf, Cardin and Masson2001; Calkins et al. Reference Calkins, Orvedahl and Featherstone2021; Nicoski et al. Reference Nicoski, O'Connor and Calkins2024). Similar force calculations can be useful to study various dynamical regimes of convection with
$Pr\neq 1$
(Calkins et al. Reference Calkins, Aurnou, Eldredge and Julien2012; Guzmán et al. Reference Guzmán, Madonia, Cheng, Ostilla-Mónico, Clercx and Kunnen2021), double-diffusive convection (Tassin et al. Reference Tassin, Gastine and Fournier2021) or geodynamo simulations (Calkins et al. Reference Calkins, Orvedahl and Featherstone2021; Mound & Davies Reference Mound and Davies2023).
Supplementary materials
Supplementary materials are available at https://doi.org/10.1017/jfm.2025.259.
Acknowledgements
This work used the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk) and ARC2, part of the high-performance computing facilities at the University of Leeds, UK. We would like to thank J.A. Nicoski and M.A. Calkins for providing the data to validate our force calculations. We gratefully acknowledge two anonymous reviewers for their suggestions that improved this work.
Funding
S.N., C.J.D. and J.E.M. are supported by Natural Environment Research Council research grant NE/W005247/1. A.T.C. and C.J.D. are supported by NE/V010867/1.
Data availability statement
The data that support the findings of this study are openly available in National Geoscience Data Centre (NGDC), British Geological Survey at https://doi.org/10.5285/6c555767-5a85-4f32-9db6-d947b1aa3387.
Declaration of interests
The authors report no conflict of interest.
Author contributions
C.J.D., J.E.M. and S.N. conceptualized and designed the research. S.N. and A.T.C. developed and conducted numerical simulations and analysis. S.N., C.J.D. and J.E.M. wrote the manuscript.
Appendix A. Table of results
A summary of the characteristics of the three new simulations performed in this study is reported in table 3. The model resolution, input parameters and selected output diagnostic quantities complement table 5 of Appendix B of Mound & Davies (Reference Mound and Davies2017), and table 12 in the appendix of L20. Here,
$N$
is the numerical resolution, which equals both the number of radial points and the maximum spherical harmonic degree and order, and
$N_{\delta i}$
and
$N_{\delta o}$
are the number of radial points within the VBLs at the inner and outer boundaries, respectively, where the VBL thicknesses are estimated from the linear intersection method as described in § 2.3. Definitions of the Ekman and modified Rayleigh numbers are given in § 2.1. The Reynolds number is defined as
$Re=U^{*}h/\nu =U=\sqrt {2\,KE/V_s}$
, where
$U$
is the non-dimensional velocity (the asterisk indicates a dimensional quantity),
$V_s$
is the shell volume, and
$KE=\iiint _{V_s}\boldsymbol {u}\boldsymbol {\cdot}\boldsymbol {u}\, {\rm d}V$
is the kinetic energy integral. Here,
$Re_{pol}$
is found by retaining only the poloidal velocity in the kinetic energy integral;
$Re_{zon}$
is found by retaining only the terms with spherical harmonic order
$m = 0$
, from the spherical harmonic expansion of the toroidal velocity in the kinetic energy integral;
$\mathcal {P}$
is the time average of the buoyancy production throughout the shell; and
$\epsilon _U$
is the time average of the viscous dissipation throughout the shell.
Summary of the three new runs at
$E=10^{-6}$
.
Appendix B. Table of critical Rayleigh numbers
The critical values of Rayleigh numbers and wavenumber at the onset of thermal convection at various Ekman numbers are listed in table 4. The critical values of Rayleigh numbers are provided for all three definitions used in this paper as defined in § 2.2.
Critical values of Rayleigh numbers and the critical wavenumbers at various Ekman numbers.
Open access
