Planetary Core‐Style Rotating Convective Flows in Paraboloidal Laboratory Experiments

Turbulent convection in a planet's outer core is simulated here using a thermally‐driven, free surface paraboloidal laboratory annulus. We show that the rapidly rotating convection dynamics in free‐surface paraboloidal annuli are similar those in planetary spherical shell geometries. Three experimental cases are carried out, respectively, at 35 revolutions per minute (rpm), 50 and 60 rpm. Thermal Rossby waves are detected in full disk thermographic images of the fluid's free surface. Ultrasonic flow velocity measurements reveal the presence of multiple azimuthal (zonal) jets, with successively more jets forming in higher rotation rate cases. The jets' cylindrical radial extent is well approximated by the Rhines scale. Over time, the zonal jets migrate to larger radial position with migration rates in good agreement with prior theoretical estimates. Our results suggest that planetary core rotating convection will be comprised of flow structures found in other turbulent geophysical fluid dynamical systems: convective turbulence dominates the small‐scale flow field, and also act to flux energy into larger‐scale, slowly evolving zonal flow structures. How the ambient magnetic fields in planetary core settings affect such turbulent flows remains an open question.

Core convection experiments in polar and spherical geometries have all, to date, been carried out in closed containers. The global no-slip conditions imparted by these containers generates flows that are massively damped by viscous effects, relative to those in planetary cores. (In contrast, a number of annular ocean-atmosphere experiments have been carried out using free upper fluid surfaces [e.g., Condie & Rhines, 1994;Read, Jacoby, et al., 2015;Smith et al., 2014].) Even for the most extreme core convection experiments Gillet et al., 2007), viscous forces are likely to be at least six orders of magnitude greater than those in planetary core settings. These viscous effects act to damp out turbulent interactions, limiting possible cascades and mean zonal flows, and removing dynamical modes that could dominate in less dissipative environments (Aurnou & Heimpel, 2004;Heimpel et al., 2022;Jones & Kuzanyan, 2009;Plumley et al., 2016;Stellmach et al., 2014).
Here we propose a novel experimental configuration to model low latitude hydrodynamic outer core convection, as exists outside the tangent cylinder ( Figure 1b). The set-up is that of a cylindrical fluid annulus with a free upper fluid surface that strongly deforms into the shape of an axisymmetric paraboloid in rapidly rotating cases. This paraboloidal annular geometry can adequately mimic the zeroth-order topographic features of a spherical shell, as shown in Section 3.3. Convection is modeled using water as the working fluid, which is cooled at the inner annular boundary and with no thermal flux across the outer annular boundary. Thermally-driven convection in water with a no flux outer boundary condition is, somewhat non-intuitively, a good analog for compositional convection in Earth's core (Buffett et al., 2000;Cardin & Olson, 1992), the dominant driver of present-day core turbulence (Driscoll & Du, 2019). Our device, shown schematically in Figure 1b and in detail in Figure 3, has been dubbed the Coreaboloid and is described further in Section 2.
The Coreaboloid is a natural extension of the differentially heated rotating annulus. Originally devised by R. Hide to simulate flows in planetary cores (Ghil et al., 2010), the classical Hide annulus has been used primarily to study atmospheric and ocean dynamics (Bastin & Read, 1997;Read, 2001; von Larcher & Egbers, 2005; Wordsworth Figure 2 compares polar images of a numerical model of rotating convection in Earth's core against a laboratory Coreaboloidal model in order to demonstrate that convection in the Coreaboloid can serve as a useful proxy for planetary core convection. Figure 2a shows a snapshot of the equatorial temperature field taken from the numerical modeling study of Gastine et al. (2016). Figure 2b shows a top view image of dye injected into the convecting Figure 1. (a) Sideview schematic of Earth's core. The core-mantle boundary bounds the molten metal, fluid outer core, which bounds the solid inner core. The so-called tangent cylinder is the imaginary axial cylinder that circumscribes the inner core equator and separates the low latitude core fluid outside the tangent cylinder (OTC) from the high latitude fluid regions inside the tangent cylinder. Polar core flow is often experimentally simulated via cylindrical fluid domains (gray cylinder). (b) Sideview schematic of the paraboloidal Coreaboloid device, by which we experimentally simulate OTC rotating convective turbulence.  Gastine et al. (2016) (Ra ⊥ /Fr = 5.5 × 10 9 , Ek = 2.1 × 10 −7 , Pr = 1 and R i /R o = 0.35). (b) Laboratory polar view of a paraboloidal convection experiment with water as the working fluid. Blue-green (red) dye has been injected in this 55 rotating per minute (rpm) case near the inner (outer) cylinder boundary (Ra ⊥ ≃ 2 × 10 10 , E ≃ 3.9 × 10 −7 ; Fr = 1.26; Pr ≃ 4, and R i /R o = 0.27). The inset shows a thermographic image of the free surface temperature field, albeit from a 45 rpm case. In all images, warmer (cooler) fluid is colored red (blue). LONNER ET AL. Coreaboloid. The red dye near the outer boundary is visually oversaturated, but the aquamarine dye in the inner half of the fluid annulus marks the fluid motions sufficiently well. The Figure 2b inset shows a snapshot of the free surface temperature, imaged from above via infrared (IR) thermography. The Coreaboloidal dye tracer and surface temperature images are both comprised of flows that are similar in zeroth-order scale and structure to those in the three-dimensional spherical shell core flow model. This shows, albeit qualitatively, that the Coreaboloid provides a good experimental system to characterize core-style OTC dynamics.
Figure 2 also shows that the temperature gradient in the Coreaboloid is reversed relative to that of a self-gravitating sphere, similar to many previous laboratory core convection experiments (e.g., Busse & Carrigan, 1974, 1976Cardin & Olson, 1992Carrigan & Busse, 1983;Sumita & Olson, 2000). In a self-gravitating sphere, gravity points radially inwards ( Figure 1a) and varies as = − (4∕3) , where is the mean material density and G is the gravitational constant. In the centrifugally-driven case, effective gravity points radially outwards, varying as = Ω 2 −̂≈ Ω 2 in the limit in which centrifugal acceleration greatly exceeds lab gravity g. To generate similar buoyancy forcing as in planets where the inner boundary fluid tends to be warmer (red) and less dense than the cooler (blue) outer fluid, the temperature gradient must be reversed in centrifugally-driven laboratory systems with the colder fluid (blue) along the inner boundary and the warmer fluid (red) along the outer boundary. Thus, the signs of the radial temperature gradients must be flipped between Figures 2a and 2b in order for them to have similar buoyancy-driven convection dynamics. Where the centers of planets must be relatively warm in order to drive thermal convection (Figure 2a), the center of a centrifugally-driven system must be made relatively cool (Figure 2b) to make the fluid unstable to convective motions. This paper introduces the Coreaboloid and demonstrates it to be a flexible new tool for simulating remote, extreme convection-driven processes in planetary cores, and other geophysical fluid systems as well. The Coreaboloidal set up takes advantage of the flexibility of the Hide Annulus and minimizes viscous drag via the free upper fluid surface. Further, the free upper surface allows for detailed global-scale thermographic measurements, comparable in spatial resolution to numerical models. The paraboloidal annular device is described in Section 2. Section 3 10.1029/2022JE007356 5 of 26 makes comparisons between paraboloids and spherical shells and then makes theoretical estimates for the flows that will arise in our device. Section 4 presents the results of three experimental cases, carried out at 35, 50, and 60 revolutions per minute (rpm). Finally, in Section 5, our laboratory results are extrapolated to core conditions, and we discuss future Coreaboloidal experiments, modifications and improvements to be made. Figure 3a shows a three-dimensional rendering of the Coreaboloid laboratory apparatus, with the IR thermographic image from the 50 rpm case overlain atop the paraboloidally-deflected free upper surface that matches that rotation rate. The rotary table is constructed of T6061 structural aluminum. The 96.5 cm diameter, double-decker rotary platform is attached to a 40.6 cm diameter thrust bearing, and is belt driven by a one horsepower AC motor. The experimental tank is placed on the upper deck, whereas all the electronics are situated on the lower deck. The upper scaffolding is constructed of MiniTec extruded aluminum beams and houses our thermographic camera and a polished aluminum plate. The plate acts as an IR reflector that creates a sufficiently long optical path such that it is possible to image the working fluid's entire free surface. Figure 3b presents a schematized top view of the experimental tank. The acrylic cylindrical sidewall is 1.27 cm thick with an inner radius of R o = 37.25 cm, and is 50 cm tall. The inner cylinder is made of 0.635 cm thick stainless steel and has an outer radius of R i = 10.2 cm and is 30 cm tall. Thus, the cylindrical annular gap has a width of R o − R i = 27.05 cm. The tank bottom is made of 1.27 cm thick polycarbonate. Two thru holes are drilled through the sidewall to accommodate two counter-facing ultrasonic Doppler velocimetry (UDV) probes at a height of 7 cm above the floor of the tank. The chord connecting these probe holes is L = 70.5 cm in length and passes within s off = 2.0 cm of R i . Figure 3c shows a snapshot of the Coreaboloid spinning at roughly 40 rpm. The centrifugally-deflected free surface of the working fluid, water, equilibrates into the shape of a paraboloid (Basta et al., 2000;Cushman-Roisin & Beckers, 2011). The volume of water is fixed here such that the non-rotating depth of the fluid layer is maintained at H NR = 17.08 cm.

The Coreaboloid Laboratory Device
The system's diagnostic capabilities include point thermometry, free surface thermography and UDV. Temperature point measurements are acquired via an Arduino Mega-XBEE 10-bit data acquisition system that connects to a radial chain of seven evenly-spaced Ametherm PANE-103395 thermistors that stretches along the tank bottom, leading to 1 Hz temperature readings with 25 mK precision. The ends of the thermistor chain, on the cylindrical boundaries at R i and R o , are used to calculate the time-averaged horizontal temperature difference across the annular fluid gap, ΔT ⊥ . This quantity is calculated by taking the average of the temperature difference at the surface (IR data) between inner and outer cylinder and the temperature difference at the bases of the inner and outer cylinders (thermistors). In addition, the five other thermistors are used in estimating the characteristic vertical temperature difference across the fluid layer, ΔT ‖ , which is the difference between top and bottom temperatures measured in the middle of the fluid gap, at γ = (s − R i )/(R o − R i ) = 0.5, and averaged in time from t = 1,800-2,400 s. Here, the ⊥ subscript denotes the s-direction perpendicular to the rotation vector and the ‖ subscript denotes the axial z-direction parallel to the rotation vector.
The free surface temperature field is mapped in space and time using an InfraTec ImageIR 8300 thermographic imaging system. The ImageIR's CCD has a spatial resolution of 640 × 512, such that each pixel corresponds to an area of approximately 1.25 mm 2 on the fluid's surface. In our temperature range, the ImageIR has a precision of 0.02 K. Calibration cases were made at each rotation rate with the water held at room temperature and with no ice in the aluminum cylinder. These isothermal cases were used to best fit the IR temperature readings, T IR (i, j), at each individual pixel (i, j) to the mean temperature measured by the thermistors, , using the empirical formula Here a(i, j) is a coefficient matrix which does not greatly deviate from unity in any of our experiments and b(i, j) is a matrix of offset temperatures with maximum values never exceeding 1 K. Expression (1) was used in processing the thermographic free surface temperature data used throughout this study, but that we must stress does not qualitatively change the temperature field or any of the following temperature field analysis.

of 26
Line of sight velocities, u DOP = (u s , u ϕ , 0), are measured along the horizontal chord of length L using a DOP 3010 UDV system from Signal Processing SA. The UDV multiplexes between two transducers placed opposite each other along the chord L shown in Figure 3b. The 16-bit data yields velocity measurements precise to 1μm/s and with accuracy within 5%. The UDV data from x = 0 to L/2 is from the Doppler probe located at x = 0 and the data x = L/2 to L is from the Doppler probe located at the x = L. We make the assumption that the radial component of the convective flow is time varying with a zero-valued temporal mean, whilst the azimuthal (zonal) component can have a non-zero time-mean. We thus project the u DOP data purely into the azimuthal, ϕ-direction by first calculating the local angle, θ DOP , between the Doppler beam and ̂ : where x is the line-of-sight distance from the transducer and s off = 2 cm is the offset between the inner cylinder and the midpoint of the Doppler chord, x = L/2, as shown in Figure 3b. The azimuthal velocity estimate is made by projecting the Doppler line-of-sight velocity onto ̂ : Averaged over time, we use Equation 3 to estimate the quasi-steady mean zonal flow (cf. Vogt et al., 2021). In future Coreaboloid experiments, a radially-oriented Doppler probe will be employed to test and quantify the assumption made here that the radial mean flow is negligible.
The rotary platform is slightly imbalanced due mainly to the off-axis mass of the ImageIR camera. This leads to the generation of a free surface wave with angular frequency matching the table's frequency of rotation. Thus, we bandpass filter the UDV data at the table's rotation frequency in order to remove this wave from the UDV data. Further, we acquire IR images at the frequency of rotation so that this wave is never sampled thermographically. Given the low Rossby numbers in our experiments (Table 2), the flows we are interrogating likely exist and evolve on sub-inertial time scales that should not be strongly affected by the faster inertial oscillatory flows that are being forced at the table's period of rotation.
Each of the three experimental cases presented here were made following the same experimental sequence, which represent the first three production runs made on the Coreaboloid. The cylindrical annulus would be filled with water to a resting (non-rotating) fluid layer depth of H NR = 17.08 cm. A heating element and mixing pump were then put into the fluid until it reached a temperature of 50° ± 1°C. (Since the outer boundary is thermally passive, the heating of the working fluid allowed us to raise the effective buoyancy forces in our experiments by roughly an order of magnitude, greatly increasing the vigor of the convective flows.) The heating element and pump were then removed from the working fluid, and a cylindrical block of ice was placed within the inner stainless steel cylinder, along with some water to thermally couple the ice block to the inner cylinder. The rotary platform and upper scaffolding were then covered in a heavy burlap sheet in order to minimize air drag on the water's free surface. The system was then spun at angular rotation velocity Ω. The fluid was given 20 min to spin-up, corresponding to approximately 12 linear spin-up times (Greenspan & Howard, 1963;Warn-Varnas et al., 1978). Thermistor, IR thermographic and UDV data was then acquired for the next 40 min. The mean temperature of the fluid layer monotonically decreased over this 40 min window due to cooling via the cold inner cylinder and via heat loss through the free upper fluid surface. Thus, the mean temperature of the top and bottom fluid surfaces varied over the data acquisition windows from 46.7°C to 39.2°C in the 35 rpm case; from 47.1°C to 42.1°C in the 50 rpm case; and from 46.7°C to 41.5°C in the 60 rpm case. However, the UDV flow velocities, shown in the results section, show that the kinetic energy in the zonal flow field was quasi-stationary over the data acquisition window. Thus, even though the mean temperature of the working fluid decreases by approximately 6°C through each experiment's data acquisition phase, we argue here that, kinematically, the data corresponds to that of an equilibrated convective state since the measured velocities in the jet cores do not significantly change.
Dimensional and nondimensional experimental parameter values are given, respectively, in Tables 1 and 2. The thermographic and UDV data sets are openly available for download (Lonner et al., 2022).

Governing Equations
Momentum conservation for a Boussinesq, rotating fluid is given by the Navier-Stokes equation: where u is the solenoidal fluid velocity vector measured with respect to the rotating frame, P is the pressure, ρ is the fluid density, =̂ is laboratory gravity, = Ω̂ is the constant angular velocity vector, and ν is the kinematic viscosity. The left two terms on the left hand side of Equation 4 describe the fluid inertia, and the remaining term is the Coriolis acceleration. On the right hand side of Equation 4, we have, from left to right, the acceleration due to pressure gradients, gravitational acceleration, centrifugal acceleration and viscous drag.
Density variations exist only in the effective gravity terms in Boussinesq fluids (Tritton, 2012) and are related directly to variations in the fluid temperature by where α is the thermal expansion coefficient of the fluid and T is the temperature relative to the mean T o . Hydrostatic balance is recovered when u = 0 Applying the Boussinesq approximation and subtracting Equation 6 from Equation 4 gives where ∕ ≡ ( ∕ + ⋅ ∇) is the material derivative, p is the dynamic pressure, and the second and third terms on the right hand side are the gravitational buoyancy and centrifugal buoyancy that drive convection in this system. The evolution of the temperature field is expressed via where κ is fluid's thermal diffusivity.

The vorticity equation arising from the curl of Equation 7 is
where ω = ∇ × u is the fluid vorticity. The terms on the right hand side of Equation 9 represent vorticity generation via, respectively, vortex stretching, baroclinic torques due to gravitational buoyancy, baroclinic torques due to centrifugal buoyancy, and viscous torques.
The ̂component of Equation 9 is the thermal wind balance Note. Fluid properties are calculated using the mean fluid temperature , averaged over the top and bottom surfaces of the fluid at the end of each data acquisition window. The Ekman drag time scale is = ℎ∕ √ Ω . Both τ Ek and the internal deformation radius δ def are calculated here using h at mid-gap (γ = 1/2). which holds for slow, quasi-steady, inviscid conditions (e.g., Vallis, 2006). In Equation 10, azimuthal baroclinic torques due to both laboratory gravitational buoyancy and centrifugal buoyancy are balanced by vortex stretching occurring via axial shearing of the azimuthal velocity field.

The axial, -component of Equation 9 is
where =̂⋅ and =̂⋅ . In rapidly rotating systems, the relative vorticity is smaller than the system's background, or planetary, vorticity, |ω/2Ω| ≪ 1, resulting in a quasi-two-dimensional ζ = ζ(s, ϕ) field at leading order and a simplified axial vorticity equation (e.g., Aubert et al., 2003;Calkins, Noir, et al., 2012), Vertically integrating the expression above from z = 0 to h yields where ( ) is the axially-averaged temperature value. In the Coreaboloid, there are two sources of axial stretching, u z (h) − u z (0). The first is due to fluid advection through topography. Since the non-planar topography in the Coreaboloidal is due only to the axisymmetric paraboloidal upper fluid surface h(s), this stretching effect can be expressed as where =̂⋅ . The second source of vertical stretching is due to Ekman pumping (EP) along no-slip, non-vertical boundaries. In the Coreaboloid, EP occurs only along the bottom boundary at h = 0 and is described via is the Ekman boundary layer thickness (Stellmach et al., 2014;Tritton, 2012). Substituting Equations 14 and 15 into Equation 13 yields the quasi-geostrophic vorticity equation Local vorticity is generated in Equation 16 by baroclinic torques that arise when centrifugal buoyancy forces act on azimuthally-varying density perturbations, Ω 2 . Such flows are, in fact, the centrifugally-driven convection (e.g., Busse & Or, 1986;Calkins et al., 2013) describes the fractional topographic stretching of axial vortices upon being displaced in radius (e.g., Lemasquerier et al., 2021;Vallis, 2006). The variation of fluid layer height stretches the planetary vorticity, 2Ω, and generates local vorticity ζ when a fluid parcel is advected in radius. This is the so-called topographic β effect. Lastly, the Ekman damping time gives the characteristic time scale over which boundary friction causes the local vorticity ζ to decay away (e.g., Brito et al., 2004;Burmann & Noir, 2018). Calculated at the center of the fluid gap, the τ Ek values are 102, 81, and 68 s for the 35, 50, and 60 rpm cases, respectively.

Non-Dimensionalization
We nondimensionalize the above equations using R o as the characteristic length, the radial temperature difference across the annulus ΔT ⊥ (where the ⊥ subscript denotes the cylindrical direction perpendicular to the rotation axis), and centrifugal free fall ⟂ ∼ √ Δ ⟂Ω 2 2 for velocity (e.g., . Time is then scaled The rescaled equations, where all variables are now non-dimensionalized, then take on the following forms.
The nondimensional Navier-Stokes equation and thermal energy equations are, respectively, The ratio of the centrifugal to gravitational buoyancy is characterized in Equation 19 by the Froude number The centrifugal Rossby number estimates the ratio of centrifugal buoyancy-driven inertial and Coriolis accelerations (e.g., Aurnou et al., 2020;, where the centrifugal Rayleigh number characterizes the buoyancy forcing, the Ekman number estimates the ratio of the viscous and Coriolis accelerations, and the Prandtl number is the ratio of the fluid's kinematic and thermal diffusivities,

= ∕
For our experiments, Ro ⊥ ≲ 0.03 (Table 2). This value provides an upper bounding Rossby number estimate since the centrifugal free fall velocity is used in constructing Ro ⊥ . Centrifugal free fall is accurate when the top and bottom bounding surfaces are planar (Hu et al., 2022;Jiang et al., 2020). When either of these surfaces are nonplanar and Ro ⊥ ≪ 1, topographic vortex stretching will act to limit the centrifugal velocity to smaller thermal wind speeds (e.g., Aubert et al., 2001;Aurnou et al., 2020).

10.1029/2022JE007356
10 of 26 The nondimensional form of the vorticity Equation 9 is The scaled thermal wind balance is and the quasigeostrophic vorticity equation is recast as Given that Ek ≪ Ro ⊥ ≪ 1 in our experiments (as well as in planets), expression (21) requires the leading order flows to be weakly varying in (e.g., Calkins, 2018). The axial vorticity dynamics will tend to be dominated by topographic β effects in the first term on the right hand side of Equation 23. Topographic β effects are indeed crucial to rapidly rotating dynamics in deep fluid shells (e.g., Bire et al., 2022;Calkins, Noir, et al., 2012;Gastine et al., 2014;Heimpel & Aurnou, 2007), whereas boundary frictional effects are argued to act weakly on large-scale geophysical flows but may still be non-negligible on local-scale convective motions Stellmach et al., 2014).
Following this section, we will use only dimensional variables, but the experimental cases will be described using both dimensional and nondimensional input and output parameters.

Fluid Layer Heights
In a spherical shell of radius R o , the axial height, h sph , varies as a function of cylindrical radius s as in the region outside the tangent cylinder (Heimpel & Aurnou, 2007). In a rotating cylindrical annulus with inner and outer radii R i and R o , respectively, hydrostatic balance Equation 6 requires that the free surface is centrifuged into a parabolic shape (Cushman-Roisin & Beckers, 2011): where g = 9.81 m 2 /s and h 0 is the height of the parabola at s = 0 and is given by in a cylindrical annulus. It may seem odd to express the free surface shape (25) in terms of the parabola's height at s = 0 when no actual fluid exists within R i in the cylindrical annulus. However, expression (25) is useful because it remains unchanged whether we are describing paraboloidal annuli or full paraboloids. One simply sets R i = 0 in Equation 26 in order to describe non-annular, fully cylindrical systems (e.g., Cabanes et al., 2017;Zhang & Afanasyev, 2021).

Topographic β Comparisons
The stretching of planetary vorticity, controlled by βu s in Equation 16, is the dominant vorticity generation process in rapidly rotating geophysical fluid systems (Vallis, 2006). Thus, it is important to compare β(s) in spheres and paraboloids to see under what conditions their profiles are comparable.
In a sphere of radius R o , where χ = s/R o is the non-dimensional cylindrical radius and β o = 2Ω/R o is our non-dimensional scale for β.
For the paraboloid, the β formulation is The parameter  is the paraboloidal deformation scale and η is the cylindrical radius normalized by  . The ratio of the outer radius to  , Γ = ∕ , characterizes the degree of rotationally-induced parabolic deformation across a given fluid annulus. When Γ < 1, gravity dominates over centrifugation and the free surface remains relatively flat; when Γ > 1, centrifugation dominates over gravity and the free surface is strongly deflected. With Γ, we can relate the two nondimensional radii via = Γ and recast (28) as which is similar in gross morphology to Equation 27. The parabolic deformation scale  = √ 0∕Ω 2 is structurally similar to the external deformation radius that arises in the free surface responses of many rotating fluid systems. In particular, the external deformation radius shows up in the shallow water solutions for geostrophic adjustment and baroclinic instability (Vallis, 2006), where it estimates the radius within which inertial effects dominate over Coriolis effects. In constrast,  contains no dynamics here and instead represents the characteristic radial distance over which the parabolic free surface hydrostatically deforms. It should also be reiterated that the β expressions (27) and (29) are structurally rather similar, although not isomorphic. This suggests that β sph (s) and β para (s) will be adequate zeroth order proxies for one another, which, in turn, implies that vortex stretching dynamics in paraboloids will share qualitative similarities to spherical systems. Figure 5a plots unsigned profiles of β/β o versus χ = s/R o for the spherical case (red curve) as well for the three Coreaboloidal experimental cases carried out here at 35 rpm (green), 50 rpm (blue) and 60 rpm (cyan). All four profiles go to zero at the origin. At s = R o , β sph diverges whereas β para remains finite. The spherical curve | ∕ | intersects all three paraboloidal curves in the fluid bulk with crossing points, respectively, located at χ ≃ 0.5, 0.75 and 0.80. Thus, the normalized β values are of the same order of magnitude across most of the annular gap (χ ≲ 0.85). Figure 5a supports the conjecture made above that the spherical and paraboloidal β profiles will be comparable in structural complexity and will have similar characteristic magnitudes (except near s = R o ).

Conductive Heat Flux Comparisons
Parallel to the above subsection, conductive heat flux profiles are compared between paraboloidal and spherical systems here. The heat flux profiles provide a base level approximation for the thermal state for a given geometry. We assume that a steady, uniform heat flux q i is conducted through an inner cylinder or sphere of radius R i , passes through the fluid layer and is then travels out through the outer cylindrical or spherical boundary at R o . Thus, we assume here that there is no free surface heat transfer in the paraboloidal configuration. Conserving this flux as it is conducted through a spherical shell requires 4 2 = 4 2 ( ) , such that where s is employed in Equation 30 in place of r in order to facilitate comparison with the paraboloidal case. Carrying out these same steps in a paraboloid leads to where ℎ = − Ω 2 2 − 2 ∕(4 ).
The purely radial conductive profile (31) very roughly approximates our current experimental set-up in which the outer boundary is actually thermally insulating whilst the free surface is not. In future studies, we will strive to more accurately describe the Coreaboloid's thermal field by accounting for both radial heat transfer and axial heat transfer through the fluid's free upper surface. Figure 5b shows the idealized spherical (30) and paraboloidal (31) heat flux profiles normalized by the heat flux through the inner boundary and plotted as a function of normalized cylindrical radius χ = s/R o . The spherical profile (red) is in good agreement with the paraboloidal profiles, and remains nested between the profiles for the 50 and 60 rpm cases over the entire χ-range. This shows that a high degree of morphological similarity exists between idealized thermal base states of the spherical and paraboloidal systems.
Taken together, Figures 5a and 5b demonstrate that the essential vortex stretching properties and the thermal conductive states are comparable between the OTC region of a spherical shell and the paraboloidal annulus. It is this basic agreement in β and in conductive heat flux profiles that underlies our contention that convective flows in paraboloidal annuli will adequately simulate turbulence within the OTC regions of rotating spherical shells (i.e., Figure 2).

Convective Flow Scales
Laboratory and centrifugal gravity will act on the density differences associated respectively with ΔT ⊥ and ΔT ‖ to generate thermal winds that exist even in the base state of the system, following Equations 10 and 22. However, it is difficult to make accurate analytical predictions of these thermal wind flow fields, for example, (Lewis & Nagata, 2004), due to the geometric complexity of the cylindro-paraboloidal fluid layer and the related complexity of the temperature fields in the Coreaboloid. This differs from spherical shells in which the base state features spherically radial temperature gradients that generate no baroclinic azimuthal torques and thus no base state thermal wind flows (u = 0).
Thermal wind flows tend to become baroclinically unstable in the form of axially-aligned, non-axisymmetric flow structures (Read, Jacoby, et al., 2015;Vallis, 2006). The internal deformation scale, δ def , estimates the cross-axial width of these instabilities, We have modified the standard δ def definition (Read, Jacoby, et al., 2015) in order to account for thermal winds from both vertical lab gravity and horizontal centrifugation in our = (1) experiments. Estimates of δ def range between 2 and 3.4 cm for our experiments.
The secular cooling from the free surface of the initially heated fluid layer will tend to drive vertical rotating convective instability (Boubnov & Golitsyn, 1986;Nakagawa & Frenzen, 1955;Ravichandran & Wettlaufer, 2020). This second instability likely co-exists alongside the baroclinically unstable flow field, and also exists in the form of axially-aligned flow structures. The near onset scale of convective instability is (Stellmach & Hansen, 2004) LONNER ET AL.

10.1029/2022JE007356
14 of 26 with estimated δ onset values ranging between 1.2 and 2.9 cm in these experiments. Thus, the δ onset estimates are comparable to δ def in our experiments. Figure 6 shows an example of a mixed baroclinically-convectively unstable experiment. The case shown was carried out using an 18 cm radius DIYnamics set up (Hill et al., 2018), rotating at 16 rpm, with a 5 cm fluid layer depth, and an approximate radial temperature difference of 10 K. The red dye was emplaced in these experiments first and is patterned by a smaller scale vertical convective instabilities driven by surface evaporation from the fluid layer, for which δ onset ≃ 0.6 cm. Ice was placed in the inner cylinder after the convective instability had already set in. Then the blue dye was emplaced just adjacent to the inner cylinder. The red dye shows the smaller scale convective instabilities in this experiment, whereas the blue dye marks the larger scale baroclinic instability, with estimated scale δ def ≃ 3.8 cm. At longer times, the two instabilities spatially cohabitate and are both essentially space filling.
The disparate values of δ onset and δ def in this desktop experiment differs significantly from our three Coreaboloid cases. The similarity in baroclinic and vertical convective scales in our Coreaboloid experiments, δ onset ≈ δ def , make them difficult to deconvolve from one another. In future experiments, we will seek to separate these scales in order to better understand the dominant mechanisms at play in the system. With baroclinic instability also referred to as sloping convection (e.g., Hide & Mason, 1975), we will simply discuss thermally-driven convective flow in following discussions, without seeking to separate the baroclinic from the vertical convective instabilities.

Thermal Rossby Waves
Rapidly rotating OTC convection tends to propagate as so-called thermal Rossby waves (e.g., Busse, 2002;Busse & Or, 1986;Calkins et al., 2013;Hindman & Jain, 2022). The thermal Rossby waves dispersion relation in the plane wave, β-plane approximation is where σ is the angular frequency, U is the characteristic azimuthal velocity, k is the azimuthal wavenumber, l is the radial wavenumber, and K 2 = k 2 + l 2 . (We do not consider radially varying β effects here or global thermal Rossby modes, following the treatment of Lemasquerier et al. (2021).) The thermal Rossby wave azimuthal phase speed is then where we have assumed that k ≃ l. Expression (36) shows that phase fronts will be carried in azimuth by time-mean zonal flows, but will also propagate as thermal Rossby waves with speed −β/[(1 + Pr)(2k 2 )]. The k −2 dependence in Equation 36 implies that larger scale structures will have a higher phase speed. Further, for U ϕ = 0, the −β dependence in Equation 36 requires that thermal Rossby wave phase fronts will propagate in the +̂-direction in spheres (prograde) and in the −̂-direction in paraboloids (retrograde). Since |β| is comparable in spheres and paraboloids ( Figure 5), the Rossby wave phase speed, | | , should be comparable in paraboloidal and comparable spherical experimental systems.

Rhines Scaling and Jet Migration
In rapidly rotating turbulence, energy can naturally be transferred to larger flow scales (e.g., Boffetta & Ecke, 2012;Chen et al., 2006;Rubio et al., 2014). This so-called inverse cascade becomes strongly anisotropic in the presence of Figure 6. Example of a dually unstable rotating annulus experiment made using an 18 cm radius DIYnamics kit (Hill et al., 2018) (Vasavada & Showman, 2005). Kinetic energy then cascades into azimuthally-oriented zonal jets, which tend to saturate in their latitudinal extent on the Rhines scale (Rhines, 1975;Vallis & Maltrud, 1993) where U is the characteristic velocity of the vorticity containing flows (Gastine et al., 2014;Gillet et al., 2007;Heimpel & Aurnou, 2007). The spatiotemporal coherence of zonal jet flows is often characterized by the ratio of length scales (Galperin et al., 2006(Galperin et al., , 2014Sukoriansky et al., 2007), where is the smallest scale at which β affects the dynamics (Vallis & Maltrud, 1993), and where the energy transfer rate through the system is approximated via ϵ ≈ U 2 /τ Ek following Cabanes et al. (2017). Substituting in these length scales into Equation 38 yields We have expressed the Ekman drag time as . In calculating the  , we have used the time and radially averaged values of u ϕ , and the radially averaged values of β and h. Further, we have set all coefficients to unity in Equation 40, similarly to Cope (2020) and Lemasquerier et al. (2021). This lowers the calculated  slightly, giving a conservative estimate of its value. The zonostrophy index expresses the available range of scales over which β-modulated jet dynamics develop, with more coherent, steady jets forming in "zonostrophic",  ≳ 2.5 environments (Galperin et al., 2010;Sukoriansky et al., 2007). For comparison,  ≲ 2 in Earth's oceans, while Jupiter's jets are estimated to have  ≈ 5 . For hydrodynamic flow in Earth's core, we use the gap width as the characteristic system scale, Ek ≈ 10 −15 , and take U ≃ 1 mm/s (Barrois et al., 2018) such that R o ≈ 3 × 10 −6 . With these estimated values, we arrive at a mid-gap estimated value of  ≈ 10 in Earth's core. Thus, (non-magnetic) planetary core dynamics exist deep in the zonostrophic regime. In our thermally-driven Coreaboloid experiments, 2.5 ≲  ≲ 2.75 (Table 2), suggesting that these experiments simulate planetary zonostrophic dynamics.
Zonal jets have been observed to migrate perpendicular to their azimuthal flow direction, traveling in latitude in shallow layer systems (Chemke & Kaspi, 2015;Cope, 2020) and in cylindrical radius in quasi-two-dimensional models of deep spherical shells (Barrois et al., 2022;Gastine, 2019). This migration is due to asymmetries in the β-modulated Reynolds stress fluxes measured relative to the centerline, or core, of the jet, as shown in the modeling study of Chemke and Kaspi (2015). Cope (2020) argued, using a quasilinear approach, that zonon dynamics controlled jet migration. However, no study to date has provided a closed form theoretical migration rate prediction. Instead, Cope (2020) empirically found the zonal jet migration velocity, V mig , to scale as via a best fit to her survey of two-dimensional zonal jet simulations. In Equation 41, μ is the characteristic dissipation time scale and the coefficient value is c ≃ 0.8. Using β para , taking μ ≈ τ Ek and setting c = 1, relation (41) can be recast in the paraboloid, using all dimensional quantities, to yield ( where the right hand side term in square brackets is nondimensional. Expressions (42) and (43) will be compared below to the jet migration observed in our 60 rpm Coreaboloidal experiment.

Results
31. This leads to flow dynamics that remain qualitatively similar to those found in large gap, slowly rotating annulus models of planetary atmospheres (e.g., Bastin & Read, 1998;Hide & Mason, 1975;Wordsworth et al., 2008). Baroclinic instability in the vicinity of R i generates relatively large-scale coherent vortices that are shed into the fluid bulk (e.g., see Movie S1 in Supporting Information S1). The zonal flow corresponds to a set of two alternating jets, with strongly prograde flow near the inner boundary and retrograde flow near the outer boundary. The mean value of the zonal flow radial profile is not located at ⟨ ( )⟩ = 0 , which we postulate, for all our experiments, is due to the effects of axisymmetric thermal wind flows. This 35 rpm case demonstrates that slowly rotating experiments in the Coreaboloid overlap with the extensive studies made in the classical Hide's annulus configuration, and that slowly rotating cases can be understood in terms of the well-established Hidean dynamical framework (Read, 2001;Read, Pérez, et al., 2015).  Figure 7b and Movie S2 in Supporting Information S1, the thermal structures are smaller in this case relative to the 35 rpm case, and distinct, long lived coherent vortices are far less prominent in the 50 rpm case. The thermal structures are also more strongly filamentary than in Figure 7a, and compare well with the thermal field near the inner radius of the three-dimensional spherical shell simulation shown in Figure 2a. The time-averaged zonal flow profile is comprised here of five alternating jets. Similar to the 35 rpm case, the time-radial mean zonal flow is non-zero, and instead is net prograde. Overall, this case compares well to models of deep rapidly rotating spherical shell convection, as well as to the rapidly rotating GT10 and GT11 experiments in Smith et al. (2014). Figure 7c and Movie S3 in Supporting Information S1 correspond to the 60 rpm case, with its strong free surface deflection such that (h o − h i )/(R o − R i ) = 0.92. The characteristic scale of thermal structures is the smallest of the three cases, as should be the case for rapidly rotating convection (e.g., Aurnou et al., 2020;Guervilly et al., 2019). The overall morphology of the thermal field compares qualitatively well with the most rapidly rotating, three-dimensional numerical simulations of planetary core convection in the literature (e.g., Gastine et al., 2016;Schaeffer et al., 2017). Relative to Figure 2a, the convection appears to be weaker near R o in Figure 7 Figure 8a shows a stacked fast Fourier transform (FFT) of all the detrended γ = 0.5 temperature data, T(γ = 0.5, ϕ, t), acquired between t = 1,800 and 2,400 s. The peak in the spectrum (black dotted vertical line) occurs at wavenumber k max /(2π) = 0.262 cm −1 , corresponding to a characteristic azimuthal thermal field length scale of 1.9 cm for an individual thermal structure, exceeding our δ def by almost a factor of four. The spectral peak in panel A is relatively broad, due to the broad band of convective scales that exist in the turbulent convection field. The green and magenta vertical lines mark the upper and lower wings of the spectral peak with k/(2π) = 0.148 and 0.356 cm −1 , respectively. These correspond to azimuthal thermal field scales of 3.4 and 1.4 cm, respectively. Our δ def ∼ δ onset ∼ 1-3 cm estimates given in 3.4.1 both agree adequately with the thermal field length scales determined via the FFT analysis in Figure 8a. Further, the finite spectral breadth in the thermally-driven flow field here differs from the sharper spectral forcing employed in mechanically driven flows of Cabanes et al. (2017) and Lemasquerier et al. (2021). Each forcing method has its strengths: the monochromatic mechanical forcing allows for cleaner analysis, whereas the broad thermal forcing is likely to more accurately represent the broadband flow fields found in planetary settings. Figure 8b is a Hovmöller diagram of the mid-gap IR free surface temperature anomaly field plotted over time from t = 2,000-2,400 s in the 60 rpm experimental case. The thermal anomalies are calculated at fixed γ = 0.5 (s = 0.636R o ) as The thermal anomaly phase characteristics travel to decreasing ϕ over time, in basic agreement with the predicted retrograde behavior for a paraboloid. We further test expression (36) by plotting the predicted TRW azimuthal phase speeds using the three k values identified in Figure 8a, with the green, black dashed and magenta lines in Figure 8b corresponding to the values predicted by Equation 36 using the low, peak and high wavenumbers from the spectral peak in Figure 8a. Noting that the Hovmöller is rather complex overall, the thermal Rossby wave phase speed predictions agree moderately well with the dominant thermal phase drifts that exist in Figure 8b. Thus, we argue that thermal Rossby waves are the dominant convective mode in rapidly rotating Coreaboloid experiments, similar to those found in models of planetary core convection (e.g., Calkins et al., 2013;Carrigan & Busse, 1983;Hori et al., 2015). Further, our IR thermographic data represents the first reported laboratory thermal field mapping of turbulent thermal Rossby wave dynamics. Figure 9 shows UDV measurements from the 35, 50, and 60 rpm experiments shown, respectively, in panels a-c. The raw velocimetry data are projected into the azimuthal direction via Equation 3, yielding the u ϕ (γ, t) data shown. Red regions in Figure 9 correspond to prograde azimuthal velocities (u ϕ > 0) and blue regions correspond to retrograde velocities (u ϕ < 0). The azimuthal velocities appear to be statistically quasi-steady in these data. Further, the zonal jet scale is not changing during migration. These two points support our hypothesis that the data describe a kinematically equilibrated stated. However, the experiment is inherently transitory in its current form and it is always possible that time variability in the thermal field affects or even controls jet migration. The inner and outer cylindrical walls will be thermostated (e.g., von Larcher & Egbers, 2005) in future modifications to the Coreaboloid in order to test this possibility.
Vertical stripping in Figure 9 is due to temporal aliasing of the ultrasonic Doppler data that cannot be removed in post. These aliasing artifacts limit the accuracy of our time-mean zonal flow values, and likely cause our mean zonal flow values to be lower in amplitude than the true values. However, they do not appear to affect the overall patterns in the mean zonal flow fields.  Figure 1. In contrast to the 35 rpm case, this data set has more time-averaged zonal jets, with dominantly prograde flow adjacent to the inner boundary, near 0.6γ where β is reaches its maximum in this case, and adjacent to the outer boundary. In contrast, strong coherency of the zonal flow field is found in the 60 rpm,  ≃ 2.75 case shown in Figure 9c, in contrast to the two more slowly rotating cases. Here seven coherent zonal jets exist across between γ = 0.2 and γ = 1 for t > 2,040 s.
The innermost jets (γ ≲ 0.7) in Figure 9c migrate outwards to larger radial position over time, comparably to Smith et al. (2014). The solid black curves in Figure 9c represent the time integration of Cope (2020)'s jet migration formula (42). In these integrations, we estimate U using the peak zonal velocity for each jet core (e.g., 1.5 mm/s for the inner jet and 2.3 mm/s for the outer jet) and take into account the variation of β as a function of the jet's evolving radial position. The radial migration of the jets in Figure 9c show that turbulent, three-dimensional laboratory experimental jet flows are agree well with the empirical scaling to Cope (2020)'s two-dimensional jets simulations. The dashed black line in Figure 9b shows the integration of Equation 42 for far less coherent retrograde jet flow it lies adjacent to. This line is dashed because we are not convinced that migrating jets actually exists in this case.
In Figure 10 we test the idea put forth in Section 2 that the flows measurements describe a quasi-stationary kinematic state. Figure 10 shows UDV azimuthal velocity measurements, u ϕ (γ), with the green curve showing a time averaged profile over the first minute and the red curve showing the time averaged profile over the last minute of the 60 rpm case's V mig path in Figure 9c. Given the inherent variability in the system, it would be difficult to decipher which profile corresponds to which point in time if no other information were given, which supports the idea that the structure of the two profiles do not fundamentally differ over this time frame. They both have similar peak prograde and retrograde jet velocity values. The peak prograde velocity increases by 7% whereas the peak retrograde velocity decreases by 4% over this nearly 10 min, 560 rotation window. Further, even though a new prograde jet exists at γ ≃ 0.3 in the red curve, the total number of jets remain comparable. The ratio of azimuthal energy, ∫ 2 , for the two profiles yields a value of 0.79. This 21% decrease over time is predominantly due to the mismatch in zonal flow fields at γ < 0.3 and with the decrease in the peak velocity of the prograde jet that migrates from γ = 0.55 to γ = 0.65. Overall, our claim that the flows analyzed approximate a kinematically equilibrated state is supported by Figure 10.  Figure 11 shows the zonal velocities averaged in time from t = 1,800-2,400 s, ⟨ ⟩ , and plotted versus gap coordinate γ. The black line in each panel is the raw ⟨ ⟩ data. The colored line each panel is the raw data corrected by a best linear fit to it across the gap, which is shown in each panel via the red dashed line. The mid-gap thermal wind velocities are estimated via Equation 10 to be of order 1 mm/s. Thus, we postulate that these linear fits approximate the thermal wind flow fields in each case, which have not been modeled here. In order to focus on the time-mean Reynolds-stress driven zonal flows, we subtract the linear fits from the time averaged zonal flow profiles, yielding the time averaged jet flows, ⟨ ⟩ . In future studies, we plan to quantify the axisymmetric thermal wind flow via axisymmetric and three-dimensional spherindrical models (Ellison et al., 2022) of our paraboloidal experiments. The difference in the velocity ranges between Figures 9 and 11 results from the aliasing in the UDV data as well as from the migration of the jets. Taking a peak time-averaged jet velocity of U ≃ 2 mm/s in these experiments yields a Rossby number of Ro = U/(2ΩR o ) ≃ 7.2 × 10 −4 in the 35 rpm case and a characteristic Reynolds number of Re = (UR o )/ν ≃ 1.1 × 10 3 . Thus, the multi-jet flows found here arise from the quasi-geostrophic, zonostrophic turbulence that exist in these Ro −1 ≈ Re ≫ 1 experiments.
The time-averaged jet flows are shown together in Figure 12a. With roughly comparable zonal velocities in each case, this panel shows that the number of zonal jets increases as a function of increasing rotation rate, which corresponds to increasing β and τ Ek . The thin vertical dotted black line denotes the λ S = Ek 1/4 h(R o ) Stewartson sidewall boundary layer thickness, where the local height of the fluid layer at R o is used here in calculating Ek. All three cases give the same δ S estimate to within 1%. The Stewartson layer thickness agrees well with the thickness of the jets nearest R o , which are likely generated by processes connected to the sidewall-attached convective modes that exist within δ S (e.g., de Wit et al., 2020;Ecke et al., 2022;Favier & Knobloch, 2020;Vogt et al., 2021). Further,  the jets are larger than the predicted instability scales, δ def and δ onset , and the measured thermal field scales, for example, by a factor of approximately five in the 50 rpm case. Figure 12b shows the radial width of each jet L J , measured between zero crossings of ⟨ ⟩ normalized by the predicted Rhines scale (37) for each jet, L Rh . These Rhines scales are calculated for each individual jet using its peak zonal jet velocity for U and corresponding β(s) value at the location of the peak zonal velocity. This L J /L Rh ratio has values ranging from ≃0.8 to 2.5 in our experiments, with a mean value of 2.05 (black dashed horizontal line). Given the approximate nature of our analysis and of the Rhines scaling theory itself, Figure 12 provides evidence that Rhines scale jets dominate the zonal flow field in our buoyancy-driven experiments. The 50 and 60 rpm cases clearly show spatiotemporally coherent, multiple zonal jet flow can develop in thermally-driven laboratory experiments, thus extending the findings of, for example, Smith et al. (2014); Zhang and Afanasyev (2016); Cabanes et al. (2017); Lemasquerier et al. (2021). Further, L J ≃ L Rh is consistent with the validation in Figure 9 of Cope (2020)'s L Rh -dependent expression (41) for jet migration rate.

Discussion
This work presents the results of the first "production runs" made on our Coreaboloid experimental device. This novel paraboloidal thermal convection set up, extended from the mechanically-driven parabolic design of Cabanes et al. (2017), differs in many ways from a planetary spherical shell. Yet the Coreaboloid contains the essential ingredients of planetary core hydrodynamics: convective turbulence, rotational forces and topographic β effects. These ingredients lead to core-style flows comprised of convection in the form thermal Rossby waves as well as the existence of coherent, migrating Rhines scale zonal jets in the most rapidly rotating case. These zonostrophic zonal flows agree well with the gross structure of those found in quasi-2D models of planetary core convection (Barrois et al., 2022;Gastine, 2019;Guervilly & Cardin, 2017). Further, their migration rates are similar to those found in the two-dimensional shallow layer models of Cope (2020). Our results suggest that Rhines scale jets will tend to develop in high  planetary fluid cores, similarly to other turbulent geophysical systems such as surface oceans (e.g., Smith et al., 2014;, sub-surface oceans (Bire et al., 2022;Soderlund, 2019), planetary atmospheres (e.g., Lemasquerier et al., 2021), and stellar convection zones (e.g., Showman et al., 2019). However, in planetary cores and stellar systems, it may be that ambient magnetic fields leads to the formation of larger magneto-Rhines scale jet flows (Tobias et al., 2007), or that other magnetohydrodynamic balances dominate altogether (e.g., Cao et al., 2018;Orvedahl et al., 2021).
It is not clear why the zonal jets are far more coherent in the 60 rpm case than the 50 rpm case given that  ≃ 2.75 in both experiments. This result implies that the zonostrophy index alone does not describe the zonal jet dynamics. For simplicity, we employed a single temperature scale ΔT ⊥ in nondimensionalizing the governing equations in Section 3.2. In doing so, we have assumed that the large-scale centrifugal and gravitational buoyancy forces are fueled by identical vertical and horizontal large scale temperature differences, leading to a simple form of the Froude number, Fr = Ω 2 R o /g which is above unity for both the 50 and 60 rpm cases in Table 2. However, our lab measurements show significant differences between ΔT ⊥ and ΔT ‖ in all three cases (Table 1). Taking this into account, we formulate a new output parameter, the modified Froude number, ΔT ⊥ Fr/ΔT ‖ , which makes use of our laboratory measurements to more accurately estimate the system-scale centrifugal versus vertical buoyancy forces than Fr. Given in the fourth column of Table 2, the value of the modified Froude number is above unity only in the 60 rpm case (Table 2). Based on this, we hypothesize that centrifugation only dominates in the 60 rpm case, leading to more coherent thermal Rossby wave dynamics that, in turn, generate more coherent zonal jet flows.
The laboratory results presented here act as a proof of concept, demonstrating that the Coreaboloid provides a rich experimental framework for exploring planetary core dynamics. It is the flexibility and configurability of the paraboloidal system that makes it a broadly useful tool. By varying the rotation rate, Ω, and fluid layer thickness, H NR , it is possible to study quasigeostrophic turbulent dynamics over a broad (β(s), U, τ Ek ) parameter space (cf. Scott & Dritschel, 2012). It has been argued throughout this paper that this paraboloidal annulus is an excellent analog for outer core flow, but in the future it can be applied to other rapidly rotating geophysical and astrophysical systems, including the deep atmospheres of gas planets, sub-surface oceans and possibly stellar convection zones (e.g., Hindman & Jain, 2022;Vasil et al., 2021). We will also consider replacing water with GaInSn as the working fluid. By including an array of rare Earth magnets below the fluid layer, it may prove possible to simulate planetary core magnetoconvection, zonal flow processes and magnetohydrodynamic wave dynamics using this free surface paraboloidal platform.

Data Availability Statement
This study's thermographic and ultrasonic Doppler velocimetry data sets are available on the Dryad repository at (Lonner et al., 2022).