Asymptotic scaling relations for rotating spherical convection with strong zonal flows

Abstract We analyse the results of direct numerical simulations of rotating convection in spherical shell geometries with stress-free boundary conditions, which develop strong zonal flows. Both the Ekman number and the Rayleigh number are varied. We find that the asymptotic theory for rapidly rotating convection can be used to predict the Ekman number dependence of each term in the governing equations, along with the convective flow speeds and the dominant length scales. Using a balance between the Reynolds stress and the viscous stress, together with the asymptotic scaling for the convective velocity, we derive an asymptotic prediction for the scaling behaviour of the zonal flow with respect to the Ekman number, which is supported by the numerical simulations. We do not find evidence of distinct asymptotic scalings for the buoyancy and viscous forces and, in agreement with previous results from asymptotic plane layer models, we find that the ratio of the viscous force to the buoyancy force increases with Rayleigh number. Thus, viscosity remains non-negligible and we do not observe a trend towards a diffusion-free scaling behaviour within the rapidly rotating regime.


Introduction
Rotating convection plays an important dynamical role in stars and planets, where it is believed to be one of the primary drivers of global scale magnetic fields (Stanley & Glatzmaier 2010;Jones 2011;Aurnou et al. 2015), and possibly gives rise to coherent large-scale flows such as zonal jets and vortices, as observed on the giant planets (Heimpel et al. 2022;Siegelman et al. 2022).Understanding the physics of turbulence driven by rotating convection remains challenging due to the vast range of spatiotemporal scales.As a result, the parameter space accessible to direct numerical simulation (DNS) and laboratory experiments remains distant from that which characterises natural convective systems.An approach that is often taken to overcome this limitation is to identify asymptotic scaling behaviour in model output so that extrapolation to the conditions of natural systems is possible (e.g.Christensen 2002;Aurnou 2007).The development of asymptotically reduced equation sets is another strategy that has been particularly helpful for improving understanding of rotating convective turbulence and dynamos in simplified planar geometries since more extreme parameter regimes can be accessed in comparison to DNS (Sprague et al. 2006;Calkins et al. 2013Calkins et al. , 2015;;Yan & Calkins 2022).In general, excellent agreement has been found between the results of plane layer asymptotic models and DNS when an overlap in parameter space is possible (Plumley et al. 2016;Yan & Calkins 2022).However, asymptotic models have not yet been developed for global spherical geometries.Towards this end, the present investigation utilises DNS in spherical shell geometries to determine the asymptotic scaling behaviour of various key quantities, including force balances, length scales, and convective and zonal flow speeds.
We consider Boussinesq convection in a rotating spherical shell with angular velocity Ω.The geometry is specified by the aspect ratio η = r i /r o , where r i is the radius of the inner sphere, and r o is the radius of the outer sphere.The fluid is forced via a temperature contrast ∆T between the inner and outer boundaries, and gravity varies linearly with radius.For this system, the convection dynamics are determined by the sizes of several non-dimensional parameters, including the Rayleigh number and Ekman number, defined by, respectively, where g o is the gravitational acceleration at the outer boundary, α is the coefficient of thermal expansion, H = r o − r i is the depth of the fluid layer, κ is the thermal diffusivity and ν is the kinematic viscosity.The most unstable state consists of convective Rossby waves that align with the rotation axis and drift in the prograde direction (Busse 1970).Asymptotic theory, valid in the limit Ek → 0, has shown that the critical Rayleigh number scales as Ra c = O(Ek −4/3 ), and the critical azimuthal (zonal) wavenumber scales as m c = O(Ek −1/3 ) (Roberts 1968;Busse 1970;Jones et al. 2000;Dormy et al. 2004).Thus, large Rayleigh numbers are needed to drive rotating convection and the subsequent motions become increasingly smaller scale as the Ekman number is reduced.
As the Rayleigh number is increased beyond critical, strong zonal flows develop in rotating spherical convection provided stress-free mechanical boundary conditions are applied on the inner and outer spherical surfaces (e.g.Aurnou & Olson 2001;Christensen 2002).These zonal flows are characterised by alternating regions of prograde and retrograde motion that are approximately invariant in the direction of the rotation axis.For a fixed value of Ek, the number of jets is controlled both by the Rayleigh number and the shell aspect ratio, η.In full sphere and small aspect ratio geometries (η ≲ 0.6), simulations typically find a single prograde jet in the equatorial region and retrograde jets at higher latitudes (Aurnou & Olson 2001;Christensen 2002;Lin & Jackson 2021).As η and Ra are increased there is a tendency for multiple jets to form at higher latitudes, leading to a banded structure that is reminiscent of the flows observed on the gas giant planets (Christensen 2001;Heimpel et al. 2005;Jones & Kuzanyan 2009;Gastine et al. 2014;Heimpel et al. 2022).The number of zonal jets that appear can be related to the Rhines length scale (Rhines 1975;Heimpel et al. 2005;Gastine et al. 2014).For sufficiently large Rayleigh numbers there is an eventual transition to a retrograde equatorial jet and prograde high latitude jets (Aurnou et al. 2007;Soderlund 2019).
Steady zonal flows are driven by Reynolds stresses and damped by global scale viscous stresses.Thus, the scaling behaviour of the zonal flow is intrinsically linked to the scaling of the correlations of the convective flows (Christensen 2002).Such correlations are not known a priori, though various scaling theories have been presented in the literature.For perfectly correlated small-scale velocity components, as relevant near the onset of convection, the zonal flow amplitude exhibits a quadratic dependence on the small-scale velocity (Aubert et al. 2001).However, Christensen (2002) found that correlations decrease with increasing supercriticality, which causes the quadratic scaling of the zonal flow to eventually break down.Later investigations have found that this loss of correlation in the small-scale velocity is a monotonically decreasing function of Ra (Showman et al. 2011;Gastine & Wicht 2012).For sufficiently large Rayleigh numbers, Christensen (2002) and Lin & Jackson (2021) find evidence that the zonal flow scaling approaches a 'diffusion-free' regime, so-called because of the lack of dependence on ν and κ, though within this regime the dynamics are no longer geostrophic since inertia and the Coriolis force are of the same order of magnitude.
Aside from their correlations, it is also important to understand the scaling of the convective flow speeds themselves.One scaling theory that is often invoked to explain convective flow speed scaling behaviour in rotating convection is the so-called Coriolis-Inertia-Archimedean (CIA) balance (e.g.Jones 2015;Aurnou et al. 2020).This balance predicts that the dominant convective length scale behaves like ℓ ∼ Ra Ek 1/3 , and the global scale Reynolds number should scale like Re = U H/ν ∼ EkRa/P r, where U is a characteristic flow speed, the reduced Rayleigh number is defined as Ra = RaEk 4/3 and the thermal Prandtl number is P r = ν/κ.King & Buffett (2013) analyzed length scales and flow speeds in a broad suite of numerical dynamo simulations in spherical geometries and concluded that viscous effects remained important in controlling these quantities.CIA theory is often contrasted with the scaling of the linear instability scale, ℓ ∼ Ek 1/3 .However, it is important to emphasise that, from the point of view of asymptotics, both the linear 'viscous' length scale and the length scale predicted by CIA theory are of the size O(Ek 1/3 ) given that Ra is an order unity asymptotic parameter.
Several previous investigations have tested these CIA scaling predictions in spherical geometries with no-slip boundary conditions (Gastine et al. 2016;Guervilly et al. 2019;Long et al. 2020), laboratory experiments in rotating cylindrical geometries (Madonia et al. 2021), as well as asymptotic models of plane layer convection with stress-free boundary conditions (Maffei et al. 2021;Oliver et al. 2023).Gastine et al. (2016) carried out a comprehensive survey of rotating spherical convection and found that the length scale for their smallest Ekman number cases (Ek = 3 × 10 −7 ) approached the Rhines scaling predicted by the CIA balance, and that the interior dissipation could also be approximated through a CIA balance.However, as we discuss in the present study, Gastine et al. (2016) did not investigate the convection and zonal flow separately.A similar approach was taken by Long et al. (2020) in which constant heat flux thermal boundary conditions were used.Guervilly et al. (2019) simulated rotating spherical convection with no-slip boundary conditions and P r = 0.01 using both a two-dimensional quasi-geostrophic model and three-dimensional DNS, and found that the length scale increases with Rayleigh number for all parameter ranges used; they find evidence of length scales and flow speeds approaching the CIA scaling predictions as the Ekman number is reduced.
The rotating convection experiments of Madonia et al. (2021) exhibited an increase of the integral length scale with Rayleigh number, but at a slower rate than that predicted by CIA theory.In asymptotic models, Maffei et al. (2021) found flow speed scaling behaviour consistent with CIA theory over a finite range of Rayleigh numbers; the deviation from CIA theory at large Rayleigh numbers was attributed to the effects of the large scale vortex (LSV) that is generated in this system.Oliver et al. (2023) also found that certain measures of the convective length scales show an increase with Ra, but at a rate that is slower than the exponent of 1/2.Importantly, however, the asymptotic models do not find a CIA force balance in the fluid interior.Instead, the buoyancy force is only comparable to the Coriolis and inertial forces within the thermal boundary layers, though viscous effects are equally important in these regions of the flow domain.Moreover, the ratio of the viscous force to the buoyancy force was found to be an increasing function of Ra, indicating that the CIA balance is never achieved in plane layer convection (Maffei et al. 2021;Oliver et al. 2023).The conclusion from these asymptotic studies is that convective length scales remain viscously controlled.This same conclusion was reached by Yan & Calkins (2022) who found similar behaviour in DNS of rapidly rotating convection driven dynamos in the plane layer geometry.An additional important finding in the studies of Madonia et al. (2021); Yan & Calkins (2022); Oliver et al. (2023) is that the viscous dissipation length scale remains approximately constant with increasing Ra -this indicates that length scale evolution in rotating convective turbulence is fundamentally different than non-rotating convective turbulence where the dissipation length decreases strongly with increasing Rayleigh number (e.g.Yan et al. 2021).We observe similar behaviour for the length scales and force balances in the spherical simulations reported in the present investigation.
The scaling behaviour of key quantities such as flow speeds and length scales is linked to the force balances in the governing equations.To our knowledge, the asymptotic scaling behaviour of the force balances (and terms in the heat equation) in spherical convection simulations have not been reported to date, though several previous dynamo studies have computed these forces over a range of parameters.For an Ekman number of Ek = 10 −4 in a spherical dynamo, Soderlund et al. (2012) noted that the Lorentz force was smaller than the inertial term and that the Lorentz force did not significantly alter the convective flow as compared to the purely hydrodynamic model.However, Soderlund et al. (2012) also noted that at smaller Ekman numbers, the Lorentz force seemed to make larger changes to the convection.In general, simulations find that the force balance in the mean equations is thermal wind to leading order, with the Lorentz force entering at higher order in the zonal component of the mean momentum equation (Aubert 2005;Calkins et al. 2021;Orvedahl et al. 2021).For the small scale convective dynamics, dynamo studies find that the zeroth order force balance is geostrophic, the first order force balance is between the ageostrophic Coriolis force, the buoyancy force and the Lorentz force, and inertial and viscous forces enter at the next order (Yadav et al. 2016) (note, however, these authors did not separate the mean and fluctuating dynamics).Here, we find a similar sequence of balances in the fluctuating momentum equation, though there does not appear to be an asymptotic difference between the buoyancy force, the viscous force and the inertial force, which is similar to plane layer rotating convection.Also in agreement with previous asymptotic studies, we find that the ratio of the viscous force to the buoyancy force is an increasing function of Ra.For the large-scale dynamics, we find that the flows are geostrophically balanced to leading order.
In this paper, we investigate the asymptotic behaviour of rapidly rotating Focus on Fluids articles must not exceed this page length convection and the associated zonal flows in a spherical shell with stress-free boundary conditions.Knowledge of such scaling behaviour is crucial for the development of asymptotic models.Overall, we find excellent agreement between the asymptotic predictions and the results of the nonlinear simulations.We develop a prediction for the asymptotic scaling of the amplitude of the zonal flow and conclude that the zonal flow must remain dependent on viscosity given that viscous stresses are the sole saturation mechanism for this component of the flow.The paper is organized as follows.In section 2, we describe the model and governing equations.We give a brief overview of the asymptotic theory in section 3 and numerical results are analysed in section 4. A discussion is provided in section 5.

Model
The governing equations consist of the conservation laws for momentum, thermal energy and mass.We non-dimensionalise these equations using the length H, the large-scale viscous diffusion time H 2 /ν and the temperature scale ∆T .With this non-dimensionalisation, the governing equations are given by 3) where u = ⟨u r , u θ , u ϕ ⟩ is the fluid velocity, T is the temperature, P is the pressure and r is radius.The 'axial' direction points in the direction of the rotation axis and the axial and radial unit vectors are denoted by z and r, respectively.In all of the simulations presented here we fix P r = 1.
The boundary conditions are impenetrable (u r = 0), stress-free and fixed temperature.We use the code Rayleigh to numerically solve the governing equations (Featherstone et al. 2022).Rayleigh is a pseudo-spectral code which uses spherical harmonics to represent data on spherical shells and Chebyshev polynomials to represent data in the radial direction.A 2/3 de-aliasing is used for both the spherical harmonics and the Chebyshev polynomials.Time-stepping is carried out using a second-order semi-implicit Crank-Nicolson method for the linear terms and a second-order Adams-Bashforth method for the nonlinear terms.The time step is chosen adaptively to maintain numerical stability.

Notation and outputs
Due to the symmetry of the model setup around the rotation axis, it is convenient to define mean and fluctuating components of some scalar quantity X relative to an azimuthal or zonal average, i.e.
We define the Reynolds number as where ⟨•⟩ denotes a volume average and (•) t denotes a time average.We further define the mean and fluctuating (convective) Reynolds numbers as, respectively, In all of the simulations, we find that the zonal component of the mean flow dominates and we therefore refer to the mean Reynolds number as the 'zonal' Reynolds number.We analyse several different length scales in this paper.Christensen & Aubert (2006) (see also Gastine et al. 2016;Long et al. 2020) define the spherical harmonic length scale as where l max is the maximum spherical harmonic degree in the simulation and E m l is the radially averaged kinetic energy density of spherical harmonic degree l and order m.Thus, the volume averaged kinetic energy density is given by (2.8) Due to the strong influence of the zonal flow in our simulations, we separate the spherical harmonic length scale into a zonal and a non-zonal component.We define these length scales as (2.9) We define the fluctuating and mean Taylor microscales as (2.10) respectively.The Taylor microscale can be considered a viscous dissipation length scale since it characterises the length scale at which viscous effects become important.
The Nusselt number is calculated according to .11)where (•) ϕ,θ,t is a shell and time average, and T c is the conductive temperature profile, which satisfies (2.12) We also define the viscous dissipation rates of the mean and fluctuating velocity fields according to (2.13) respectively.

Theory
Here we provide arguments for the scaling behaviour of various quantities.We use the term asymptotic to mean rotationally constrained motions in which the Ekman number and Rossby number, Ro = U/ (2ΩH), are both small relative to unity, i.e. (Ro, Ek) ≪ 1.The dimensional flow speed U characterises the magnitude of fluctuating velocity field.We expect many aspects of the asymptotic theory for rotating convection in spherical geometries presented in Jones et al. (2000) and Dormy et al. (2004) to hold here, though some of these scalings must be modified to account for nonlinear terms in the governing equations.In particular, the scaling of the convective flow speeds and temperature perturbation need to be reduced by a factor of Ek 1/3 relative to Jones et al. (2000) and Dormy et al. (2004), though this difference does not influence the leading order force balance in the fluctuating momentum equation when zonal flows are weak.
The scaling of the large-scale zonal flow depends on the small-scale convective velocity, so we first consider theoretical scaling laws for the convective velocity and the corresponding convective length scale.Such scaling laws can be found by examining the fluctuating momentum equation and the fluctuating heat equation, which are given by, respectively, We find that the terms in brackets can be large compared to some terms due to the large amplitude of the zonal flow.However, summing the terms in brackets leads to results smaller than the individual terms in the brackets, which physically means that the Lagrangian time derivative is smaller than the Eulerian time derivative.We will therefore consider the sum of the terms in brackets rather than each individually.It is well known from linear theory that the Coriolis force and pressure gradient force are dominant terms in the limit (Ro, Ek) → 0. We will see below that the zonal flow can also modify the leading order force balance.In order to study the first-order effects, we eliminate the pressure gradient by taking the curl of the fluctuating momentum equation.This operation yields Assuming length scales of order one in the z-direction, length scales of order one for azimuthally averaged terms, and length scales of order ℓ otherwise, the vorticity equation and heat equation can be approximately written as where factors of order one have been dropped.We now follow Aurnou et al. (2020) and assume a balance between u ′ T and u ′ T ′ ℓ −1 in the temperature equation.Using T = O(1) yields T ′ ∼ ℓ.Plugging this relation in for T ′ in the momentum equation and assuming a CIA balance where the fluctuating-fluctuating advection term is used yields which can be solved for the convective velocity and length scale to give We can rewrite these expressions in terms of the reduced Rayleigh number as Here again we note that both the CIA and viscous length scale have the same Ek 1/3 dependence.One of the points that we stress in the present study is that all length scales in this system are viscously selected to leading order, with order one variations away from this viscous scale since Ra = O(1) (e.g.Yan & Calkins 2022;Oliver et al. 2023).
We can understand why these different balances produce the same Ekman dependence by making a different set of assumptions than the assumptions used for a CIA balance.The first assumption we make is that the Ekman dependence for any term can be written in terms of a power law.We will assume that the Ekman dependence of the convective flow speeds, the length scale, and the fluctuating temperature can be written as The second assumption is that ratios of certain terms in the momentum and heat equations do not change when the Ekman number is reduced.We will assume that the advection of fluctuating velocity by fluctuating velocity, viscosity, and buoyancy follow the same Ekman number scaling in the fluctuating momentum equation.In the heat equation we will assume that conduction and advection of the mean temperature by the fluctuating velocity follow the same Ekman number scaling.These assumptions produce the system of equations given by (3.9) Note that we have assumed the mean temperature and the length scale of the mean temperature do not depend on the Ekman number.Solving this system of equations yields x u = −1/3, x ℓ = 1/3, and x T = 1/3, which implies that u ′ = O Ek −1/3 , ℓ = O Ek 1/3 , and T ′ = O Ek 1/3 .Note that the Ekman number dependence derived here is the same as the Ekman number dependence derived using the CIA balance written in terms of the reduced Rayleigh number.Therefore, we can get the Ek 1/3 scaling for the length scale by assuming that various terms follow the same Ekman number scaling without actually assuming any balances a priori.These scalings are equivalent to the scalings used to derive the asymptotic model of rotating convection in a Cartesian geometry (e.g.(Sprague et al. 2006)).
An asymptotic constraint on the amplitude of the zonal flow can now be obtained upon examination of the mean momentum equation.The mean momentum equation is given by The zonal component is then (3.11) where the square brackets and corresponding subscript are used for brevity.
The above equation allows for a straightforward interpretation of the zonal flow dynamics.Time dependence and advection of the zonal flow by mean meridional flows are captured by the first two terms on the left-hand side.The term u ′ • ∇u ′ ϕ includes the divergence of the Reynolds stresses and acts as the primary source of the zonal flow.The zonal component of the mean Coriolis force can be written as [ z × u] ϕ = cos θ u θ + sin θ u r such that only meridional circulation appears in this term.Averaging equation (3.11) along the z-direction eliminates the Coriolis term due to continuity, and averaging equation (3.11) in time eliminates the time derivative term.Therefore, averaging equation (3.11) along the z-direction and in time leaves only the two advective terms and the viscous term.For a careful derivation of the balance between the averaged advection and viscosity in spherical coordinates, see Dietrich et al. (2017).Other studies have found that the mean-mean advection term is small (e.g.Dietrich et al. 2017), so there must be a balance between the fluctuating advection term and the viscous term.Letting (•) ϕ,z,t denote an average over ϕ, z, and time, the preceding argument implies that (3.12) The above balance is expected to hold for all values of Ek and Ra; the zonal flow is therefore intrinsically dependent on viscosity and we should not expect its scaling behaviour to be 'diffusion-free'.We note that a similar balance holds for mean flows in planar geometries (e.g.Nicoski et al. 2022).Finally, since averaged quantities vary on order one length scales, this balance suggests where C R represents the correlation of the fluctuating velocity components.We might expect that this correlation gets weaker as the reduced Rayleigh number is increased and the flow becomes less constrained by rotation.Christensen (2002) confirmed this weakening of the correlation as the Rayleigh number is increased by directly calculating the correlation and by noting that the zonal flow amplitude only scales as the square of the small-scale velocity for Rayleigh numbers near the onset of convection.However, the correlation cannot be greater than one and we would not expect the correlation to get closer to zero as Ek is decreased with Ra fixed.We might therefore anticipate that C R does not depend on the Ekman number for small values of the Ekman number, in which case we can use equation (3.13) to find a scaling relationship for the zonal flow with Ekman number.The asymptotic analysis predicts that the convective velocity scales as u ′ = O Ek −1/3 , which implies that (3.14) thus indicating that the zonal flow is intrinsically dependent on viscosity, albeit in an asymptotic sense.For time-dependence in the zonal flow, we require that the time derivative is comparable to the viscous force so that which, along with equation (3.14), indicates that the zonal flow time scale is O(1) in our non-dimensional, large-scale viscous diffusion units.Thus, the zonal flow varies on a large-scale viscous diffusion time.This property is one of the reasons that computations of rotating spherical convection with stress free boundary conditions are so demanding -very long integration is necessary to saturate the amplitude of the zonal flow and reach a statistically stationary state.
The asymptotic constraint on the zonal flow amplitude allows for additional insight into the force balance by which it is constrained.In the limit Ek → 0, the Rayleigh number must scale as Ra = O Ek −4/3 to generate convection (Roberts 1968).Along with the fact that the magnitude of the mean temperature is independent of the Ekman number, this indicates that the mean buoyancy force scales as O Ek −4/3 .The radial and co-latitudinal components of the mean Coriolis force both contain u ϕ , thus indicating that these components scale as O Ek −5/3 , which is larger than the mean buoyancy force by a factor of O Ek −1/3 .Thus, the zonal flow is geostrophically balanced to leading order, i.e.
It is informative to compare the scaling behaviour of zonal flows that are geostrophically balanced with a zonal flow that is in thermal wind balance (e.g.Calkins et al. 2021).An order of magnitude estimate for the scaling of the thermal wind component of the zonal flow can be obtained if we balance the mean Coriolis force with the mean buoyancy force, Using the definition of the reduced Rayleigh number this becomes And since Ra = O(1) this implies u tw ϕ = O(Ek −1/3 ), which is of the same order as the convective flow speeds.Thus, zonal flows that satisfy a thermal wind balance are substantially weaker than those that are geostrophic.Interestingly, the above scaling also provides an estimate for the scaling of the thermal wind with Rayleigh number and represents a 'diffusion-free' scaling in the sense that it indicates that the thermal wind does not depend on either ν or κ.The thermal wind scaling seems to be loosely consistent with the zonal flows present in the dynamos of (Calkins et al. 2021); when a magnetic field is present the Lorentz force strongly damps the geostrophic component of the zonal flow.

Numerical Results
4.1.Overview To test the theoretical arguments given in the previous section, we perform a suite of direct numerical simulations of convection across a range of Ekman number and Rayleigh number.We consider the aspect ratios η = 0.35 and η = 0.7.The smallest Ekman numbers that were simulated are Ek = 10 −6 and Ek = 3 × 10 −5 for the η = 0.35 and η = 0.7 aspect ratios, respectively.Details of the simulations are contained in Tables 1-2 in the appendix.
The qualitative nature of the zonal flows was similar for much of the parameter space covered, and representative cases for both shell thicknesses are shown in figure 1.The zonal flows we observe in this study are similar to zonal flows observed in previous works.These zonal flows are characterised by a nearly invariant structure in the axial (z) direction, and consist of a single prograde jet at the equator with retrograde jets at higher latitudes.A small number of the thin shell (η = 0.7) cases at low Ekman number and high Rayleigh number developed high latitude jets, as shown in figure 1(c).A subset of our cases exhibit relaxation oscillations in which the convection mainly occurs during short bursts; this behaviour was also observed in previous work (e.g.Christensen 2002).Movie 1 from the supplementary material shows the radial velocity in the equatorial plane over the course of one relaxation oscillation with the convective Reynolds number shown for reference.From this movie, we see that during times of weak convection, the convection is strongest near the inner boundary in the equatorial plane.However, during times of strong convection, the convection fills the whole region in the equatorial plane.
We find that using Ra, as opposed to the supercriticality measure, Ra/Ra c , results in improved collapse of our data when comparing with asymptotic predictions; Christensen (2002) also found that using Ra improved the collapse of some data.This effect likely arises from the slow rate of convergence of the critical Rayleigh number to the predicted asymptotic scaling of Ra c ∼ Ek −4/3 (e.g.Dormy et al. 2004;Barik et al. 2023).

Flow speeds
Global rms values of both the fluctuating and mean velocity are computed to determine their scaling behaviour with respect to the Ekman number.Note that the mean velocity is dominated by the zonal flow (i.e. the ϕ-component of the mean velocity).Figure 2(a) shows the fluctuating Reynolds number for both the η = 0.35 cases and the η = 0.7 cases.Figure 2(b) shows the asymptotically rescaled fluctuating Reynolds number, i.e.Re c = Ek 1/3 Re c , for the two different aspect ratios.We find that the rescaled data is order unity and collapses onto a single curve, which supports the u ′ = O Ek −1/3 asymptotic scaling for the fluctuating velocity.However, we note that the Ekman number scaling of the convective Reynolds number might be time dependent.If we calculated the convective Reynolds number using data only during the convective peaks of the relaxation oscillations, we would obtain a steeper scaling closer to Ek −1/2 .This suggests that the time series for the convective velocity becomes more strongly peaked at lower Ekman number.Note that we expect deviation from this asymptotic scaling behaviour as the system loses rotational constraint; this deviation is particularly noticeable in figure 2(b) for the high Rayleigh number aspect ratio, which suggests these cases may be scaling as Re c ∼ Ra.However, this scaling behaviour may be localised in Ra space.For sufficiently small Ekman number and large Rayleigh number, the compensated plot for Re c Ra −3/2 collapses the data well, though the scaling appears slightly weaker than Ra 3/2 which suggests that the convective Reynolds number scales approximately as Re c ∼ Ra 3/2 in this regime.
Figure 3(a) shows the mean Reynolds number as a function of Ra, and figure 3(b) shows the corresponding asymptotically rescaled mean Reynolds number, Re z = Ek 2/3 Re z .As mentioned previously, the mean flow is dominated by the zonal component in all of our simulations.While there is clearly some spread in the rescaled data for the thick shell cases, there is an indication that the data collapses to an asymptotic state as Ek → 0.Moreover, the rescaled values are order unity.There appears to be better collapse for the thin shell cases, indicating that the fluid depth may play an important role.Also note the excellent collapse for the three thin shell data points near Ra ≈ 60 -the two lower Ekman number cases of these three develop prograde high latitude jets as shown in figure 1(c).Taken together, the data seems to support that the zonal flow scales as Ek −2/3 in the rapidly rotating regime.We note that Gastine et al. (2014) found an empirical scaling for the zonal Rossby number of Ro zon ∼ Ra 0.6 Ek 0.99 , which, when converted to a Reynolds number and using Ra ∼ Ek −4/3 gives Re z ∼ Ek −0.81 .This result broadly agrees with the scaling of Re z ∼ Ek −2/3 derived in this paper.
While the balance between viscosity and Reynolds stresses predicts a zonal flow scaling of Ek −2/3 , this balance is unable to predict how the zonal flow scales with the reduced Rayleigh number due to the fact that the correlation of the fluctuating velocity components is an unknown function of Ra.Thus, we make an empirical fit of Re z with respect to the reduced Rayleigh number, which is shown in figure 3(b).We find that a line does a good job of fitting our small Ekman number cases.However, our larger Ekman cases at η = 0.35 are run to larger values of the reduced Rayleigh number and are not well fit by a line at larger values of Ra.One possibility for this effect is that these larger Ekman number cases at large Rayleigh number are no longer in the rapidly rotating regime and therefore follow a different trend.Another possibility is that the affine dependence we have used to fit the data only holds near the onset of convection and that the zonal Reynolds number follows a different trend for large enough values of the reduced Rayleigh number, even as the Ekman number is decreased.
It is also interesting to consider the relationship between the zonal and convective Reynolds numbers.Equation (3.13) suggests that the relevant quantity to consider is Re z /Re 2 c , which we show in figure 4. We predicted that this ratio is independent of the Ekman number, which we see is a good approximation for our range of parameters.Christensen (2002) noted that the relation Re z ∼ Re 2 c holds only when the correlation between the fluctuating velocity components is constant, which occurs near the onset of convection.Anelastic simulations of rotating spherical convection also find that Re z /Re 2 c is nearly constant near the onset of convection, and that Re z /Re 2 c decreases with increasing Ra for sufficiently large Rayleigh number (Gastine & Wicht 2012).Our simulations are also consistent with this behaviour, showing that Re z /Re 2 c is approximately constant up to Ra ≈ 15.For Ra ≳ 15, we observe that Re z /Re 2 c decreases with increasing Ra.This trend of Re z /Re 2 c decreasing as Ra is increased does not show a systematic dependence on the Ekman number, since all data points appear to collapse well.Because Re z /Re 2 c is related to the correlation of the convective velocity components, this suggests that the correlation of the convective velocity components decreases as Ra is increased, even at asymptotically small Ekman numbers.

Length scales
In this section, we compute length scales for both the mean flow and the smallscale convection.To provide a physical picture of how these length scales vary with Ekman number, we show snapshots of the radial velocity in the equatorial plane for three different Ekman number cases with Ra ≈ 50 in figure 5.As expected, the typical length scale of the convection decreases with decreasing Ekman number.
One way to quantitatively study how the length scale varies with Ekman number is to consider how the kinetic energy power spectrum varies with Ekman number.We define the sum of the kinetic energy power spectrum over the order m and the normalization of the kinetic energy power spectrum summed over the order m as (4.1) Due to the influence of the zonal flow, we again only sum over m ⩾ 1 modes to remove the influence of the zonal flow.Figure 6 shows how these kinetic energy spectra vary with Ekman number.The critical azimuthal wavenumber m c is shown for reference, and m c takes the values 6, 8, 11, 16, and 23 for the cases of Ek = 3×10 −4 , Ek = 10 −4 , Ek = 3×10 −5 , Ek = 10 −5 , and Ek = 3×10 −6 , respectively.These values were calculated using the eigenvalue solver Kore (Barik et al. 2023).We see that the smaller Ekman number cases have more power in each mode and extend to higher values of l in comparison to larger Ekman number cases.We attempt to collapse this data by multiplying the degree l by Ek 1/3 , the expected scaling for the length scale, and multiplying the amplitude of the data by Ek 1/3 .We choose the amplitude scaling such that the area under the kinetic energy power spectrum, the volume averaged kinetic energy density, scales as Ek −2/3 , which is consistent with the scaling of our convective Reynolds number.This collapsed data for the kinetic energy power spectrum is shown in figure 6(b).We see that the collapse of the data is reasonable for values of l greater than m c , but not as good for values of l less than m c .This discrepancy might suggest that large length scales are following a different scaling than anticipated.
We investigate the behaviour of the kinetic energy power spectrum for the small l modes in figure 6(c,d).The power spectrum for these small l modes appears to shift more slowly to the right with Ekman number than the large l modes.Several papers on the linear onset of convection predict a longer length scale in the cylindrical radial direction than in the azimuthal direction.Dormy et al.(c) Ek = 3 × 10 −6 , Ra = 1.2 × 10 9 ( Ra ∼ 51.9).
(2004) predicts a Ek 2/9 length scale for cases with differential heating, and a Ek 1/6 length scale for cases with certain boundary conditions and heating conditions.We show the collapse of our data for the predicted Ek 2/9 length scale in figure 6(d), which does a reasonable job of collapsing our data.However, we do not have a large enough range in Ekman number to discern between a Ek 1/6 and a Ek 2/9 scaling.
Figure 7 shows the kinetic energy spectra for the η = 0.35, Ek = 10 −5 cases averaged in time and radius.Figure 7(b) shows the spherical harmonic spectra normalized to have area one.We find that the kinetic energy contained in high l values increases slightly with increasing Ra, but that the overall shape of the spectra does not change significantly.This behaviour indicates that small-scale convection becomes more prominent as the Rayleigh number is increased.
The computed length scales are all shown in figure 8.The spherical harmonic length scale is shown in panels (a) and (b), and the Taylor microscale is shown in panels (c) and (d); the asymptotically rescaled lengths are shown in panels (b) and (d).The spherical harmonic length scale is observed to decrease rapidly for Ra ≲ 50, then levels off for larger values of Ra; not surprisingly this behaviour is consistent with the trends observed in the kinetic energy spectra.There may be a trend suggesting that this length scale converges to a nearly constant value at the largest values of Ra, though this behaviour may be occurring outside of the rapidly rotating regime.Both the decreasing and constant dependence of the length scale on Rayleigh number observed here is in contrast to the prediction made by the CIA balance where the length scale increases as the Rayleigh number increases.Figure 8(b) shows the asymptotically rescaled spherical harmonic length scale.As with the flow speeds, this rescaled quantity is order unity and we find significantly less scatter in the data when viewed in this rescaled coordinate, suggesting that this particular length scale does approximately scale as Ek 1/3 .However, we also note that a subtle Ekman number dependence still exists in the rescaled length scale; larger Ekman number cases seem to level off at smaller rescaled length scales in comparison to the smaller Ekman number cases.This behaviour might indicate that the system is either converging slowly toward the O(Ek 1/3 ) length scale with decreasing Ekman number, or that the asymptotic dependence of this length scale is slightly weaker than Ek 1/3 .Our data for the kinetic energy power spectrum indicates that a second length scale may be present at small l values which follows a different Ekman number scaling.Thus, the spherical harmonic length scale may have a composite asymptotic dependence on Ek since it may be capturing both the Ek 2/9 and Ek 1/3 asymptotic scalings.Figure 8(c,d) shows the Taylor microscale, where we see that it follows a very similar trend in comparison to the spherical harmonic length scale.The collapse of the rescaled quantity indicates that the viscous dissipation length scale is consistent with a Ek 1/3 dependence across the entire range of parameters used.The Taylor microscale appears to decrease with Ra over a wider range of Ra than the spherical harmonic length scale, which becomes approximately constant with Ra.If we interpret the spherical harmonic length scale as the energy containing length scale (similar to an integral length scale), then our computed length scales indicate that the scale separation present in these simulations is relatively small across the entire range of parameters since ℓ ′ sh and ℓ ′ tm are not too dissimilar in value.For instance, at Ra ≈ 50, ℓ ′ sh ≈ 4 and ℓ ′ tm ≈ 0.9.We note that the trends observed in both length scales contrast sharply with computed length scales in non-rotating Rayleigh-Bénard convection, where both the energy containing length scale and the dissipation length scale decrease rapidly with increasing Rayleigh number (e.g.Yan et al. 2021).
In order to study the behaviour of the Ek 2/9 length scale found in the power spectrum, we define a length scale based on the degree l peak where the kinetic energy power spectrum (again with the m = 0 mode removed) reaches a maximum.The peak length scale is then given by ℓ ′ peak = π/l peak .We fit a 30th degree polynomial to our kinetic energy power spectrum in order to smooth l peak , although we obtained similar results when we fit our data with a 15th degree polynomial.Figure 9(a) shows the peak length scale and figure 9(c) shows the peak length scale rescaled by Ek −2/9 .We note that while the peak length scale approximately scales as Ek 2/9 , we obtain a somewhat better fit for Ek 1/6 as shown in figure 9(d).Regardless, this peak length scale seems to scale differently than ℓ ′ sh , which suggests that ℓ ′ peak does in fact correspond to a different length scale than the small-scale convective length scale.We also note that the length scale ℓ ′ peak increases with increasing Rayleigh number, in contrast to ℓ ′ sh and ℓ ′ tm , which either decrease or remain constant as the Rayleigh number is increased.Figure 9(b) shows the peak length scale as a function of the Rossby number characterising the convective flow speeds, Ro = EkRe c , where it can be seen that ℓ ′ peak approximately scales as Ro 1/4 , which is in contrast to the Ro 1/2 scaling that is predicted from the CIA balance.To ensure the validity of our results with polynomial interpolation, we also estimated the peak length scale by replacing E m l in equation 2.9 with (E m l ) 10 , which weights the length scale more towards the peak.Using this estimate of the peak length scale, we found qualitatively similar results to what is shown in figure 9. We now consider length scales of the zonal flow, which are shown in figure 10.Unlike the convective length scales, we expect these mean length scales to be order unity and approximately independent of the Ekman number.The mean length scale reaches a maximum at approximately Ra ∼ 25 and, for larger values of the reduced Rayleigh number, decays slowly over the range of Rayleigh numbers used in this study.We also computed a Taylor microscale for the mean flow, as shown in 10(b); although this length scale is slightly smaller than the spherical harmonic length scale, it shows a similar behaviour with Ek and Ra in comparison to the spherical harmonic length scale.

Heat equation balances
The asymptotic scaling behaviour of terms in both the momentum and heat equations are linked.For this purpose, we study the scaling of various terms in both equations beginning with the heat equation in the present section.Figure 11 shows the rms of the fluctuating temperature for all of the thick shell cases.We observe a systematic decrease in the magnitude of the fluctuating temperature as the Ekman number is reduced.Figure 11(b) shows that the fluctuating temperature is well described by the relationship T ′ = O(Ek 1/3 ) for fixed Ra, as predicted by asymptotic theory.A line gives a decent approximation for how the fluctuating temperature varies with Ra, in contrast to the CIA balance, which predicts a reduced Rayleigh number dependence of T ′ ∼ Ra 1/2 .Figure 11(c) shows the compensated fluctuating temperature using the reduced Rayleigh number dependence predicted by the CIA balance: Ra 1/2 .We see that the data is not very well described by this Rayleigh number dependence; for none of the Ekman numbers does there exist a large range of Ra where the compensated data is horizontal.We also note that for the Ek = 3 × 10 −4 cases, the largest Ra cases exhibit a slight decrease in the rms value of T ′ as compared to lower Ra cases, so T ′ is not strictly increasing with Ra in our data.
Figure 12 shows how the various terms from the fluctuating heat equation depend on Ekman number and reduced Rayleigh number.Figure 12(d,e,f) shows how well the asymptotic prediction of the scaling collapses the data.The diffusion The symbols are the same as defined in figure 2. Thermal terms term and the advection of the mean temperature by the fluctuating velocity term are well collapsed by the predicted Ek −1/3 scaling, although the advection of the fluctuating temperature by the fluctuating velocity is not as well collapsed by the asymptotic scaling.The data suggests that a slightly stronger dependence on the Ekman number is present since the rescaled data from the smaller Ekman number cases lies above that of the larger Ekman number cases.However, we note that the magnitude of the rescaled terms are all comparable, thus providing support for the asymptotic theory.
Figure 13 shows how the various terms from the fluctuating heat equation vary with the reduced Rayleigh number for two fixed values of the Ekman number.The large magnitude of the time derivative term is a consequence of the zonal flow; advection by the zonal flow is not balanced and causes large accelerations.Therefore, the sum of the advection by the zonal flow and the time derivative is smaller by comparison and this sum is approximately balanced by the term u ′ • ∇T ′ .Another interesting observation from figure 13 is that there is a near balance of the diffusion term and the advection of the mean temperature term for a wide range of Ra.This balance was also noted by Maffei et al. (2021).

Force balances
In this section, we investigate the scaling behaviour of the forces present in the radial component of the fluctuating momentum equation (similar scaling was observed in the other two components).Figure 14 shows the time-average and rms over the entire domain of the viscous force, buoyancy force, and the fluctuating advective term.For the fluctuating advective term, we remove the spherically symmetric l = 0 mode as this mode is not dynamically relevant.From the arguments in the theory section, we expect these three terms to all scale as Ek −1 in our non-dimensional units.We test these scalings by multiplying the data in figure 14(a,b,c) by Ek, which is shown in figure 14(d,e,f).We see that the predicted Ekman number scalings are consistent with the data.Similar to the  Forces heat equation terms in the previous section, the collapse for the advective term is not as good as the collapse for the other terms.This difference either indicates that the advective term is converging more slowly than the other terms or that the advective term is converging to a slightly different scaling than predicted by asymptotic theory.We also note that the scaling of the advective term is time- dependent; removing the intervals of time in which strong convective bursts are occurring in the time series produces a scaling for the fluctuating advective term that more closely follows the predicted Ek −1 scaling.This effect suggests that the relaxation oscillations are playing some role in the scaling, which may be due to the larger Rossby numbers that occur during these times.The Coriolis force is not shown since its scaling follows from the scaling of the fluctuating velocity given in figure 2 -in our units it scales as Ek −4/3 .These scalings suggest that, for a fixed value of the reduced Rayleigh number, the ratio of the viscous force, buoyancy force, and fluctuating advection term remain approximately fixed as the Ekman number is decreased.On the other hand, the ratio of either the viscous force, the buoyancy force, or the fluctuating advective term to the Coriolis force will scale as Ek 1/3 .While we have only shown the radial component of the fluctuating force balance, we note that the θ and ϕ components of the fluctuating force balance show similar behaviour.We also tested removing the thermal boundary layers, and did not observe a qualitative change in any of the trends shown.
Some terms from the momentum equation follow a stronger Ekman number scaling due to the presence of the zonal flow, and these terms are shown in figure 15.The zonal flow strongly advects the fluctuating velocity, which is unbalanced and causes large accelerations in the small-scale fluid structures.Because the zonal flow scales as Ek −2/3 , the fluctuating velocity scales as Ek −1/3 , and the length scale of the fluctuating velocity scales as Ek 1/3 , we expect the advection by the zonal flow to scale as Ek −4/3 , which we also expect for the scaling of the time derivative.Figure 15(d,e) shows the collapse of the mean advection term and the time derivative term for this scaling.We see that both the advection by the zonal flow and the time derivative follow a Ek −4/3 scaling, which we note is the same scaling followed by the Coriolis force.Therefore, the advection by the zonal flow and the time derivative appear at leading order asymptotically, and the sum of these two terms is somewhat smaller than either individually, as shown in figure 15(c,f) and in figure 16.
Figure 16 shows how the different terms in the fluctuating radial momentum equation depend on the reduced Rayleigh number for two Ekman numbers.For both Ek = 3 × 10 −4 and Ek = 10 −5 , the buoyancy force is about an order of magnitude larger than viscosity or advection at small Ra, and both viscosity and advection grow more rapidly with Rayleigh number as compared to the buoyancy force.For Ek = 3×10 −4 , advection becomes larger than buoyancy near Ra ≈ 140 and continues growing larger than buoyancy.For Ek = 10 −5 , we do not reach large enough Rayleigh numbers for advection to become as large as buoyancy.Therefore, we do not observe a CIA balance in our simulations since the buoyancy force and advection scale differently with Rayleigh number for the parameter space we surveyed.
In figure 17 we plot the ratio of the viscous force to the buoyancy force, and the ratio of the fluctuating advective term to the buoyancy force.We find that the ratio of the viscous force to buoyancy is an increasing function of the Rayleigh number, which suggests that viscosity does not become negligible at more extreme parameters.This behaviour is in contrast to the diffusion-free scaling that is used in many previous studies, which assume that viscosity is negligible.Therefore, the trends we observe in our data do not support viscosity becoming negligible, although it is possible that this trend changes at higher values of the reduced Rayleigh number than we achieved in this study.

Heat flow and dissipation
Figure 18 shows the Nusselt number for all of the η = 0.35 cases from this paper with two scalings shown for reference.The solid line shows the predicted scaling for diffusion-free heat transport, N u ∼ Ra 3/2 , and the dashed line shows the scaling N u ∼ Ra.It appears that our data has not quite yet reached the predicted N u ∼ Ra 3/2 scaling, although some of our Ek = 3 × 10 −6 cases appear to be close to this scaling.This might suggest that the N u ∼ Ra 3/2 scaling can be reached at lower Ekman number and higher Rayleigh number.
Figure 19(a) shows the dissipation by both the convective flow and the mean flow.The dissipation by the convective flow is greater than the dissipation by the mean flow for all of the cases we studied.Figure 19(b) shows that the dissipation follows the expected scaling of Ek −4/3 , and that the convective dissipation follows the scaling Ek −4/3 slightly better than the mean flow.Figure 20 shows the ratio of the dissipation by the mean flow to the dissipation by the convective flow.We see that for all Ekman numbers, this ratio reaches a maximum when Ra ∼ 30, and that for higher values of the reduced Rayleigh number the ratio of mean to convective dissipation decreases.This trend might suggest that for large values of the reduced Rayleigh number, the dissipation is mainly controlled by the convection.There is also a trend where the maximum value of the ratio of the mean dissipation to the convective dissipation increases as the Ekman number is decreased.However, this trend is rather weak, with only about a factor of two change in the ratio while the change in Ekman number was two orders of magnitude.

Discussion
A systematic investigation of rotating convection in a spherical shell with stressfree boundary conditions was carried out.The scaling behaviour of various quantities was investigated for varying Rayleigh number and Ekman number, and two different aspect ratios were employed.An emphasis was placed on characterising the asymptotic nature of the system as the Ekman number is reduced.Here we have utilised the known scalings obtained from both the linear theory of rotating spherical convection (Jones et al. 2000;Dormy et al. 2004) and the closely related nonlinear reduced model for rotating Rayleigh-Bénard convection (e.g.Sprague et al. 2006).Overall we find good agreement between the asymptotic scalings and the DNS.
In general, we find that the asymptotic nature of the system can be demonstrated when plotting asymptotically rescaled quantities as a function of the reduced Rayleigh number, Ra.Using this approach, we find that the convective flow speeds, as measured by the large-scale Reynolds number, scale as Re c = O(Ek −1/3 ) for both thin and thick shell simulations.This scaling can be deduced from the governing equations by requiring the convection to be geostrophically balanced to leading order, with the buoyancy force entering the next asymptotic order.As discussed by previous work, the amplitude of the zonal flow is limited solely by large-scale viscous diffusion; thus, we do not anticipate a diffusion-free scaling for this component of the velocity field at any combination of Ekman or Rayleigh number in the rotationally constrained regime.By balancing the large-scale viscous diffusion of zonal momentum with the Reynolds stresses, and noting that Re c = O(Ek −1/3 ), we deduce that the zonal flow speeds scale as Re z = O(Ek −2/3 ).The numerical data supports a trend toward this scaling as Ek → 0.Moreover, we find that because of this diffusion-Reynolds-stress-balance, the correlations in the convective flow field are independent of the Ekman number.As noted in previous work (Christensen 2002;Gastine & Wicht 2012), these correlations decrease rapidly with increasing Rayleigh number (approximately as ∼ Ra −3/2 ).
Kinetic energy spectra exhibit a trend with increasing Rayleigh number in which energy builds at nearly all spatial scales, with no clear preference for scales significantly larger than the linear convective instability, as indicated by a Ek 1/3 or Ek 2/9 scaling.Using the length scale calculation approach of Christensen & Aubert (2006) and Gastine et al. (2016), we find that the dominant convective length scale decreases with increasing Ra and scales approximately as Ek 1/3 .At the largest accessible Rayleigh numbers, but also the largest Ekman numbers, this length scale appears to saturate.The Taylor microscale also scales as Ek 1/3 and shows a very similar trend in comparison to the spherical-harmonic-based length scale; after initially decreasing with increasing Ra it tended to level off at the largest accessible Rayleigh numbers.Moreover, both of these length scales remain comparable in value across the entire range of parameters that were simulated, and the collapse of the asymptotically rescaled length scales suggests that this trend continues for more extreme parameter values.The largest length scales present in the fluctuating kinetic energy spectra scale weaker than the primary convective instability and may be consistent with the Ek 2/9 radial scale from linear theory (Dormy et al. 2004), although a Ek 1/6 scaling is also consistent with out data.This Ek 2/9 length scale increases with increasing Rayleigh number, and approximately follows the scaling ℓ ′ peak ∼ Ro 1/4 c .Our results suggest that there are multiple length scales present in rotating spherical convection, and that these length scales exhibit different dependencies on the input parameters.Nevertheless, all of these length scales appear to depend on the Ekman number and are therefore viscously controlled to some degree.
Several previous investigations have utilized no-slip boundary conditions that yield different scaling behaviour of the length scales in comparison to the present study.Gastine et al. (2016) computed length scales over a broad range of parameters with η = 0.6.These authors used a similar definition of the sphericalharmonic-based length scale (equation (2.7)), but also included m = 0 in their calculations.In general, they found non-monotonic behaviour in their values of ℓ sh , though they did find an increase in the convective length scale provided Ra was sufficiently large and Ek was sufficiently small.We have found that including the m = 0 component in the length scale calculation causes the length scale to increase with increasing Rayleigh number over some range of Ra, suggesting it is necessary to remove this component when deducing length scale behaviour of the small scale convection.Using constant heat flux thermal boundary conditions, Long et al. (2020) also used the full spectrum version of equation (2.7) (i.e. they included m = 0) and found an increase in ℓ sh at sufficiently small Ekman number.
The mechanical boundary conditions that are employed in the model might play an important role in the scaling behaviour of the convective length scales.No-slip boundary conditions result in the formation of Ekman layers, which tend to significantly damp the large scale flows (Soward 1977;Christensen 2002).One might expect that the strong radial shear generated by the zonal flow limits the size to which convective structures can grow, therefore yielding distinct scaling behaviour in comparison to simulations in which the zonal flow is suppressed.It is also possible that a transition to a regime in which the convective length scale increases occurs in a parameter regime outside of that used in the presented study.However, our data shows no evidence toward such a transition.Another possibility for the observed differences in the length scale could be that the m = 0 component of the flow is influencing the length scale for some of the no-slip papers.While the zonal flow is smaller for no-slip cases, zonal flows can still be present and could influence l sh .It is also worth noting that Gastine & Aurnou (2023) found that heat transport for convection in rotating spherical shells is strongly dependent on latitude.This suggests that other quantities, such as the length scale and Reynolds number, may also depend strongly on latitude.In this case, it may be illuminating to study the Rayleigh number dependence of the length scale as a function of latitude rather than as a globally averaged quantity.
The asymptotic scaling of various terms in the fluctuating momentum and heat equations was investigated.In particular, for the fluctuating momentum equation, we found that the viscous force and the buoyancy force scale as Ek −1 , and that the fluctuating-fluctuating advective term scales approximately as Ek −1 , or perhaps slightly stronger than Ek −1 (with our large scale viscous diffusion nondimensionalisation), whereas the Coriolis force scales as Ek −4/3 .These scalings again agree with asymptotic theory up to a small difference in the advective term.We find that viscosity does not become smaller compared to buoyancy as the reduced Rayleigh number is increased, so we do not observe a trend in which viscosity becomes negligible.We also do not observe a balance between buoyancy and advection for any of our simulations; the advective term is always growing more rapidly with Rayleigh number than the buoyancy force.Indeed, for the largest Ekman number considered for η = 0.35 (Ek = 3 × 10 −4 ), the magnitude of inertia surpasses that of the buoyancy force, though this behaviour is outside the regime of rapidly rotating convection.This finding suggests that, at least for the parameter space observed in this study, there is no CIA balance that occurs.
The force balances that were computed for our spherical simulations are similar to the results found by Guzmán et al. (2021) for rapidly rotating convection with no-slip boundaries in Cartesian coordinates.In particular, Guzmán et al. (2021) found that for large Rayleigh numbers, the viscous force was approximately as large as the buoyancy force and that the inertial term would be larger than buoyancy.Simulations of the asymptotic model for the plane layer also show that the ratio of the viscous force to the buoyancy force approaches unity in the high Rayleigh number regime (Maffei et al. 2021;Oliver et al. 2023).This similarity indicates that the Rayleigh number dependence of the forces may not depend too sensitively on the geometry of the system.Some dynamo simulations also show similar scalings compared to the hydrodynamic cases.Yan & Calkins (2022) carried out numerical simulations of a dynamo in a plane layer and also found that viscosity and buoyancy follow the same Ekman number dependence with no indication of a regime where viscosity becomes small compared to buoyancy.On the other hand, Yadav et al. (2016) carried out numerical simulations of spherical dynamo cases with no-slip boundaries and found that the ratio of viscosity to buoyancy decreased as the Ekman number was decreased at constant Ra/Ra c .Some other papers studied the spectral decomposition of the forces for spherical dynamo cases and also found that viscosity was much smaller than buoyancy at small values of the Ekman number (e.g.Schwaiger et al. 2019Schwaiger et al. , 2021)).These dynamo simulations where viscosity seems to become small compared to buoyancy at small values of the Ekman number could be following a different asymptotic scaling than we observe in the hydrodynamic cases in this paper, although dynamo plane layer simulations seem to follow similar scaling laws as the hydrodynamic model (Yan & Calkins 2022).
Funding.The authors gratefully acknowledge funding from the National Science Foundation through grants EAR-1945270 (JAN, MAC) 2: Summary of the cases with η = 0.7.Nr is the number of radial grid points used in the simulation and lmax is the maximum spherical harmonic degree used in the simulation.The standard deviation in time of given quantities is shown after the '±'.

Figure 2 :
Figure 2: Reynolds number characterising the flow speeds of the fluctuating (convective) velocity versus the reduced Rayleigh number, Ra: (a) the convective Reynolds number Rec; (b) the rescaled convective Reynolds number Rec = Ek 1/3 Rec; (c) the compensated convective Reynolds number Rec Ra −1 ; (d) the compensated convective Reynolds number Rec Ra −3/2 .The filled symbols represent η = 0.35 cases and the hollow symbols represent η = 0.7 cases.

Figure 3 :Figure 4 :
Figure 3: Reynolds number characterising the flow speeds of the mean (zonal) velocity field versus Ra: (a) the zonal Reynolds number Rez; (b) the rescaled zonal Reynolds number Rez = Ek 2/3 Rez.The solid line is a least squares fit ofthe data for the η = 0.35, Ek = 10 −5 cases to a line and is given by Rez = (0.022Ra − 0.19)Ek −2/3 ; the dashed line is a least squares fit of all the η = 0.7 data to a line and is given by Rez = (0.031Ra − 0.016)Ek −2/3 .The symbols are the same as defined in figure 2.

Figure 6 :
Figure 6: Time-and radially-averaged spherical harmonic kinetic energy spectra.Cases with η = 0.35 and similar values of the reduced Rayleigh number are shown.Circles denote the critical azimuthal wavenumber (mc) at the onset of convection for each parameter set and squares denote the degree l derived from the spherical harmonic length scale ℓ ′ sh = π/l.(a) Spherical harmonic spectra.(b) Rescaled spherical harmonic spectra.(c) Spherical harmonic spectra at small l.(d) Rescaled spherical harmonic spectra at small l.

Figure 7 :
Figure 7: Time-and radially-averaged spherical harmonic spectra summed.Select cases from η = 0.35, Ek = 10 −5 are shown.Circles denote the critical azimuthal wavenumber at the onset of convection and squares denote the spherical harmonic degree calculated from ℓ ′ sh , as in figure 6.(a) Spherical harmonic spectra.(b) Spherical harmonic spectra normalized to have an area of one.

Figure 8 :
Figure 8: Convective length scales for the η = 0.35 cases: (a), (b) spherical harmonic length scale; (c), (d) Taylor microscale.Plots (b) and (d) show the asymptotically rescaled length scales.The symbols are the same as defined in figure 2.

Figure 9 :Figure 10 :Figure 11 :
Figure 9: Length scale calculated from the degree l peak where the fluctuating kinetic energy power spectrum peaks: (a) peak length scale; (b) peak length scale versus the Rossby number characterising the convective flow speeds (two Rossby number scalings are shown for reference); (c) peak length scale rescaled by Ek −2/9 ; (d) peak length scale rescaled by Ek −1/6 .

Figure 14 :
Figure 14: Time-averaged volume rms of various terms from the radial component of the fluctuating momentum equation: (a) viscous force; (b) buoyancy force; (c) fluctuating advection term; (d) rescaled viscous force; (e) rescaled buoyancy force; (f) rescaled fluctuating advection term.The symbols are the same as defined in figure 2.

Figure 15 :
Figure 15: Time-averaged volume rms of various terms from the radial component of the fluctuating momentum equation: (a) advection by the zonal flow; (b) time derivative (inertia); (c) the sum of the time derivative and advection by the zonal flow; (d) rescaled advection by the zonal flow; (e) rescaled time derivative; (f) rescaled sum of the time derivative and advection by the zonal flow.The symbols are the same as defined in figure 2.

Figure 17 :
Figure 17: Ratios of forces in the fluctuating radial momentum equation where Fv = ∇ 2 u ′ r , Fa = [u ′ • ∇u ′ ] r , and F b = (Ra/P r)(r/r0)T ′ : (a) ratio of the viscous force to the buoyancy force; (b) ratio of the fluctuating advective term to the buoyancy force.The symbols are the same as defined in figure 2.

Figure 18 :
Figure 18: Time-averaged Nusselt number calculated at the outer boundary: (a) Nusselt number; (b) compensated Nusselt number.The symbols are the same as defined in figure 2.

Figure 19 :Figure 20 :
Figure 19: Viscous dissipation rates for the convective flow and the zonal flow: (a) dissipation rate; (b) rescaled dissipation rate.The filled symbols denote small-scale dissipation rate (ε ′ ) and hollow symbols denote large-scale dissipation rate (ε).