Predicting the Rossby number in convective experiments

The Rossby number is a crucial parameter describing the degree of rotational constraint on the convective dynamics in stars and planets. However, it is not an input to computational models of convection but must be measured ex post facto. Here, we report the discovery of a new quantity, the Predictive Rossby number, which is both tightly correlated with the Rossby number and specified in terms of common inputs to numerical models. The Predictive Rossby number can be specified independent of Rayleigh number, allowing suites of numerical solutions to separate the degree of rotational constraint from the strength of the driving of convection. We examine the scaling of convective transport in terms of the Nusselt number and the degree of turbulence in terms of the Reynolds number of the flow, and we find scaling laws nearly identical to those in nonrotational convection at low Rossby number when the Predictive Rossby number is held constant. Finally, we describe the boundary layers as a function of increasing turbulence at constant Rossby number.


INTRODUCTION
Rotation influences the dynamics of convective flows in stellar and planetary atmospheres. Many studies on the fundamental nature of rotating convection in both laboratory and numerical settings have provided great insight into the properties of convection in both the rapidly rotating regime and the transition to the rotationally unconstrained regime (King et al. 2009;Zhong et al. 2009;Schmitz & Tilgner 2009;King et al. 2012;Julien et al. 2012;King et al. 2013;Ecke & Niemela 2014;Stellmach et al. 2014;Cheng et al. 2015;Gastine et al. 2016) The scaling behavior of heat transport, the nature of convective flow structures, and the importance of boundary layer-bulk interactions in driving dynamics are well known. Yet, we do not know of any simple procedure for predicting the magnitude of vortical flow gradients purely from experimental control parameters, such as bulk rotation rate and thermal input.
In the astrophysical context, many studies of rotating convection have investigated questions inspired by the solar dynamo (Glatzmaier & Gilman 1982;Busse 2002;Brown et al. 2008Brown et al. , 2010Brown et al. , 2011Augustson et al. 2012;Guerrero et al. 2013;Käpylä et al. 2014). Even when these simulations nominally rotate at the solar rate, they frequently produce distinctly different behaviors than the true Sun, such as anti-solar differential rotation profiles (Gastine et al. 2014;Brun et al. 2017). It seems that these differences occur because the simulations produce less rotationally constrained states than the Sun. The influence of rotation results from the local shear gradients, and these are not direct input parameters. Recent simulations predict significant rotational influence in the deep solar interior, which can drastically affect flows throughout the solar convection zone (Featherstone & Hindman 2016;Greer et al. 2016). In the planetary context, the balance between magnetic and rotational forces likely leads to the observed differences between ice giant and gas giant dynamos in our solar system (Soderlund et al. 2015). The work of Aurnou & King (2017) demonstrates the importance of studying a dynamical regime with the proper balance between Lorentz, Coriolis, and inertial forces when modeling astrophysical objects such as planetary dynamos.
In short, simulations must achieve the proper rotational balance if they are to explain the behavior of astrophysical objects. In Boussinesq studies, rotational constraint is often measured by comparing dynamical and thermal boundary layers or deviation in heat transport from the non-rotating state (King et al. 2012;Julien et al. 2012;King et al. 2013). Such measurements are not available for astrophysical objects, where the degree of rotational influence is best assessed by the ratio between nonlinear advection magnitude and the linear Coriolis accelerations. The Rossby number is the standard measure of this ratio, where Ω denotes the bulk rotation vector. Many proxies for the dynamical Rossby number exist that are based solely on input parameters, most notably the convective Rossby number. However, all proxies produce imperfect predictions for the true dynamically relevant quantity.
In this letter, we demonstrate an emperical method of predicting the output Rossby number of convection in a simple stratified system.
In Anders & Brown (2017) (hereafter AB17), we studied non-rotating compressible convection without magnetic fields in polytropic atmospheres. In this work, we extend AB17 to rotationally-influenced, f -plane atmospheres (e.g., Brummell et al. 1996Brummell et al. , 1998Calkins et al. 2015). We determine how the input parameters we studied previously, which controlled the Mach and Reynolds numbers of the evolved flows, couple with the Taylor number (Ta, Julien et al. 1996), which sets the magnitude of the rotational vector.
In section 2, we describe our experiment and paths through parameter space. In section 3, we present the results of our experiments and in section 4 we offer concluding remarks.

EXPERIMENT
We study fully compressible, stratified convection under precisely the same atmospheric model as in AB17, but here we have included rotation. We study polytropic atmospheres with n ρ = 3 density scale heights and a superadiabatic excess of = 10 −4 such that flows are at low Mach number. We study a domain in which the gravity, g = −gẑ, and rotational vector, Ω = Ωẑ, are antiparallel (as in e.g., Julien et al. 1996;Brummell et al. 1996).
We evolve the velocity (u), temperature (T ), and log density (ln ρ) according to the Fully Compressible Navier-Stokes equations in the same form presented in AB17, with the addition of the Coriolis term, 2Ω × u, to the left-hand side of the momentum equation. We impose impenetrable, stress-free, fixed-temperature boundary conditions at the top and bottom of the domain.
We specify the kinematic viscosity (ν), thermal diffusivity (χ), and strength of rotation (Ω) at the top of the domain by choosing the Rayleigh number (Ra), Prandtl number (Pr), and Taylor number (Ta), where L z is the depth of the domain as defined in AB17, ∆S ∝ n ρ is the specific entropy difference between the top and bottom of the atmosphere, and the specific heat at constant pressure is c P = γ/(γ − 1) with γ = 5/3. Throughout this work we set Pr = 1. The Taylor number relates to the often-quoted Ekman number by the equality Ek ≡ Ta −1/2 . Due to stratification, Ra and Ta both grow with depth as (Ra,Ta) ∝ ρ 2 (see AB17). We nondimensionalize our atmospheres at the top of the domain, and so all values of Ra and Ta quoted in this work are the minimal value of Ra and Ta in the domain at z = L z . For direct comparison to Boussinesq studies, past work has found that the value of Ra at the atmospheric midplane (z = L z /2) varies minimally with increasing stratification (Unno et al. 1960). For the atmospheres presented in this work, midplane Ra and Ta values are larger than reported top-of-atmosphere values by a factor of ∼70, and values at the bottom of the atmosphere are larger by ∼400.
When Ta is large, the wavenumber of convective onset increases according to k crit ∝ Ta 1/6 (Chandrasekhar 1961;Calkins et al. 2015). We study horizontallyperiodic, 3D Cartesian domains with extents of x, y = [0, 4(2π/k crit )] and z = [0, L z ]. At large values of Ta, these domains are tall and skinny, as in Stellmach et al. (2014). We evolve our simulations using the Dedalus 1 pseudospectral framework, and our numerical methods are identical to those presented in AB17. The supplemental materials of this paper include a .tar file which contains the code used to perform the simulations in this work, and this tarball is also published online in a Zenodo repository .
The critical value of Ra at which rapidly rotating convection onsets also depends on Ta (see the black line in figure 1a), roughly according to Ra crit ∼ Ta 2/3 (Chandrasekhar 1961; Calkins et al. 2015). Even taking account of linear theory, the dependence of the evolved nonlinear fluid flows on the input parameters makes predicting the rotational constraint very challenging. We will explore three paths through Ra-Ta space: (Ro p ) 2 Pr 1/2 Ta 3/4 (III). (3) Paths on constraint I are at constant supercriticality, S ≡ Ra/Ra crit (Ta) (blue dash-dot line in figure 1a). Paths on constraint II (green dashed line in figure 1a) are at a constant value of the classic convective Rossby number, which has provided (e.g., Julien et al. 1996;Brummell et al. 1996) a common proxy for the degree of rotational constraint. This parameter measures the importance of buoyancy relative to rotation without involving dissipation. Paths on constraint III (e.g., orange solid lines in figure 1a) set constant a ratio which we call the "Predictive Rossby number," Unlike paths through parameter space which hold Ro c constant, paths with constant Ro p feel changes in diffusivities but not the depth of the domain. To our knowledge, these paths have not been reported in the literature, although the importance of Ra/Ta 3/4 = Ra Ek 3/2 has been independently found by King et al. (2012) using a boundary layer analysis. We compare our results to their theory in Section 4. In this work, we primarily study three values of Ro p . These values are shown in Fig. 1a and Table 1. Table 1 lists the values of (Ra crit , Ta crit ) for each value of Ro p , and also the maximum value of (Ra, Ta) studied in this work for each path. We additionally walked two pathways at constant supercriticality (constraint I, S = {2, 3}) and three pathways at constant convective Rossby number (constraint II, Ro c = {1, 0.3, 0.1}). Full details on all cases are provided in Appendix A and the supplemental materials. 0.60 (10 4.88 , 10 7.10 ) (10 9.09 , 10 12.72 ) 0.96 (10 2.44 , 10 3.30 ) (10 8.58 , 10 11.49 ) 1.58 (10 1.39 , 10 1.33 ) (10 7.14 , 10 8.99 ) Note-Values of the critical Ra and Ta for each Rop track are reported, as well as the maximal values of Ra and Ta studied on each track. All values reported are for the top of the atmosphere. A fuller set of simulations are reported in Table 2 with midplane Ra and Ta values as well.

RESULTS
In our stratified domains, for Ta ≥ 10 5 , a bestfit to results from a linear stability analysis provides Ra crit (Ta) = 1.459Ta 2/3 and k crit (Ta) = 0.414Ta 1/6 for direct onset of convection. In figure 1a, the value of Ra crit (Ta) is shown. Sample paths for each criterion in equation 3 through this parameter space are also shown. In this work, we often find it instructive to use one critical Ra for an entire Ro p path. This Ra crit is determined by the intersection of the onset curve and Ro p path (indicated by the orange circles in figure 1a, and quoted in Table 1). In the high Ta regime, we find that Ra crit = 18.5Ro −16 p . In figure 1b, we display the evolution of Ro with increasing Ra along various paths through parameter space. We find that Ro increases on constant Ro c paths, decreases on constant S paths, and remains roughly constant along constant Ro p paths. In figure 1c, the value of Ro is shown simultaneously as a function of Ro p and Ro c for all experiments conducted in this study. We find a general power-law of the form Ro = CRo α c Ro β p . In the rotationally-dominated regime where Ro < 0.2 and Re ⊥ > 5 (see Eqn. 6), we find α = −0.02, and Ro can be said to be a function of Ro p alone. Under this assumption, we report a scaling of Ro = (0.148 ± 0.003)Ro 3.34±0.07 p . In the less rotationally dominated regime of Ro > 0.2 and Re ⊥ > 5, we find {C, α, β} = {0.2, −0.19, 1.5}.
In figure 2, sample snapshots of the evolved entropy field in the x-y plane near the top and at the middle of the domain are shown. In the left column, flows are at Ro ∼ 1 and resemble the classic granular structure of nonrotating convection (see e.g., figure 2 in AB17), where strong narrow downflow lanes punctuate broad upwellings. The narrow downflows at the top organize themselves into intense coherent structures at the mid- plane, and at the midplane the downflows have much stronger entropy fluctuations than the broad and slower upflows.
Ro decreases from left to right into the rotationally constrained regime. As Ro decreases, the narrow downflow lanes begin to disappear and the flows at midplane become more symmetric. In the rotationally constrained regime (Ro ∼ 0.03), the convective structures are distinctly different. Here we observe dynamically persistent, warm upflow columns surrounded by bulk weak downflow regions. At the midplane, the upflow columns have substantially higher entropy perturbations than the surrounding weak downflows which sheathe them, and the locations of the columns are tightly correlated with their positions at the top of the domain. These quasitwo-dimensional dynamics are similar to those seen in rapidly rotating Rayleigh-Bénard convection (e.g., Stellmach et al. 2014). The select cases displayed in figure 2 each have an evolved volume-averaged Re ⊥ ≈ 32 (defined below in equation 6).
We measure the Nusselt number (Nu), which quantifies heat transport in a convective solution, as defined in AB17.
In figure 3a, we plot Nu as a function of Ra/Ra crit at fixed Ro p . We find that Nu ∝ {Ra 0.29±0.01 , Ra 0.29±0.01 , Ra 0.24 } for Ro p = {0.6, 0.957, 1.58}. In the regime of Ro 0.1, these scaling laws are indistinguishable from a classic Ra 2/7 power law scaling, which is observed in nonrotating Rayleigh-Bénard and stratified convection (Ahlers et al. 2009, AB17). Our results seem consistent with the stress-free, rotating Rayleigh-Bénard convection results of Schmitz & Tilgner (2009), whose re-arranged Eqn. 7 returns a best-fit of Nu ∝ Ra 0.26 at fixed Ro p . Their work primarily spans the transition regime between rotationally constrained and unconstrained convection, and so it is perhaps not surprising that their power law is a blend of our rotationally-constrained Ra 2/7 power law and the fairly rotationally unconstrained Ra 0.24 at Ro p = 1.58.
Flows are distinctly different parallel to and perpendicular from the rotation vector, which aligns with gravity and stratification. We measure two forms of the RMS Reynolds number, where the length scale in Re ⊥ is the wavelength of convective onset, and is related to the horizontal extent of our domain (see section 2). From our work in AB17, we expect the RMS velocity to scale as |u| ∝ √ ∆S. By definition, ν ∝ Ra/(Pr ∆S), and L z is a constant set by the stratification while k crit ∝ Ta 1/6 . Along paths of constant Ro p , we thus expect Re ∝ Ra 1/2 and Re ⊥ ∝ Ra 5/18 when Pr is held constant.   In figure 3b, we plot Re and Re ⊥ as a function of Ra/Ra crit at fixed Ro p . We find that Re ∝ {Ra 0.44±0.01 , Ra 0.45±0.01 , Ra 0.44 } and Re ⊥ ∝ {Ra 0.22±0.01 , Ra 0.23±0.01 , Ra 0.21 } for Ro p = {0.6, 0.957, 1.58}. These scalings are similar to but slightly weaker than our predictions in all cases.
However, the scaling of Re ∝ Ra 0.45 , is once again a power law observed frequently in nonrotating convection (Ahlers et al. 2009, AB17). We also observe that Re ⊥ collapses for each Ro p track, while Re experiences an offset to larger values as Ro p shrinks. The offset in Re is unsurprising, because more rotationally constrained flows result in smaller boundary layers relative to the vertical extent of our stratified domain. The horizontal extent of our domain scales with the strength of rotation, and so regardless of Ro p , flows perpendicular to the rotational and buoyant direction are comparably turbulent at the same Ra/Ra crit . We find Re ⊥ and Re are, respectively, good proxies for the horizontal and perpendicular resolution required to resolve an experiment.  ior from low Ra (yellow) to high Ra (purple) is denoted by the color of the profile. As Ra increases at a constant value of Ro p , both the thermal (σ s ) and dynamical (Ro) boundary layers become thinner. We measure the thickness of the thermal boundary layer (δ s ) at the top of the domain by finding the location of the first maxima of σ s away from the boundary. We measure the thickness of the Ro boundary layer (δ Ro ) in the same manner. In figure 4e, we plot δ Ro /δ s , the ratio of the sizes of these two boundary layers. As anticipated, the dynamical boundary layer (δ Ro ) becomes relatively thinner with respect to the thermal boundary layer (δ s ) as Ro and Ro p decrease. However, the precise scaling of this boundary layer ratio with Ro p and Ra is unclear, and we cannot immediately compare these ratios to similar measures from the Rayleigh-Bénard convection literature, such as Fig. 5 of King et al. (2013). They measure the dynamical boundary layer thickness as the peak location of the horizontal velocities, but our horizontal velocities are subject to stress-free boundary conditions, and we find that the maxima of horizontal velocities occur precisely at the boundaries. In figure 4f, we plot δ s in units of the density scale height at the top of the atmosphere, and we plot vertical lines when this crosses 1. We find no systematic change in behavior when δ s is smaller than the local density scale height.

DISCUSSION
In this letter, we studied low-Mach-number, stratified, compressible convection under the influence of rotation. We examined three paths through Ra-Ta space, and showed that the newly-defined Predictive Rossby number, Ro p = Ra/(Pr 1/2 Ta 3/4 ), determines the value of the evolved Rossby number.
Shockingly, along these constant Ro p pathways, particularly when Ro 0.1, we find Nu ∝ Ra 2/7 and Re ∝ Ra 0.45 . These scalings are indistinguishable from the scalings of Re and Nu with Ra in non-rotating Boussinesq convection (Ahlers et al. 2009). Julien et al. (2012) theorized that in the rapidly rotating asymptotic limit, (Nu − 1) ∝ (Ra 3/2 /Ta) = (Ra/Ra crit (Ta)) 3/2 . Thus, at fixed Ta, a very sharp Ra 3/2 scaling law is expected. At a fixed Ta = 10 14 , Stellmach et al. (2014) found that the Ra 3/2 scaling described the results of stress-free DNS in Boussinesq cylinders very well. Gastine et al. (2016) studied Boussinesq convection in spherical shells with no-slip boundaries, and also found good agreement with the theory of Julien et al. (2012) for various Ra at Ta ≥ 10 10 .
Here, when we run simulations at fixed Ro p , the value of Ta is coupled to the value of Ra, and both increase simultaneously. Recasting the scaling of Julien et al. (2012) into this perspective, we find (Nu − 1) ∝ Ra 3/2 /Ta = Ro 8/3 p Ra 1/6 ∝ (Ra/Ra crit ) 1/6 , where in this final result we use the Ra crit value of the whole Ro p path, such as those specified in Table 1. This Ra 1/6 scaling is much weaker than the Ra 2/7 law we find here. We leave it to future work to explain this discrepancy between Boussinesq theory and our observed Nu vs. Ra scaling.
In this work, we experimentally arrived at the Ra/Ta 3/4 = Ra Ek 3/2 scaling in Ro p , but this relationship was independently discovered by King et al. (2012). Arguing that the thermal boundary layers should scale as δ S ∝ Ra −1/3 and rotational Ekman boundary layers should scale as δ Ro ∝ Ta −1/4 = Ek 1/2 , they expect these boundary layers to be equal in size when Ra/Ta 3/4 ∼ 1. They demonstrate that when 2 Ra/Ta 3/4 20 flows are in the transitional regime, and for Ra/Ta 3/4 2, flows are rotationally constrained. We remind the reader that Boussinesq values of Ra and Ta are not the same as their values in our stratified domains here, as diffusivities change with depth (see section 2). Taking into account this change with depth, our simulations fall in King et al. (2012)'s rotationally constrained (Ro p = 0.6) and near-constrained transitional regime (Ro p = {0.957, 1.58}). The measured values of Ro in Fig. 1b and the observed dynamics in Fig. 2 agree with this interpretation.
We note briefly that the scaling Ra ∝ Ta 3/4 is very similar to another theorized boundary between fully rotationally constrained convection and partially constrained convection predicted in Boussinesq theory, of Ra ∝ Ta 4/5 (Julien et al. 2012;Gastine et al. 2016). This Ta 4/5 scaling also arises through arguments of geostrophic balance in the boundary layers, and is a steeper scaling than the Ta 3/4 scaling present in Ro p . This suggests that at sufficiently low Ro p , a suite of simulations across many orders of magnitude of Ra will not only have the same volume-averaged value of Ro (as in Fig. 1b), but will also maintain proper force balances within the boundary layers.
Our results suggest that by choosing the desired value of Ro p , experimenters can select the degree of rotational constraint present in their simulations. We find that Ro ∝ Ro 3.34±0.07 p , which is within 2σ of the estimate in King et al. (2013), who although defining Ro very differently from our vorticity-based definition here, find Ro ∝ Ro 3.84±0.28 p . We note briefly that they claim that the value of Ro is strongly dependent upon the Prandtl number studied, and that low Ro can be achieved at high Pr without achieving a rotationally constrained flow. We studied only Pr = 1 here, and leave it to future work to determine if the scaling of Ro p ∝ Pr −1/4 is the correct scaling to predict the evolved Rossby number.
Despite the added complexity of stratification and despite our using stress-free rather than no-slip boundaries, the boundary layer scaling arguments put forth in King et al. (2012) seem to hold up in our systems. This is reminiscent of what we found in AB17, in which convection in stratified domains, regardless of Mach number, produced boundary-layer dominated scaling laws of Nu that were nearly identical to the scaling laws found in Boussinesq Rayleigh-Bénard convection.
We close by noting that once Ro p is chosen such that a convective system has the same Rossby number as an astrophysical object of choice, it is straightforward to increase the turbulent nature of simulations by increasing Ra, just as in the non-rotating case. Although all the results reported here are for a Cartesian geometry with antiparallel gravity and rotation, preliminary 3D spherical simulations suggest that Ro p also specifies Ro in more complex geometries (Brown et al. 2019 in prep).
We thank Jon Aurnou and the anonymous referee who independently directed us to the work of King et al. (2012) Table 2. The simulation at minimum (Ra, Ta) and maximum (Ra, Ta) for each of the Ro p , Ro c , and S paths in Fig. 1b are shown. This information for the displayed simulations and all other simulations in this work is included as a .csv file in the supplemental materials and is published online in a Zenodo repository ).