3D Simulations of Convective Entrainment in Gas Giants: Rotation and Decreasing Luminosity as Barriers to Mixing

TL;DR

This study uses 3D hydrodynamic simulations to show that decreasing luminosity and rotation can inhibit convective mixing in gas giants, potentially explaining the persistence of fuzzy cores with compositional gradients.

Contribution

It demonstrates that declining luminosity alone, especially when combined with rotation, can significantly suppress convective mixing in gas giants, a novel insight into planetary interior evolution.

Findings

Declining luminosity can prevent convective mixing in gas giants.

Rotation slows but does not entirely stop mixing; combined effects are significant.

Fuzzy cores may persist over long timescales due to suppressed convection.

Abstract

Observations from Juno and Cassini suggest that Jupiter and Saturn may possess fuzzy cores -- central regions where the abundance of heavy elements varies smoothly with depth. Such gradients pose a longstanding puzzle for models of planetary evolution and formation, which predict that vigorous convection would homogenize the interior of gas giants within the first $\sim 10^6$--$10^8~\mathrm{years}$ of cooling. Previous 3D simulations and analytic predictions for the propagation of a convection zone into a stable region have demonstrated that the rapid rotation of gas giants can significantly slow convective mixing, but not enough to stop it. Another piece of the puzzle is luminosity. Gas giants cool as they age, and with that comes a declining heat flux over time. Recent ideas suggest that when this declining luminosity is combined with rotational effects, convection may stall. We explore this possibility using 3D hydrodynamic simulations that include both rotation and a surface cooling flux that decreases as $1/t$. Our results demonstrate that, even without rotation, a declining luminosity can suppress mixing sufficiently to preserve an initial compositional gradient in the deep interior of gas giants. If confirmed by more realistic simulations, this may help to explain the long-term survival of fuzzy cores.

Figures

Equations

$$\frac{1}{2}gh\Delta \rho \frac{dh}{dt}\approx \frac{\gamma }{2}{\rho }_0{U}_c^3,$$

$$U_{\mathrm{NR}}\approx {\left(\frac{\alpha g{F}_{T}h}{{\rho }_0{c}_P}\right)}^{1/3},U_{\mathrm{R}}\approx {\left(\frac{\alpha g{F}_T}{{\rho }_0{c}_P}\right)}^{2/5}{\left(\frac{h}{2\Omega }\right)}^{1/5}.$$

$$\Delta \rho ={\rho }_0(\beta \Delta C-\alpha \Delta T),$$

$$h_{\mathrm{NR}}(t)={\left[h_0^2+4\gamma H^{2}Nt\right]}^{1/2},$$

$$h_{\mathrm{R}}(t)\approx {\left[h_0^{12/5}+\frac{24\gamma }{5}{H}^{12/5}{\left(\frac{N}{2\Omega }\right)}^{3/5}(Nt)\right]}^{5/12}.$$

$$h_{\mathrm{NR}}(t)\approx {\left[h_0^2+4\gamma H_0^2(N{t}_0)\ln (t/t_0)\right]}^{1/2},$$

$$h_{\mathrm{R}}(t)\approx {\left[h_0^{12/5}+24\gamma H_0^{12/5}(N{t}_0){\left(\frac{N}{2\Omega }\right)}^{3/5}\left(1-\frac{t_0^{1/5}}{t^{1/5}}\right)\right]}^{5/12},$$

$$\displaystyle\nabla\cdot{\mbox{{\boldmath$u$}}}=0\,,$$

$$\displaystyle\dfrac{\partial{\mbox{{\boldmath$u$}}}}{\partial t}+{\mbox{{\boldmath$u$}}}\cdot\nabla{\mbox{{\boldmath$u$}}}+\dfrac{\rm Pr}{\rm Ek}{\mbox{{\boldmath$\hat{z}$}}}\times{\mbox{{\boldmath$u$}}}=-\nabla P+\mathrm{Ra}\mathrm{Pr}\left(R_{\rho}C-T\right){\mbox{{\boldmath$\hat{z}$}}}$$

$$\displaystyle+\mathrm{Pr}\nabla^{2}{\mbox{{\boldmath$u$}}}\,,$$

$$\displaystyle\dfrac{\partial C}{\partial t}+{\mbox{{\boldmath$u$}}}\cdot\nabla C=\tau\nabla^{2}C\,,$$

$$\displaystyle\dfrac{\partial T}{\partial t}+{\mbox{{\boldmath$u$}}}\cdot\nabla T=\nabla^{2}T\,,$$

$$\mathrm{Ra}=\frac{\alpha g|{\partial }_z{T}_0|H^4}{{\kappa }_T\nu },\mathrm{Ek}=\frac{\nu }{2\Omega H^2},$$

$$R_{\rho }=\frac{\beta |{\partial }_z{C}_0|}{\alpha |{\partial }_z{T}_0|},\mathrm{Pr}=\frac{\nu }{{\kappa }_T},\tau =\frac{{\kappa }_C}{{\kappa }_T},$$

Full text

Planetary cores (1247); Planetary interior (1248); Solar system gas giant planets (1191); Hydrodynamics (1963); Hydrodynamical simulations (767); Convective zones (301)

1 Introduction

Modern interior models of Jupiter and Saturn, which fit observational constraints from gravitational fields (for both planets) and ring seismology (for Saturn), suggest inhomogeneous interiors, some of them predicting extended composition gradients (e.g., Fuller, 2014; Wahl et al., 2017; Debras & Chabrier, 2019; Mankovich & Fuller, 2021; Nettelmann et al., 2021; Militzer et al., 2022; Miguel et al., 2022; Idini & Stevenson, 2022; Helled et al., 2022; Dewberry, 2023; Howard et al., 2023). Currently, no consensus on a formation pathway to these current-day structures has emerged (Liu et al., 2019; Müller et al., 2020; Stevenson et al., 2022; Bodenheimer et al., 2025). These findings caused significant interest in gaining a better understanding of the mixing processes inside giant planets, particularly how convection–driven by cooling from the planet’s surface–may lead to the dispersal and/or homogenization of their cores.

In the past few years, several studies have explored this topic using 1D thermal evolution models (e.g., Müller et al., 2020; Knierim & Helled, 2024; Tejada Arevalo et al., 2025). A key result from these works is that if a planet forms hot–with a high internal entropy–vigorous convective mixing can take place on short timescales ($\sim$ 1 Myr), in some cases nearly homogenizing the entire planet’s interior. These simulations, however, lack a detailed model of how composition and energy are transported at convective boundaries. In an attempt to better model convective boundary mixing, Fuentes & Cumming (2020) used 2D hydrodynamic simulations of a local Cartesian box to investigate the propagation of a convection zone into a stably-stratified region by a composition gradient. Their experiments, which aimed to model convective mixing between the planet’s core and the envelope, showed that composition gradients are ineffective barriers against moving convective zones. Similar results were found by Anders et al. (2022) in the context of core-convection in stars.

Another way the planet’s interior could avoid complete mixing is through the formation of a semiconvective staircase beneath the outer convective envelope (e.g., Radko, 2003; Leconte, J. & Chabrier, G., 2012; Leconte & Chabrier, 2013). In this scenario, the competition between thermal and compositional gradients in the deep interior triggers double-diffusive instabilities. The nonlinear outcome is a remarkable form of self-organization, a series of convective layers separated by thin, stably stratified interfaces (see, e.g., Garaud, 2018, for a review). However, this mechanism has been challenged by hydrodynamical simulations, showing that staircases rapidly mix into a single convective layer (e.g., Wood et al., 2013; Moll et al., 2017; Fuentes et al., 2022; Tulekeyev et al., 2024).

Rapid rotation, inherently present in giant planets, is well-known to reduce the efficiency of convective heat transport when the Rossby number $\mathrm{Ro}$ — the ratio of the rotation period to the convective turnover time — is much smaller than one (e.g., Julien et al., 1996). Building on this, Fuentes et al. (2023) used classical scaling arguments for convective velocities in rapidly rotating flows (e.g., Stevenson, 1979; Barker et al., 2014; Aurnou et al., 2020) together with 3D simulations in boxes to revisit the problem of a convection zone advancing into a stably stratified region. They found that for typical values of the Rossby number in gas giants ($\mathrm{Ro}\sim 10^{-6}$; Guillot et al. 2004), the mixing time of primordial composition gradients can be delayed by more than two orders of magnitude, suggesting that rotation plays an important role to preserve inhomogeneous interiors. Along similar lines, Fuentes et al. (2024), using 3D local Cartesian simulations, and Fuentes (2025), using simulations in the full sphere, demonstrated that rotation also slows the merger of semiconvective layers. Yet, in all cases, the simulations ultimately evolved toward a fully convective state, i.e., mixing is slowed, not stopped.

Seeking a mechanism to fully suppress convective mixing, Hindman & Fuentes (2023) made an analytic model for the propagation of a convective interface. They suggested that convection in a gas giant can stall before fully mixing the interior, not because of a strong stabilizing composition gradient, but because the system runs out of energy. As the luminosity declines, convective motions weaken. This, in combination with rotation further suppressing the ability to entrain deeper material, causes the growing convective zone to stall, leaving behind a partially unmixed interior. While promising, this theory has yet to be tested with hydrodynamical simulations. In this paper, we aim to test these ideas. In Section 2 we summarize the analytic predictions of Hindman & Fuentes (2023). In Section 3 we present our numerical experiments and results. We conclude in Section 4 with a summary and discussion of how our results apply to gas giants.

2 Mixing-Length Scalings and Convective Entrainment

In this section, we summarize the mixing-length scalings for convective velocities in both nonrotating and rotating flows, and discuss how these scalings are applied to model the propagation of a convection zone into a stably stratified compositional gradient.

When the boundary between the convection zone and the adjacent stable region is thin, energetic arguments predict the speed at which the convection zone propagates. The entrainment hypothesis (e.g., Linden, 1975; Fernando, 1987; Molemaker & Dijkstra, 1997) states that the gravitational potential energy required to mix denser fluid with the overlying light fluid within the convection zone must come from the kinetic energy of the convective motions; explicitly, the rate of potential energy increase due to the upward transport of heavy material is proportional to the vertical kinetic energy flux carried by the convection. Such a hypothesis leads to a differential equation for the depth of the convection zone, $h$,

[TABLE]

where $g$ is the acceleration due to gravity, $\rho_{0}$ is the mass density within the convection zone, $\Delta\rho$ is the density jump across the thin interface, $U_{c}$ is a typical convective flow speed, and the proportionality constant $\gamma$ (known as the mixing efficiency) is of order unity for turbulent flows (Fuentes & Cumming, 2020). We emphasize that $\Delta\rho$ arises because composition in the convection zone is well mixed, creating a contrast in density with the underlying region which has a smooth composition gradient.

This entrainment equation can be solved for the thickness of the convection zone $h$ as a function of time, given prescriptions for the density jump across the convective interface $\Delta\rho$ and convective speed $U_{c}$. In the non-rotating and rapidly-rotating limits, typical convective speeds ($U_{\rm NR}$ and $U_{\rm R}$, respectively) can be estimated using mixing-length theory and turbulent, diffusion-free scaling laws for convective heat transfer (e.g., Barker et al., 2014; Aurnou et al., 2020; Fuentes et al., 2023)

[TABLE]

In these expressions, $F_{T}$ is the heat flux applied at the top boundary, $\alpha=-(\partial\ln\rho/\partial T)_{C}$ is the coefficient of thermal expansion at fixed composition, $c_{P}$ is the specific heat capacity at constant pressure, and $\Omega$ is the rotation rate.

On other hand, the density jump across the interface is determined by the corresponding jumps in composition and temperature,

[TABLE]

where $\beta=(\partial\ln\rho/\partial C)_{T}$ is the coefficient of compositional contraction at fixed temperature, and both $\Delta C$ and $\Delta T$ can be obtained by applying conservation of mass and energy within a layer of size $h$ before and after convective mixing. Since a sharp stably-stratified interface is dominated by the compositional jump across111For example, Fuentes & Cumming (2020) found that the $\Delta T$ term contributes a $\approx 20$–$30$% correction to $\Delta\rho$ at low Prandtl number (see their Fig. 5)., Hindman & Fuentes (2023) assumed $\Delta\rho\approx\rho_{0}\beta\Delta C$. For a constant driving heat flux and an initially linear density stratification, integration of Equations (1)–(2) in the limit of no rotation and neglecting heat diffusion across the interface yields the well-known result that the convective layer grows with a square root of time dependence

[TABLE]

where $h_{0}$ is an integration constant, $N^{2}=g\beta|dC_{0}/dz|$ is the initial buoyancy frequency and $H=(g\alpha F_{T}/\rho c_{P}N^{3})^{1/2}$ is an overshooting length scale that characterizes the depth to which convective plumes penetrate locally into the stable region (Molemaker & Dijkstra, 1997; Fuentes & Cumming, 2020). However, in the presence of rapid rotation, the growth of the convective layer is slowed, following a power law with an index of 5/12 (Fuentes et al., 2023)

[TABLE]

Adapting the model of Fuentes et al. (2023) to include a more realistic, time-dependent boundary condition for the heat flux, $F_{T}=F_{0}(t_{0}/t)$, Hindman & Fuentes (2023) developed an analytic model and showed that in rapidly rotating flows, the growth of the convective layer eventually stalls, whereas in the absence of rotation, the layer continues to deepen over time

[TABLE]

where $H_{0}=(g\alpha F_{0}/\rho c_{P}N^{3})^{1/2}$. While the predictions in Equations (4) and (5) have been investigated and validated through 3D numerical simulations (Fuentes et al., 2023), the predictions corresponding to the time-dependent flux in Equations (6) and (7) have not yet been verified by numerical simulations. The stalled solution is particularly relevant to the internal structure of gas giants, as it implies the preservation of primordial composition gradients beneath the convective front. In this work, we verify these predictions using 3D hydrodynamical simulations.

3 Numerical Simulations

We present four 3D simulations of penetrative convection, both with and without rotation, in a plane-parallel domain with equal horizontal and vertical extent $H$. We follow the same numerical setup and use the same code as described in Fuentes et al. (2023). The lateral boundaries are periodic, while the top and bottom are impermeable and stress-free, with zero compositional flux. The fluid is initially stabilized by a linear compositional gradient. Cooling is imposed at the top boundary by fixing the thermal gradient with a negative heat flux, either held constant at $F_{0}$, or allowed to decay with time as $F_{0}~(t_{\star}/(t_{\star}+t))$, where we introduced $t_{\star}\sim 5\times 10^{-3}\ll t$ to avoid divergence of the flux at $t=0$. In both cases, the effect is to destabilize the upper layers until convection sets in. As the system evolves, convection progressively erodes the compositional gradient underneath, establishing a well-mixed region that deepens with time. We investigate how the growth rate of convective zone changes by with the imposed heat flux and rotation.

In our simulations, as in the predictions above, we assume that fluid motions occur on scales much smaller than a density scale height. This justifies neglecting curvature effects and adopting the Boussinesq approximation (Spiegel & Veronis, 1960). Although the density in giant planets can vary by many orders of magnitude, this approximation fully captures nonlinear mixing near the convective boundary, which is the focus of our study. In the rotating cases, we consider uniform rotation throughout the fluid, and assume that gravity and rotation are aligned and point in the z direction, {\mbox{{\boldmath\Omega}}}=\Omega\hat{{\mbox{{\boldmathz}}}}, and {\mbox{{\boldmathg}}}=-g\hat{{\mbox{{\boldmathz}}}}, i.e., the simulation is at a polar latitude. The density of the fluid depends on both composition and temperature fluctuations. In what follows, all results are presented in dimensionless form, with time normalized to the thermal diffusion time across the box, $\tau_{T}=H^{2}/\kappa_{T}$, where $\kappa_{T}$ is the thermal diffusivity, and length normalized to the box height $H$. For more details on the numerical setup, nondimensionalization, and code, we refer the reader to the Appendix.

Figure 1 shows 3D snapshots of the $z$ component of the fluid vorticity, \hat{{\mbox{{\boldmathz}}}}\cdot(\nabla\times{\mbox{{\boldmathu}}}). Since both non-rotating and rotating cases with fixed-flux, boundary-driven convection zones have already been presented in Fuentes et al. (2023), we focus here on new simulations where convection is driven by a gradually declining heat flux. Snapshots are shown for both the non-rotating (panel a) and rotating (panel b) cases at the same time, $t\approx 8\times 10^{-2}$, when the flux heat has decreased to about $10\%$ of its initial value. As expected, in the absence of rotation, the convection zone has entrained and mixed more fluid than in the rotating case, since the kinetic energy flux available for mixing is larger. The flow is also isotropic, with length scales that are approximately equal in all directions. In contrast, in the rotating case, the convective flow is anisotropic, forming vortices and columns whose horizontal length scales are much smaller than the depth of the convection zone.

Figure 2 shows the temporal evolution of the horizontally-averaged composition. The dashed curves show the extent of the convection zone as measured from density profiles. For constant heat flux, convection mixes the entire box as expected. When the flux decreases with time, the convection zone stalls even without rotation. This departs from the predictions in Section 2, where stalling was expected only in the rotating case.

To test the analytic predictions, we first verify that both scalings for the convective velocity in Equation (2) are consistent with our numerical experiments. To do so, we compute the average vertical component of the rms velocity across the convection zone.

In the rotating cases, a large-scale vortex develops over time, enhancing the horizontal velocities, often by an order of magnitude compared to the vertical flow. The analytical model does not capture this vortex, including only the vertical velocity. However, we find that it reproduces the numerical results well even in the presence of the vortex, suggesting that the vortex does not affect the entrainment process and that the entrainment is set by the strength of the vertical velocity.

As shown in Figure 3(a), the scaling of the convective velocity with the heat flux and the thickness of the convective zone is well supported by our simulations, except at early times, when the convection zone is still shallow and not fully developed.

Figure 3(b) shows the measured size of the convection zone compared with the analytic predictions. Overall, the simulations are consistent with the theory in that, for rotating cases, the rate of entrainment and the growth of the convective layer are significantly reduced compared to the nonrotating cases. For the cases with constant flux, there is excellent agreement with Equations (4)–(5). In contrast, for the cases with a decreasing flux, we observe noticeable deviations from Equations (6) and (7), specifically the measured $h(t)$ curves become much flatter than predicted at late times. Preliminary analyses suggest that these discrepancies arise from a time-dependent mixing efficiency in the entrainment equation. In cases with declining flux, the efficiency decreases asymptotically to 0, whereas it remains approximately constant and of order unity in cases with constant flux. A complete characterization of the mixing efficiency requires a thorough analysis of the heat and composition transport across the interface, which is beyond the scope of this paper. Nevertheless, we emphasize that the most important result, the stall of the convection zone, is confirmed by our hydrodynamical simulations. Furthermore, even the nonrotating case seems to stall.

4 Discussion and Conclusion

In this work, we explored how the deep interiors of gas giants may resist convective mixing when the planet cools over time with a decreasing luminosity. Using 3D hydrodynamical simulations, we tested recent analytical arguments by Hindman & Fuentes (2023) (see also Equations 4–7), which suggest that convection can stall at a finite depth because rotation and a declining heat flux act together to inhibit overturning motions. In this picture, the convective envelope fails to reach the core, potentially preserving a primordial interior that remains largely untouched since formation.

As shown in Figure 3, there is excellent agreement between the simulations and the theoretical predictions for the growth rate of the convective layer in cases with constant driving flux, both with and without rotation (Equations 4–5). However, when the driving flux decays with time, the models in Equations (6) and (7) tend to overpredict the extent of the convection zone. The simulations show that the convection zone eventually stalls in these declining flux cases, even without rotation. The size of the stalled convection zone is different in the cases with and without rotation, however. In terms of the size of the box, the convection zone size is $h\approx 0.61$ in the non-rotating case, whereas with rotation it is only $h\approx 0.38$, nearly half as deep, showing that rotation limits the extent even though the qualitative behavior is similar.

Whereas the analytic model in Equation (1) (from which Equations (6) and (7) are derived) assumes a constant mixing efficiency $\gamma$, the simulations show that $\gamma$ decreases significantly in the dwindling-flux cases, leading to a stalled convection zone even in the absence of rotation. This reduction in $\gamma$ likely reflects a change in the physical mechanism of entrainment at low fluxes, where viscous dissipation may compete with or even overwhelm buoyant driving, although further work will be needed to clarify the origin of this behavior.

Our findings could have important implications for 1D evolution codes (e.g., Vazan et al., 2018; Müller et al., 2020; Tejada Arevalo et al., 2025; Sur et al., 2024, 2025). While these models account for a declining luminosity over time, if the initial entropy is high enough (a hot start), the outer convection zone can reach deep into the planet, sometimes nearly all the way through. In contrast, for low-entropy initial conditions (cold starts), the outer convection zone stalls. However, we emphasize that such models rely on mixing-length theory and lack a physically motivated prescription for turbulent entrainment across convective boundaries. Incorporating entrainment laws such as the ones in this work could significantly improve the predictions of 1D evolution models. A first step in this direction was taken by Scott et al. (2021) in the context of core convection in stars, where they showed that including an entrainment-based mixing scheme helps to reproduce key observational features, such as the mass dependence of the main-sequence width in the HR diagram. Similar improvements could be expected in planetary evolution models with compositional gradients, likely improving constraints on the mass and size of fuzzy cores.

Our results may also offer new insight into the energetics of core erosion. Previous estimates by Helled et al. (2022), and more recently by Fuentes et al. (2025) suggest that if the energy released during the planet’s cooling history is available to drive mixing (a scenario they called the “global story”), then even a compact or fuzzy core could be fully mixed into the outer convective envelope. In contrast, if mixing is restricted to the thermal energy available near the core-envelope interface (the “local story”), then the energy is insufficient to fully erode any core. These two scenarios, however, neglect the role of rotation in modulating convective mixing. Our findings suggest that rapid rotation, combined with a declining cooling flux, can substantially reduce the mixing efficiency of convection, potentially inhibiting core erosion even in the global story.

In light of current observations of Jupiter and Saturn, there is no consensus on the extent of the well-mixed convective region, which determines the size of the fuzzy core (see Figure 2 in Helled et al., 2022). For example, the models of Militzer et al. (2022), constructed from Jupiter’s zonal gravity data, and those in Idini & Stevenson (2022), based on tidal interactions between Io and Jupiter, constrain the extent of the fuzzy core up to $0.7R$, leaving the outer $30\%$ of the radius well-mixed and vigorously convective. In contrast, the models of Nettelmann et al. (2021) truncate the core at $\sim 0.4R$ for both Jupiter and Saturn. Although the analytic models of Hindman & Fuentes (2023) predict that the convection zone is confined to the outer $30\%$ of Jupiter and $20\%$ of Saturn, we stress that both those models and our simulations remain too limited to provide a reliable estimate of the convection zone’s depth.

We encourage that future work should relax some of the simplifications adopted in this work. Our simulations neglected density stratification, as we adopted the Boussinesq approximation, which formally applies when the vertical extent of the fluid is small compared to a density scale height. In real planets, however, convection is driven by surface cooling, while the region of interest (the core) lies many scale heights below. Therefore, it is not clear whether the entrainment rates seen in Boussinesq models extend to more realistic, stratified simulations. Including density stratification could modify mixing through several ways. The energy required to dredge up heavy material from depth is different, and the nature of convective mixing itself changes: Boussinesq convection homogenizes density, while stratified convection restores an adiabatic gradient. Stratification also enhances asymmetries between upflows and downflows, altering the net kinetic energy flux available for mixing. Therefore, hydrodynamical models that go beyond Boussinesq would be of great interest.

In addition, our simulations were performed in a Cartesian box where gravity and rotation are aligned, which is appropriate near the poles. Real planets are spherical, so gravity decreases radially with depth and the rotation vector is aligned with gravity only at the poles. Fuentes (2025) investigated semiconvective mixing in the full sphere, and showed that the latitude-dependent misalignment of gravity and rotation, combined with the radial variation of gravity, alters the buoyancy flux required for mixing. A non-constant gravity profile changes the potential energy cost of lifting heavy elements upward and therefore modifies the entrainment rate. They also found that these effects produce anisotropic mixing. In particular, Figures 6 and 7 of Fuentes (2025) show stronger heat and composition fluxes near the equator than at higher latitudes, favoring cylindrical rather than spherical mixing and raising the possibility of a non-spherical fuzzy core. The global simulations on the sphere also permit the development of large-scale flows such as meridional circulation and zonal flows that are inhibited in Cartesian domains.

J.R.F. is supported by the Sherman Fairchild Postdoctoral Fellowship at Caltech, and NASA Solar System Workings grant 80NSSC24K0927. A.C. and S.Z. acknowledge support from NSERC Discovery Grant RGPIN-2023-03620. S.Z. was supported by a Max Stern Recruitment Fellowship from McGill University.

Appendix A Fluid equations and numerical methods

We express the fluid quantities as the sum of a linear hydrostatic background (denoted by the subscript 0) and a dynamic perturbation to the background (denoted by primes), e.g., the total temperature and composition are expressed as $T=T_{0}(z)+T^{\prime}$, and $C=C_{0}(z)+C^{\prime}$, respectively. The density perturbations satisfy $\rho^{\prime}/\rho_{0}\ll 1$, and are related to $T^{\prime}$ and $C^{\prime}$ through $\rho^{\prime}=\rho_{0}(\beta C^{\prime}-\alpha T^{\prime})$, as demanded by the Boussinesq approximation (Spiegel & Veronis, 1960). Here, $\beta$ and $\alpha$ are the coefficients of compositional and thermal contraction/expansion (both assumed positive constants), respectively. Before presenting the fluid equations, we non-dimensionalize them using $[T]=|\partial_{z}T_{0}|H$, $[C]=|\partial_{z}C_{0}|H$ as units of temperature and composition. We use the domain’s depth $H$ as the unit of length, and the thermal diffusion time $\tau_{T}=H^{2}/\kappa_{T}$ (where $\kappa_{T}$ is the thermal diffusivity) as the unit of time. By this choice, a unit of pressure corresponds to $[P]=\rho_{0}(\kappa_{T}/H)^{2}$. The dimensionless equations are

[TABLE]

where D/Dt=\partial_{t}+{\mbox{{\boldmathu}}}\cdot{\mbox{{\boldmath\nabla}}}, and $u$ is the velocity field.

There are 5 dimensionless numbers that characterize the evolution of the flow. These are the Rayleigh number Ra, Ekman number Ek, stability ratio $R_{\rho}$, Prandtl number Pr, and diffusivity ratio $\tau$, which are defined respectively as

[TABLE]

where $\Omega$ is the rotation rate, $\nu$ is the kinematic viscosity, and $\kappa_{C}$ is the compositional diffusivity. $R_{\rho}$ is the ratio of the stabilizing and destabilizing effect of the compositional and thermal buoyancy.

We initialize the fluid with uniform temperature profile, $T_{0}(z)=1$, and a linear composition gradient, $C_{0}(z)=1-z$. This setup ensures that the initial stratification is stable against both convection and double-diffusive instabilities, regardless of the value of the input stability ratio $R_{\rho}$. We impose zero compositional flux at both the top and bottom boundaries, $\partial_{z}C|_{z=0,1}=0$. Similarly, the thermal boundary conditions consist of zero heat flux at the bottom, $\partial_{z}T|_{z=0}=0$), and a time-dependent (or constant) heat flux at the top, $\partial_{z}T|_{z=1}=f(t)$, where $f(t)=-1$ for the cases of constant cooling flux, and $f(t)=-t_{\star}/(t+t_{\star})$, with $t_{\star}\sim 5\times 10^{-3}$, for the cases with declining flux. The boundary conditions for the velocity are impenetrable and stress free at both boundaries (\hat{z}\cdot{\mbox{{\boldmathu}}}=\hat{x}\cdot\partial_{z}{\mbox{{\boldmathu}}}=\hat{y}\cdot\partial_{z}{\mbox{{\boldmathu}}}=0 at $z=0,1$).

All the simulations in this work use the following input parameters: Rayleigh number $\mathrm{Ra}=2\times 10^{10}$, stability ratio $R_{\rho}=0.1$, Prandtl number $\mathrm{Pr}=0.1$, diffusivity ratio $\tau=0.1$, and Ekman number $\mathrm{Ek}=3\times 10^{-6}$. For the non-rotating cases, we set the Coriolis force to zero by taking $\mathrm{Ek}=\infty$. In comparison, under Jovian interior conditions, typical values are $\mathrm{Ra}=\mathcal{O}(10^{43})$, $\mathrm{Ek}=\mathcal{O}(10^{-17})$, and $\tau\approx\mathrm{Pr}=\mathcal{O}(10^{-2})$ (see, e.g. Guillot et al., 2004; French et al., 2012; Dhouib et al., 2024). The value of the stability ratio $R_{\rho}$ remains uncertain due to limited knowledge of the thermal and compositional gradients deep inside gas giants, but we expect it to vary from $R_{\rho}\lesssim 1$ in convective regions, to $R_{\rho}\gg 1$ in the core. Nevertheless, by selecting parameters such that $\mathrm{Ra}\gg 1$, $\mathrm{Ek}\ll 1$, $\mathrm{Pr}<1$, and $\mathrm{Le^{-1}<1}$, we ensure that our simulations lie within the same qualitative dynamical regime as that of gas giant interiors.

We time-evolve Equations (A1)–(A4) using the Dedalus pseudospectral solver (Burns et al., 2020) version 3, using timestepper SBDF2 (Wang & Ruuth, 2008) and CFL safety factor of 0.2. All variables are represented using a Chebyshev series with 384 terms for $z\in[0,1]$ and Fourier series with 384 terms in the periodic $x$ and $y$ directions, where $x,~y\in[0,1]$. We use the 3/2-dealiasing rule in all directions, so that nonlinearities are calculated in physical space on a $576^{3}$ grid. To start our simulations, we add random-noise temperature perturbations sampled from a normal distribution with a magnitude of $10^{-5}$ compared to the initial temperature field.