{'fraction_non_alphanumeric': 0.02960830312077913, 'fraction_numerical': 0.01945451529584607, 'mean_word_length': 4.693288364249579, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Compositional Convection in the Deep Interior of Uranus

Uranus and Neptune share properties that are distinct from the other giant planets in the solar system, but they are also distinct from one another, particularly in their relative internal heat flux. Not only does Neptune emit about ten times the amount of heat that emitted by Uranus, the relative amount of emitted heat to the energy they absorb from the sun also differs greatly, being comparable at Uranus and the largest of all giant planets at Neptune. As a result, it is questionable whether thermal convection occurs within the interior of Uranus. However, the presence of an intrinsic magnetic field implies that interior fluid motions must exist. Here, we consider compositional convection driven by the release of hydrogen associated with the formation of large organic networks or diamond precipitation in the deep interior. We test this hypotheses using a set of numerical rotating convection models where the convective driving is varied between thermal and compositional sources and is sufficiently vigorous to not be strongly constrained by rotation. In most cases, we find ice-giant-like zonal flows develop, with three bands characterized by a retrograde equatorial jet and prograde jets at higher latitudes. Large-scale circulation cells also develop and lead to heat and mass fluxes that tend to exhibit local maxima along the equatorial plane. This similarity between convective flows driven by thermal and compositional buoyancy therefore predict Uranus and Neptune to have similar interior dynamics despite Uranus' minimal internal heat flow and may thus explain why both ice giants have comparable magnetic fields.


INTRODUCTION
Uranus and Neptune share a number of similarities that have lead them to be considered twin planets, often referred to collectively as the ice giants, differentiating them from the gas giant planets, Jupiter and Saturn. Uranus and Neptune have similar masses (14.5 M ⊕ and 17.1 M ⊕ ), radii (4.0 R ⊕ and 3.9 R ⊕ ), rotational periods (1.012 × 10 −4 s −1 and 1.083 × 10 −4 s −1 ), atmospheric compositions (dominated by hydrogen/helium and enriched in heavier elements), three-jet zonal winds, non-dipole dominated magnetic fields, and measured gravitational moments (Jacobson 2009(Jacobson , 2014Williams 2015). However, despite these similarities, there are important differences between the two planets that may give insights into the evolution of planetary systems. For example, Uranus has a lower mean density than Neptune (1200 kg m −3 compared to 1600 kg m −3 ) and a uniquely extreme axial tilt of 97.8 • , compared to the more Earth-like tilt of 28.3 • of Neptune. Though both Uranus and Neptune have colder effective temperatures than the other giant planets, Neptune emits the greatest amount of heat relative to that received through insolation while Uranus has the lowest among the giant planets (Pearl et al. 1990;Pearl & Conrath 1991). As will be described below, we hypothesize that compositional convection-rather than thermal-in Uranus' deep interior may explain the dynamical similarities of the ice giant planets.
In this paper, we first discuss the background and scientific rationale for studying different drivers of deep convection in ice giant interiors. Section 2 further discusses the dynamical regimes of convective motions that may be expected in giant planet interiors. Section 3 then details the numerical model of rotating convection in a spherical shell used to test thermal and/or compositional buoyancy as the origin of flows in the deep interior of Uranus. The velocity and thermodynamic fields found in our numerical experiments are described in Section 4, and the underlying processes are discussed in Section 5. Section 6 extrapolates these results to ice giant interiors.
1.1. Internal structure, composition, and heat flow The distribution of mass in planetary interiors is measured by gravitational moments, which are non-unique with respect to both the radial dependence of density as well as composition and temperature. Traditionally, internal structures of Uranus and Neptune are modeled with three distinct layers: a solid rocky core, a mantle consisting primarily of molecular ices (e.g., water, methane, and ammonia), and an atmosphere consisting of molecular hydrogen and helium and other gases such as methane and ammonia (e.g., Hubbard et al. 1991;Podolak et al. 1991;Nettelmann et al. 2013). The ratio of rock to ice remains an open question, and the measured gravitational moments, J 2 and J 4 , may be explained by either an ice-dominated or a rock-dominated composition (e.g., Teanby et al. 2020). There may also be a continuous variation of density and concentration of heavy elements through the interior rather than distinct phase transitions separating layers (Fortney et al. 2010;Helled et al. 2010a;Movshovitz & Fortney 2020;Vazan & Helled 2020).
Nonetheless, the interiors of Uranus and Neptune are thought to contain a mixture of materials containing carbon, nitrogen, and oxygen (Atreya et al. 2020). Hydrocarbons, primarily methane, are present in the atmospheres of the ice giants, forming both clouds and hazes in the upper atmosphere, which obscure deeper regions of the atmosphere, making deep mixing ratios of carbon to hydrogen difficult to measure (Lindal et al. 1987;Baines et al. 1995). Measurements of the abundances of atmospheric methane show that both planets are substantially enriched in carbon, approaching 100 times solar ratios of carbon-to-hydrogen (Atreya et al. 2020). At high pressures and temperatures, methane has been shown to dissociate into diamond (carbon) and hydrogen, rather than forming a superionic solid (Ross 1981;Kraus et al. 2017). Simulations indicate that mixtures of methane, ammonia, and water also enter such a dissociative state (Chau et al. 2011;Naumova et al. 2020).
While atmospheric abundances of methane and other carbon-rich molecules can be measured through ground based observations, the atmospheric abundances of ammonia and water poorly known (Mousis et al. 2020). Water has been studied under planetary interior conditions both experimentally and theoretically, predicting that both electrically conductive ionic fluids and superionic solids may occur within the ice giants (French et al. 2016;Millot et al. 2018). Ammonia and mixtures of ammonia and water exhibit similar states of matter (Cavazzoni et al. 1999;Ravasio et al. 2021;Bethkenhagen et al. 2017a). Figure 1 shows the pressure-temperature diagram for several models of the interiors of Uranus and Neptune and relevant chemical states that may occur in their interiors. "Synthetic Uranus" approximates the interior composition by a ratio of hydrogen to carbon to nitrogen to oxygen (H:C:N:O) of 28:4:1:7, which has been shown to exist in different phases of matter and to become electrically conductive with increasing pressure and temperature (Nellis et al. 1988;Chau et al. 2011;Guarguaglini et al. 2019). Ab initio simulations predict that synthetic Uranus enters what is referred to as a reticulating phase, where C-C and C-N bonds become energetically favorable and form networks within the material (Chau et al. 2011). In a similar mixture of water, ammonia, and methane, carbon may demix from the solution to form diamonds (Naumova et al. 2020). In both cases, hydrogen may be released as a product of the reactions. This release of light components may then be a potential source of compositional buoyancy, similar to the release of light elements as Earth's inner core solidifies.
Besides composition, the thermal structures of the ice giants are similarly not well constrained, which are reflected by the inclusion of both adiabatic and non-adiabatic interior models in Figure 1 (Nettelmann et al. 2013;Podolak et al. 1991;Stixrude et al. 2020). Magnetic fields were discovered at Uranus and Neptune by the Voyager 2 mission (Ness et al. 1986(Ness et al. , 1989. The existence of planetary magnetic fields can be explained by dynamo theory, which states that a planetary magnetic field is sustained against diffusion by induction caused by complex motion within an electrically conductive fluid in the planet's deeper interior (Lühr et al. 2018). For this reason, we can infer the existence of convective motion in, at least a part of, the interiors of the ice giants. However, adiabatic models of Uranus and Neptune that assume fully convective interiors throughout their history are unable to correctly predict the thermal properties of either planet in the present, suggesting that the planetary interiors are only partially convective or that their thermal behaviors are time-dependent (Scheibe et al. 2019  Temperature, K U1 (Nettelmann et al. 2013) N1 (Nettelmann et al. 2013) N2b (Nettelmann et al. 2013) Icy Uranus (Bethkenhagen et al. 2017) Frozen core (Stixrude et al. 2020) Conduction, cold (Podolak et al 2019) Conduction, hot (Podolak et al 2019) Reticulating phase (Chau et al. 2011) Diamond forming (Naumova et al. 2021) No diamond forming (Naumova et al. 2021) Figure 1. Pressure-temperature diagrams for both materials thought to be present in the ice giant interiors as well as specific internal structure models of Uranus and Neptune that feature both fully convective models (Nettelmann et al. 2013;Bethkenhagen et al. 2017a) and models that include significant non-adiabatic regions (Podolak et al. 2019;Stixrude et al. 2020). The green shaded area represents the chemically active reticulating phase of synthetic Uranus, though the boundaries of this region are approximate (Chau et al. 2011). The region fades with increasing temperature, as the extent in temperature of this region of phase space are also unknown. Diamonds represent the states of a mixture consisting of water, methane, and ammonia in a ratio of 56.5:32.5:11, where filled symbols denote where demixing into diamonds occurs and empty symbols denote where this demixing does not occur (Naumova et al. 2020).
Thermal properties of the giant planets are summarized in Table 1. The majority of energy emitted by Neptune originates from the interior, compared to only about 6% for Uranus (Pearl et al. 1990;Pearl & Conrath 1991). The uncertainty in Voyager 2 measurements is quite large, and it is therefore possible that no internal heat contribution was measured at Uranus. The thermal properties of planets have traditionally been assumed to be constant in time over decadal timescales; however, measurements made by Cassini of Saturn revealed a measurable decrease in emitted power from the period of 2005 to 2009, which have raised questions as to the temporal variation of thermal properties of the other giant planets (Li et al. 2010). Furthermore, the instruments aboard the Cassini spacecraft were able to make measurements of Jupiter that revealed that more heat is emitted than had been previously measured by the Pioneer or Voyager missions (Li et al. 2018). In contrast, the temporal variability of Uranus' and Neptune's thermal properties are unknown because subsequent measurements have not yet been made.
The low luminosity of Uranus may be accounted for in several possible ways, such as a continuous compositional gradient within the interior or the existence of a thermal boundary layer that could prevent convective transport of heat to the surface (Podolak et al. 1991;Nettelmann et al. 2016;Vazan & Helled 2019). Condensation of heavier molecules, such as methane or water, out of hydrogen dominated atmospheres has also been suggested as a mechanism that may inhibit thermal convection (Guillot 1995;Cavalié et al. 2017;Friedson & Gonzales 2017;Leconte et al. 2017;Markham & Stevenson 2021). If convection is inhibited in both Uranus and Neptune, but condensation only occurs in the atmosphere of Neptune, then latent heat associated with this process may warm the atmosphere and prevent it from cooling to a temperature similar to Uranus (Kurosaki & Ikoma 2017). Additionally, the difference in heating may be caused by demixing of water and hydrogen in the deep interior that occurs in Neptune, but not in Uranus (Bailey & Stevenson 2021). Regions of the deep interior may experience heat transport primarily by conduction, which also predicts a smaller amount of heat reaching the surface (Podolak et al. 2019). The low luminosity of Uranus is also consistent with a frozen-core, which is able to sequester primordial heat as the freezing process continues by inhibiting convection throughout the deep interior of the planet (Stixrude et al. 2020). Here, we hypothesize that Uranus' deep Table 1. Thermal properties and estimates of non-dimensional parameters of the giant planets. Values for Jupiter are reported in Li et al. (2018), for Uranus in Pearl et al. (1990), and for Neptune in Pearl & Conrath (1991). Approximate values of viscosity and thermal diffusivity are based on first principles calculations in French et al. (2012) and French & Nettelmann (2019). Uranus and Neptune are assumed to consist of water, therefore other properties are estimated using the equation of state presented in Mazevet et al. (2019). Thermal Rayleigh numbers are estimated following Soderlund (2019 interior is unstable convection that is driven either partly or entirely by compositional buoyancy because of the light element release at depth and likely non-adiabaticity of the planet's interior. This could differentiate Uranus from Neptune, whose significant internal heat flux suggests thermally driven convection.

Velocity fields
As shown in Figure 2, Uranus and Neptune have surface winds that are characterized by three zonal (east-west) jets, with retrograde (westward) winds at the equator and prograde (eastward) winds at higher latitudes (Mousis et al. 2018). In addition to physical units, velocities are expressed in terms of the he Rossby number, Ro, the ratio of inertial to Coriolis forces: where u is velocity, Ω is rate of rotation of system, and r o is outer radius. The Rossby number is dependent on length scale, and we choose r o so that these Ro values serve as a lower bound; this length scale will be revisited in later sections. Equatorial winds on Neptune are more intense, reaching 400 m s −1 compared to ∼ 100 m s −1 on Uranus; the extent of these winds is also wider on Neptune, extending between ±20 • on Uranus compared to ±50 • on Neptune. However, these measurements are calculated on the basis of the rotational periods measured by Voyager 2, and modified rotational periods can result in similar winds of ∼ 200 m s −1 at both Uranus and Neptune (Helled et al. 2010b). The prograde winds tend to be similar in intensity, ∼250 m s −1 on Uranus and ∼300 m s −1 on Neptune, but these prograde bands tend to be broader on Uranus than on Neptune. It is unknown how deep these winds extend into the interior, although gravity measurements and Ohmic dissipation estimates suggest they are confined to about the top 1-2 thousand km in both planets (Kaspi et al. 2013;Soyuer et al. 2020). Zonal winds may develop from a number of possible physical processes that typically rely on being either driven by insolation or by internal heat. Jupiter-like or ice giant-like winds can develop in a thin weather layer through a turbulent inverse energy cascade (Cho & Polvani 1996) or moist convection (Lian & Showman 2009). Deep-convection can similarly generate zonal flows that have a prograde equatorial jet with alternating jets at higher latitudes (i.e. jupiter-like) in systems strongly constrained by rotation or three jets with retrograde equatorial flow (i.e. ice giant-like) when the system is not rotationally dominated (e.g., Busse 1976;Aurnou et al. 2007;Soderlund et al. 2013;Gastine et al. 2013b). While these latter studies have focused on thermal buoyancy, we hypothesize that compositional convection in the deep interior may also be a viable mechanism to generate Uranus' zonal flows.

REGIMES OF ROTATING CONVECTION
The properties of convective motion are dependent on the relative influence of rotation and the physical characteristics of the fluid, and varying these parameters can result in different regimes of rotating convection (Gastine et al. 2016;Cheng et al. 2018Cheng et al. , 2020. The relative importance of rotation in the system can be quantified using the Ekman number, Winds are expressed both in terms of physical units (m s −1 ) and non-dimensionally in terms of the Rossby number, Ro = u/roΩ. Measurements were made by Voyager 2 and Earth-based observations (Lindal et al. 1987;Limaye & Sromovsky 1991;Sromovsky et al. 2015Sromovsky et al. , 2001. Empirical fits are provided by Sromovsky et al. (2015) and French et al. (1998) for Uranus and Neptune, respectively. Grey lines mark the equator (dash-dot) and the zero wind velocity (dashed).
while the properties of thermal transport can be quantified using the Prandtl and Rayleigh numbers: Here, D is length scale, ν is kinematic viscosity, κ is thermal diffusivity, α is the coefficient of thermal expansion, g is gravitational acceleration, and ∆T is the superadiabatic temperature contrast that drives convection. Convective heat transfer efficiency of the system is characterized by the ratio of total to conductive heat flux, known as the Nusselt number: where q is heat flux, ρ is density, and c p is specific heat capacity. Trends in N u behavior as a function of Ra, P r, and Ek are often used to identify different convective regimes. Cheng et al. (2018) describe six regimes, which are defined by the strength of convective forcing: sub-critical, columnar, plumes, geostrophic turbulence, unbalanced boundary layer, and non-rotating heat transfer. In rotating systems, the onset of convection occurs once the Rayleigh number reaches a critical value, Ra C ∼ Ek −4/3 . As the convective forcing is increased, the dynamics of the system and properties of the heat transfer will be determined by the influence of rotation. When buoyant forcing is weak, rotation is able to impose structure in the fluid that is constant on surfaces parallel to the axis of rotation. As the convective forcing is increased, the organization imposed by rotation begins to break down in the so-called unbalanced boundary layer regime and then ultimately the system behaves as if it is non-rotating.
The thermal convective regime for a planetary system can then be predicted using estimates for Ek, P r, and Ra. While Ek and P r are relatively straight forward, the Rayleigh number and ∆T in particular is more challenging. Following Soderlund (2019), the superadiabatic temperature contrast can be obtained by using N u − Ra scaling laws for rapidly rotating (∆T RR ) and non-rotating heat transfer (∆T N R ) as upper and lower bounds for a given heat flux q: ∆T RR = 2.1 Ω 4 q 2 κD ρ 2 c 2 p να 3 g 3 1/5 .
Table 1 contains estimates of the non-dimensional parameters that describe the physical properties of the giant planets, assuming a water-dominated interior composition. At high pressures and temperatures relevant to the ice giant interiors, water is predicted to have a kinematic viscosity of ν ∼ 3 × 10 −8 to 10 −7 m 2 s −1 and Prandtl numbers between 0.1 and 1 (French & Nettelmann 2019). If we allow the length scale to be between 30% and 70% the planetary radius (Helled et al. 2010a;Nettelmann et al. 2013), the Ekman numbers are between 10 −18 and 10 −17 . The thermal expansion coefficient can be estimated as α ∼ 10 −4 K −1 , and the value of ρc p is estimated to be between 8.7 × 10 6 and 2.6 × 10 7 J m 3 K −1 (French & Nettelmann 2019). The rates of rotation and gravitational field are assumed to match planetary values of Ω ∼ 10 −4 s −1 and g ∼ 10 m s −2 . Estimates for the properties of the interior of Jupiter are estimated based on the equations of state of hydrogen-helium mixtures, given a viscosity of ν ∼ ×10 −7 to 10 −6 m 2 s −1 and a length scale between 30% and 50% the radius of Jupiter results in an Ekman number between 10 −19 and 10 −18 (French et al. 2012;Wahl et al. 2017). To estimate the Rayleigh number, we let ρc p ≈ 3 × 10 8 J m 3 K −1 and allow the Prandtl number to be between 0.01 and 0.1 ). Using these properties, the range of internal heat fluxes from Table 1, and the temperature contrasts from equations 6 and 7, we estimate Rayleigh numbers of Ra ∼ 10 27 to 10 31 for Neptune based on the mean value of heat in Table 1 and Ra 10 31 for Uranus taking the possibility of zero internal heat source into account. Estimates of the Rayleigh and Ekman numbers for the giant planets are shown in the context of these regimes of rotating convection in Figure 3b. Uranus and Neptune fit within the weakly-rotating unbalanced boundary layer regime and the transitional regime between unbalanced boundary layers and rotationally dominated convection. In contrast, Jupiter appears to be confined mostly to the transitional regime between rotationally dominated and weakly rotating regimes. When the Rayleigh numbers are calculated on the basis of the mean heat flux, Jupiter is confined to this transitional region, and only enters into the unbalanced boundary layer regime as larger values of the heat flux within one standard error of the mean are considered. However, these only take thermal convection into account, and any compositional contributions to convection in these planets is unconstrained by observations, meaning that either Neptune or Uranus may have a higher total Rayleigh number, which would lift them further towards the non-rotating regime. This lack of a strong rotational constraint will guide the choice of input parameters for the models presented in this work. Choices for our models are shown in Figure 3a. At high Ekman number, these regimes are less well differentiated, so we chose parameters in the non-rotating regime to reflect this weak rotational constraint.

Model
In this study, we investigate the effect of varying the type of convective driving on the resulting dynamics of the system. Given the significantly different thermal behaviors displayed between Uranus and Neptune, we will consider Ra GT U ∼ Ek −8/5 P r 3/5 Ra C ∼ Ek −4/3 Dashed: P r = 1 Dotted: Sc = 10 Dash-dot: P r = 0.1 Dash-dot-dot: P r = 0.01 Figure 3. Convective regime diagram for a) the models in this work and b) estimates for the giant planets found in Table 1, represented by the shaded boxes for Uranus (U), Neptune (N), and Jupiter (J). In a), the Rayleigh number of the double diffusive case is expressed as the total Rayleigh number. The grey region indicates Rayleigh numbers that cannot sustain convection at a given Ekman number (RaC ). The purple lines and shaded bands represent transitions between rotationally dominated regimes of convection and the unbalanced boundary layer regime (RaGT U ), and red lines and shaded bands represent transition between unbalanced boundary layers and non-rotating convection (RaUNR).  cases that are purely thermal, purely compositional, and driven by both thermal and compositional buoyancy (i.e. double-diffusive convection, albeit with both unstable gradients). Figure 4 shows schematics of the internal structures assumed in our models. Convection is assumed to occur in the ionic ocean and gaseous envelope, while the solid core that lies underneath is unable to move on a timescale relevant to dynamics of the system. While the existence of magnetic fields in the ice giants implies convective motion, using scaling law arguments, the Lorentz force is estimated to be several orders of magnitude weaker than the Coriolis force in ice giant interiors (Soderlund et al. 2015). As a result, the magnetic fields are not expected to modify the convective flows significantly, which justifies our use of hydrodynamic models.
We model thermo-compositional convection of a Boussinesq fluid confined within a spherical shell that rotates at a constant rate about theẑ-axis. The fluid's behavior is governed by a system of dimensionless equations that describe the time evolution of velocity u, temperature T , and concentration of light components dissolved in the fluid ξ under the assumption of incompressibility: The gravitational field g is assumed to be proportional to the radius: g = g o r/r o , where g o is gravitational acceleration at the outer boundary. The equations are non-dimensionalized using shell thickness D as the characteristic length scale of the system and viscous diffusion time τ ν = D 2 /ν as the characteristic time scale, while temperature and the concentration are expressed in terms of the contrast between the inner and outer boundaries, ∆T and ∆ξ, respectively. In addition to the Ekman, Prandtl and Rayleigh number, the system is defined by three other non-dimensional input parameters: Schmidt number Sc, compositional Rayleigh number Ra ξ , and shell geometry χ: Inner boundary radius Outer boundary radius .
Here, κ ξ is diffusivity of light components, and r i and r o are the inner and outer shell radii, respectively. We assume a Boussinesq fluid, where fluid properties are assumed constant except for density anomalies associated with temperature and composition perturbations. Density is equal to a mean value plus a time varying part, and governed by the equation of state: where α ξ is the coefficient of compositional expansion. Both inner and outer boundaries are stress-free, meaning that viscous stresses vanish at these boundaries as appropriate for the nearly inviscid interiors of giant planets. The values of temperature and composition are fixed at both boundaries, with T i = ξ i = 1 and/or T o = ξ o = 0, which allows both of these thermodynamic properties to drive convection.
As shown in Table 2, the Ekman numbers in our models are set to Ek = 10 −3 (cases T1, D, and C1) and Ek = 3 × 10 −4 (cases T2 and C2); while these values are large compared to planetary estimates, they reflect our available computational resources given the strong degree of supercriticality and this assumption will be revisited in Section 6. In thermal convection cases, the Prandtl number is chosen to be equal to one following French & Nettelmann (2019) who predict a range of Prandtl numbers between P r ∼ 0.1 − 10 for ice giant interiors, with most values being roughly equal to 1. In geophysical and planetary environments, the ratio of heat and mass diffusivity can be very small (e.g. Monville et al. 2019), but we choose Sc = 10 so that the diffusivities differ by at least one order of magnitude while keeping in mind the limitations of computational resources. Finally, we choose a value for χ = 0.35, which is consistent with the models of Soderlund et al. (2013) and motivated by the transition to chemically active phases near ∼ 100 GPa (e.g., χ ≈ 0.35 for a uniform density ofρ ≈ 1000 kg m −3 ).
The Rayleigh numbers were chosen so that our models would have convective turbulence that is three-dimensional with limited vertical stiffness as argued in Section 2. This is achieved by considering the most extreme scaling for the regime transition between the unbalanced boundary layer-non-rotating heat transfer regime identified by Cheng et al. (2018): Ra U N R ∼ Ek −2 P r, which can be rewritten in terms of the convective Rossby number: Ro C ∼ 1. This parameter represents the ratio of buoyancy to Coriolis forces: Convection is weakly constrained by rotation when Ro C 1 and the dynamics are rotationally dominated when Ro C 1 (Gastine et al. 2016). In our cases, we will choose Ro C to be equal to √ 2 following Soderlund et al. (2013). In the double diffusive case, the Rayleigh numbers are chosen so that the coefficient of the thermal and compositional buoyancy term in Equation 8 are equal. The resulting Rayleigh numbers range in supercriticality from Ra/Ra C ∼ 20 for case T1 to nearly Ra ξ /Ra C ∼ 400 for case C2 (Table 2), where Ra C are the critical Rayleigh numbers for the onset of convection calculated using the Spherical INertia-Gravity Eigenmodes (SINGE) 1 code (Vidal & Schaeffer 2015).
We next consider these dimensionless input parameters in terms of the associated dimensional properties. The outer radius, shell thickness, and rotation rate are fixed to r o ≈ 2.5 × 10 7 m, D ≈ 1.6 × 10 7 m and Ω ≈ 10 −4 s −1 in all models. Given the choices of Ekman number, this implies a kinematic viscosity of ν = [3 × 10 7 , 8 × 10 6 ] m 2 s −1 for Ek = [10 −3 , 3×10 −4 ]. Setting P r = 1 implies that the thermal diffusivity is equal to the viscosity: κ = [3×10 7 , 8×10 6 ] m 2 s −1 for Ek = [10 −3 , 3 × 10 −4 ]. Conversely, Sc = 10 implies that the mass diffusivity is an order of magnitude smaller: κ ξ = [2.6 × 10 6 , 8 × 10 5 ] m 2 s −1 for Ek = [10 −3 , 3 × 10 −4 ]. If we assume that the fluid in the interior is water, then the thermal expansion coefficient is on the order of α ∼ 10 −4 K −1 ) and the temperature contrast is approximately 390 K in both thermal cases and 195 K in the double diffusive case. In order to estimate the composition contrast, we first calculate the density of a mixture of hydrogen and water using the linear mixing approximation (Bethkenhagen et al. 2017b). Interpreting ξ as the mass fraction of hydrogen, the density of the mixture is given by where ρ M is the density of the mixture, ρ H is the density of hydrogen and ρ W is the density of water, taken at identical pressure and temperature conditions. For small values of ξ, we can approximate the expansion coefficient as It is reasonable to assume that in the planetary interiors, light elements would infuse slowly into the convective region of the interior. For this reason, the concentration of the fluid most enriched in light elements near the inner boundary will be sufficiently small that this approximation is applicable and ∆ξ 1. Based on equations of state for water and hydrogen, we find that typical densities that coexist in hydrogen and water are ρ H ∼ 600 kg m −3 and ρ W ∼ 2, 500 kg m −3 , respectively Mazevet et al. 2019). This gives α ξ ∼ 3.2, which corresponds to ∆ξ ∼ 0.02 kg kg −1 for both compositional cases and ∆ξ ∼ 0.01 kg kg −1 for the double diffusive case.
Many of these physical values (e.g., diffusivities, thermal and compositional anomalies) are necessarily much larger than those expected for the ice giants. Considering kinematic viscosity as an example, ν is more than 10 orders of magnitude too large, although exploiting a larger turbulent eddy viscosity interpretation would reduce this disparity. However, given the limited resolution of Uranus and Neptune observations, we assume kinematic similarity exists between the large-scale flows of ice giants interiors and our models. It is also broadly instructive to understand how varying the source of convective driving effects the dynamics in rotating, spherical systems, and how these differences may manifest in nature. Table 2. Control and time-averaged derived parameters for each case. The control parameters are defined in Equations 2 -4, 12, and 13. Spatial resolution is given by the number of radial grid points Nr and maximum spectral harmonic degree nmax. The duration of time-averaging is given in terms of the viscous dissipation time scale. Three derived parameters that characterize the velocity field include: the Reynolds number, Re, the local Rossby number Ro l , and characteristic spherical harmonic degree, n, defined in Equation 19. The Nusselt number, N u, and Sherwood number, Sh, characterize the thermodynamic fields and are defined in Equation 5 and 20, respectively. Coefficients of variation (CV) are provided for the temperature, composition, and the axial component of absolute angular momentum which are calculated using Equation 24. Equations 8 -10 are solved numerically using the pseudospectral fluid dynamics code MagIC 2 compiled using the SHTns library 3 of spherical harmonic transforms to improve performance (Wicht 2002;Breuer et al. 2010;Schaeffer 2013). MagIC takes advantage of Equation 11 by decomposing the velocity into two scalar potentials: a poloidal potential and toroidal potential. Non-linear terms of the model equations are calculated on a grid in spherical coordinates (r, θ, φ) of N R × N φ /2 × N φ , where N R is the number of radial points and N φ is the number of longitudinal grid points. Radial grid points are distributed over Chebyshev nodes, which allows for better resolution near the inner and outer boundary in addition to allowing for efficient transforms between the grid and spectral representation. Latitudinal grid points placed on Gauss-Legendre quadrature nodes and longitudinal grid points are distributed uniformly. All quantities have a spectral representation, with angular dependence in terms of spherical harmonics, Y m n (θ, φ) and radial dependence in terms of Chebyshev polynomials, with the spectral representation being truncated at a maximum spherical harmonic degree n max and Chebyshev degree N R − 2. Table 2 contains the resolutions used each case. The higher Ekman number cases were initialized as solidly rotating bodies (u = 0), with perturbations to the temperature and/or composition seeded as white noise. The lower Ekman number cases were initialized using equilibrated states of the corresponding higher Ekman number cases. Once a statistically consistent state is achieved, the solutions are time-averaged over at least one viscous diffusion time to obtain mean quantity values. Figure 5 shows the relative percentages of the kinetic energy contained in different components of the flow. The axisymmetric poloidal and axisymmetric toroidal components of the kinetic energy correspond to the zonal flows and meridional circulations, respectively, while the non-axisymmetric components correspond to convective motions. The zonal flows tend to dominate the dynamics of these systems in the lower Ek cases (T2 and C2), while non-symmetric components dominate in the higher Ek cases (T1, D and, C1). The zonal flows, in addition to being axisymmetric, are mostly equatorially symmetric (see Figure 6). In all cases, meridional circulation represents less than 2% of the total

Convection
Meridional circulation Zonal flow Figure 5. Bar plot of the percentage of kinetic energy that is represented by axisymmetric zonal flows, axisymmetric meridional circulations, and the non-axisymmetric convection. The time-averaged value of the total non-dimensional kinetic energy for each case is given above the bar corresponding to that case.
kinetic energy, meaning that these systems may be interpreted as consisting primarily of zonal flows and convection. Comparing T2 and C2, while the absolute kinetic energy associated with convection is larger in T2 (2.1 × 10 5 in T1 compared to 9.8 × 10 4 in C2) , the fractional contribution of convection increases with the decrease in diffusivity. However, comparing T1 and C1, the same changes do not occur, with similar percentages for convective and zonal flow kinetic energy relative to the total despite different absolute kinetic energies. The zonal flow thus tends to decrease in total fraction of kinetic energy moving from T2 and C2, and moving from T1 to D, indicating that the amplitude of the zonal winds is influenced by the thermodynamics of the system. This trend does not continue when moving from D to C1, which behaves differently from the others in several regards as described further below. The characteristic spherical harmonic degree of the flow can be calculated from the kinetic energy spectra following Christensen & Aubert (2006) where u n is the velocity associated with all spherical harmonics of degree n. With the notable exception of C1, the cases tend to have similarn values for a given Ekman number, differing by about 2% between C2 and T2 and T1 and D. Case C1, however, is a clear outlier withn more than twice that of T1 and D. When the Rossby number (defined in Equation 1) is calculated on the basis of the length scale Dπ/n as described in Equation 19, it is known as the local Rossby number, and better reflects the ratio of inertia to the Coriolis force in the force balance of the system (Christensen & Aubert 2006). The local Rossby numbers, tabulated in Table 2, span the range 0.25 ≤ Ro ≤ 0.43 and are roughly similar in cases with the same mode of convective driving. The largest values occur in T2 and T1, while the smallest values appear in C2 and C1, with D being intermediate between the thermal and compositional cases. The Coriolis force is thus expected to exceed inertia on convective length scales, but not by a significant amount. Figure 6 shows the time-and azimuthally-averaged meridional circulations and zonal flows. All cases have threebanded zonal flows, with four of the five cases having qualitatively ice-giant-like flows along the outer boundary, i.e. a retrograde jet at low latitudes and prograde jets at higher latitudes. As the Ekman number is decreased for the thermal convection cases (T1 → T2), equivalent to decreasing the fluid viscosity, the retrograde flows intensify by 17%. In contrast, the zonal jets reverse directions between the two compositional convection cases (C1 → C2), with prograde equatorial flows and retrograde higher latitude jets developing in the case of C1. While the zonal flows tend to have relatively small variations along the axialẑ direction at large cylindrical distances, more significant vertical variations develop near and within the tangent cylinder. As convective driving is varied (T → C), there is a decrease in the flow intensity across the system, by 9% between T2 and C2 and by 120% between T1 and C1 (which is made especially large given the reversal). There is also a change in the characteristics of vertical variations of these flows, with a particularly intense prograde belt in the deep interior of T1 and T2, which does not appear in D or C2. In C1, in both the polar region (within the tangent cylinder) and far from the axis of rotation near the outer boundary, the zonal flow is constant along nearly vertical surfaces.
Meridional circulations contain comparable fractions of kinetic energy across all of the models and similarly all feature, to some extent, fluid upwelling from the lower boundary near the equator ( Figure 6). However, the circulation patterns also have notable differences between them. In each case, mid-latitude circulation cells are present, but the number and characteristics of the circulation cells change as both Ekman number and convective driving are varied. T2 has two circulation cells outside the tangent cylinder in each hemisphere, whereas single cells develop with faster radial velocities appear in the case of T1. Inside the tangent cylinder, additional circulation cells appear near the inner boundary as the Ekman number is decreased from T1 to T2. Cases T1 and D behave quite similarly, although flows are weaker in the latter. Compared with T2, C2 has two additional circulation cells that form outside of the tangent cylinder at the mid latitudes. Decreasing the Ekman number in the compositional convection cases (C1 → C2, the extent of the polar cell increases and an additional cell appears near the axis of rotation. In varying the convective driving from D to C1, this polar cell diminishes in its extent and the mid-latitude cell crosses into the tangent cylinder. Figure 7 shows slices in the equatorial plane of radial and zonal velocity and thermodynamic fields. Similar trends appear in the intensity of velocities as those seen in Figure 6. The flow becomes more disorganized in the lower Ekman number cases, where there are more small scale features in the components of the velocity field. In the cases with retrograde equatorial winds, one hemisphere is thermally/compositionally active, while the other is quiet, which leads to differences in the structure of the velocity field. Convective motions associated with the radial velocity are more intense and show more small scale features in the active hemisphere. As the fluid flows from the quiet to the active hemisphere, a region of subsumption forms, causing the colder (light element depleted) fluid to flow back toward the inner boundary. Along with this, the zonal flows in the active hemisphere become more disturbed, with a smoother transition from prograde to retrograde flow from the inner to outer boundary in the quiet hemisphere. In the lower Ekman cases, the equatorial jet wraps more completely around the circumference of the system.
In the case of C1, the properties of the flow are significantly different. Prograde zonal flows appear to encircle the entire system and the system no longer appears to be divided into an active and quiet hemisphere, with radial upwelling and downwelling occurring across all longitudes. For this reason, the zonal flows appear to be disturbed at all longitudes as well. Figure 8 plots typical three-dimensional isosurfaces of the axial component of vorticity for each of the cases, illustrating the three-dimensional nature of the velocity fields that has developed. When rotation dominates the system, the Taylor-Proudman theorem dictates that the properties of the fluid should be organized into columnar structures (i.e. invariant in z). In all cases except for C1, the flows that develop are not strongly organized by rotation as constructed a priori. The distribution of isosurfaces is similar to the features that develop in the temperature and composition fields in Figure 7c, f, i, l, and o. That is, in cases with active and quiet hemispheres, the vorticity is mostly concentrated in the active hemisphere. Again, C1 appears to be different, with features distributed more evenly across all longitudes and that tend to exhibit some alignment in the axial direction. Figure 9 plots meridional slices of time-and azimuthally-averaged temperature and/or composition for each case. Heat and mass are transported by advection and diffusion as described in Equations 9 and 10. Temperature (composition) profiles in which heat (light elements) is transported only by diffusion can be predicted using Laplace's equation, which will have a particular gradient and the property that these solutions are radially symmetric. In these cases, the presence of turbulent motion will act to improve the ability of heat and mass to become more well mixed in the bulk fluid, which would decrease the steepness of the interior temperature and/or concentration gradient compared to that predicted by diffusion alone. When controlled for viscosity, greater diffusivity of heat or mass leads to a steeper gradient near the boundaries, but a shallower gradient throughout the bulk of the fluid. When controlled for convective driving, though the effect is less pronounced, the gradient in the bulk of the fluid becomes shallower. Because this is

C1
-0.28 -0.21 -0.14 -0.07 0.00 0.07 0.14 0.21 0.28 equivalent to reducing the diffusivity of the thermodynamic field, it effectively causes the field to adhere less closely to the diffusive solution. In addition, the latitudinal dependence of this gradient reflects the differing dynamics within the convective region. The greatest degree of mixing occurs along the equatorial plane and mid-latitudes, whereas the gradient tends to become steeper near the polar regions. This results from the fluid motions in the equatorial plane, which tend to be more intense that those that form within the polar regions. Each case has equatorial upwellings, with multiple cells forming outside the tangent cylinder, which will mix the fluid. Most of the poloidal kinetic energy is equatorially symmetric, which implies a local extremum along the equatorial plane. The Nusselt number, N u, and its compositional analog, the Sherwood number, Sh, are the non-dimensional heat and compositional gradients across the outer boundary, respectively (Incropera et al. 2007). The Sherwood number is calculated similarly to the Nusselt number in Equation 5: whereξ is the mass flux of light elements, measured in kg m −2 s −1 . Figure 10a contains time series of these quantities, and time-averaged values are also tabulated in Table 2. They tend to be consistent across the duration of simulations, with most cases being distinct from one another. As the Ekman number is decreased for a fixed Ro c , the Nusselt Values of the isosurfaces are a) |ωz| = 1.14 × 10 4 for T2, b) |ωz| = 6.12 × 10 3 for C2, c) |ωz| = 2.28 × 10 3 for T1, d) |ωz| = 1.81 × 10 4 for D, and e) |ωz| = 1.28 × 10 3 for C1. The grey sphere in the interior represents the inner boundary, and the black circles indicate the equator and a great circle facing the viewer on the outer boundary.
and Sherwood numbers both increase due to the stronger convective motions that result from the lower viscosity and stronger driving force. Comparing T1 to C1 and T2 to C2, there is also an increase in the convective transport, as reducing the Prandtl or Schmidt number decreases the relative importance of diffusive heat or mass transport. While D has a comparable, yet smaller, Nusselt number relative to T1, it has a larger Sherwood number than C1; this is likely a consequence of the faster flows that develop in D compared to C1, resulting in more convective transport. The range of Nusselt and Sherwood numbers does not have a distinct trend with respect to the convective driving; however, because compositional fields are associated with a higher flux than the temperature fields, the relative variability decreases with the decrease in thermal (mass) diffusivity. Figure 10b plots heat and mass flux as a function of latitude. Local maxima occur near the equator in all cases except C1. The cases show two trends in the flux near the poles. T2, C1, and C2 have polar maxima, while D (both temperature and composition) and T1 do not, though additional local maxima do appear at higher latitudes. In all cases except C1 that has a global minimum at the equator, local minima appear at mid-to high-latitudes. Because these measure the relative importance of convective heat and mass transfer, they depend on the structure of the velocity field. The effect of decreasing viscosity tends to be opposite when controlled for convective driving. Comparing T2 to T1, the peak at the equator is suppressed, with the global maximum appearing at the poles. The opposite is true in comparing C2 to C1, where the variation with latitude is smaller in C2 as peak near the equator has been raised in comparison to C1. Similarly, the effects of varying convective driving do not have consistent effects, when controlled for viscosity. Furthermore, we find that the contrast in heat (mass) flux as a function of latitude becomes smaller (i.e. the curve is flatter) as Péclet number increases.

Angular momentum
Properties of the velocity field are related to angular momentum, and conservation of absolute angular momentum is important to determining the dynamics of the system. Absolute angular momentum is the total angular momentum in the non-rotating reference frame resulting from solid body rotation and differential rotation of the fluid. Radial transport of angular momentum outward by individual fluid parcels can result in decreasing zonal flow, if angular momentum for the parcel is conserved, in the same way that a spinning skater extending their arms is able to reduce their rate of rotation (Gilman & Foukal 1979). Retrograde equatorial winds are hypothesized to form when absolute momentum becomes mixed, either completely or partially, in the fluid contained within a spherical shell (Aurnou et al. 2007;Gastine et al. 2013a). Because the system rotates around theẑ-axis, we focus on the component of absolute angular momentum in this direction, h z . We normalize this quantity by the value of solid-body angular momentum at a distance of r o from the axis of rotation: where s = r sin θ is the perpendicular distance from the axis and Ro φ is the Rossby number of the zonal flow. Figure 11a plots h z as a function of latitude at the outer boundary for all models as well as end-member scenarios representing spatially homogeneous angular momentum and solid-body rotation. For homogeneous angular momentum, the profile is a vertical line. In regions where the curve is flatter, the gradient is less steep, and therefore we interpret these regions as being more well mixed than where the curve changes more rapidly with latitude (Aurnou et al. 2007;Soderlund et al. 2013). The thermal cases have quite efficient mixing of angular momentum in their equatorial region, whereas the corresponding compositional cases tend to have less efficient mixing, resulting in a more recurved shape. Figure 11b plots the corresponding zonal flows, including the zonal flow that results from homogenized angular momentum. In this homogenized case, we expect the formation of not only strong equatorial flows, but also nonphysical infinitely fast prograde flows approaching the poles. In the case of T1 and T2, relatively fast equatorial flows develop, as well as a parabolic shape to the curve around the equator. In all cases, the flows near the poles are damped out, resulting in local maxima (i.e. prograde bands) at higher latitudes. However, in the compositional and double-diffusive case, the equatorial flows are significantly diminished compared to the thermal cases. In C1, rather than seeing well mixed angular momentum, the prograde flow near the equator does not vary strongly with latitude. Thus, despite identical Ro C values, a range of angular momentum behaviors are observed, ranging from quite well mixed near the equator in T1 and T2 to poorly mixed and closer to solid body rotation in C1. These mixing behaviors are approached quantitatively in the next section.

Degree of mixing
A method for quantifying the quality of mixing in a system is the coefficient of variation (CV), which is equal to the standard deviation divided by the mean of a given quantity, which measures the intensity of segregation of the solute (Kukukova et al. 2009). The mean, µ ψ , and variance, σ 2 ψ , of a field ψ are volume weighted. In the case of a Boussinesq fluid, this choice is equivalent to mass weighted averages.
By construction, a CV of zero represents a field that has been fully homogenized. We are interested in the mixing of temperature, composition, and the axial component of absolute angular momentum. For each case, time-averaged CV quantities are obtained by averaging CV values calculated for snapshots of the three-dimensional velocity, temperature, and composition fields with the numbers of snapshots used given in Table 2. To give context to these results, we also calculate the CV for three special cases: the temperature or composition transported by diffusion alone, the angular momentum of a rotating solid body, and for a fluid that moves as a rotating solid body inside the tangent cylinder but is fully mixed outside the tangent cylinder (CV T C ).
In the case of heat transport by diffusion, the system reaches steady state when ∇ 2 T = 0 with boundary conditions T (r i ) = 1 and T (r o ) = 0 such that the temperature profile will only be a function of normalized radius for a given shell geometry: where r i = χ/(1 − χ) and r o = 1/(1 − χ) (with equivalent arguments for the compositional field). As a result, the mixing properties of these diffusive profiles are: (28) Figure 12a shows a plot of CV T for all possible values of χ. For χ = 0.35, the CV T is equal to 1.02, which establishes an upper bound on cases where fluid motions cause the heat/light elements to be mixed. Angular momentum of a rotating solid body can be obtained by setting u φ equal to zero in Equation 21, leading to As shown in Figure 12a, the CV ranges from 0.65 for the full sphere (χ = 0) to 0.45 as the spherical shell becomes infinitesimally thin (χ → 1). When χ = 0.35, CV is equal to 0.62. For relatively thick shells (χ 0.4), the CV deviates from the full sphere value proportionally to χ 3 : For thinner shells (χ 0.4), the deviation from the infinitely thin shell value is proportional to (χ − 1) 2 : As shown in Figure 12b, decreasing the Ekman number increases the efficiency of both heat (mass) and angular momentum mixing due to the increased intensity of the flows that mix these fields. This is supported by the general downward trend with increasing Reynolds number seen in Figure 12d. This trend is also true when controlled for convective driving by comparing T1 to T2 and C1 to C2.
Varying the convective driving has the opposite effect in angular momentum as it does in heat (mass), as shown in Figure 12c. The reduction of CV from the diffusive end-member value of 1.02 demonstrates the role that convective transport plays in homogenizing the fluids. Figure 12f shows that while there is not a distinct trend in the angular momentum mixing, increasing the Nusselt (Sherwood) number is associated with more efficient mixing of heat (mass). The Péclet number is the ratio of advective to diffusive transport and given by P e T = ReP r in the case of heat and P e ξ = ReSc in the case of light elements. As seen in Figure 12e, there is similarly no distinct trend in the angular momentum mixing, but an increase in the Reynolds (Péclet) number is associated with a decrease in the CV of heat (mass). Despite generally weaker flows in the compositional and double diffusive cases compared to the thermal cases, the relative importance of advective heat (mass) transport to diffusion increases with decrease in diffusivity of mass compared to heat. Lines connect markers from cases with similar convective driving (i.e. T1 to T2 and C1 to C2). In a), d) and f), square markers indicate thermal (P r = 1) cases, triangles indicate double diffusion, and starts indicate compositional cases (Sc = 10). In b) and e), circular makers indicate Ek = 10 −3 cases and diamonds indicate Ek = 3 × 10 −4 cases. Gold markers indicate absolute angular momentum, filled blue barkers indicate temperature, and empty blue markers composition. In b), e) and f), the same CV hz for the double diffusive case is plotted at both P r = 1 and Sc = 10; however, CVT is plotted only at P r = 1 and CV ξ is plotted only at Sc = 10. Horizontal lines correspond to heat (mass) transport by diffusion (blue dotted), solid body rotation (gold dashed), and tangent cylinder models (gold dashdot). c) Coefficients of variation as a function of χ, where the vertical grey line marks χ = 0.35. The diamonds marking the intersection of these curves with the vertical line represent the values of the horizontal lines in the other subfigures.
This can be understood through a feedback process, by which the reduction of diffusivity allows the system to achieve the same degree of heat (mass) mixing with less intense fluid motion, and that a more thoroughly mixed temperature (composition) field causes the influence of buoyancy on the dynamics of the system to become reduced. We can quantify this one way by considering the total non-dimensional buoyant power in each case (Yadav et al. 2013): Comparing T2 and C2, we find that T2 has a buoyant power of 1.9 × 10 9 compared to 5.8 × 10 8 in C2 (only 30% of T2). A similar decrease occurs from T1 to D to C1, which have buoyant powers of 8.9 × 10 7 , 5.2 × 10 7 and 2.2 × 10 7 , respectively. This means that as the convective driving is varied to less diffusive modes of convection, the buoyant force is contributing less to the energetics of the system, and therefore the strength of the flow is weakened. Despite this, efficient mixing of mass is maintained.

Force balance
The dominant force balance determines the evolution of the dynamics in the system. The force terms of Equation (8), f , can be expressed in terms of components in spherical coordinates, The root-mean-square (RMS) value as a function of spherical harmonic degree n may be time-and volume-averaged as where ∆t is the duration of the time-averaging. The global value may then be obtained by summing over all spherical harmonic degrees: Figure 13 plots of the time-averaged magnitude of forces in the system throughout the full shell. At small n, the most dominant forces are the pressure gradient and Coriolis in all cases; however, different forces become prominent at smaller length scales. In T1, the ageostrophic Coriolis term exceeds the magnitude of the Coriolis force at n 11, which indicates that a substantial amount of the Coriolis force is anti-parallel to the pressure gradient. As n → n max , the force balance becomes dominated by viscosity. In C1, there is still a leading order geostrophic balance at length scales larger than the characteristic wavenumber, with pressure and Coriolis force peaks at n = 2 and 4 that substantially exceed the n = 1 values in contrast to the other cases. As n → n max , while the viscosity again becomes the greatest force in the force balance, the compositional buoyancy remains relatively strong because smaller scale features of the composition field persist. Case D maintains properties similar to both T1 and C1, passing through configurations like both cases with increasing n. The strength of the compositional buoyancy, though it begins at small n with a low value, persists to higher spherical harmonic degrees, crossing over the thermal buoyancy at n ≈ 24. Viscosity crosses inertia at n ≈ 62. At this point, the dominant force balance resembles T1 more, but compositional buoyancy crosses inertia at n > 100.
In T2 and C2, there are peaks in pressure at n = 1, 3, and 4 and in the Coriolis force at n = 2 and 4. Above n = 5, pressure continues to dominate the force balance of T2. Inertia crosses over the Coriolis force at n = 10. The viscosity only eventually becomes comparable to the pressure and inertial forces as n → n max . In C2, inertia crosses over the Coriolis force at n ≈ 35, however as n increases, the Coriolis force eventually becomes the least significant force. At n ≈ 80, viscosity crosses over the pressure, where it becomes the dominant force in the system.
In each case, inertia crosses over the Coriolis force as n is increased. We can attempt to find where this crossing occurs by assuming that the Coriolis term is similar in strength to inertia at the length scale CI : which can be recast as Because our models tend to have Ro ∼ 0.1, this crossing is expected to occur at moderate values of n consistent with Figure 13. The crossover length can be obtained from the spectra by comparing the magnitude of the Coriolis force and inertia. Because we expect a linear relationship between the Rossby number and this length scale, we are able to obtain an empirical linear fit of the form CI = aRo between these quantities using the least-squares method. The fit we obtain is These values are plotted in Figure 13 as the vertical dashed lines. This prediction is more successful in the thermal cases, where these lines pass almost exactly through the crossover, than in the compositional and double diffusive cases. However, in the non-thermal cases, the crossover occur at degrees where other forces have become significant. This is especially true in C1, where all the forces are approximately equal to one another. The magnitude of the force terms in Equation 8 are dependent on the choice of characteristic time scale, which depends on the choice of viscosity. However, if instead we choose rotational time as the characteristic time, we are able to compare the magnitudes of the forces between each case, assuming that the rate of rotation is the same between them. This is equivalent to considering the dimensional forces to be normalized by DΩ 2 Values for the integrated forces, as the viscosity and convective driving is varied are tabulated in Table 3. First, the magnitudes of buoyancy are less effected by decreasing Ek than by varying the convective driving. Comparing the compositional cases to the thermal cases, the compositional buoyancy falls off less rapidly with increasing n than thermal buoyancy, but thermal buoyancy has a higher peak value than compositional buoyancy. This has the effect of reducing the integrated buoyant force in the compositional cases. In T2, the integrated buoyant force is 0.14, compared to 0.07 in C2. In T1, the integrated buoyant force is 0.15, compared to 0.09 in D and C1. When this is taken with the discussion of the buoyant power in Section 5.2, as the heat or mass become more thoroughly mixed, the influence of the buoyancy is reduced.
As the system undergoes weaker forcing, the velocity field decreases in intensity. Figure 14 contains plots of the globally integrated forces against the input parameters. The convective Rossby number is intended to be a proxy for the ratio between buoyancy and the Coriolis force, with Ro C ∼ 1 often delineating the transition between eastward and westward equatorial flows (Aurnou et al. 2007;Gastine et al. 2013b). As the Ekman number is decreased, the total buoyancy is roughly consistent to slightly decreased for a given mode of convection. However, for constant Ekman number, total buoyancy decreases moving from thermal to compositional convection. While there is no change between T1 and T2, but a notable decrease in the total Coriolis force between C1 and C2. Furthermore, in cases T1 and T2 and the thermal component of D, buoyancy is never greater than either the Coriolis force nor inertia at any spherical harmonic degree. In cases C1, C2, and D, the compositional buoyancy force becomes greater than Coriolis at large values of n. The ratios of the integrated buoyancy to Coriolis force are tabulated in Table 3. In the case of C1, the volume integrated buoyancy compared to Coriolis is the strongest among the different cases with a force ratio of ∼ 0.5 compared to 0.2 for the other cases. Comparing the properties of each system in Figure 7 shows that the dynamics of T1, T2, D, and C2 tend to be dominated by lower order terms, while the relative  importance of higher order terms is reflected in the spectrum for C1, and the peaks that form in the Coriolis, pressure and buoyancy terms at n = 2 and n = 4.

Gravitational moments
Processes occurring within planetary interiors are not directly observable and, therefore, their properties must be inferred from observable properties such as their gravitational fields. Anomalies in temperature and composition give rise to buoyancy in the system create anomalies in the density. These anomalies are not necessarily radially symmetric, like the hydrostatic background, meaning that the gravitational field that results from them may likewise not have a radial symmetry.
The gravitational field surrounding a planet is irrotational and can therefore be represented in terms of a scalar potential and multipole expansion. The strength of each term is determined by the distribution of mass inside the planet and quantified using the gravitational moment, J n , M r n o J n = 2π 0 π 0 ro 0 ρ(r, θ, φ)r n+2 sin θP n (cos θ)drdθdφ, where P n is the Legendre polynomial of degree n. When m = 0, the moments are called the zonal gravitational moments. These moments can be further divided into two components, arising from the flattening that is associated with solid body rotation and the dynamics of the system, ∆J n . Solid body rotation causes flattening that is symmetric with respect to the equatorial plane and, therefore, the gravitational moments associated with odd values of n can only result from the dynamics of the system. In contrast, even moments reflect both planetary structure and dynamics.
Here, we calculate the density anomaly associated with dynamics in our models, noting there are two effects that are not represented: the flattening associated with the rotation of the system and the increase in background density with depth. The first two even zonal gravitational moments, J 2 and J 4 , were measured in Uranus and Neptune by Voyager 2 (Stone & Miner 1986. Gravitational science has remained one of the primary of objectives for more recent scientific missions to the giant planets, with higher order gravitational moments being measured to J 10 in Jupiter by Juno and to J 12 in Saturn by Cassini-Huygens (Debras & Chabrier 2019;Galanti et al. 2019). Mission concepts to the ice giants typically seek to constrain the low degree gravitational moments up to at least J 6 (e.g., Hofstadter et al. 2017;Jarmak et al. 2020;Rymer et al. 2020).
Equation 15 relates the density anomaly to the temperature and composition fields, from which we can derive an equation for the gravitational moment resulting from the dynamics of the fluid: where ρ c is density of the core region, and f T = Ek 2 Ra/P r and f ξ = Ek 2 Ra ξ /Sc are the fractional contributions of thermal and compositional to buoyancy, respectively. In our calculations, we assume that ρ c ≈ 10ρ as an order-of-magnitude approximation, which is similar to the ratio of densities in other interior models (e.g., Hubbard et al. 1991;Nettelmann et al. 2013). The quantity R 3 Ω/GM is approximately 0.03 in ice giant planets. This estimate is obtained by considering reasonable choices of constant density for the fluid and core region to approximate the interior of Uranus, assuming that a distinct transition between these layers exists. For example, if the density of the fluid is approximately equal to 1000 kg m −3 , then the density of the core would be approximately 7300 kg m −3 so that the mean density of the entire system matches that of Uranus. Assuming that the density of the fluid is approximately 500 kg m −3 , the core density would be 18000 kg m −3 (36 times larger) so that the mean density matches Uranus. Decreasing the choice of density in the fluid would require a likely un-physically large density in the core while a fully uniform interior would result in a ratio of one, with a choice intermediate between these extremes more likely to approximate the physical reality of the planetary interiors. Figure 15 plots the time-averaged gravity spectrum up to spherical harmonic degree 10. Though we have ignored flattening due to rotation, gravitational moments associated with a uniform spheroid with flattening similar to the ice giants has been included to give context to the dynamical gravitational moments calculated from the simulation results. Calculations for J 2n are provided by Hofmeister et al. (2018): where e is the ellipticity of the spheroid. Similar to the coefficients of variation, gravitational moments are calculated from snapshots of the temperature and composition fields, then indexed temporally to create a time series. For each radial level, the spectral representation of temperature and composition are determined, which are then integrated to determine the values of ∆J n by reconstructing the spherical harmonic representation of the field and relying on the orthogonality of these functions (Holmes & Featherstone 2002;Wieczorek & Meschede 2018). From this time series, the statistical properties of the gravitational moments can be calculated. At low values of n, this solid body rotation produces a stronger gravitational signature than the dynamics in the system, but also falls off more rapidly. As a result, the dynamics produce a stronger signal at n 6 when Ek = 10 −3 (Figure 15b) and n 8 when Ek = 3 × 10 −4 (Figure 15a). Some insights can be gained into these results by considering the features of the temperature and composition fields in Figure 9. There is a reduction in the mean odd moments as the Ekman number is decreased. Comparing T1 to T2 and C1 to C2, the equatorial bulge in temperature and composition are larger in T1 and C1. That being said, these moments are composites of many snapshots in time, and we see that this bulge is, in most of the cases, associated with the non-symmetric features from the right column of Figure 7. Comparing T1 (Figure 7i) to T2 (Figure 7c), the temperature anomalies extend higher in the spherical shell, where the gravitational moments will be more sensitive. Furthermore, due to the three-dimensional nature of these flows, equatorial asymmetries will develop that are probed by the gravitational moments. In the calculated zonal gravitational moments, solid lines connect the even series, and dotted lines connect the odd series. The grey dotted line represents the even gravitational moments that result from a rotating, uniform spheroid with flattening that approximates the ice giants (Hofmeister et al. 2018). Measured gravitational moments, J2 and J4 are also included for Uranus and Neptune, though the values are very similar and overlap one another (Jacobson 2009(Jacobson , 2014Williams 2015).
In the case of ∆J 2 and ∆J 4 , the variance of the temperature and composition fields and the resulting uncertainty of the gravitational moments in our models is only moderate, but the standard deviation of the gravitational moments increases with n. The time series of the gravitational moments are stationary, however, meaning that the values presented in Figure 15 do not change over the course of the simulation. In case T1 and T2, only values of ∆J 2 and ∆J 4 were statistically non-zero (determined by a students-t test), but in the case of D only, ∆J 2 was non-zero. In cases C1 and C2, the harmonics ∆J 2 , ∆J 4 , ∆J 6 , and ∆J 10 were statistically non-zero. For odd zonal gravitational moments, values of the calculated harmonics tend to alternate between positive and negative values, such that the time-average is not statistically different from zero because the non-equatorially symmetric features of the temperature and composition tend to alternate between higher concentrations in the northern and southern hemisphere over time.

DISCUSSION AND CONCLUSIONS
We performed a suite of rotating convection models driven by thermal and/or compositional buoyancy sources at different Ekman numbers for a fixed convective Rossby number to probe hypothesized processes in ice giant planets. Most of the cases develop similar flow behaviors, including the formation of large overturning circulation cells and three zonal jets, but there is more diversity in the thermal and compositional behaviors. Equatorial peaks in thermal and mass fluxes appear in T1, D, and C2, with diminished equatorial flux in the case of T2 and C1 (see Figures 6 and  10). Case C1 (Ek = 10 −3 , Sc = 10) exhibits significantly different behavior. While this case has three zonal jets, their directions are reversed with prograde equatorial flows and retrograde bands at high latitudes; it similarly lacks an equatorial peak in mass flux that is likely caused by the disruption of circulation cells in the region near the equator. Reducing viscosity (C1 → C2 and T 1 → T 2) tends to improve the efficiency of absolute angular momentum mixing, which produces and/or amplifies a retrograde equatorial jet (Figure 11). Reducing the diffusivity of the thermodynamic field (T 1 → C1 and T 2 → C2) tended to produce zonal flows that adhered more closely to solid-body rotation, i.e. with weaker flow speeds ( Figure 11).

Zonal flows and mixing
Physical properties of the fluid in the interior of planets, such as viscosity and thermodynamic diffusivities, are more extreme than those used in our models, but we can try to predict the effect of scaling these parameters to planetary values based on how well the angular momentum and thermodynamics might be mixed as discussed in Section 5.2. The main determinant of how well mixed the thermodynamics became in our models was the choice of diffusivity ( Figure 12c). Planetary values of thermodynamic diffusivities are much smaller than those assumed in our models, likely even relative to the viscosity, and therefore the ice giant interiors may be more thoroughly mixed in terms of the temperature and composition. Reducing the viscosity causes the flow to become more turbulent, which allows the fluid to become more well mixed (Figure 12b). Continuing to reduce the viscosity to planetary values while maintaining convection characterized by 3D turbulence should result in increasingly well mixed systems, eventually exhibiting nearly uniform thermodynamic fields in the bulk and zonal winds that approach what is predicted by fully homogenized absolute angular momentum. In the case of χ = 0.35, full homogenization of angular momentum results in extreme equatorial jet speeds of 1500 m s −1 , and realistic speeds are not predicted for any spherical shell geometry if full homogenization is assumed. Restriction of angular momentum homogenization to occur only outside of the tangent cylinder still results in westward zonal speeds of 1300 m s −1 . The equatorial wind speeds of Uranus and Neptune can be be reproduced by fully-homogenized angular momentum mixing for χ = 0.99 and χ = 0.87, respectively, with thicker domains possible if only partial mixing were to occur or if full mixing of some different initial value of total angular momentum are assumed.
In contrast, decreasing the thermodynamic diffusivity has an antagonistic effect (Figure 12c) on the absolute angular momentum mixing, which reduces the intensity of zonal winds throughout the system. Because of the complex interactions between these effects, it is conceivable that with a sufficiently low thermal or mass diffusivity, that the system could be nearly perfectly homogenized in heat or mass, but still maintain a relatively poor state of angular momentum mixture, preventing the formation of winds as intense as those predicted by perfect angular momentum mixing.
These velocity fields emerge from dynamics that are, for the most part, intended to represent convection that is only weakly influenced by the effects of rotation. Global geostrophic balance is present in neither T1 nor T2, nor is it present in D or C2 (Table 3), with this type of balance present only at the largest length scales of the system. However, this is not the case in C1, where the integrated pressure gradient and Coriolis force are roughly balanced. Furthermore, the convective Rossby number failed to predict the direction of the zonal jet and the ratio of integrated total buoyancy to the Coriolis in case C1.
The density anomalies used to calculate the gravitational moments give rise to buoyancy that drives the flows within the system, meaning that we would expect similarities in the flows to be reflected in similarities between the gravitational moments of the system. Decreasing the Ekman number in the thermal cases does not result in significantly different zonal flows as shown in Figure 6, and despite the fact that there is a reversal of direction, the flows in the compositional cases have similar morphology. The magnitude of J 2 and J 4 tend to be similar as Ekman number is decreased as well as when the convective driving is varied. However, there are sign changes evident when comparing J 2 between the thermal and compositional cases and for J 4 in the case of T2 and C2. These gravitational moments are those that arise only from the dynamics of the fluid contained within the shell. If the effects of a spheroidal geometry can be ignored, then these flows would modify the gravitational signature of an oblate spheroid in a way such that C2 is more representative of the ice giants than T2. It would be useful to study the effects of decreasing the diffusivity of heat and light elements further to learn if this trend continues.

Light element and heat fluxes
The patterns of light element fluxes in our compositional and double diffusive cases may also be compared against the latitudinal variations in composition that have been observed in the ice giants. In particular, methane and hydrogen sulfide appear to have greater relative abundances in polar regions, although these observations are based on infrared observations of the ice giants that probe the top few bars of the atmosphere (e.g., Moses et al. 2020). Analogously to how sea water becomes brackish from the influx of freshwater near the mouth of a river, the measured deficiency of methane relative to other atmospheric gases may result from an influx of molecular hydrogen. We assume that there are no additional chemical reactions occurring in the deep atmosphere as the light elements are transported from the interior.
While the flux of light elements is large in the polar regions in the case of C1, this feature appears to be caused by the disruption of equatorial circulation cells that reappear as the viscosity is decreased. In the case of C2, local maxima in light element flux appear in the polar regions, but are still accompanied by strong equatorial upwellings that transport mass upward from the interior. The injection of a greater amount of light elements could affect the composition of the atmosphere, and while C1 may be consistent with the observed methane depletion in the polar regions of the ice giants, the equatorial peak in C2 implies that the composition of the atmosphere near the equator should be similarly affected. In the double diffusive case, the equatorial maxima also disappear, and the mass flux becomes dominated by an equatorial peak. For this reason, none of our cases are able to account for both the dynamical properties and the atmospheric compositions in a self consistent manner. Furthermore, based on trends with the Péclet number seen in our simulations, we may expect the mass flux to become more uniform with latitude when approaching planetary values of both the Reynolds number (which is very large, assuming a turbulent interior) and Schmidt number. Given this dissimilarity and the simplicity of our chemistry assumptions, additional observations and determination of high pressure-temperature material properties will be necessary to better understand the chemistry of the deep atmosphere.
While the internal heat flux of Neptune has local maxima at the equator and at the poles, the heat flux of Uranus appears to be less variable with latitude (Pearl et al. 1990;Pearl & Conrath 1991). The simulated heat flux patterns tend to have significant latitudinal variability with distinctive peaks at the poles and a secondary peak at the equator (T2), a distinctive peak and the equator and secondary peaks at high latitudes (T1 and D). Though Neptune's equatorial peak is larger than either polar peak, none of these cases are particularly reminiscent of the observed flux pattern. The change in simulated heat flux pattern as the Ekman number is reduced, however, indicates that further work is needed to determine how these results scale toward more realistic parameters. Similarly, none of the cases produce results that appear similar to the heat flux of Uranus. The mass flux of C2 is the least variable of all the cases and bears the greatest similarity to the heat flux of Neptune, suggesting that similar parameters (e.g., P r = 10) or perhaps a double diffusive case with a small relative contribution of thermal driving (e.g. f T ≤ 25%) may produce a heat flux with less latitudinal variability.
Tropospheric processes are not considered in our models, but the time scales involved in atmospheric transport are likely to be much faster than in interior processes. In the context of shallow atmospheric dynamics, mass and heat flux from the deep interior may be potentially significant forcing terms and this interior-atmosphere coupling is another exciting direction for future inquiry.

Density stratification
We assume that the fluid in our models is incompressible, which neglects any background density stratification, and that the outer boundary represents the 1 bar surface level in order to better facilitate comparisons with observations. Note that an alternative interpretation is that the outer boundary corresponds to the top of the ionic ocean, in which case the link to observations is obscured.
Interior models suggest that the density contrast across the ionic ocean is expected to be relatively moderate with much more rapid density changes closer to the surface (e.g., Helled et al. 2020). The density contrast between the inner and outer boundaries can be quantified by the number of density scale heights, N = log(ρ i /ρ o ). Model U1 of Nettelmann et al. (2013), for example, has a rapid decrease with increasing radius above r = 0.9r o . If we approximate the outer density as ρ o ≈ 100 kg m −3 at r ≈ 0.95r o and the inner density as ρ i ≈ 4000 kg m −3 at r ≈ 0.3r o , then about 3.6 scale heights exist across the shell. However, if the convective region of interest is confined to the ionic ocean, then the increase in density across this layer is not predicted to be very large, perhaps only representing a single scale height (Nettelmann et al. 2013).
The effect of anelasticity on the dynamics of the system is reflected primarily in the overturning circulations, which become concentrated in a layer near the surface; while the speed of zonal flows has been found to increase in the nonrotating regime of convection with increased density stratification, the morphology of the zonal flows are somewhat insensitive to a change in density (Gastine et al. 2012(Gastine et al. , 2013a. The density anomaly that gives rise to the zonal gravitational moments is most dependent upon zonal winds, and so they should also not be strongly dependent on density stratification. This not the case for non-axisymmetric gravitational moments, however, so their properties could be changed as the circulation is modified by density stratification. This should be the subject of future work.

Magnetic fields
Although numerous hypotheses exist for mechanisms that generate non-dipole-dominated magnetic fields and how the magnetic field morphology might be predicted via scaling laws (see Soderlund & Stanley (2020) for a recent review), we focus here on the local Rossby number because Ro l > 0.1 predicts the occurrence of multipolar magnetic fields in Boussinesq dynamo models reasonably well (e.g., Christensen 2010; Gastine et al. 2012). As discussed in Section 4.1, the local Rossby numbers in our simulations have Ro l 0.3; we, therefore, expect multipolar magnetic fields would be generated by these convective flows, not unlike those in the ice giants, which will be investigated in a follow-on paper.

Summary
We have investigated whether the properties of Uranus (particularly in comparison to Neptune) are consistent with deep convection that is driven by compositional, rather than thermal, gradients. While the relatively large internal heat flux of Neptune may be explained by primarily thermal convection, the low heat flux of Uranus may suggest a primarily chemical source of convection, which we have hypothesized results from the properties of mixtures of carbon-containing and oxygen-and nitrogen-containing compounds that comprise the planets and lead to the release of buoyant light elements in the deep interior. In chemically driven cases, thermal contrasts across the interior are replaced, either partially or totally, by compositional contrasts, which may produce detectable effects such as changes to relative atmospheric abundances. In most cases, zonal flows akin to those observed at Uranus and Neptune are generated, although the flow speeds are considerably stronger. In addition, buoyancy drives large convective cells that contribute to the transport of either heat or light elements to the surface, featuring larger heat fluxes in near the equator and in the polar regions, but suggesting that mass flux may become less variable with latitude than heat flux. Furthermore, convection in the interior provides a means for sustaining a dynamo, and the dynamic similarities between thermal and compositional suggest that these magnetic fields should be non-dipolar, which may explain the similarities of the magnetic fields of Uranus and Neptune, despite thermal differences.
Current observations of atmospheric composition (particularly its vertical structure) are limited to the top few bars of the atmosphere, but future in situ measurements, including those made by an atmospheric probe and gravity science, could improve our understanding of deep structure and dynamical processes, including the ice-to-rock ratio and mass fluxes predicted by our compositional convection models. Furthermore, the decrease of diffusivity of mass, relative to heat, results in more locally concentrated parcels of fluid enriched in light elements that may result in detectable gravitational features that reflect the interior dynamics at spherical harmonic degrees above n ≥ 6, which may be revealed by improved gravity science at the ice giants.