3D Simulations of Semiconvection in Spheres: Turbulent Mixing and Layer Formation

We consider a sphere of radius $R$ filled with a stably-stratified fluid comprised of two components, one light and the other relatively heavy (e.g., hydrogen and heavy elements). The concentration of the heavy component is represented by $C$. We adopt the Boussinesq approximation (Spiegel & Veronis, 1960), i.e., we express the fluid density $\rho$ with a linear relationship between the temperature and composition, $\rho=\rho_{0}(1-\alpha T+\beta C)$, where $\rho_{0}$ is a fiducial density and $\beta$ and $\alpha$ are the coefficients of compositional and thermal contraction/expansion (both assumed positive constants), respectively. Since the Boussinesq approximation is valid for subsonic flows and spatial scales smaller than a density scale height, our simulations should not be interpreted as models of entire stars or planets. Instead, they are more appropriate for the inner core of a gas giant, where the density varies by a factor of 2–4 within 20%–25% of the radius (see, e.g., Helled & Howard, 2024). In gas giant interiors, the density scale height is comparable to the planetary radius, making the Boussinesq approximation reasonable for a small core. For stars, this assumption is more questionable, and we discuss its limitations later in the paper.

We initialize the temperature and composition profiles to increase quadratically with depth, i.e., $T_{0}(r)=\Delta T_{0}(1-r^{2}/R^{2})$ and $C_{0}(r)=\Delta C_{0}(1-r^{2}/R^{2})$. By this choice the initial density ratio is constant across the sphere ($R_{\rho}=\beta\Delta C_{0}/\alpha\Delta T_{0}$).

We present the fluid equations in nondimensional form, using the radius of the sphere $R$ as the unit of length, and the thermal diffusion time $R^{2}/\kappa_{T}$ as the unit of time. For the units of temperature and composition, We use $[T]=\Delta T_{0}$ and $[C]=\Delta C_{0}$, respectively. A unit of pressure corresponds to $[P]=\rho_{0}(\kappa_{T}/R)^{2}$. Under this normalization, the surface of the sphere is located at $r=1$. The dimensionless equations are

where $u$ is the velocity field, and $g(r)=r$ is the dimensionless gravity. To be consistent with the initial profiles, and to ensure a stationary solution, we fix the gradient of temperature and composition at the surface of the sphere to match the initial conditions, so that $\partial_{r}T|_{r=1}=\partial_{r}C|_{r=1}=-2$ . The boundary conditions for the velocity are impenetrable and stress free, $u_{r}=\partial_{r}(u_{\theta}/r)=\partial_{r}(u_{\phi}/r)=0$, where $u_{r}$, $u_{\theta}$ and $u_{\phi}$, the radial, polar, and azimuthal components of the velocity, respectively.

The set of equations (3)–(6) is governed by five dimensionless numbers: the Prandtl number, the diffusivity ratio, and density ratio (defined in Equation 1), as well as the Rayleigh number and Ekman number, which are defined respectively as

where $g_{0}$ is the gravity at $r=1$, and $\Omega$ is the rotation rate.

All the simulations were run with the same Rayleigh number, Prandtl number, and diffusivity ratio (${\rm Ra}=10^{9}$, ${\rm Pr}=0.1$, and $\tau=0.1$). Under this choice, double-diffusive instabilities are expected to arise for $1<R_{\rho}<5.5$. Since we are interested in the layered regime, we varied the stability of the fluid by changing the density ratio $R_{\rho}$ from 1.25 up to 2. We set the Coriolis force to be identically zero in the non-rotating cases ($\mathrm{Ek}=\infty$), while for the rotating cases we vary the Ekman number between $\mathrm{Ek}\approx 2.6\times 10^{-6}$–$10^{-5}$. Real astrophysical conditions, where $\mathrm{Ra}\sim 10^{30},\mathrm{Ek}\sim 10^{-20},~{}\mathrm{Pr}\sim\tau\sim 10^{-7}$–$10^{-2}$ (Jermyn et al., 2022) are inaccessible for current computational capabilities. However, by fixing $\mathrm{Ra}\gg 1$, $\mathrm{Ek}\ll 1$, $\mathrm{Pr}<1$, and $\tau<1$, we ensure that the simulations are qualitatively in the same dynamical regime as stars and giant planets.

We time-evolve equations (3)–(6) using the Dedalus pseudospectral solver (Burns et al., 2020) version 3. The variables are represented in spherical harmonics for the angular directions and Jacobi polynomials for the radial direction. The number of radial, latitudinal, and longitudinal coefficients in all the simulations are $(N_{r},N_{\theta},N_{\phi})=(256,384,768)$, respectively. For time-stepping, we use a second order semi-implicit BDF scheme (SBDF2, Wang & Ruuth, 2008), where the linear and nonlinear terms are treated implicitly and explicitly, respectively. To ensure numerical stability, the size of the time steps is set by the Courant–Friedrichs–Lewy (CFL) condition, using a safety factor of 0.2 (based on trial and error). To prevent aliasing errors, I apply the “3/2 rule” in all directions when evaluating nonlinear terms. To start the simulations, we add random-noise temperature perturbations sampled from a normal distribution with a magnitude of $10^{-5}$ compared to the initial temperature field.

We begin with a reference case with stability parameter $R_{\rho}=1.25$. This case is the closest to overturning convection. We find that the fluid becomes rapidly unstable to semiconvection everywhere. There is a short phase where the system is in the regime of homogeneous oscillatory convection, with weak deviations from the initial density profile. Between $t\approx 0.03$–0.04, a few convective layers form spontaneously, with separating interfaces that are highly distorted due to overshooting motions. Since gravity (and in turn, the buoyancy fluxes) increases towards the surface, the system is dominated by a vigorous outer convection zone that propagates inwards, engulfing and mixing everything on its way (Figure 1a–d). To clearly see the step-like structures of the layered regime of semiconvection (convective staircases), we show radial profiles of the density in Figure 1e). Instead of performing a full spherical average, we measure the profiles at particular angles $\theta$ and $\phi$ because fluctuations due to overshooting motions cause an artificial interface broadening when horizontally-averaging over all $\theta$ and $\phi$. Although the profiles are noisy, they are consistent with the snapshot in panel (b), where 3 convective layers are clearly seen at $t\approx 0.036$. Finally, at $t\approx 0.12$, the outer convective layer reaches the center of the sphere and the entire fluid becomes fully-mixed.

Simulations with larger density ratios ($R_{\rho}=1.5,~{}1.75,2$) exhibit qualitatively similar behavior: multiple convective layers form in the fluid (Figure 2), which subsequently merge and mix until the sphere becomes fully convective. However, the stability of the diffusive interfaces between the layers increases with larger $R_{\rho}$, as the stabilizing effect of the compositional gradient becomes stronger. Consequently, the timescale for the sphere to become fully mixed increases with $R_{\rho}$, with $t_{\rm{mix}}\approx 0.12,~{}0.175,~{}0.23,~{}0.32$ for $R_{\rho}=1.25,~{}1.5,~{}1.75,~{}2$, respectively.

We are interested in the radial fluxes of heat and composition through semiconvection. These fluxes are usually measured in terms of thermal and compositional Nusselt numbers, which are defined as the ratio of the total flux to the diffusive flux

where the notation $\langle\cdot\rangle$ represents a volume average over the entire sphere.

Figure 3(a) shows time series of the thermal and compositional Nusselt numbers during the evolution of the simulations. In all cases, the Nusselt numbers increase during different phases characterized by different slopes of the curves with time. First, there is an exponential increase where the semiconvective instability grows and saturates once nonlinear effects become important. Once semiconvection fully establishes and layers form and evolve, both numbers increase over time but a smaller rate. Later on (starting from the times denoted by the vertical lines in Figure 3a), the increase in the turbulent transport becomes steeper with time as the fluxes become dominated by the more vigorous outer convection zone. Note that due to the reduced overall stability of the flow at lower $R_{\rho}$, the contribution from convective fluxes and consequently the magnitude of $\mathrm{Nu}$, is higher for smaller $R_{\rho}$. Finally, as the entire sphere transitions to a fully convective and well-mixed state, the Nusselt numbers reach a peak and settle to a constant value, independent of $R_{\rho}$, but controlled by the imposed fluxes at the surface of the sphere which are the same for all the simulations.

There exist different models for the relationship between the thermal and compositional Nusselt numbers. For example, in the analytical model of Spruit (2013), the Nusselt numbers are related via

where $q$ is a constant of order unity.

On the other hand, the numerical simulations of Rosenblum et al. (2011) and Wood et al. (2013) propose

respectively, where the prefactors $A\approx 0.1$ and $B\approx 0.03$ are weakly dependent on $R_{\rho}$, $\mathrm{Pr}$, $\tau$, or $\mathrm{Ra}$.

We tested the three models against measurements of the Nusselt numbers from the simulation results. When averaging over the semiconvective phase (i.e. the time span between the end of the initial exponential grow and the vertical lines in Figure 3(a)), we find that none of the models give the exact measured values from the simulations. However, Spruit’s theoretical prediction (Equation 10) gives a better fit to the data (Figure 3(b)). Note that Rosenblum’s prediction (Equation 11) cannot be ruled out, as it provides the correct scaling with the density ratio $R_{\rho}$, and adjusting the prefactor to $\approx 0.33$ also gives a good fit.

Another important quantity in recent studies of semiconvection is the “buoyancy flux ratio” $\gamma^{-1}$, defined as the ratio of the buoyancy flux associated with compositional transport to the equivalent buoyancy flux associated with heat transport (see, e.g., Garaud, 2018, 2021)

As the sphere becomes mixed over time, the gravitational potential energy released from the thermal stratification must exceed the potential energy required to mix the compositional stratification, implying $\gamma^{-1}<1$. Figure 4 shows that $\gamma^{-1}$ remains below unity these simulations. Initially, during the homogeneous phase (before layer formation), $\gamma^{-1}$ decreases as a function of $R_{\rho}$. Over time, $\gamma^{-1}$ increases and reaches a peak, although no clear trend emerges for how its maximum value varies with $R_{\rho}$. Previous studies have proposed that $\gamma^{-1}=\tau^{1/2}$ in experiments of semiconvection in salty water (e.g., Turner, 1965; Linden, 1975), or $\gamma^{-1}\sim\tau(\mathrm{Pr}+1)/(\mathrm{Pr}+\tau)$ in numerical simulations of semiconvection in triply periodic boxes (e.g., Moll et al., 2017). For $\tau=\mathrm{Pr}=0.1$ these predictions yield $\gamma^{-1}\approx 0.32$ and $0.55$, respectively. In comparison, our results range from 0.3 to 0.4 during the semiconvective regime. Once the sphere becomes fully mixed, both the composition and heat flux reach a statistically steady state, varying linearly with radius. At this stage, the volume-averaged value of $\gamma^{-1}$ approaches the steady-state solution, $\tau R_{\rho}$ (horizontal lines in Figure 4). On average, we find that our results are broadly consistent with those from non-rotating simulations in Cartesian boxes (see, e.g., Garaud, 2021, for a review).

The evolution of the rotating runs is qualitatively similar to that of the nonrotating cases, i.e., once the semiconvective instability is triggered, a few convective layers form, which grow and merge until the entire sphere is well-mixed and fully convective. However, there are some significant differences with respect to the nonrotating counterpart. First, the thermal and compositional transport measured through the Nusselt numbers and the buoyancy flux ratio are much smaller than for the nonrotating case (see Figure 5). Consequently, the sphere becomes fully mixed at later times, with the difference respect to the nonrotating case being greater as the Rossby number decreases ($t_{\mathrm{mix}}\approx 0.4,~{}0.3,~{}0.26$ for $\mathrm{Ro_{c}}=0.25,~{}0.5,~{}1$, respectively). Further, unlike the nonrotating cases, once the sphere becomes fully convective, the Nusselt numbers are not the same between the different simulations, despite that $R_{\rho}$ and the surface gradients are equal in all the simulations. The explanation for this discrepancy is the weaker convective transport and the well-known enhancement of the thermal gradient in rotating convection, which becomes more pronounced as the Rossby number decreases (see, e.g., Barker et al., 2014). This implies that, compared to the non-rotating case, a larger fraction of the heat flux is carried by conduction rather than convection as rotation increases. Overall, rotation exhibits a stabilizing effect on convective fluxes, similar to that of increasing the stability ratio $R_{\rho}$ in nonrotating simulations. Note since we run simulations at fixed $R_{\rho}=1.25$, we cannot properly test the analytical predictions in Equations 11–12. However, we found that the ratio between the Nusselt numbers varies between 2.5 and 3.5, i.e., close to that of the nonrotating simulation at the same density ratio.

Another difference is that rotation significantly affects the spatial scales of the convective flow. In the nonrotating simulations, the flow exhibits approximately isotropic length scales, whereas in the rotating cases, the horizontal length scale (perpendicular to the rotation vector) becomes much smaller than the vertical scale as rotation increases. This behavior is consistent with the Taylor-Proudman theorem, which predicts the formation of vortices and thin columns aligned with the rotation axis (e.g., Proudman, 1916; Taylor, 1917). This difference is clearly seen for the case at $\mathrm{Ro_{c}}=0.25$, at mid and high latitudes (Figure 6c).

To illustrate how the fluid mixes over time, we show in Figure 6(d) meridional planes of the density field for $\mathrm{Ro_{c}}=0.5$. We selected this case because it has an intermediate level of rotational influence, and evolves in a similar way to $\mathrm{Ro_{c}}=0.25$, whereas the case $\mathrm{Ro_{c}}=1$ behaves more similar to the nonrotating cases. Initially, the flow is dominated by elongated structures aligned with the rotation axis, which subsequently merge to form 3–4 convective layers. The lack of isotropic density redistribution, which would be characteristic of spherical mixing, further emphasizes the role of rotational effects in organizing the flow. For example, the extent of mixing in the equatorial regions is clearly larger than at the poles. This indicates that mixing is governed by cylindrical geometry, shaped by the interaction between buoyancy and Coriolis forces.

We also observe the formation of prograde zonal flows at the surface, which alternate and reverse to retrograde and then back to prograde as a function of radius. The radial shear between the flows significantly distorts the interfaces between the layers. A comprehensive study exploring the fluid’s transport properties, flow morphology, and zonal flow characteristics as a function of $\mathrm{Ra}$, $R_{\rho}$, $\mathrm{Pr}$, and $\tau$ will be presented in a future publication.

Given the significant differences between convective flows in the polar regions and those at lower latitudes near the equator, we quantify these variations by measuring local Nusselt numbers as a function of latitude. Since we aim to compare with results from previous studies of thermal convection in rotating shells, we show azimuthally averaged profiles of the Nusselt numbers at $r\approx 0.9$, i.e., in the outer convective layer where the fluxes and the velocities are larger. To ensure the flux measurements are captured early in the semiconvective phase and to further reduce noise, we time-average the profiles over $\Delta t\approx 0.01$ centered around $t\approx 0.05$.

Figure 7 shows that both compositional and heat transport in the radial direction are more effective near the equator while becoming less efficient at higher latitudes. This aligns well with results from previous studies focusing on thermal convection alone (see, e.g., Figure 3c in Wang et al., 2021). Recently, Gastine & Aurnou (2023) showed that differences in heat transport as a function of latitude depend strongly on the supercriticality of the flow, defined as $\mathcal{R}=\mathrm{Ra}\mathrm{Ek}^{4/3}$. For $\mathcal{R}\sim 1$–$10$, convective transport is expected to be more efficient at the equator than at the poles, leading to larger Nusselt numbers at the equator. However, for larger supercriticalities ($\mathcal{R}>10$), heat transport at the poles becomes increasingly efficient as $\mathcal{R}$ increases, reaching a point where the equator and poles transport heat equally efficiently for $\mathcal{R}>100$, when rotational effects become less influential. The Nusselt numbers in Figure 7 are measured at times when the fluid is not fully convective, leading to $\mathcal{R}\sim 10$–$100$. This means that our results align well with the regime of more efficient convection at the equator.

In this work, we have studied the onset and evolution of semiconvection in astrophysical fluids, i.e., in fluids of small Prandtl number $\mathrm{Pr}=\nu/\kappa_{T}$, and small diffusivity ratio $\tau=\kappa_{C}/\kappa_{T}$. Unlike previous work which was limited to conducting simulations on Cartesian boxes, or spherical shells imposing artificial boundary conditions on the inner surface and potentially altering the dynamics, we have solved the Boussinesq fluid equations for the entire sphere, including $r=0$, for the first time, and including the Coriolis force. For our choice of $\mathrm{Pr}=\tau=0.1$, we varied the density ratio $R_{\rho}$ between 1.25 and 2 for the non-rotating runs, and we fixed $R_{\rho}=1.25$ while varying the convective Rossby number $\mathrm{Ro_{c}}$ between 0.25 and 1 in the rotating cases.

In all cases, we observed an initial phase of weak, instability-driven turbulence, followed by the formation of a few convective layers in the fluid (see Figures 1, 2, and 6). These layers quickly begin to merge until the sphere becomes fully mixed, with the mixing time increasing with $R_{\rho}$ due to the enhanced stability of the composition gradient, and also increasing for smaller convective Rossby numbers $\mathrm{Ro_{c}}$ due to the much smaller heat and compositional transport hampered by rotation (see measurements of the radial Nusselt numbers in Figures 3(a) and 5(a)). We observe a similar trend during the semiconvective phase (i.e. before the sphere becomes fully convective and well mixed) when comparing the buoyancy flux ratio $\gamma^{-1}$ for different $R_{\rho}$ and $\mathrm{Ro_{c}}$ (see Figures 4 and 5b). Interestingly, when testing our measurements of the radial Nusselt numbers during the semiconvective phase against theoretical predictions from previous work (Rosenblum et al., 2011; Spruit, 2013; Wood et al., 2013), we found that Spruit’s theory (Equation 10) gives a better agreement with the data. However, Rosenblum’s prediction (Equation 11) gives the right scaling with the density ratio ($[\mathrm{Nu}_{C}-1]/[\mathrm{Nu}_{T}-1]\propto 1/R_{\rho}$), matching the data with a different prefactor equal to 0.33. Therefore, more work is needed to better distinguish between the two models, especially if the prefactor depends on the boundary conditions or geometry of the system.

A noticeable difference between the nonrotating and rotating runs is the morphology of the convective flow. In the nonrotating case, convection exhibits nearly isotropic spatial scales, whereas in the rapidly rotating cases, the flows become anisotropic, with longer axial scales compared to horizontal scales (see Figures 2 and 6). These differences are also reflected in the radial transport as a function of colatitude (Figure 7). Near the poles, the radial component of the convective flow is weaker due to the columnar morphology of the fluid, which favors horizontal fluid motion. In contrast, radial transport is more efficient at the equator because the Coriolis force does not significant affect radial motions.

Semiconvection in our non-rotating simulations in the full sphere behaves qualitatively similar to previous work on non-rotating Cartesian boxes, with only minor differences with respect to the critical density ratio at which layers form. For example, in the suite of simulations by Mirouh et al. (2012), for $\mathrm{Pr}=\tau=0.1$, layer formation occurs when $R_{\rho}\approx 1.5$. In these simulations, we find layer formation even for $R_{\rho}=2$. The fluxes measured through the buoyancy flux ratio are similar. For example, Mirouh et al. (2012) measured $\gamma^{-1}\approx 0.3$–0.5 for the same values of $R_{\rho}$, Pr, and $\tau$ considered in this work (see their Table 5). We measured $\gamma^{-1}\approx$ 0.3–0.4 during the semiconvective phase. A more significant difference is in the Nusselt numbers, where Mirouh et al. (2012) measured $\mathrm{Nu}_{C}\sim 10$, $\mathrm{Nu}_{T}\sim 1$–5, while we measured $\mathrm{Nu}_{C}\sim 100$, $\mathrm{Nu}_{T}\sim 10$. These differences are not surprising given the significant differences in simulation setups. Our simulations were conducted on a full spherical domain, with gravity and fluxes increasing toward the surface. In contrast, the simulations in Mirouh et al. (2012) were performed in Cartesian boxes with constant gravity, triply periodic boundary conditions, and no external flux forcing beyond the initial profiles. However, we emphasize that the formation and evolution of the system are qualitatively similar, with no significant differences.

Comparison with previous work that includes rotation is more difficult because very little is known about rotating semiconvection. Blies et al. (2014) investigated semiconvection in rotating spherical shells, but it remains unclear whether they observed layer formation. The Rayleigh number used in their simulations was $\sim 10^{7}$, which corresponds to a domain size of about $50d$, where $d\approx(\kappa\nu/\alpha g|dT/dz|)^{1/4}$ is the characteristic scale of double-diffusive eddies. Typically, in simulations of semiconvection, the first layers that form are never thinner than 30$d$–50$d$ (Garaud, 2021). So, it is likely that Blies et al. (2014) did not find layers because of the smaller Rayleigh number (our simulations have $\mathrm{Ra}=10^{9}$, which corresponds to $\sim 180d$, so that a few layers are expected to form). Nevertheless, consistent with the findings of this study, they reported that rotation reduces both heat and compositional fluxes. When testing the theoretical predictions described earlier, they also found that Spruit’s theory provides a better approximation to the data.

An intriguing suite of simulations was conducted by Moll & Garaud (2017) on Cartesian boxes. Similar to the results presented here, they demonstrated that rotation consistently diminishes thermal and compositional transport across the system. Interestingly, in some of their simulations, where the flow was strongly constrained by rotation (low Rossby numbers), layer formation was suppressed, and large-scale vortices formed instead, enhancing compositional transport. However, their results appear to be highly sensitive to the size of the boxes adopted in the simulations, likely due to their use of triply periodic boundary conditions. More recently, Fuentes et al. (2024) conducted 3D simulations on Cartesian boxes but using fixed flux boundary conditions. They demonstrated that rotation reduces the kinetic energy flux available for mixing, and, as consequence, rotation prolongs the lifetime of semiconvective layers.

Finally, the latitudinal dependence of the fluxes in our simulations aligns well with previous studies of thermal convection in rapidly rotating shells (e.g., Wang et al., 2021; Gastine & Aurnou, 2023).

We have made a number of approximations that must be relaxed before making conclusions for stars and planets. First, we limited the simulations to fluids of $\mathrm{Pr}=\tau=0.1$. These values are at the middle of values expected to occur in planetary interiors, where the Prandlt number and diffusivity ratio may extend down to $\sim 0.01$ (e.g., Stevenson & Salpeter, 1977; French et al., 2012). In stellar interiors, even lower values are expected. For example, in non-degenerate regions of MS and RGB stars, $\mathrm{Pr}\sim 10^{-6}$ and $\tau\sim 10^{-7}$ (e.g., Garaud et al., 2015). These parameter values are beyond current computational capabilities.

Second, we adopted the Boussinesq approximation, which restricts the vertical scale to be much smaller than a density scale height. Incorporating density stratification would change the mixing rate of the layers through a combination of two effects. First, the potential energy needed for mixing would decrease. While composition is uniformly mixed within a convective layer in both cases, stratified convection adjusts the density to the adiabatic gradient, whereas Boussinesq convection maintains a constant density, which requires more energy. Second, stratification enhances the asymmetry between upflows and downflows (Meakin & Arnett, 2007), thereby increasing the net kinetic energy flux available for mixing. At present, we cannot speculate on the net impact of these effects, since, to the best of our knowledge, there are no simulations of semiconvection incorporating density stratification.

Third, we imposed fixed temperature and composition gradients at the surface of the sphere, which, despite that do not affect the layer formation process, are certainly not realistic boundary conditions for stars or planets. For example, in giant planets, the surface flux decreases over time, and there are no composition fluxes at the surface. In stars, internal heat generation and compositional changes from nuclear burning play a significant role, both of which are not included in these simulations.

Improvements in numerical modeling, focusing on adding density stratification, and the use of more suitable boundary conditions, would be of great interest to inform 1D evolution models and interpret the rich set of observations of stars and planets.