Canonical models of geophysical and astrophysical flows: Turbulent convection experiments in liquid metals

© 2015 by the authors; licensee MDPI, Basel, Switzerland. Planets and stars are often capable of generating their own magnetic fields. This occurs through dynamo processes occurring via turbulent convective stirring of their respective molten metal-rich cores and plasma-based convection zones. Present-day numerical models of planetary and stellar dynamo action are not carried out using fluids properties that mimic the essential properties of liquid metals and plasmas (e.g., using fluids with thermal Prandtl numbers Pr < 1 and magnetic Prandtl numbers Pm ≪ 1). Metal dynamo simulations should become possible, though, within the next decade. In order then to understand the turbulent convection phenomena occurring in geophysical or astrophysical fluids and next-generation numerical models thereof, we present here canonical, end-member examples of thermally-driven convection in liquid gallium, first with no magnetic field or rotation present, then with the inclusion of a background magnetic field and then in a rotating system (without an imposed magnetic field). In doing so, we demonstrate the essential behaviors of convecting liquid metals that are necessary for building, as well as benchmarking, accurate, robust models of magnetohydrodynamic processes in Pm ≪ Pr < 1 geophysical and astrophysical systems. Our study results also show strong agreement between laboratory and numerical experiments, demonstrating that high resolution numerical simulations can be made capable of modeling the liquid metal convective turbulence needed in accurate next-generation dynamo models.


Introduction
System-scale magnetic fields are observed to develop in galaxies, stars, planets, moons and even asteroids [1][2][3][4]. In all these systems, magnetohydrodynamic dynamo action occurs via the conversion of the kinetic energy of flowing electrically-conductive fluids into magnetic field energy [3,5]. Generated by the motion of the electrically-conductive fluid, dynamo fields typically evolve on length and time scales related to the underlying flows [6].  [7]. The field is dominated by the axial dipole component, but also has rich higher order structure; (b) Secular variation of the 2005 radial CMB field, measured in units of µT /yr [7]. The radial field component can change by a few percent per year; (c) Schematic of Earth's interior with the silicate mantle overlying the iron-rich core, all rotating with angular velocity Ωẑ. Convective flows in the outer core are likely the key driver of the geomagnetic field. (Convection image adapted from [8].) For example, Earth's geomagnetic field is shown in Figure 1a [7]. The field is predominantly dipolar, with magnetic flux emerging from Earth's southern hemisphere and returning via the northern hemisphere. The large-scale field cannot be generated in the silicate mantle because temperatures there exceed the Curie temperature of its magnetic minerals at relatively shallow depths, too shallow to generate the 0.1 mT fields observed on the planet's surface. Furthermore, the fields vary significantly in time (Figure 1b). Such time variations cannot be explained by permanent magnetization of planetary materials. Thus, this strong, spatiotemporally complex magnetic field is best explained via continual magnetohydrodynamic (MHD) induction occurring through the action of turbulent liquid metal flows that exist in Earth's iron-rich molten outer core (Figure 1c). Similarly, large-scale flows in the sun drive its famous quasi-periodic 22-year magnetic field reversal cycle [9].
However, it is unlikely that there is a simple, direct mapping between flows and magnetic fields in natural dynamo systems. This is predominantly due to the particular material properties of the metallized fluids in which these dynamos are generated.
In galaxies and stars, the electrically-conducting fluid is typically a plasma; in planets and asteroids, the electrically-conducting fluid likely resides within a liquid metal core (Figure 1a). The (inverse) thermal conductivity of the fluid is often characterized non-dimensionally via the Prandtl number, where ν is the viscous diffusivity, κ = k/(ρC p ) is the thermal diffusivity, and ρ and C p are, respectively, the fluid's density and heat capacity. In liquid metals and astrophysical plasmas, the high thermal conductivity of the fluid yields typical P r values in the range 10 −2 to 10 −7 . Similarly, the electrical conductivity of the fluid is characterized by the magnetic Prandtl number, where η = 1/(µ o σ) is the magnetic diffusivity, µ o is the permeability of free space and σ is the electrical conductivity. In liquid metals and astrophysical plasmas, the high magnetic diffusivity of the fluid yields typical P m values in the range 10 −5 to 10 −9 . The ordering of the values of the different diffusivities, has interesting physical implications for convection-driven MHD processes in geophysical and astrophysical fluids. This ordering implies that the characteristic length and time scales of the magnetic fields will tend to exceed those of the thermal structures, which will again exceed the characteristic scales of the underlying flow fields, e.g., [10]). Indeed, energy is expected to be injected at small-scales via buoyancy-driven convection processes in natural dynamo systems ( [11][12][13]; Figure 1c), and is likely to then cascade upscale to larger scales [14], where we hypothesize that dynamo action likely occurs more efficiently [15] and where large-scale magnetic structures are likely formed. Thus, a detailed understanding of turbulent convection processes in metals is essential for the development of broadly accurate models of geophysical and astrophysical dynamo processes. Unfortunately, relatively little information exists to describe convective turbulence in liquid metals. Even less information exists concerning turbulent convection in liquid metals in the presence of magnetic fields and rotation, as occurs in geophysical and astrophysical settings.
To address this knowledge gap, we have developed a mixed laboratory-numerical platform to investigate the essential properties of small-scale rotating magnetoturbulence in convecting liquid metals. In studies of all types of convecting fluids, an idealized, simplified version called Rayleigh-Beńard convection (RBC) is the canonical, reference system [16]. In RBC, fluid is contained between two co-planar, impermeable boundaries, whose unit normal vectors are aligned with gravity. Typically, the boundaries provide some spatially-uniform thermal forcing, such that an adverse vertical temperature gradient is imposed across the fluid, with warmer fluid underlying cooler fluid. Because the local fluid density is inversely proportional to the temperature, the system releases potential energy when the warm (lower density) fluid becomes convectively unstable to upward motions and the cool (denser) fluid to descending motions.
In order to overcome the magnetic diffusivity of liquid metals and plasmas, generating a dynamo requires either a vast amount of fluid or an extremely strong forcing, e.g., [17]. Instead of studying the generation of a dynamo magnetic field from a turbulent flow, we propose to take the inverse approach and instead investigate the effects of imposed magnetic fields and system rotation on a given RBC system. Thus, we consider such RBC systems to represent idealized parcels of a planet's or star's dynamo generation region in which dynamo action does not locally occur, and instead, the magnetic field, assumed to be generated on larger scales, is imposed.
The manuscript is organized as follows. In Section 2, we describe our laboratory experimental device and associated numerical modeling framework. In Section 3, we present the results of a mixed laboratory-numerical investigation of Raleigh-Bénard convection (RBC; non-magnetic, non-rotating) in a cylindrical cavity, using the liquid metal gallium (P r 0.025; P m 10 −6 ) as the working fluid. (See Table 1 for the thermophysical properties of Ga.) In the parameter ranges accessible to our experiments, turbulent convection in liquid gallium is dominated by large-scale inertially-dominated flywheel-style flows. In Section 4, we show how the imposition of a vertical magnetic field constrains the flow, lessening the inertial effects dominant in the RBC system. In Section 5, we remove the magnetic field and instead rotate the system. Rotation also constrains the fluid dynamics, but mechanically instead of electromagnetically. Rotation leads to oscillatory weakly-nonlinear convective motions and differing heat transfer behaviors than in the magnetoconvection experiments. Finally, in Section 6, we summarize our laboratory-numerical findings and discuss their relevance to our understanding of natural dynamo systems and present-day numerical dynamo models. Table 1. Thermophysical properties of liquid gallium. Estimates from [19,20].

Laboratory and Numerical Methods
We have developed a new coupled laboratory-numerical tool, well-suited to investigate convective turbulence in liquid metals. Here, we will describe each component of this novel system: first the laboratory experimental device and then the numerical modeling framework.

Laboratory Experimental Framework
In the laboratory, we have conducted thermal convection experiments using liquid gallium as the working fluid. Our convection tank is a right cylindrical chamber with an outer diameter of 20 cm (see Figure 2a,b). The device can accommodate sidewalls of varying heights. In this study, we have utilized two tank heights, one with a H = 5 cm high stainless steel sidewall (3.0 mm wall thickness) and the other with an H = 10 cm high acrylic sidewall (6.35 mm wall thickness). The 5 cm tank has an aspect ratio Γ 19.6/5.0 = 3.92 and the 10 cm tank has an aspect ratio Γ 18.73/10.0 = 1.87, where Γ = D/H is the ratio of the diameter of the fluid layer D versus height H. To minimize heat losses, we thermally insulate the system with a 20 cm layer of thermal blanketing that surrounds a 10 cm thick inner layer of closed-cell foam insulation. In addition, several thermocouples are inserted within these layers in order to quantify heat losses.
As shown in Figure 2a, the bottom thermal block consists of a 1 cm thick copper plate; and the top thermal block consists of a 2.5 cm thick copper plate. In order to determine the full temperature difference across the fluid layer, ∆Θ, the top and bottom thermal blocks each contain 6 lid thermistors, homogeneously distributed in azimuth, and located 1 mm above and below the fluid layer, respectively. In addition the top thermal block has 5 thru-holes ( Figure 2b) allowing us to insert thermistors of various lengths into the fluid interior. In the 5 cm high tank, the internal thermistors are located at depths of: 4 cm (T 1 ), 1cm (T 2 ), 2.5 cm (T 3 ), 1 cm (T 4 ) and 1.25 cm (T 5 ). In the 10 cm high tank (not shown in Figure 1), the internal thermistor depths are: 5 cm (T 1 ), 2.5 cm (T 2 ), 5 cm (T 3 ), 1.5 cm (T 4 ) and 2.5 cm (T 5 ). The internal thermistors have depth tolerances of ±0.5 mm. T 3 is always located at the center point of the tank. The other thermistors allow us to measure the 3D structure of the temperature field.
The fluid must be heated from below and cooled from above to drive thermal convection in the fluid layer. Thus, the device has a 19 cm diameter, non-inductively wound, electrical resistance heater, capable of delivering 4900 W of heating power, that is attached to the base of the bottom thermal block (Figure 2a). The heat is removed by a water-based heat exchanger that is situated above the top thermal block. The attached cooling block is a right cylinder of aluminum into which has been cut a double wound spiral flow channel. This double wound spiral acts to minimize temperature gradients in the block [18]. The temperature of the water is thermostated by a precision chiller capable of removing up to 10 kW of heating power while maintaining the bath temperature between a 5 o C and 35 o C, stable to within ±25 mK.
Temperature data is acquired via a digital data acquisition system capable of simultaneously sampling up to 32 thermistor channels, typically at a sampling rate of 10 Hz. In addition, 32 more channels are available to measure temperatures via K-type thermocouples, which we use to evaluate thermal losses and to monitor ambient temperature in the room. Figure 2. Schematic of the right cylindrical laboratory and numerical convection set-ups used in this study. (a) Side view of the experimental convection device. This illustration shows the device with the 5 cm high stainless steel sidewall in place. Lid thermistors measure the temperature difference across the fluid layer, ∆Θ; internal thermistors, T 1 through T 5 , measure local (dimensional) temperatures, T , inside the fluid layer; (b) Plan view (upward facing perspective) of the laboratory experiment's top thermal block; (c) Illustration of a meridian slice of the mesh used in the numerical simulations. In this schematic, the mesh shows only one-tenth of the actual numerical resolution. Points P 1 through P 5 denote temperature measurement locations of the numerical (non-dimensional) temperature field Θ; (d) Plan view of the numerical domain. Note that the laboratory and numerical internal temperature measurements coincide in spatial location, and thus can be compared precisely.
In laboratory convection systems, it is important to ascertain the thermal properties of the boundaries, in order to be comparable with theoretical models, cf. [21]. To evaluate the thermal boundary conditions, we determine the value of the Biot number, Bi, for the top and bottom thermal blocks. The Biot number, which is the ratio of the effective thermal conductance of the fluid layer and the boundary material, is defined as: where K ga = 31 W/(m K) is the thermal conductivity of gallium, K cu = 401 W/(m K) is the thermal conductivity of the copper in the thermal blocks, H bl the thickness of the thermal blocks and H ga is the thickness of the gallium layer. For the 5 cm tank, we estimate the maximum values of the Biot number to be Bi = 0.064 for the bottom block and Bi = 0.092 for the top block. For the 10 cm high tank, we estimate maximum values of Bi = 0.046 for the bottom block and a Bi = 0.093 for the top block. For Biot number values less than 0.1, the thermal boundaries typically remain isothermal to within a few percent [22]. Our experimental Bi values fall within this range (Bi 0.1), such that both thermal boundaries should be treated as being isothermal. These arguments are verified by our laboratory magnetoconvection and rotating convection experiments, in which convection first onsets in accordance with linear predictions for isothermal boundaries. This was also found to be the case in a preceding study of rotating convection in a larger volume tank of gallium [20].
The entire device described above is situated on a rotary table that is belt-driven by a 10 N.m brushless servomotor. The device can be rotated between 0.1 and 2π rad/s (i.e., 0.5 to 60 revolutions per minute), steady to within ±0.1%.
A 35 cm bore solenoidal electromagnet is situated in the (non-rotating) lab frame and can be lowered on jack screws around the convection tank. The solenoid can supply a vertical, imposed magnetic field of up to 0.1325 T that is uniform to within ±0.5% over the volume of our largest gallium convection tank, which is approximately 20 cm 3 .
Thus, this laboratory experimental device allows us to investigate turbulent convection in liquid metal subject to magnetic forces, rotational forces as well as their combined effects, as are all relevant to our understanding of natural systems. For further device details, we refer the reader to [23].

Numerical Modeling Framework
In close association with our laboratory experiments, we have carried out numerical simulations using a modified version of the SFEMaNS code, which has been used previously in studies of mechanically-forced magnetohydrodynamic flows in liquid metals [24][25][26]. SFEMaNS is a hybrid code that makes use of Fourier series in the azimuthal direction. In meridian planes, it uses nodal quadratic Lagrange finite elements for the velocity field and nodal linear Lagrange finite elements for the pressure field. The discretization of time is made by means of a semi-implicit, second-order backward finite difference method (BDF2). This code is efficiently parallelized in Fourier space and on meridian sections (see Appendix 1). The domain decomposition is made using METIS [27], and the parallelization is made using MPI and PETSC libraries [28,29]. The code's unstructured meridional mesh, shown in Figure 2c, allows SFEMaNS to be used to simulate MHD flows in any axisymmetric geometry.
The key modification we have made to the code is the inclusion of a thermal field solver, allowing for thermally-driven MHD simulations of Boussinesq fluids in which the density is constant except in the buoyancy term [30,31]. This new version of the code, called SFEMaNS-T, has already been used successfully to simulate convection-driven, planetary dynamo action in a spherical shell geometry in the community benchmarking exercise of [32].
In this study, we carry out numerical calculations in right cylindrical geometries (see Figure 2c,d). The numerical tank geometries are nearly identical to our laboratory experimental cases, with numerical aspect ratio values of Γ = 2 and 4. The points P 1 through P 5 (shown in Figure 2c,d) represents numerical probes that are directly equivalent to the laboratory thermistors T 1 through T 5 (shown in Figure 2a,b), allowing point-to-point comparisons of the temperatures in each system.
Within the fluid layer, denoted as domain I in Figure 2c,d, SFEMaNS-T simultaneously evolves the velocity, thermal and magnetic fields by solving the Boussinesq forms of the Navier-Stokes momentum conservation equation in a reference frame rotating at constant angular velocity Ωẑ; the energy conservation equation; and the magnetic induction equation, derived from Maxwell's equations under the MHD approximation in which the velocities of the electrically-conducting fluid are assumed to be far less than the speed of light. For detailed derivations of these equations see [33][34][35].
The non-dimensional governing equations that SFEMaNS-T iteratively solves are subject to the the solenoidal conditions ∇·u = 0 and ∇·B = 0 for the non-dimensional velocity and magnetic fields, u and B, respectively. The non-dimensional dynamic pressure is denoted p and the non-dimensional temperature is given by Θ.
In the momentum Equation (4), the left hand side terms give the Lagrangian acceleration of the fluid. The terms on the right hand side are, respectively, the pressure gradient force (per unit volume), the buoyancy force, the viscous force, the Lorentz force and the Coriolis force. The centrifugal force is absorbed inside the pressure term. Lorentz forces arise in this via a uniform, axial magnetic field B oẑ that can be imposed. Similarly, the system can axially rotate at a steady angular velocity Ωẑ, leading to Coriolis forces.
In the energy Equation (5), the Lagrangian time derivative is on the left hand side and thermal diffusion is on the right hand side. As there are no internal buoyancy sources in this system, the thermal buoyancy forcing must be maintained by the thermal boundary conditions.
In the magnetic induction equation (6), the magnetic field is generated by the first term on the right hand side, the "stretching" term, and it is diffused by the second term on the right hand side. Flows can generate magnetic fields, but in doing so, they generate Lorentz forces that act back on the velocity field such that system's total energy is conserved.
Non-dimensional control parameters are grouped in front of the various terms of Equations (4)- (6). To arrive at these groups, we have non-dimensionalized the governing equations in the following way. Velocities are scaled using the free-fall velocity estimate for an undamped thermally buoyant parcel that traverses the fluid layer, where α is the thermal expansivity of the fluid, g is gravity oriented in the verticalẑ-direction, ∆T is the vertical temperature difference across the fluid layer, and H is the vertical height of the fluid layer.
Typically, U f is the highest velocity that can be reached in a convection system, where all the fluid's buoyancy is transferred into inertia, u · ∇u = αg∆T . This choice of velocity scale is useful in analyzing convecting liquid metals, in which strongly inertial flows tend to develop [48]. Length is scaled by H. After subtracting out the background hydrostatic pressure field, the residual, dynamic pressure p is scaled by ρU 2 f /H, where ρ is the mass density of the fluid. The magnetic field is scaled by the strength of the vertically-imposed background field B o .
We have previously defined the two Prandtl numbers that non-dimensionally describe the fluid properties. Three other non-dimensional control parameters arise in Equations (4)- (6). The Rayleigh number estimates the effective buoyancy effects relative to diffusion and is defined as the ratio of the thermo-viscous diffusion time scales, τ κ ∼ H 2 /κ and τ ν ∼ H 2 /ν, and the squared free-fall time, τ f , When Ra is large, buoyancy effects dominate over diffusion. The Chandrasekhar number, is the ratio of the magneto-viscous diffusion time scales and the squared Alfven speed crossing time, where τ η ∼ H 2 /η and τ A ∼ H/U A . Here, U A = B 2 o /(ρµ o ) is the Alfven wave speed, where µ o is the permeability of free space. Large Q implies that the Lorentz forces dominate the viscous forces in the flow. Lastly, the Ekman number is defined as which is the ratio of the characteristic rotation time scale, τ Ω ∼ 1/(2Ω), and the viscous diffusion time scale. In this case, small values of E imply that the rotational effects dominate over viscous effects. The Prandtl numbers defined in (1) and (2) can be re-cast as P r = τ κ /τ ν and P m = τ η /τ ν , and the free-fall velocity estimate can be re-cast as U f ∼ (Ra/P r) 1/2 (H/τ ν ) = (Ra/P r) 1/2 (ν/H).
For the velocity field, non-penetrative, non-slip conditions are applied on all the boundaries of fluid layer, which is denoted as domain I in Figure 2c. The non-dimensional temperature Θ is scaled by the vertical temperature difference ∆T . The top and bottom boundary conditions are isothermal with values of Θ(z = 0.5) = 0 and Θ(z = −0.5) = 1 , respectively. The cylindrical sidewalls are thermally insulating: In order to correctly solve for the magnetic field, the fluid region is surrounded by a "vacuum" region, denoted as domain II in Figure 2c,d. In the vacuum, no electrical currents can exist and, thus, the magnetic field satisfies a potential field solution. The continuity of the magnetic field between domains I and II is weakly imposed on the interfaces Σ 2 , Σ 4 and Σ 5 in Figure 2c,d through an interior penalty method, e.g., [24]. The magnetic field is set to zero on the outer boundary of domain II. Our finite element solution requires no special conditions to be imposed along the axis of rotation.
We set the initial temperature field to the diffusive solution: where the non-dimensional grid is always defined such that z = ±0.5. Based on these definitions, we also define the temperature anomaly ζ defined as: which gives the temperature perturbation on a given z-level with respect to the diffusive solution. For the velocity field, we apply a small amplitude, random initial perturbation on all Fourier modes.

Convective Free-Fall Parameters
The control parameters in Equations (4)-(6) are the same as those used to describe the laboratory experiments. However, P r, Q and E are all ratios relative to diffusive time-scales, whereas we have non-dimensionalized the equations using the (strongly turbulent) free-fall velocity and time scales, U f and τ f . To better account for this choice of free-fall non-dimensionalization, we define the following "convective" parameters.
The convective Reynolds number is defined as When Re c is large (and constraining magnetic and rotation forces are comparatively small), then fluid inertia dominates viscous diffusion of momentum and the convective flows are turbulent [36]. The convective Peclet number is defined as Large values of P e c imply that thermal advection processes dominate thermal diffusion effects. In such cases, thermal turbulence tends to mix and often isothermalize the interior of the convective fluid layer, e.g., [16]. The convective magnetic Reynolds number is defined as For small Rm c , Ohmic dissipation dominates magnetic induction such that the magnetic field never greatly differs from the imposed field B oẑ ; for large Rm c , induction can dominate and self-sustaining dynamo action can become possible [5,33]. The convective interaction parameter is This parameter estimates the ratio of Lorentz and inertial forces in the momentum equation. For large N c , Lorentz forces dominate fluid inertia [33,34]. The convective Rossby number predicts the ratio of fluid inertia and Coriolis force and is defined as For small Ro C , convective turbulence is likely to be constrained by the effects of rotation, e.g., [37]. With these free-fall-based parameters, the governing equations are then re-cast as: In the following sections, we will present our results typically in terms of the standard parameters Ra, Q and E. We will typically use the standard parameters when describing physical phenomena occurring near to the onset of convection, where molecular diffusive effects are often relevant. In contrast, we will usually interpret the physics using the convective Reynolds number based parameters (16)- (20) in fully-developed, turbulent cases in which the flow speeds may be approaching U f and advection processes overwhelm diffusive transport phenomena, cf. Table 2.

Compilation of Heat Transfer Measurements
The behavior of thermal convection systems is often characterized by measurements of the efficiency of global heat transfer across the fluid layer. This heat transfer efficiency is parameterized by the non-dimensional Nusselt number, N u, which is the ratio of total heat transfer to that transported only by thermal conduction, where q = q cond + q conv is the total heat flux, which is the sum of the conductive heat flux q cond = K∆T /H and the convective heat flux q conv = ρC p u z Θ. Here, K is the fluid's thermal conductivity in and the vertical component of fluid velocity is u z . For weak buoyancy forcing (i.e., sufficiently low Ra), the fluid layer is stable and perturbations in the flow and temperature field will decay to zero. Thus, the vertical velocities will tend to zero and the heat flux must diffuse across the fluid layer via thermal conduction alone, such that N u = 1. Then at the so-called critical Rayleigh number, Ra crit , the fluid in the layer becomes gravitationally unstable and convective motions onset. Warmer fluid parcels will tend, on average, to rise through the fluid layer and colder fluid parcels will tend to fall, releasing gravitational potential energy stored in the initial configuration of the density field in which least dense fluid lies at the base of the fluid layer. These convective motions advect heat across the fluid layer and the Nusselt number will tend to increase above unity. In order to measure N u in our laboratory experiments, we obtain measurements of the total temperature difference, ∆T , across the thermistors embedded in the top and bottom thermal blocks, and the heat flux, q supplied by the basal heat pad. The material properties of the liquid gallium, such as k, are calculated as in [20]. Figure 3 is a compilation of all our heat transfer measurements. The non-dimensional buoyancy forcing, Ra, is plotted on the x-axis versus the heat transfer efficiency, N u, on the y-axis. Rayleigh-Bénard convection experiments are demarcated by open square symbols. Magnetoconvection (MC) experiments, carried out at fixed Chandrasekhar number Q = 3.5 × 10 3 , are marked by diamond shaped symbols. Lastly, rotating convection (RC) experiments, all made using a fixed Ekman number of E = 2 × 10 −5 , are represented by circular symbols. Open symbols denote laboratory experimental results and red-filled symbols correspond to numerical simulation results. The colors of the symbol lines express the fluid layer geometry: black lines correspond to Γ 4 cases; blue corresponds to Γ 2 cases; and green corresponds to Γ 1 cases. (The Γ 1 rotating convection laboratory data is from the study of [20].) The three different suites of experiments, RBC, MC and RC, all show significant differences in their heat transfer behaviors. The lowest value of N u in the RBC data set is near two. Exceeding N u = 1, this implies that the Rayleigh numbers investigated in this set of experiments all exceed the critical Rayleigh number for RBC. In addition, the RBC data appears to display a single power law behavior, suggesting that we are interrogating a single behavioral regime of liquid metal convection. This scaling behavior is discussed in detail in Section 4. In contrast, the N u-Ra data shows that the zeroth order effect of either an imposed magnetic field or an imposed system rotation is to stabilize the layer of liquid metal against convective motions. In the MC cases, the critical Rayleigh number is just below Ra = 3 × 10 4 . In the RC cases, convection onsets closer to Ra = 1.5 × 10 5 . After convection onsets in the MC and RC systems, the heat transfer behaviors are more complex than in the RBC experiments. Just after the onset of convection, the MC N u data rises sharply with Ra and has a concave downward curvature until, at a sufficiently high Rayleigh number of about 10 6 , the MC data merges with the RBC data. In contrast, the RC data has a very shallow slope just after convection onsets. This gives way to a steeper trend at close to Ra = 10 6 , which then intersects with the RBC data in the vicinity of Ra = 10 7 .
The data in Figure 3 shows that the MC and RC return to the RBC trend at sufficiently high Ra. Qualitatively, this behavior is unsurprising: MC systems behave similarly to RBC systems when N c is sufficiently low and RC systems behave similarly to RBC systems when Ro c is sufficiently high. Thus, the RBC heat transfer data provides an approximate upper bounding curve the MC and RC heat transfer data sets. This further implies that knowledge of RBC convective scaling behavior is necessary to understand the behaviors of the MC and RC systems. A similar result exists for rotating convective heat transfer in moderate P r fluids, such as water [13]. Thus, the import of understanding RBC as the essential basis for more complicated convection systems holds broadly.
In the following three sections, we will consider in greater detail the behaviors of Rayleigh-Bénard convection in liquid metal, magnetoconvection and then rotating convection.

Essential Theory
Without the effects of magnetic fields or background rotation, the governing equations for RBC simplify to subject still to the solenoidal ∇ · u = 0 condition for Boussinesq flow. The critical Rayleigh number at which RBC onsets is independent of P r and occurs via steady (not oscillatory) motions at Ra crit = 658 (Ra crit = 1708) for impenetrable, stress-free (non-slip) boundary conditions [31,33,34,38].
In strongly forced (high Ra) convection systems, buoyancy-driven flows develop with relatively high convective Reynolds numbers, Re c = Ra/P r. For Re c values above roughly 10 3 , multi-scale, three-dimensional (3D), turbulent motions will tend to develop in RBC systems. However, in liquid metals, the low value of P r means that the P e c = √ Ra P r may remain low even in relatively high Re c settings, e.g., [39]. Because of this, strong convective flows in liquid metals (e.g., Re 10 3 ; P e 10 3 ) will often be unable to strongly mix the thermal field. Instead, such flows will tend then to create large amplitude, large scale temperature perturbations in the fluid bulk (e.g., see Figure 1c in [40]). Vigorous small-scale flows in liquid metals will be largely unable to perturb the thermal field.
In RBC, convective heat transfer is parameterized via scaling laws predicting the value of the Nusselt number, N u, as function of Ra and P r: where α and β are the powers of the scaling laws.
In moderate Prandtl number convection (P r 1), studies often find the heat transfer is controlled by boundary layer physics resulting in α values between 2/7 and 1/3 [13,16,36,41,42]. The α 1/3 law is predicted to arise under the following when vigorous convective mixing creates a nearly isothermal fluid bulk (P e c 1). Then the temperature gradients in the fluid become concentrated into the quasi-static regions of thickness δ that exist adjacent to the top and bottom boundaries (where impenetrability requires u z → 0). In these boundary layers, the heat transfer is dominantly diffusive and the temperature drop, δT , across each of these regions must scale as δT 1/2∆T . Then the Nusselt number can be re-cast as If it is assumed that the thermal boundary layers operate independently of one another at large N u, then the heat transfer should be independent of the fluid layer depth H, e.g., [43]. For fixed P r, dimensional analysis then requires N u ∼ Ra 1/3 [43][44][45]. This result is indeed found in experiments carried out at Ra 10 10 ; N u 100. At lower Ra values, it is typically found that RBC heat transfer follows a N u ∼ Ra 2/7 scaling. It was first postulated by Shraiman & Siggia (1991) [46] that this α 2/7 scaling occurs when the boundary layers are subjected to a large-scale shear flow. They argue that this shear is typically produced by large-scale circulations, making the heat transfer dependent on the container geometry. In this α 2/7 regime, the heat transfer is still controlled by the boundary layer physics, but the shear modifies the heat transfer efficiency (Recent, two-dimensional convection simulations by Goluskin et al. (2014) [47] suggest that such shear inhibition effects may be capable of lowering the heat transfer coefficient well below 2/7.) In sharp contrast, boundary layer processes do not necessarily control the heat flux in low P r fluids. Instead, heat transfer is controlled by large scale flows in the fluid bulk. Strong bulk flows naturally arise in low P r fluids because inertial effects can dominate the system even at the onset of convection [48]. Various laboratory and numerical studies have shown that convection in liquid metals or low P r fluids tends then to be dominated by inertial flywheels [49]. In an inertial flywheel, the fluid motions cascade upward in scale, nearly to the size of the container, and approach the free-fall velocity of the system in which fluid buoyancy is transferred directly into inertia. However, such high Re c , moderate P e c flows are unable to isothermalize the fluid bulk, as in moderate P r fluids. Instead, these high velocity, large scale flows are able to generate container-scale gradients in the perturbation temperature field [40,50].
It is hypothesized that in these low P r, inertial flows heat transfer will follow an α = β = 1/4 heat transfer law [20,49]. This scaling behavior arises from a balance between thermal advection by the inertial flywheel flows and thermal diffusion across the quasi-static boundary layers: Substituting in the free-fall velocity relation (7) and using equation (28) leads to Note that low P r inertial flywheel flows are predicted to transfer heat less efficiently (α 1/4) than convection occurring in moderate P r fluids (α 1/3). Heuristically, this difference in α values arises because the rate at which heat is transferred across the fluid layer in inertial systems depends on both the free fall time τ f necessary to advect heat across the fluid layer and the time needed to diffuse heat across the thermal boundary layers. In contrast, in moderate P r convection, the rate of heat transfer is higher because it depends only upon the time needed to diffuse heat across the thermal boundary layers as heat is shunted directly across the isothermal interior fluid.
However, the inertial heat transfer scaling (30) is not firmly established in the liquid metals literature. In fact, experiments in liquid metals have yielded a wide range of results providing α values from 1/5 up to 1/3 [18,20,[51][52][53][54]. This broad range of possible α values shows that heat transfer in liquid metals is not as well understood as it is in moderate P r fluids [16].
In our RBC experiments presented below, we show that a α 1/4 scaling law best fits the heat transfer data over the range 2 × 10 4 Ra 3 × 10 7 . Further, we demonstrate that large-amplitude temperature variations arise on the scale of the various inertial flywheel morphologies that develop in our experiments. Irrespective of the morphology of the flywheel flows, we find a the mid-plane thermal anomaly to remain near to 60% of the vertical temperature difference across the fluid layer. Table 3. Input and output parameters in our suite of laboratory Rayleigh-Bénard convection (RBC) experiments made using liquid gallium. The fluid layer depth is H. The basal heating power is P , which defines the heat flux q = P/A with A = π(D/2) 2 2.8 × 10 −2 m 2 denoting the fluid layer area. The mean temperature of the fluid is T and the vertical temperature difference is ∆T . The fluid layer aspect ratio is Γ = D/H. The parameters P r, Ra and N u are defined in Equations (1), (8) and (24), respectively. The rightmost column show the associated numerical RBC cases.

Experimental Results
Along with twenty-seven laboratory experiments (Table 3), we have carried out four fully three-dimensional numerical simulations, RBC 1 n through RBC 4 n made using SFEMaNS-T (Table 4). The input parameters in these numerical cases correspond closely to the four laboratory experimental cases demarcated RBC 1 l through RBC 4 l in the left hand column of Table 3. Table 4. Input and output parameters for the numerical RBC cases. The Prandtl number has a fixed value of P r = 0.025. The number of meridional grid points is n p , the number of meridional subdomains is n S , and the number of azimuthal Fourier modes is n F .  In Figure 4, we plot the N u-Ra data from all laboratory and numerical experiments. Hollow symbols represent laboratory data; red-filled symbols represent numerical simulation results. Square symbols connote a 5 cm high fluid layer, corresponding to a Γ 4 tank geometry; diamond-shaped symbols connote a 10 cm high fluid layer, corresponding to a Γ 2 tank geometry; circular symbols connote a 20 cm high fluid layer, correspond to a Γ 1 tank geometry.
The combined laboratory-numerical data set, over the range 2 × 10 4 Ra 3 × 10 7 , is best fit by the power-law: The data match the predicted heat transfer law (30) describing inertial convection in low P r fluid, N u = 1/2(RaP r) (1/4) . In fact, this fit agrees with the theoretical prediction to within 14% of the value of the coefficient and to within ±0.4% of the predicted value of the exponent. This quantitative agreement with theory demonstrates that inertial convection transfers heat via the predicted 1/4 scaling law in the range of parameters investigated here. This finding appears theoretically sensible based on the values of Re c and P e c accessed in our experimental range. For instance, inertial flywheel flows should be expected to develop as 860 Re c 3.6 × 10 4 in our RBC experiments. In contrast, the convective Peclet number range is far more moderate, covering 20 P e c 830. In Γ = 1/2 RBC experiments in liquid mercury, Glazier et al. (1999) [41] found that the 1/4 scaling law broke down to a 2/7 scaling at higher Ra values, typically exceeding ∼ 10 9 .
We also find good quantitative agreement between the local temperature measurements in the laboratory and numerical experiments. For example, Figure 5a shows the non-dimensional temperature time series, Θ, acquired on thermistors T 3 (black line) and T 5 (grey line) in the RBC 2 l laboratory experiment, here with temperatures normalized by ∆T . The corresponding Θ time series are also plotted on the point probes P 3 (red line) and P 5 (green line) from the corresponding numerical simulation RBC 2 n . There is good agreement between mean values at locations 3 and 5. The low frequency behaviors of the time series are also comparable between the two system. High frequency temporal variations in the laboratory data are not present in the numerical simulations with variances in the temperature measurements differing by up to 50%.
The strong agreement we find between all the laboratory-numerical internal mean temperature measurements implies that the same large-scale flow structures develop in our numerical simulations and in each associated laboratory experiment. Thus, the numerical simulations provide a means to indirectly visualize the large-scale flows that develop in the optically opaque liquid gallium laboratory experiments.
The agreements of laboratory-numerical point-by-point temperature measurements, as well as global heat transfer measurements, serve as a benchmark and provide cross-validation of the two approaches. Figure 5b further underscores the laboratory-numerical cross-validation. In this panel, we plot the time-averaged, non-dimensional temperature on each thermistor (i.e., Θ(T i ) t for the i th thermistor), normalized by the time-averaged, non-dimensional temperature of the corresponding numerical point-probe ( Θ(P i ) t ). The values of this ratio are close to unity (within ±5%) on all five thermistor/point-probe pairs, in each of the four RBC comparison cases. For such pointwise agreement to occur, the large-scale thermal and flow structures must be similar across all our laboratory-numerical comparison cases, further validating our coupled approach. We find that inertial flows develop in all four numerical simulations. The convective flows reach speeds near the free-fall estimate. This is shown in Figure 6a, which plots the typical maximum velocity amplitude, |U z |, measured on the mid-plane (z = 0) and averaged over time t for the four RBC SFEMaNS-T simulations ( Table 4). In all the simulations, the ratio of |U z | and U f falls well within 10% of unity, verifying that these RBC experiments are in the inertially-dominated convection regime. Because P r 1, these strongly inertial flows exist at only moderate P e c , and thus might be expected to be inefficient at advecting temperature anomalies across the fluid bulk. Instead, large-scale, large-amplitude thermal anomalies are generated in the bulk (cf. [40,50]). Figure 6b shows the time-averaged, maximum horizontal thermal anomaly on the mid-plane, |δζ|, normalized by the imposed vertical temperature difference ∆Θ = 1, for each of the four RBC SFEMaNS-T simulations. The value of this ratio remains fixed near to 0.6, a significant fraction of the imposed vertical temperature difference. Figure 7 shows the results from snapshots of numerical simulations carried out in Γ = 4 fluid layer geometries, with case RBC 1 n results displayed in the left hand column and case RBC 2 n results displayed in the right hand column. The top row of images, Figure 7a,c, show snapshots of the temperature anomaly field on a meridional slice through the fluid layer. Red color contours represent anomalously warm fluid, while blue contours represent anomalously cool fluid. In addition, three-dimensional streamlines in the vicinity of this meridional plane are also rendered. The color of the streamlines demarcates the local vertical flow velocities, with red denoting upward flow, and blue denoting downward flow.  The streamlines in Figure 7b show that three inertial convection cells span the fluid layer in case RBC 1 n . Upward flow velocities correspond to warm fluid being advected towards the top boundary; downward flow velocities advect cold fluid downward. This azimuthal wavenumber m 3 thermal and flow pattern is in good agreement with the temperature profiles shown in Figure 7a. We find the horizontal temperature differences on the mid-plane of the fluid reach values of 60% of the imposed vertical temperature difference, verifying the |δζ| value for the RBC 1 n case shown in Figure 6b. The oblique view in Figure 7c shows that the convection pattern is not dominantly axisymmetric in this case. The scatter in the dash ied snapshot profiles in Figure 7a is produced by precession and low-frequency variations of the non-axisymmetric convection pattern displayed in Figure 7c.
In sharp contrast to case RBC 1 n , the more strongly forced case RBC 2 n develops a predominantly axisymmetric flow pattern with a large-scale upwelling along the centerline (r = 0) of the fluid layer and downward flows along the container sidewalls. As shown in Figure 7d, this large-scale axisymmetric, inertial circulation (with power dominantly at m = 0 and subdominantly at m = 4) generates similarly strong horizontal temperature gradients on the mid-plane of the fluid, again, approaching 60% of the imposed vertical temperature gradient. The decrease in dominant horizontal mode number from m = 3 in case RBC 1 n to the axisymmetric mode in case RBC 2 n is reminiscent of the inferred change in mode number in the Γ = 4.5 liquid sodium convection study of [53].   Figure 8a,b shows that a complex, non-axisymmetric flow structure, dominantly characterized by a mode number m = 2 pattern, develops in case RBC 3 n . Case RBC 4 n , which has the highest Ra value and relatively low aspect ratio of Γ = 2, yields a classical large-scale circulation (LSC), with a single non-axisymmetric m = 1 inertial cell dominating the fluid domain [55][56][57]. The meridional slice in Figure 8c shows upwelling, warm fluid on one side of the container and downwelling, cool fluid on the opposite side. The oblique, snapshot view in Figure 8d shows that the LSC is spatially well-constrained, circulating about a relatively well-defined axis. This m = 1 flow pattern generates an approximately monotonic mid-plane temperature profile, and, yet again, generates maximum mid-plane temperature anomalies of |δζ| 0.6 ( Figure 6b).
The large-scale structures forming in our experiments are produced by a so-called inverse cascade process, in which turbulent convection excites small-scale modes that drive larger-scale inertial modes [49,58]. These inertial modes can approach the size of the fluid domain (Figure 8c). We have carried out lower resolution numerical test cases that do not contain small-scale modes and find that these cases fail to generate large-scale structures. This implies that upscale fluxes are important to the formation of large-scale structures in our convection experiments. Such cascades can also occur in sufficiently turbulent convection in moderate P r fluids, including water [16] and air [59]. However, these moderate P r LSC flows generate relatively large P e c values and, thus, tend to thermally mix and isothermalize the bulk fluid. These P r 1 LSC flows, then, do not generate such strong mid plane temperature anomalies, as are found here as well as in other comparable low P r convection studies [40,50].
Our strongly inertial low P r flows generate rather moderate values of P e c , and thus are expected to be inefficient at generating small-scale temperature anomalies in the fluid. In contrast, though, we find that large-scale flywheels are extremely well-suited at generating significant large-scale temperature anomalies throughout the convecting layer of liquid metal. Heuristically, these large-scale flows advect the large-scale temperature field with little small-scale thermal mixing effects. Low P e c inertial flywheels act, effectively, to rotate the imposed background temperature gradient around the flywheel's horizontal rotation axis, thereby generating large thermal anomalies in the bulk which exceed those found in lower thermal conductivity (higher P r) fluids.
The strong thermal anomalies that arise in liquid metal RBC experiments produce remarkably clean temperature time series. Figure 9 shows temperature time series recorded on the central thermistor T 3 in the left hand column and associated Fourier spectra in the right hand column. The top row, labeled a, shows the results from the laboratory RBC  l (c,d) laboratory experiments, carried out in the Γ 2 tank, and at Ra = 5.6 × 10 5 and Ra = 2.8 × 10 6 , respectively. Data collection starts once statistically-steady conditions have been reached. The spectral peak at f /f LSC 1.5 in panel b was first predicted using the RBC 3 n SFEMaNS case results that showed that the flow spent significant time in the n = 2 two-roll state.
In the Fourier spectra plots, frequency is normalized by the predicted overturn frequency of a container scale (n = 1) LSC, which we estimate in the following way. The typical time for a parcel, traveling at near the free-fall velocity, to make a vertical traverse of the the container is τ f = H/U f . Roughly maintaining that speed, the time for a parcel to circulate around a given large-scale convection cell is where n characterizes the number of convection cells spanning the tank diameter D = ΓH. The inertial overturn frequency of a cell at a given n value is then The results of associated numerical case RBC 4 n predicts that the dominant frequency in the RBC 4 l case should be characterized by that of an (n = 1) LSC. In contrast, the RBC 3 n simulation vacillated between the two roll, n = 2 state shown in Figure 8a,b, and a single roll n = 1 LSC state. Thus, these Γ = 2 numerical models predict that the spectrum in the RBC 4 l case should have a single dominant peak at f /f LSC 1, whereas the spectrum in the RBC 3 l should contain power at f /f LSC 1 for the n = 1 LSC flow as well as at f /f LSC 3/2 for the n = 2 state.
The peaks of the Fourier spectra in Figure 9b,d agree well with predictions based on our visualizations of the corresponding SFEMaNS-T cases RBC 3 n and RBC 4 n shown in Figure 8. The fact that the SFEMaNS-T results successfully predict the laboratory thermistor Fourier spectral content likely provides the strongest cross-validation of our coupled laboratory-numerical modeling environment.
Further, it implies that we can interpret strong peaks in laboratory Fourier spectra in terms of the morphological structure of the flywheels that develop in inertially-dominated convection cases.
In the RBC cases presented here, we have shown broad agreement between theoretical predictions, laboratory measurements and numerical modeling results. However, the fixed strength of |δζ| 0.6 found here is particularly surprising. One-dimensional scaling theories fail to predict a fixed value for the large-scale temperature anomalies on the mid-plane. For example, modifying the low P r scaling arguments of [60] to make use of equation (30), we predict that the mid-plane anomalies will decrease approximately as |δζ| ∼ (RaP r) −1/4 ∼ N u −1 . However, this behavior is clearly not verified by our SFEMaNS-T simulation results, nor those of [40,50]. This |δζ| 0.6 result must, however, break down in the limits of either small P r → 0, where the temperature field must asymptote back to the diffusive profile, or large P e c , where the bulk fluid becomes isothermalized. Advanced theoretical descriptions of low P r, moderate P e c heat transfer dynamics should aim to describe this apparently unique feature of Rayleigh-Bénard convection in liquid metals.

Essential Theory
We have carried out a laboratory-numerical suite of non-rotating, magnetoconvection (MC) experiments, using vertically-imposed magnetic fields with Chandrasekhar number Q 3.5 × 10 3 (see Figure 10 and Tables 5 and 6 The governing equations for the MC system are subject to the solenoidal velocity and magnetic field conditions. In the numerical MC simulations (Table 6), Equation (35) through (37) are solved, as described in detail in [24,25]. To be clear, the simplifications discussed below are presented to elucidate the physics, but we do not employ them in the present numerical simulations. To understand the effects of the magnetic field, we first estimate the convective magnetic Reynolds number, Rm c . At the highest Ra reached in our experiments, Rm c 10 −2 . For Rm c 1, the imposed magnetic field is not significantly induced or advected by the flows in our experiments. The imposed field in a Rm c 1 fluid remains effectively invariant; induced fields are negligible and the magnetic diffusion time-scale drops out of the problem. In this low Rm c limit, it is not necessary then to solve equation (37). Instead, it is possible estimate the Lorentz forces by calculating the induced currents using Ampere's law in an electrically conductive fluid, where V is the electrical potential and σ = 1/(µ o η) is the electrical conductivity. In the "quasi-static approximation", electrical potential effects are neglected in estimating J. Formally, this approximation holds only in 2D and in axisymmetric problems (cf. [34]), and, thus, does should not accurately apply to fully-3D magnetoconvection. However, as we show below, it predicts the MC regime boundaries well in our low Rm c experiments.   Table 6. SFEMaNS-T magnetoconvection (MC) simulation parameters. The Chandrasekhar, Prandtl and magnetic Prandtl numbers have fixed values of Q = 3.5 × 10 3 , P r = 0.025 and P m = 10 −4 , respectively. Applying the quasi-static approximation, the Lorentz force becomes Under this approximation, the Lorentz force acts only on u ⊥ , the horizontal velocity component perpendicular to B oẑ . Further, the Lorentz force acts to oppose these motions, providing a magnetic drag on the horizontal flows that acts with a characteristic drag time-scale τ drag = (σB 2 o /ρ) −1 . Since the low Rm c Lorentz force acts like an anisotropic magnetic viscosity, e.g., [19,34,60], in the limit of high Q, the onset of convection is no longer sensitive to the fluid's molecular viscosity ν. Instead, the critical Rayleigh number for magnetoconvection in an asymptotically strong vertical magnetic field is found to grow in proportion to the strength of the Lorentz force [34,38]: This predicts that MC will onset at Ra Q crit 3.45 × 10 4 in our Q = 3.5 × 10 3 experiments. The strong damping of horizontal flows acts to decrease the horizontal scale of overturning convection cells at the onset of magnetoconvection. Thus, for Q 1, the plane layer theory predicts that the onset width of MC cellular flows, non-dimensionalized using the fluid layer depth H, is [38] Thus, the horizontal length scale of overturning convective motions (i.e., half of a traditional wavelength), is expected to be Q 0.43 in our experiments. Using this estimate derived for an infinite horizontal layer, we estimate that roughly n Γ/ Q 10 steady cellular structures will span a diameter in our Γ 4 cylindrical experiments near convective onset. Similarly, we estimate an azimuthal wavenumber m = πΓ/(2 Q ) 15.
Once magnetoconvection begins to occur, the convective heat transfer will be hampered by damping effects of the magnetic fields. In order then to provide a rough prediction of the transition between the Lorentz-dominated MC regime and the undamped, inertial-dominated regime, we non-dimensionalize equation (39) and substitute it into equation (35) in place of N c (∇×B)×B. This leads to the quasi-static equations for the MC system: The non-dimensional parameter N is the interaction parameter in the quasi-static, low Rm c limit and is the ratio of the free-fall and Lorentz drag time scales, The value of N decreases below unity in our experiments in cases with Ra 3 × 10 5 . These low Rm c arguments predict that the Lorentz forces will restrict the convective motions and the associated convective heat transfer at Ra 3 × 10 5 , whereas at higher Ra, the MC flows will be effectively undamped and RBC-style inertial heat transfer will occur.
A number of MC investigations have considered low Rm c heat transfer scaling laws in the magnetically-controlled N 1 regime [19,[60][61][62]. With data at only one value of Q, we will put forward qualitative descriptions of MC heat transfer here. More detailed arguments will require data covering broader ranges of both Ra and Q.

Experimental Results
A total of 21 MC experiments have been carried out at Q 3.5 × 10 3 , 16 laboratory and 5 numerical cases (see Tables 5 and 6). The laboratory and numerical experimental set-up are nearly identical, with the largest difference existing between the value of the laboratory magnetic Prandtl number, P m 1.5 × 10 −6 , and the numerical magnetic Prandtl number, P m = 10 −4 . The N u-Ra-Q data is plotted in Figure 10, following the same conventions as in Figure 3. In all but the highest Ra MC cases, the heat transfer is suppressed relative to RBC cases at comparable Ra values. In fact, convective heat transfer is completely stifled such that N u 1 for Rayleigh number values that are below the predicted critical value Ra Q crit 3.5 × 10 4 . For Ra > Ra Q crit , the heat transfer rises sharply with Ra, but with a negative curvature to the N u-Ra trend (cf. [60]). The behavior of the MC heat transfer data merges with that of the best fit RBC trend at around Ra ∼ 10 6 . This merging occurs relatively close to the predicted N = 1 transition at Ra 3 × 10 5 . Figure 11 is constructed identically to Figure 5. Figure 11a shows local time series measurements made on thermistors T 3 and T 5 from case MC 2 l as well as on the associated measurements from temperature point probes P 3 and P 5 from case MC 2 n . The time series show similar mean values for laboratory and numerical measurements, with higher variance in the laboratory thermistor data. Long period vacillations exist in the temperature signals that greatly exceed any time scales found in the RBC data in Figure 5. These vacillations are due to changes in the slowly evolving and re-organizing planform of these quasi-steady, strongly magnetically damped cases.   Figure 11b, the ratio of these mean temperatures is close to unity on all five probes, to within ±5%. This implies, similar to our RBC experiments, that our MC simulations are generating large-scale temperature fields that agree well with the laboratory realizations. This further strengthens our argument that we can use the SFEMaNS-T code as a tool to visualize the large-scale flows that develop in our optically-opaque liquid gallium experiments. Figure 12a plots the time-averaged maximum mid-plane velocity amplitude |U z | measured in the SFEMaNS-T MC calculations, and normalized by the free fall velocity U f . In the RBC experiments, all the |U z |/U f values are close to unity. In contrast, here we see that MC 2 n , carried out at N = 2.18, yields |U z |/U f 0.40. This velocity adequately agrees with the scaling velocity predicted by balancing buoyancy and Lorentz drag terms U/U f ∼ N −1 0.46. In contrast, the velocities in the two high Ra cases, MC 3 n and RBC 4 n , cannot be estimated with this Lorentz damping argument since the |U z | values are close to U f . They are too strongly inertial for this two term balance to be applicable.  Figure 12b shows the time-averaged, maximum horizontal thermal anomaly on the mid plane, |δζ|, normalized by the imposed vertical temperature difference, ∆Θ = 1, in each of the three MC SFEMaNS-T simulations that are supercritical (Ra > Ra Q crit ; N u > 1). The value of |δζ| in the MC 2 n case is reduced by roughly 1/3 relative to the MC 3 n and MC 4 n cases. In contrast, the inertial flows that develop in the MC 3 n and MC 4 n cases both produce thermal anomalies in good agreement with the RBC results, |δζ| 0.6. Figure 13, constructed in parallel to Figure 7, shows visualizations from MC simulations carried out in Γ = 4 fluid layers, with case MC 2 n results displayed in the left hand column and case MC 4 n results displayed in the right hand column.
The velocity structures in Figure 13b shows that a spoke-like pattern of convection cells develop in case M C 2 n . The azimuthal wavenumber measured near the outer radius of the container is roughly m 5. This value is 1/3 of our estimate based upon planar linear MC theory, suggesting that the finite cylindrical container significantly effects the mode selection process. Case MC 4 n , in contrast, develops a predominantly axisymmetric flow pattern with a large-scale upwelling along the centerline and downward flows along the container sidewalls. This flow is similar to that of RBC 2 n in Figure 7d. However, in this MC case, the flow is more strongly axisymmetrized. Although the large-scale flows that develop in MC 4 n and RBC 2 n are roughly comparable, the smoother smaller-scale flow fields in MC 4 n , suggest that magnetic damping effects are still dynamically relevant at N = 1.30. Figure 14, constructed parallel to Figure 9, shows laboratory temperature time series from central thermistor T 3 in the left hand column and the corresponding Fourier spectra in the right hand column. The top row (a,b) shows results from MC 1 l ; the middle row (c,d) shows results from MC 2 l ; the bottom row (e,f) shows results from the case defined in Table 5 with Ra = 2.96 × 10 5 and N = 1.04. In the MC 1 l and MC 2 l cases, both with N > 2, the Fourier spectra contain power at frequencies far lower than the inertial frequency, f LSC , estimated using Equation (34) with n = 1. In contrast, in the N = 1.04 case, the spectrum strongly peaks at f f LSC . Based on the comparison of RBC4 l and RBC 4 n results presented here (cf. Figures 8c and 9d), we infer here that the strong spectral peak at f /f LSC 1 in Figure 14f corresponds to the existence of a container-scale LSC with convection velocities approaching the free-fall value. Thus, the regime transition from magnetically-damped MC flows to inertially-dominated MC flows occurs at N 1 in our low Rm c experiments, in agreement with the predictions made under the quasi-static approximation.  Table 5). The spectral plots show that the magnetically damped cases with N > 2 have little power at f = f LSC . However, in the Ra = 2.96 × 10 5 case, the magnetic damping no longer dominates and the spectrum peaks at the inertial frequency f = f LSC .

Essential Theory
The third canonical liquid metal convection system considered here is one in which there is a background rotation and no magnetic field (see Figure 15 and Tables 7 and 8). In comparison to numerous studies of RC in moderate Prandtl number fluids [11,13,14,20,23,[63][64][65][66][67][68][69][70], only a rather limited number of rotating convection investigations have been made in low Prandtl number fluids [18][19][20]51,71]. The governing equations for rotating convection are along with the solenoidal velocity field condition ∇ · u = 0. Centrifugal effects are not taken into account here, as the external Froude number F r = Ω 2 D/2g 4 × 10 −3 1 in our laboratory RC experiments (cf. [72]).  i.e., the angular velocity Ω = 0.84 rad/s and the rotation period P Ω = 7.5 s. The height of the layer is fixed at H = 10 cm. (The Γ 1 data shown in Figure 15 is from [20].)  In comparing the MC and RC equation sets, we note that N = (Re c /Q) −1 ∼ (buoyancy/Lorentz) −1 is the inverse magnetic analog to the convective Rossby number Ro c = Re c E ∼ buoyancy/Coriolis. However, significant dynamical differences exist between the two systems. Magnetic fields act to damp turbulent motions, directly removing kinetic energy from the flow field. In contrast, rotation acts to dynamically constrain the fluid motions, as discussed below, but without directly injecting or removing energy from the system. Thus, inertial effects are expected to be more prevalent in liquid metal RC systems than in comparable MC systems.
In the limit of small Ro c (even for large Re c [11]), the pressure term primarily balances the Coriolis force in equation (45). The leading order momentum balance then becomes Taking the curl of the above equation leads to an equation for the vorticity of the fluid, ω = ∇ × u, in the limit of rapid rotation: ∂u ∂z = 0.
Expression (48), called the Taylor-Proudman theorem (TPT), requires rapidly rotating fluid motions (E 1; Ro c 1) to be two-dimensional with no variation in fluid velocities along the direction of the rotation axis. The rapid system rotation imparts a massive vorticity to every fluid element, 2Ω. In the limit of weak inertia, buoyancy and viscous effects, there exists no way to strongly torque on a fluid element. Thus, any fluid motions must be z-invariant in order to induce no strong components of vorticity in planes perpendicular to rotation vector Ωẑ, e.g., [30,73].
In order to initiate convection in a rapidly rotating fluid layer, TPT must be overcome. In fluids with P r > 0.68, convection onsets via steady motions [38]. (This is true of the RBC and MC as well, independent of P r.) In order to overcome this dynamical constraint, steady rotating convection (S) onsets via thin columns that are of the axial height of the fluid layer and have a characteristic horizontal length scale of [74] S = 2π 4 1/6 E 1/3 H 2.4E 1/3 H.
a result verified in a number of studies [11][12][13]70,75,76]. It is the small horizontal scale of the columnar motions that allows for viscous forces that are capable of balancing the Coriolis force, and thereby relaxing TPT, e.g., [75,77]. For E 1, the steady form of rotating convective is predicted to onset in an infinite plane layer at a critical Rayleigh number of [38,74,78] Thus, convection becomes harder to initiate the more the fluid is constrained by the effects of rapid rotation, with the critical temperature gradient increasing as ∆T crit ∝ Ω 4/3 . In liquid metals, inertial effects are greatly enhanced relative to those arising in convection in higher P r fluids, as demonstrated by the strongly inertial flows that dominate our liquid gallium RBC experiments. This suggests an alternative way to break TPT. The Coriolis force can supply the restoring force that supports so-called inertial oscillations and waves in rotating fluids [79][80][81][82]. These inertial flows allow for an oscillatory-style of rotating columnar-style convection, which is typically more easily excited in low P r fluids than the steady form of convection. In the low E limit and P r < 0.68, the critical Rayleigh number for oscillatory rotating convection (O) in an infinite plane layer is [38,74] Ra O crit = 3 2π 4 P r 4 1 + P r Thus, oscillatory convection onsets before steady convection with Ra O crit 2P r 4/3 Ra S crit as P r → 0. The horizontal onset length-scale in the oscillatory case is O = 2π 4 1/6 1 + P r P r which simplifies to O 2.4(E/P r) 1/3 H for P r 1.
These oscillatory rotating convection structures are predicted to be wider than the steady onset structures by P r −1/3 . Stability analysis in a plane layer predicts then that rotating convection will onset in liquid metals via relatively wide oscillatory column-like convection cells. Narrower scale, steady convection structures may develop once the Rayleigh number exceeds Ra S crit . In our liquid gallium, Γ 2, E = 2 × 10 −5 RC experiments, oscillatory convection is predicted to onset with a horizontal scale of O H 2.25 cm, which are approximately three times the predicted width, S H 0.75 cm, of steady onset structures. Further, this O estimate suggests that n 8 oscillatory structures will span a typical tank diameter.
The infinite layer, E 1 prediction for the critical oscillation frequency, f O , normalized by the system's rotation frequency, f Ω = Ω/2π, is [38,74] In the low P r limit, this yields  [83] for rotating convection in a finite cylinder with non-slip boundary conditions. In a finite container yet another mode of rotating convection can develop in which the convective modes are attached to the sidewall and precess in the direction opposite to that of the system rotation, e.g., [83][84][85][86][87]. These sidewall modes exponentially decay inward towards the center of the tank. For E 1, the azimuthal onset length scale of the wall modes (W ) varies approximately in proportion to the aspect ratio of the container, W ∼ ΓH, and the retrograde precession frequency scales as f W /f Ω ∼ E/P r, which differs from the oscillatory convective frequencies as f W ∼ (E/P r) 2/3 f O . Further, critical Rayleigh number estimates for wall mode convection are not a function of P r [85]: In moderate P r fluids and for E 1, wall modes develop before steady convection modes fill the bulk fluid since Ra W crit /Ra S crit 3.7E 1/3 . In liquid metals, however, wall modes are not necessarily the preferred onset behavior for rotating convection. Comparing (51) and (56)  in low P r cases for which P r E 1/4 . Thus, in low P r fluids in which P r E 1/4 , it is expected that oscillatory convective motions will fill the fluid bulk before wall modes develop. In the suite of RC experiments presented below, the typical value of the Prandtl number is 0.025 and the typical Ekman number value is 2 × 10 −5 . Since P r is roughly 1/3 of E 1/4 in these experiments, we predict that RC will onset as oscillatory convection at Ra O crit 2.3 × 10 5 . Wall modes are likely to develop above Ra W crit 1.6 × 10 6 , and steady rotating convection is predicted to develop at Ra S crit 1.6 × 10 7 . Figure 15 shows laboratory heat transfer data for rotating convection at E 2 × 10 −5 in liquid gallium. The convection onsets very near the predicted critical Rayleigh number for bulk oscillatory convection. The heat transfer rises relatively shallowly over the next decade in Ra. Qualitatively similar shallow heat transfer behavior is found in non-rotating RBC in liquid metals in the vicinity of onset [52,88]. This shallow heat transfer scaling differs qualitatively from the MC results, which show a rather steep heat transfer scaling just beyond convective onset. Steep heat transfer scalings at onset have been found in recent studies of P r 1 rotating convection, e.g., [13,23,66,70].

Experimental Results
The critical Ra value for the onset of wall modes is reached at 10Ra O crit . We do not detect any wall modes directly in our laboratory experiments, likely because our interior thermistors are not located near enough to the sidewalls. However, there is a sharp increase in the N u-Ra scaling at around this point. It is possible that the inclusion of the wall modes helps to destabilize the flow field, leading to stronger interior mixing and more efficient heat transfer. However, our present laboratory data provides no definitive way to test this assertion.
Plane layer theory predicts that steady RC will develop in our experiments at Ra 10 2 Ra O crit . However, the RC heat transfer intersects, and merges with, the RBC best fit trend at nearly the same Ra value, similarly to [20]. This is unfortunate: it is not possible at E = 2 × 10 −5 to determine the convection behavior of liquid metal RC in the geophysically-relevant regime where Ra > Ra S crit and N u N u RBC . Figure 16. (a) Temperature time series acquired on thermistors T 3 in the RC 1 l laboratory experiment (with temperature normalized by ∆T ) and the corresponding time series of the non-dimensional temperatures measured on the point probe P 3 in numerical simulation RC 1 n ; (b) corresponding Fourier spectra. Figure 16a shows local time series measurements made on thermistors T 3 from case RC 2 l as well as the associated measurements on the temperature point probe P 3 from case RC 2 n . The time series show excellent agreement between laboratory and numerical measurements, both featuring comparable oscillations in the thermal field. The corresponding Fourier spectra, shown in Figure 16b, shows quantitative agreement in the peak power. This demonstrates that our combined laboratory-numerical tool is also well-suited for the investigation of rotating convective behaviors in liquid metals. Figure 17 shows visualizations from numerical cases RC 2 n (a-c) and RC 3 n (d,e). Figure 17a,b show streamlines that are color-coded by the local z-vorticity component of the flow. Figure 17c shows close-ups of the boundary layer flow fields as well as an axial rendering of the streamlines of an oscillatory convection column. Figure 17d,e show that the columnar flows are not stable in the Ro c = 0.198, RC 3 n case. The flow field is complex in this case, but Figure 17d suggests that a significant component of the flow appears to exist in the form of an n = 2 double roll large-scale flow, with maximum equatorial plane velocities of 0.88 U f . Figure 18 shows central thermistor T 3 data from cases RC 1 l , RC 2 l and RC 3 l on successive rows. The left hand panels show temperature time series, with time normalized by the system's rotation period, P Ω = 2π/Ω. The right hand panels show Fourier spectra, with the frequency normalized by the systems's rotation frequency, f Ω = 1/P Ω . The black dashed vertical lines estimate the inertial frequency of an n = 1 LSC normalized by the rotation frequency which yields f LSC /f Ω 2.2Ro c for the Γ = 1.87 laboratory cases shown in this figure. The green dashed lines denote the onset frequency, f /f Ω = 0.38, for oscillatory columns in plane layers [38] and finite cylinders [83] subject to no-slip mechanical boundary conditions. The time series data in Figure 18a,c show that the temperature signals are oscillatory in RC 1 l and RC 2 l . The RC 1 l spectrum peaks near f /f Ω 0.3, which is roughly 20% below the predicted onset frequency for planar oscillatory convection subject to no-slip conditions and 33% below the onset prediction for free-slip conditions. The RC 2 l spectrum has significant power at f /f Ω 0.3, but contains even more power at f /f Ω 0.38, in good agreement with the estimates for rotating convection in the presence of no-slip boundaries [38,83].
The RC 3 l case is qualitatively different. The n = 1 free-fall frequency is f LSC /f Ω = 0.43, which exceeds the no-slip oscillatory onset frequency prediction, f /f Ω = 0.38. Thus, the characteristic free-fall time scale is shorter than the oscillatory convection time scale, implying that the rotational dynamics should be subdominant to non-rotating inertial dynamics. If we assume an n = 2 flow somewhat akin to Figure 17d, we can use the more general version of Equation (57), f n /f Ω = (2π n Ro c )/(n + Γ) to predict an n = 2 free-fall frequency of 0.42. No dominant spectral peak was found at f /f Ω = 0.42. However, the broad band spectrum still implies that the oscillatory columnar flow regime has broken down at Ro c 0.1.

Discussion
Present-day models of planetary core (and stellar convection zone) dynamics are not yet capable of simulating the material properties of metalized fluids [89]. Thus, it is important to study convection in realistic P r < 1, P m 1 analog fluids, as has been done here. This allows us to correctly capture the canonical behaviors of a convecting volume of liquid metal, subject either to strong inertial, magnetic or rotational effects. The results presented here also provide an opportunity for modelers to test whether their codes can accurately simulate liquid metal convection dynamics. Such a benchmarking exercise may be the first step towards building realistic liquid metal geophysical and astrophysical convective dynamo models.

Present Findings
In general, we have found excellent agreement between theory, laboratory experiments and numerical simulations. In RBC experiments, an inertial heat transfer scaling regime developed in our laboratory-numerical experiments. Large-scale flywheel flows generated velocities that reached within 10% of the free-fall velocity prediction. These flows also generate significant horizontal temperature gradients in the fluid interior that appear to saturate at approximately 60% of the vertically imposed temperature gradient, a result unexplained by present theoretical models. In addition, we have successfully used the morphologies of the inertial flywheels found in our numerical RBC simulations to predict the dominant peaks in the laboratory RBC temperature spectra. The predictive capabilities of our numerical models as well as the high degree of coherence between numerical and laboratory experiments gives us great confidence in our results.
Laboratory-numerical agreement was found again in our magnetoconvection experiments. We have shown that our low magnetic Reynolds number MC experimental results are well described by the predictions of quasi-static theory. Heat transfer and flow velocities were found to be magnetically-damped for interaction parameter values N 1. Further, at N 1, the thermal spectra become strongly peaked at the overturn frequency of a single large-scale structure, with flow occurring near the free-fall velocity estimate U f .
The suite of rotating liquid metal convection experiments produced the most complex results. Rotating convection occurred in the form of oscillatory columns at Ra 2.5Ra O crit (i.e., RC 2 l and RC 2 n ). These inertial columns become unstable by Ro c 0.1. In low P r fluids such as metals, the range over which these columns are stable comprises less than a decade in Ra at E = 2 × 10 −5 , cf. [13,70].

Future Directions
In future efforts, it will be important to use fluids with the correct material properties, such as those of metals. In addition, for applications to rather extreme geophysical and astrophysical dynamo systems, it will also be important to carry out experiments both at Ra Ra crit and at Ra Ra T rans , where Ra T rans denotes the transition Rayleigh number where the system reverts to non-magnetic, non-rotating RBC convection behavior. This will allow for the study of the geophysically and astrophysically relevant regimes of high N and low Ro c flows that are also strongly turbulent such that Re c 1. In our MC experiments, convection was found to onsets at Ra Q crit = π 2 Q in excellent agreement with theory. The convection was found to transition back to the effectively non-magnetic RBC behavior at roughly N = 1, which corresponds to a transition Ra value of Ra Q T rans Q 2 P r. Then the supercritical range of Ra that is available to investigate turbulent MC dynamics is of order Ra Q T rans /Ra Q crit Q P r/π 2 which has a value close to 10 in our experiments. This is likely too small a range to determine robust turbulent scalings. In fact, it is unlikely that turbulent flow can become fully-developed over such a limited range of supercriticality (even in P r 1 metals). This implies that significantly higher Q experiments are necessary to study the geophysically and astrophysically relevant regime of high N (and high N c ) flows that are also strongly turbulent such that Re c is large.
In the next generation of MC laboratory-numerical simulations, we challenge modelers to attempt to reach Chandrasekhar numbers of Q ∼ 10 7 . For this Q value, the critical Rayleigh number is Ra Q crit 10 8 and the transitional Rayleigh number is Ra Q T rans 2.5 × 10 12 , allowing for over 4 orders of magnitude in Ra-space between onset and transition. Carrying out experiments at Ra 10 10 would allow investigation of MC with N = 15.8 and Re c 6 × 10 5 . It will be of great interest to determine whether such flows are laminarized by the effects of the magnetic field or whether turbulent magnetohydrodynamic cascades are able to develop [90,91], as are likely to be relevant in geo-and astrophysical dynamo systems [15,92]. In which has a value of close to 50 in our RC experiments. Although this seems an adequate range of supercriticality, we are unable to investigate, at Ro c 1, how the dynamics change once steady rotating convection modes develop (Ra Ra S crit ). This is problematical, as it appears that Ro c 1 and Ra Ra S crit in planetary [93] -and possibly stellar [94]-interiors. Following the MC considerations above, we propose a future liquid metal RC study at E = 10 −7 . At this E value, oscillatory convection will onset at Ra O crit 3 × 10 8 and the transitional value is approximated at Ra Ω trans ∼ 5 × 10 11 . Carrying out RC experiments at Ra 10 10 will yield results at Ro c ∼ 0.063 (= 1/15.8). Based on the numerical findings of [14,70] in moderate P r fluids, these parameters may allow for the development of inverse turbulent cascades of convective energy, and the formation of large-scale vortical flows in low Ro c liquid metal flows.
Lastly, magnetic and rotational forces can be comparable in strength in geo-and astrophysical settings [5,15,95]. Thus, future studies must also seek to provide a detailed understanding of convection-driven, liquid metal rotating magnetoturbulence for Re c 1 and N −1 ∼ Ro c 1. This is indeed a rich topic, one for which our coupled laboratory-numerical experimental platform is well suited and for which our present results provide a solid baseline.

Author Contributions
A.R. and J.M.A. designed research; A.R. and G.F. performed laboratory research; A.R. and J.L.G. carried out numerical research; A.R., G.F. and J.M.A. analyzed data; A.R. and J.M.A. wrote the paper.

Appendix 1: SFEMaNS-T Performance Results
SFEMaNS makes use of Fourier series in the azimuthal direction and finite elements in the meridional plane (see Section 2.2, and [25] for a detailed discussion of the code). Any 3D axisymmetric domain is decomposed azimuthally into n F Fourier modes, which are solved independently on n F cores. Each Fourier mode is projected onto a meridional finite element grid containing n p nodes. Scaling information is reported in [25]. The results of [25] also showed, for a given computing cluster, that the parallel efficiency was limited by inter-core communication issues. Thus, an additional spatial domain decomposition has been implemented using the METIS libraries [27] to address this limitation. In this 2D domain decomposition, the meridional grid mesh is split into n S subdomains, such that n p /n S is the number of nodes of each subdomain processed by the cluster. On Stampede, the optimal value of n p /n S is order 10 4 . Figure A1 schematically illustrates this double parallelization. Figure A1. Schematic of the double parallelization used in the SFEMaNS-T code. The solution is represented by a Fourier expansion with n F modes in azimuth. The meridian grid mesh is approximated via finite elements, and is decomposed into n S subdomains. The total number of cores employed in a given case is n c = n F × n S . This schematic shows the RBC 4 n domain geometry, but for visual clarity only 1/50th of the meridional resolution is represented. Figure A2 shows the results of a weak parallel scaling test, made using the the most extreme thermal convection case performed in this study, RBC 4 n . This case uses a meridianal mesh made up of n p = 46, 565 nodes. Using the TACC Stampede cluster, we determined that the communication time was optimal when the n p nodes are decomposed into n S = 4 subdomains. In this scaling test, we have fixed the number of subdomains at n S = 4, and varied the number of Fourier modes from n F = 16 to 512, such that the total number of cores n c = n F × n S ranges from 16 to 2048. The parallel efficiency is defined here as the ratio of the computing time for the reference case, τ (n c = 64), and the computing time for each case performed in the scaling test, τ (n c ). In all cases, we have carried out 1000 time steps. The computing time did not include the time for initialization and output operations. This weak scaling tests shows that SFEMaNS-T maintains an efficiency of 82% up to 2048 cores. Figure A2. Weak scaling test of the RBC 4 n case (Table 4) using a fixed meridional grid mesh composed of n p = 46, 565 nodes, and varying number of Fourier modes, n F . The parallel efficiency is measured by taking the ratio of the computing time required to perform 1000 iterations for the referential case, n c = n F ×n S = 64 cores, and the computing time required to perform 1000 iterations for each n c case performed in the test.

Conflicts of Interest
The authors declare no conflict of interest.