4.1. The layered phase of evolution

In this section, we study the transient layered phase in more detail, quantifying the dependence of the layers on the system parameters and considering their long-term fate. Although double-diffusive layers are well documented in periodic boxes (e.g. Rosenblum et al. Reference Rosenblum, Garaud, Traxler and Stellmach2011; Mirouh et al. Reference Mirouh, Garaud, Stellmach, Traxler and Wood2012; Wood et al. Reference Wood, Garaud and Stellmach2013), significantly less work has been done in a spherical geometry. According to the widely accepted theory of Radko (Reference Radko2003), the layering instability depends only on the density ratio $R_\rho$ , and a semi-convective fluid will develop density layers if the ratio of the composition to temperature fluxes $\gamma =F_C/F_T$ decreases as a function of the density ratio. This theory has been confirmed by numerous works on simulations in both fingering and SC regimes (e.g. Radko Reference Radko2003; Wood et al. Reference Wood, Garaud and Stellmach2013). However, this is a linear theory, based on local uniform gradient background states, without the influence of boundaries. The simulations shown in figure 2 use the same value of $R_\rho$ , but do not both produce layers, showing that in the bounded geometry that we consider, layering must also depend on the stratification or forcing.

Recently, Fuentes (Reference Fuentes2025) conducted simulations of SC in full spheres, finding that the SC instability developed into a series of convective layers, which were gradually eroded by a strong outer layer growing into the interior. In these simulations, faster rotation rates resulted in longer-lasting layers, with the outer convective layer growing towards the centre more slowly. This runs contrary to the findings of Moll & Garaud (Reference Moll and Garaud2017), who demonstrated that fast rotation suppresses layer formation. If the modified Taylor number

(4.7) \begin{equation} \textit{Ta}^*=\frac {4{\textit{Pr}}^2}{\textit{Ek}^2\textit{Ra}_T}\gg 1, \end{equation}

then the rotation is so strong that it affects even the small-scale double-diffusive instability. Fuentes et al. (Reference Fuentes, Hindman, Fraser and Anders2024) commented that this layer suppression was due to the periodic box simulations of Moll & Garaud (Reference Moll and Garaud2017) using a domain too small for layers to form. Such a domain-size constraint has been documented by Mirouh et al. (Reference Mirouh, Garaud, Stellmach, Traxler and Wood2012) and Wood et al. (Reference Wood, Garaud and Stellmach2013), who found that a domain of depth $\gtrsim 50d$ was necessary for layer formation, where

(4.8) \begin{equation} d=\left (\frac {\kappa _T\nu }{\alpha _Tg\left (\partial T_0/\partial r\right )}\right )^{1/4}, \end{equation}

represents the characteristic scale for the SC instability. By redefining the Rayleigh number in terms of the thermal gradient, rather than temperature difference,

(4.9) \begin{equation} \textit{Ra}_T=\frac {\alpha _Tg\partial T_0/\partial r R_o^4}{\kappa _T\nu } = \frac {N^2_{0,T}R_0^4}{\kappa _T\nu }, \end{equation}

we note that $\textit{Ra}_T^{1/4}(\Delta R/R_o)$ provides a measure of the shell thickness normalised by $d$ ; hence, the condition of Mirouh et al. (Reference Mirouh, Garaud, Stellmach, Traxler and Wood2012), Wood et al. (Reference Wood, Garaud and Stellmach2013) translates into a requirement that $\textit{Ra}_T\gt (50R_o/\Delta R)^4$ for layer formation. While a domain-size constraint is a potential reason for layers to be suppressed at lower $\textit{Ta}^*$ , the high- $\textit{Ta}^*$ simulations of Moll & Garaud (Reference Moll and Garaud2017) suppress the $\gamma$ -instability itself.

To study the effects of both parameters, we discuss simulations along two transects: with $\textit{Ra}_C$ fixed while $R_\rho$ varies, and vice versa. In both transects, the strength of the forcing ( $\textit{Ra}_T$ ) varies; however, the fixed- $R_\rho$ transect runs approximately parallel to the onset of SC, while the fixed- $\textit{Ra}_C$ transect spans the width from the onset of convection ( $R_\rho =1$ ) to past the onset of SC $R_\rho =({\textit{Pr}}+1)/({\textit{Pr}}+1/{\textit{Le}})$ .

4.1.1. Effect of $R_\rho$ on layer formation

To measure the effect of $R_\rho$ on layering, we show simulations with $\textit{Ra}_C=2.33\times 10^8$ fixed, and $1\leqslant R_\rho \leqslant 3.3$ , for ${\textit{Pr}}=0.3$ , $\textit{Sc}=3$ , $\Delta R=0.5$ , $\textit{Ek}=10^{-4}$ . This transect runs across the SC region, from the onset of top-heavy convection, to the stable region. From the temperature and composition equations (2.4)–(2.5), and the definitions of the background fields (2.7), we define the local radial temperature and composition fluxes as

(4.10) \begin{eqnarray} F_T^{\textit{local}} &=& \frac {1}{{\textit{Pr}}}\frac {\partial }{\partial r}\left (T_0+T\right ) + \left (T+T_0\right )u_r, \end{eqnarray}

(4.11) \begin{eqnarray} F_C^{\textit{local}} &=& \frac {1}{\textit{Sc}}\frac {\partial }{\partial r}\left (C_0+C\right ) + \left (C+C_0\right )u_r. \end{eqnarray}

Upon averaging over a spherical surface, comparison with (4.2)–(4.3) shows that

(4.12) \begin{eqnarray} \langle F_T^{\textit{local}}\rangle _S = \frac {\langle \partial T_0/\partial r\rangle _S}{{\textit{Pr}}} {\textit{Nu}}(r,t), \end{eqnarray}

(4.13) \begin{eqnarray} \langle F_C^{\textit{local}}\rangle _S = \frac {\langle \partial C_0/\partial r\rangle _S}{\textit{Sc}} \textit{Sh}(r,t). \end{eqnarray}

Hence the flux ratio may be written as

(4.14) \begin{equation} \gamma =\frac {F_C}{F_T}=\frac {\textit{Sh}}{{\textit{Nu}}}\frac {R_\rho }{{\textit{Le}}} \end{equation}.

In a pure conductive state (i.e. $\boldsymbol{u}=\boldsymbol{0}$ ), the Nusselt and Sherwood numbers are equal to unity, and hence $\gamma =R_\rho /{\textit{Le}}$ . Recall from § 3 that layering is predicted by the $\gamma$ -instability where $\partial \gamma /\partial R_\rho \lt 0$ .

As layering is a secondary instability that acts on the initial SC response, $\gamma$ must be calculated immediately after the saturation of the initial instability; measuring the fluxes after the development of layers results in spurious trends in $\gamma (R_\rho )$ due to the different evolution of the layered compared with non-layered simulations. Figure 3(a) shows how $\overline {\gamma }$ (averaged over a short period just after the saturation of the initial instability) varies as a function of $R_\rho$ across a transect of simulations. As $R_\rho$ increases from $1$ , $\overline {\gamma }$ decreases, reaching a minimum at $R_\rho \approx 1.5$ . Past this critical point, $\gamma$ increases almost linearly with $R_\rho$ , becoming very close to the conductive limit $\overline {\gamma }\sim R_\rho /{\textit{Le}}$ . This curve resembles that found in periodic box simulations by Mirouh et al. (Reference Mirouh, Garaud, Stellmach, Traxler and Wood2012) (although with quantitatively different values, likely due to the effect of rotation). Layer formation is predicted for $R_\rho \lt R_{\rho ,\textit{layer}}\approx 1.5$ . The colour of the points denotes the value of $\overline {{\textit{Re}}}$ , showing that as $R_\rho$ increases, the flow becomes weaker, and a significant drop-off in $\overline {{\textit{Re}}}$ occurs past the limit of layering. Note that the two simulations with the largest values of $R_\rho$ are past the onset of rotating SC, $R_{\rho ,{\textit{onset}}}=\textit{Ra}_C/\textit{Ra}_{T,\textit{crit}}(\textit{Ra}_C)\approx 3.17$ (i.e. are stable to both convection and SC). In these simulations, mass transfer occurs only via conduction, so they lie exactly on the conductive limit.

Figure 3.

(a) Flux ratio time averaged over a short period just after the initial instability has saturated $\overline {\gamma }$ as a function of $R_\rho$ for simulations along the transect $\textit{Ra}_C=5\times 10^9$ , $1\leqslant R_\rho \leqslant 3.3$ , with ${\textit{Pr}}=0.3$ , $\textit{Sc}=3$ , $\Delta R=0.5$ , $\textit{Ek}=10^{-4}$ . The values for $\overline {\gamma }$ are calculated (by (4.14)) during the growth phase of each simulation. The vertical dashed line marks the limit of rotating SC, $R_{\rho ,{\textit{onset}}}=3.17$ . Points are coloured by the value of $\overline {{\textit{Re}}}$ after saturation. The dashed blue line shows the limit of pure conduction $R_\rho /{\textit{Le}}$ . (b) Radial profiles of the total axisymmetric ( $m=0$ ) composition field $C_0+C_{\textit{axi}}$ , for the two simulations labelled ‘I’ and ‘II’ in (a). (c,d) Snapshots of the axisymmetric part of the composition field for the same two simulations, showing both the total field and perturbation; (e, f) snapshots of the radial and azimuthal parts of the velocity field. Snapshots are shown immediately after the initial appearance of layers. Black lines in (c–f) show the transect along which the profiles in (b) are calculated.

Figure 3. Long description

Panel A: A scatter plot shows the flux ratio time averaged over a short period just after the initial instability has saturated as a function of R rho. The vertical axis is labeled with a dimensionless parameter, and the horizontal axis is labeled R rho. Points are colored by the value of Re after saturation. The dashed blue line indicates the limit of pure conduction. Panel B: A line graph displays radial profiles of the total axisymmetric composition field C plus C axi for two simulations labeled I and II. The vertical axis is labeled C plus C axi, and the horizontal axis is labeled r. Panel C: A circular heatmap shows the axisymmetric part of the composition field for simulation I, with the total field and perturbation. The color scale ranges from 0 to 1 times 10^9. Panel D: Another circular heatmap shows the axisymmetric part of the composition field for simulation II, with the total field and perturbation. The color scale ranges from 0 to 1 times 10^7. Panel E: A circular heatmap displays the radial and azimuthal parts of the velocity field for simulation I. The color scale ranges from -2000 to 2000. Panel F: A circular heatmap shows the radial and azimuthal parts of the velocity field for simulation II. The color scale ranges from -1000 to 1000. Black lines in panels C and F indicate the transect along which the profiles in panel B are calculated.

Figure 3(b) shows radial profiles of the axisymmetric part of the composition field $C_0+C_{\textit{axi}}$ , for two simulations (one either side of $R_{\rho ,\textit{layer}}$ ). Figure 3(c–f) shows snapshots of the composition and velocity fields for the same simulations. We show the composition, rather than total density, as fast thermal diffusion (small $\textit{Pr}$ ) means that the temperature profile is rather smeared, making it more difficult to see structures in $\rho$ than in $C$ . In simulation I, where $R_\rho \lt R_{\rho ,\textit{layer}}\approx 1.5$ (figure 3 c,e), there is a clear layered structure, with three bands of relatively constant composition, separated by sharp interfaces. This structure is also visible in the radial velocity, particularly near the poles, where three different convective bands can be seen. In simulation II, by contrast, no layers are visible, and the total composition field shows little variation from the background conductive profile $C_0(r)$ . Instead, the SC instability results in a small-amplitude linear response that does not develop into layers. The weak velocities seen in figure 3(f) result in very small convective fluxes, so the flux ratio is very close to the conductive limit $\overline {\gamma }\sim R_\rho /{\textit{Le}}$ . These two snapshots are representative for all simulations with $R_\rho \lt 1.5$ and $R_\rho \gt 1.5$ respectively, so there is good agreement with the $\gamma$ -instability theory, with layers forming where $\partial \overline {\gamma }/\partial R_\rho \lt 0$ .

Note that the small departure from the conductive profile for simulation II in figure 3(b) does not mean that SC is not occurring; a linear stability analysis (see Appendix A) shows that at this value of $\textit{Ra}_C$ , the SC limit is $R_{\rho ,{\textit{onset}}}=3.17$ . For values in the range $R_{\rho ,\textit{layer}}\lesssim R_\rho \lt R_{\rho ,{\textit{onset}}}$ , layers do not develop, and instead the long-time state of the field is a small-amplitude weakly nonlinear response, with a Reynolds number (and velocities) two orders of magnitude smaller than in the layered simulation. This matches the results of Moll et al. (Reference Moll, Garaud and Stellmach2016), who found that turbulent fluxes remain small for $R_\rho \gt R_{\rho ,layer}$ .

Figure 4.

(a) Dependence of the initial number of layers $N_l$ on ${\textit{EkRa}}_T$ , for simulations along transects $R_\rho =1.2$ and $ \textit{Ek}=10^{-4}$ and $10^{-5}$ , with ${\textit{Pr}}=0.3$ , $\textit{Sc}=3$ and $\Delta R=0.5$ . Circles/crosses represent simulations where the layers eventually decay to a state with a SSL overlaying a convective zone; squares/plusses represent cases where the final state is a well-mixed region over the majority of the domain. (b–m) Snapshots of the simulations from the $\textit{Ek}=10^{-4}$ transect labelled ‘I’–‘IV’ in (a), with each column coming from a single simulation, just after layers have formed. (b–e) Axisymmetric part of the perturbation composition field $C_{\textit{axi}}$ . (f–i) Radial profiles of the total axisymmetric composition field $C_{tot}=C_{\textit{axi}}+C_0(r)$ along the white lines in (b–e), coloured by the value of $C_{\textit{axi}}$ . Dotted black lines show the conductive profile $C_0(r)$ , while the solid black line shows the final profile of $C_{tot}(r)$ after the end of the layering phase. (j–m) Axisymmetric part of the azimuthal and radial velocities.

4.1.2. Effect of $\textit{Ra}_T$ on layer formation

As discussed previously, layering in our simulations depends not only on $R_\rho$ , but on the magnitude of the thermal forcing $\textit{Ra}_T$ as well, as demonstrated by figure 2. Fuentes et al. (Reference Fuentes, Hindman, Fraser and Anders2024) suggest that this may be due to geometric constraints, with the instability scale being too large with respect to the domain depth. To investigate the dependence on the stratification, we consider simulations with fixed $R_\rho =1.2$ , and vary $\textit{Ra}_C$ and $\textit{Ra}_T$ . Thus, the magnitude of the stratification and forcing varies, but the distance from the onset is relatively constant.

Figure 4 shows how the number of layers $N_l$ depends on $\textit{Ra}_T$ along two transects with $R_\rho =1.2$ for $\textit{Ek}=10^{-4}$ and $10^{-5}$ , alongside snapshots of simulations with zero, one, two and three layers. The data from the two transects collapse onto a single trend in terms of ${\textit{EkRa}}_T$ (a parameter that is independent of viscosity). Mixed layers are identified from the composition profiles (figure 4 f–i) as broad regions where the total composition gradient is notably lower than in the conductive profile $C_0(r)$ ; where more than one layer exists, they are separated by relatively narrow, high gradient ‘interfaces’. Figure 4(j–m) shows the zonal and radial velocities. As discussed in § 4.1.1, the layers can also be observed as distinct convective bands in the profiles of $u_r$ , although this is significantly less clear than in the composition fields. Over time, the number of layers reduces through a series of mergers; however, the time scale for this process is much longer than that for the initial layer formation. As such, the initial number of layers $N_l$ can easily be identified from the snapshots.

To understand this trend in the number of layers, we consider the parameters in terms of the characteristic length $d$ and the modified Taylor number $\textit{Ta}^*$ . The condition of Mirouh et al. (Reference Mirouh, Garaud, Stellmach, Traxler and Wood2012) suggests that layering requires $\textit{Ra}_T\gt (50R_o/\Delta R)^4=10^8$ , however, for the $\textit{Ek}=10^{-4}$ transect we see layering at $\textit{Ra}_T=10^7$ , suggesting that the addition of rotation has changed the necessary length scale. According to the $\textit{Ta}^*\gg 1$ condition (4.7), the critical value of $\textit{Ra}_T$ would scale with $\textit{Ek}^2$ ; however, figure 4 shows that the number of layers actually depends on the viscosity-free quantity ${\textit{EkRa}}_T$ , with layering for ${\textit{EkRa}}_T\gtrsim 1000$ . We note that our simulations show layering up to $\textit{Ta}^*=9$ , but appear to all be in the small- $\textit{Ta}^*$ regime; the difference from the $\textit{Ta}^*=1$ cutoff of Moll & Garaud (Reference Moll and Garaud2017) is likely due to a scaling factor in the definition (4.7).

Figure 4(a) shows the onset of layering at ${\textit{EkRa}}_T\approx 1000$ , with transitions to two and three layers at $6000$ and $12\,000$ , respectively. Based on this, we identify the empirical scaling $h_{\textit{layer}}$ for layer formation

(4.15) \begin{equation} h_{\textit{layer}}\approx \frac {5}{({\textit{EkRa}}_T)^{1/3}}, \end{equation}

and note that the number of layers initially present in the domain will be the integer part of $R_o\Delta R/h_{\textit{layer}}$ .

Note the difference between the zero-layer state in figure 4(b) and that seen in figure 3(d). Here, layering is possible by the $\gamma$ -instability, so there is a full nonlinear response with convective motions visible in the velocity field, and noticeable density structure visible, but no coherent layering throughout the domain. By contrast, in figure 3(d), the system is stable to layering, with only a small-amplitude weakly nonlinear response. The composition field is rather homogeneous, with no discernible departure from the background conductive state in figure 3(d).

4.2. Long-term evolution of staircases

It is well documented that, in simulations, density staircases gradually erode through layer merger events, until only a single mixed layer remains in the domain (e.g. Wood et al. Reference Wood, Garaud and Stellmach2013; Fuentes et al. Reference Fuentes, Hindman, Fraser and Anders2024; Pružina Reference Pružina2025). In figure 4(f–i), the solid black lines represent the final profile of $C_{tot}$ after layers have decayed, showing that in our simulations, there are two possible final states: low values of $\textit{Ra}_T$ lead to a structure comprising an inner convective region overlain by a wide SSL, while higher values of $\textit{Ra}_T$ produce a single well-mixed region filling the domain, with narrow boundary layers on either side.

Figure 5 shows the long-term evolution of the radial composition profile for two simulations: one leading to each final state. In figure 5(a), for the (low- $\textit{Ra}_T$ simulation labelled ‘II’ in figure 4), the initially smooth background profile gradually sharpens, leading to a ‘one-layer’ state, with a single well-mixed layer in the interior, overlain by a stably stratified region. Over time, the mixed layer increases in size, but eventually stabilises far from the outer boundary. This state is remarkably stable for long times. By contrast, in figure 5(b), (simulation ‘IV’ in figure 4), the initial profile develops into a stack of three mixed layers separated by sharper interfaces, which merge until a single mixed layer takes up the vast majority of the domain, between narrow stably stratified boundary layers.

Figure 5.

Profiles of $C_{tot}(r)$ throughout the layered phase of two simulations, for (a) $\textit{Ra}_T = 1.15\times 10^7$ , (b) $\textit{Ra}_T = 4.2\times 10^8$ . Profiles are shown at a series of times, which can be read at the top left of each profile. The horizontal ( $C$ ) scale for each profile is the same due to the fixed- $C$ boundary conditions. (c) Time spent in layered state $t_{{layer}}$ , for simulations in the final convective zone regime.

Figure 5. Long description

Panel A: A line graph shows profiles of a variable throughout the layered phase of a simulation for a specific parameter set. The vertical axis represents the variable r, ranging from 0.6 to 1.0. The horizontal axis represents EkRaTt, ranging from 0 to 140. Multiple lines indicate profiles at different times, which can be read at the top left of each profile. Panel B: Another line graph shows profiles of a variable throughout the layered phase of a simulation for a different parameter set. The vertical axis represents the variable r, ranging from 0.6 to 1.0. The horizontal axis represents EkRaTt, ranging from 0 to 100. Multiple lines indicate profiles at different times, which can be read at the top left of each profile. Panel C: A scatter plot with a line graph shows the time spent in the layered state for simulations in the final convective zone regime. The vertical axis represents t_layer on a logarithmic scale from 10^-2 to 10^-1. The horizontal axis represents EkRaT on a logarithmic scale from 10^3 to 10^5. Blue crosses and orange plus signs represent different values of Ek, with a legend indicating Ek values of 10^-4 and 10^-5. A black line represents the trend (EkRaT)^-1.

For higher values of $\textit{Ra}_T$ , the behaviour of the layers is similar to that seen in other studies, with a gradual erosion of the layered state leading ultimately to a well-mixed state throughout the domain (e.g. Fuentes Reference Fuentes2025). However, the final state seen for lower values of $\textit{Ra}_T$ is rather different. This suggests a competition between rotational and convective processes, where the tendency of a rotating fluid to form columnar structures prevents the erosion of the stably stratified outer layer. For these low- $\textit{Ra}_T$ simulations, the layers never fully merge – rather, the wide convective layer overlain by an SSL survives without further evolution of its size – meaning that it can impact the dynamics on long time scales.

Figure 5(c) shows the survival time of the layered state, for simulations resulting in a convective zone throughout the entire domain (it is difficult to define an end of the layered phase when the final state contains a wide SSL). Larger values of ${\textit{EkRa}}_T$ result in a shorter layered phase, with an approximate scaling $t_{{layer}}\sim ({\textit{EkRa}}_T)^{-1}$ , independent of the viscosity. This agrees with the results of Fuentes et al. (Reference Fuentes, Hindman, Fraser and Anders2024), who found that higher rotation rates relative to the stratification produced longer-lasting layers, as rotational effects resisted the convection processes that cause layers to merge.

Figure 6.

(a) Dependence of the height $h_{\textit{SSL}}$ of the outer SSL on $\widetilde {{Ra}}_T$ ; (b) variation of the Nusselt number at the outer boundary as a function of the rescaled SSL height $({\textit{Pr}}/\textit{Ek})^{1/2}h_{\textit{SSL}}$ ; (c) dependence of the time-averaged Nusselt number $\overline {{\textit{Nu}}}$ on $\widetilde {{Ra}}_T$ . Points are coloured by Ekman number $\textit{Ek}$ , with different values of $R_\rho$ , $\textit{Pr}$ and $\Delta R$ represented by different symbols. Grey points in panel (a) show the width $h_{\textit{in}}$ of the inner boundary layer. The black dashed line in (a) and (c) marks a change in scaling for $h$ , corresponding with a jump in $\overline {{\textit{Nu}}}$ . The red dashed lines in all panels correspond to a transition as $(Pr/\textit{Ek})^{1/2}h_{\textit{SSL}}$ passes through unity: in (b), the horizontal line shows the value of $\textit{Ek}^{1/3}\overline {{\textit{Nu}}}$ at this point; in (c), the red shaded band shows this value of $\overline {{\textit{Nu}}}$ for the range of $\textit{Ek}$ surveyed, and the vertical line shows the value $\widetilde {{Ra}}_T\approx 300$ at which these lines cross the corresponding transects. In panel (a), the vertical red line shows that this value of $\widetilde {{Ra}}_T=300$ corresponds to a change of scaling for $h_{\textit{SSL}}(\widetilde {{Ra}}_T)$ .

To quantify the transition between the two regimes, we define the boundaries of the convective zone as

(4.16) \begin{eqnarray} R_1 &=& \min \{r | N^2(r)\lt 0\}, \end{eqnarray}

(4.17) \begin{eqnarray} R_2 &=& \max \{r | N^2(r)\lt 0\}. \end{eqnarray}

We identify the height of the inner and outer SSLs as

(4.18) \begin{equation} h_{\textit{in}}=R_1-R_i,\quad h_{\textit{SSL}}=R_o-R_2. \end{equation}

Figure 6(a) shows how $h_{\textit{SSL}}$ varies as a function of the thermal forcing, for a range of simulations with $10^{-6}\leqslant \textit{Ek}\leqslant 10^{-4}$ . The majority of these simulations take $R_\rho =1.2$ , ${\textit{Pr}}=0.3$ and $\Delta R=0.5$ , but these values were varied in a small number of simulations, while keeping ${\textit{Le}}=10$ constant. The data collapse relatively well onto a single curve with $h_{\textit{SSL}}$ controlled by the ‘reduced Rayleigh number’

(4.19) \begin{equation} \widetilde {{Ra}}_T = \textit{Ek}^{4/3}\textit{Ra}_T, \end{equation}

in which $\textit{Ra}_T$ is rescaled by the asymptotic scaling for the onset of rotating convection (Dormy et al. Reference Dormy, Soward, Jones, Jault and Cardin2004). This parameter measures the distance from the rotating convection onset, with smaller values of $\widetilde {{Ra}}_T$ being closer to the onset, and therefore more strongly rotationally constrained. It is also possible to fit the data to other rescaled parameters; for example, there is a similarly good collapse in terms of the control parameter $\textit{Ek}^{3/2}\textit{Ra}_T{\textit{Pr}}^{-1}$ , which does not depend on length. We note that a similar parameter $\textit{Ek}^{3/2}\textit{Ra}_T$ is suggested by King et al. (Reference King, Stellmach and Aurnou2012) to control the transition between rotating and non-rotating convection, for ${\textit{Pr}}\geqslant 1$ .

The grey points in figure 6(a) show the value of $h_{\textit{in}}$ , which generally takes values significantly smaller than $h_{\textit{SSL}}$ , but follows a similar trend. Three regimes are clearly visible. For $\widetilde {{Ra}}_T\lt 50$ , $h_{\textit{SSL}}/\Delta R$ is rather large, taking a value between $0.2$ and $0.7$ , and scales rather weakly, with $h_{\textit{SSL}}/\Delta R\sim \widetilde {{Ra}}_T^{-1/3}$ . Above $\widetilde {{Ra}}_T\gtrsim 300$ , the same $-1/3$ power law is seen. Between these regimes, there is a transition region where $h_{\textit{SSL}}$ varies quickly as a function of $\widetilde {{Ra}}_T$ , with significant spread in the data meaning that a dedicated study would be necessary to further classify the trend.

These transitions can be explained in terms of the dependence of $\textit{Nu}$ on $h_{\textit{SSL}}$ . In Rayleigh Bénard convection, the thermal boundary layer thickness is conventionally defined as $\delta _T=H/(2{\textit{Nu}})$ , where $H$ is the depth of the domain (see Belmonte et al. Reference Belmonte, Tilgner and Libchaber1994). Figure 6(b) shows that, in our simulations, $\overline {{\textit{Nu}}}\sim h_{\textit{SSL}}$ for $h_{\textit{SSL}}\lesssim (\textit{Ek}/{\textit{Pr}})^{1/2}$ . This cutoff corresponds to the length above which the Coriolis force acts faster than thermal diffusion. At this point, $\textit{Ek}^{1/3}\overline {{\textit{Nu}}}\approx 0.15$ , which is shown as a red shaded band on figure 6(c) for the range of $\textit{Ek}$ surveyed. The data fall in this band at $\widetilde {{Ra}}_T\approx 300$ ; the same value at which the kink is visible in figure 6(a). At these small scales, rotation has only a small effect compared with diffusion, so a scaling independent of $\textit{Ek}$ may be expected. In figure 6(a), the shift between the $\textit{Ek}=10^{-4}$ and $10^{-5}$ transects for $\widetilde {{Ra}}_T\gt 300$ is consistent with the rotation-free scaling $h_{\textit{SSL}}\sim \textit{Ra}_T^{-1/3}$ , matching that expected for a boundary layer in thermal convection (e.g. King et al. Reference King, Stellmach and Aurnou2012).

When the SSL is thicker, figure 6(b) shows a shallower slope $\overline {{\textit{Nu}}}\sim h_{\textit{SSL}}^{-1/2}$ . The transition at $\widetilde {{Ra}}_T\approx 50$ does not have a clear signature in figure 6(b), with the $\overline {{\textit{Nu}}}(h_{\textit{SSL}})$ dependence remaining the same. On the other hand, figure 6(c) shows a sharp jump of more than an order of magnitude in the value of $\overline {{\textit{Nu}}}-1$ at $\widetilde {{Ra}}_T=50$ , which can be accounted for as a transition to turbulence, occurring at the same point as the first transition in figure 6(a).

We note that a boundary layer interpretation is also possible for $\widetilde {{Ra}}_T\lt 50$ , with the scaling $h_{\textit{SSL}}\sim \widetilde {{Ra}}_T^{-1/3}$ matching that obtained for a boundary layer in rotating thermal convection (King et al. Reference King, Stellmach and Aurnou2012). This suggests the SSL obtains the maximum size of a thermally unstable layer before the onset of convection. After the initial SC instability the system stabilises into the largest possible non-convecting region, with a convective zone beneath.

To summarise, the SSL generally scales as a thermal boundary layer. When $\widetilde {{Ra}}_T\lt 50$ , the flow is non-turbulent, $h_{\textit{SSL}}$ follows the scaling $\widetilde {{Ra}}_T^{-1/3}$ for rotating thermal convection. For $\widetilde {{Ra}}_T\gt 300$ , the thermal forcing outweighs rotation, and the SSL takes a scaling for non-rotation convection $\textit{Ra}_T^{-1/3}$ . In between these limits, there is a transition range where the flow is both turbulent and influenced by rotation, and $h_{\textit{SSL}}$ quickly transitions between the non-rotating and rotating boundary layers.

Figure 7.

Time-averaged Reynolds number $\overline {{\textit{Re}}}$ plotted as a function of $\widetilde {{Ra}}_T$ for simulations with $R_\rho =1.2$ , ${\textit{Pr}}=0.3$ , $\textit{Sc}=3$ , $\Delta R = 0.5$ , for (a) $\textit{Ek}=10^{-4}$ , (b) $10^{-5}$ and (c) $10^{-6}$ . Blue points show the total Reynolds number, with orange and green showing the jet- and convection-based Reynolds numbers $\overline {{\textit{Re}}}_{\textit{jet}}$ and $\overline {{\textit{Re}}}_{{conv}}$ , respectively. Red circles mark simulations that are shown in detail in figure 8. Red and black dashed lines mark the positions of the transitions seen in figure 6, and green lines show the scaling ${\textit{Re}}_{{conv}}\sim \widetilde {{Ra}}_T$ . Black dashed line shows $\textit{Ro}=0.1$ . (d) Mixing efficiency $\overline {\chi }$ calculated in the long-term saturated state, plotted for three values of $\textit{Ek}$ .

Figure 7. Long description

Panel A: A scatter plot shows the relationship between the Reynolds number (Re) and RaT for an Ekman number (Ek) of 10^-4. The x-axis represents RaT on a logarithmic scale, and the y-axis represents Re on a logarithmic scale. Blue points indicate the total Reynolds number, orange crosses show the jet-based Reynolds number, and green plus signs represent the convection-based Reynolds number. A green line shows the scaling Re RaT. Black and red dashed lines mark specific transitions. Panel B: Another scatter plot for Ek = 10^-5, with similar axes and markers. Red circles highlight specific simulations. Panel C: A scatter plot for Ek = 10^-6, with similar axes and markers. Panel D: A scatter plot shows the mixing efficiency as a function of RaT for three values of Ek. The x-axis represents RaT on a logarithmic scale, and the y-axis represents the mixing efficiency on a logarithmic scale. Different colors represent different Ekman numbers.

4.3. Behaviour of the saturated semi-convection field

We now consider the dynamics beyond the layered phase, when the system has saturated to a statistically steady state. In § 4.2 we noted two regimes, with the post-layered state consisting of either a single convective region, or a convective zone overlain by a SSL with a mostly zonal flow. To understand the signature of these regimes on the global properties of the flow, we define Reynolds numbers based on the zonal–toroidal ( $u_{\textit{ZT}}$ ), and non-zonal–poloidal ( $u_{\textit{NZP}}$ ) components of the velocity (from the poloidal–toroidal decomposition used by XSHELLS). The radial velocities are poloidal only and the zonal jets are toroidal only. Hence, the zonal–toroidal and non-zonal–poloidal components can be used as proxies for the strength of zonal jets and convective motions, respectively,

(4.20) \begin{equation} {\textit{Re}}_{\textit{jet}} = \frac {\langle u_{\textit{ZT}}\rangle _{\textit{rms}}R_o}{\nu },\quad {\textit{Re}}_{{conv}} = \frac {\langle u_{\textit{NZP}}\rangle _{\textit{rms}}R_o}{\nu }. \end{equation}

Figure 7(a–c) shows the dependence of $\overline {{\textit{Re}}}$ , $\overline {{\textit{Re}}}_{\textit{jet}}$ and $\overline {{\textit{Re}}}_{{conv}}$ on $\widetilde {{Ra}}_T$ along the transect $R_\rho =1.2$ , for $\textit{Ek}=10^{-4}$ , $10^{-5}$ and $10^{-6}$ . In general, as the forcing increases, the flow strength also increases. In § 4.2 we noted two visible transitions. These transitions are also visible in figure 7, with a sudden jump in the value of $\overline {{\textit{Re}}}$ at $\widetilde {{Ra}}_T = 50$ , and a saturation of the value of $\overline {{\textit{Re}}}_{\textit{jet}}$ at $\widetilde {{Ra}}_T\approx 300$ . For the majority of simulations, $\overline {{\textit{Re}}}_{\textit{jet}}\gg \overline {{\textit{Re}}}_{{conv}}$ , with convective motions being more dominant only for the $\textit{Ek}=10^{-4}$ transect where $\widetilde {{Ra}}_T\gtrsim 300$ . Here, the outer SSL is thinner than $(\textit{Ek}/{\textit{Pr}})^{1/2}$ (see figure 6). At fixed $\textit{Ek}$ , there is a relatively good fit to the scaling ${\textit{Re}}_{{conv}}\sim \widetilde {{Ra}}_T$ , particularly for the smaller values of $\textit{Ek}$ . Hereafter, we will refer to the regime with $\widetilde {{Ra}}_T\lt 300$ as ‘jet dominated’, and $\widetilde {{Ra}}_T\gt 300$ as ‘convection dominated’.

Figure 7(d) shows the mixing efficiency in the convective zone, calculated (following Vidal et al. Reference Vidal, Cébron, Schaeffer and Hollerbach2018) as

(4.21) \begin{equation} \chi =1-\frac {\int _{R_1}^{R_2}\langle N^2\rangle _{S(r)} \mathrm{d}r}{\int _{R_1}^{R_2}\langle N^2_0\rangle _{S(r)} \mathrm{d}r}, \end{equation}

where $\langle \boldsymbol{\cdot }\rangle _{S(r)}$ is the mean over the spherical surface at radius $r$ , and $R_1$ and $R_2$ are defined by (4.16)–(4.17). This quantifies how well mixed the convective zone is, with $\chi =0$ when there is no departure from the background density profile, and $\chi =1$ when the convective zone is fully mixed to a constant density. Figure 7(d) shows that, in the jet-dominated regime, $\overline {\chi }$ takes small values, but increases as a function of $\widetilde {{Ra}}_T$ . By contrast, $\overline {\chi }$ is large and relatively constant in the convective regime. This means that not only is the convective zone larger for values of $\widetilde {{Ra}}_T\gt 50$ , but the convection itself is also more efficient for mixing.

To quantify the effect of rotation on the flow, it is common to define the Rossby number

(4.22) \begin{equation} \textit{Ro}=\textit{Ek}{\textit{Re}}. \end{equation}

When $\textit{Ro}\lesssim 0.1$ , flows are considered to be rotationally constrained; above this limit, rotational effects are much less important. In figures 7(a) and 7(b), $\textit{Ro}=0.1$ is marked with horizontal dotted lines; we see that for both values of $\textit{Ek}$ , as the forcing increases, ${\textit{Re}}_{\textit{jet}}$ remains close to $\textit{Re}$ until this cutoff is reached, but for $\textit{Ro}\gt 0.1$ , the jet-based Reynolds number reaches a plateau. This cutoff corresponds relatively well with the jet–convection transition already identified at $\widetilde {{Ra}}_T\approx 300$ . We note that the other transition, at $\widetilde {{Ra}}_T\approx 50$ , occurs for $\textit{Ro}=0.03$ . For the smallest value of $\textit{Ek}$ surveyed, all simulations are strongly rotationally constrained, with a maximum value of $\textit{Ro}\approx 0.01$ .

Figure 7 does not demonstrate a unified behaviour across different values of $\textit{Ek}$ and the whole range of $\widetilde {{Ra}}_T$ ; we shall discuss how data collapse onto a single master curve later in § 5.

The transitions at $\widetilde {{Ra}}_T=50$ and $300$ are most obvious in the $\textit{Ek}=10^{-4}$ transect, although are also present in the simulations with $\textit{Ek}=10^{-5}$ . To look deeper into the nature of these transitions, figure 8 shows snapshots of the velocity and temperature fields, for four simulations labelled by their position along the $\textit{Ek}=10^{-5}$ transect (cf. figure 7 b).

For $\widetilde {{Ra}}_T\lt 50$ (labelled ‘4’ in figure 7 b), the flow is generally weak and laminar, with a prograde equatorial jet, retrograde flow at higher latitudes and some wave-like structures along the shear layer. The flow here is strongly rotationally constrained, with $\textit{Ro}=0.006$ . The density perturbation is also rather weak, with a stable stratification retained throughout almost the entire domain, except a small patch near the equator. For simulation 9, with $\widetilde {{Ra}}_T\gt 50$ , the flow is turbulent, with figures 8(b), 8(f) and 8(j) showing much more intense flows and greater disruption of the temperature field than in simulation 4. The effect of rotation is still very strong, with $\textit{Ro}=0.03$ . Figure 8(n) shows that, for simulation 9, the zonal flow has a three-jet structure, with prograde flow at the poles and near the equator, and retrograde flow in between. The density field is separated into two clear layers; an inner convective zone with $N^2_{\textit{axi}}\lt 0$ , and an outer SSL with $N^2_{\textit{axi}}\gt 0$ .

For simulation 12, for a higher value of $\widetilde {{Ra}}_T$ , but still beneath the jet–convection transition, the only major difference from simulation 9 is that the outer SSL is significantly narrower.

Finally, for $\widetilde {{Ra}}_T\gt 300$ , simulation 17 shows very strong convective flows. Figure 8(l) shows that the zonal jets are rather broken up compared with the other simulations, and figure 8(p) shows that the zonal flow is on average weaker and the shear layers more smeared. Here, the effect of rotation is much weaker than the other simulations, with $\textit{Ro}=0.11$ . The profile of $N^2_{\textit{axi}}$ shows that the convective zone fills the entire domain, with no SSL visible anywhere. Note that, while the profile of $N^2$ appears rather noisy, the colour bar has been chosen to highlight changes of sign (to clearly show the border of the convective and stable layers). In a well-mixed state, $N^2$ is dominated by small-scale fluctuations around zero, resulting in the noisy profile visible in figure 8(p).

Figure 8.

Snapshots of the four simulations marked by red circles in figure 7(b). Each column shows results from a single simulation (labelled at the top). Each row shows, respectively, equatorial slices of the temperature perturbation $T$ , radial velocity $u_r$ , azimuthal velocity $u_\phi$ and meridional slices of the axisymmetric ( $m=0$ ) part of the azimuthal velocity $u_\phi ^{\textit{axi}}$ and axisymmetric squared Brunt–Väisälä frequency $N^2_{\textit{axi}}$ .

Figure 9.

Quantities describing the strength of convection for simulations along the transect $R_\rho =1.2$ : (a) time-averaged Nusselt and Sherwood numbers on the outer boundary $r=R_o$ , as a function of $\widetilde {{Ra}}_T$ . (b) Time-averaged thermal, compositional and total convective power. Vertical dashed lines denote the positions of the transitions identified in figure 6. (c) Nusselt and Sherwood numbers plotted against convective Reynolds number ${\textit{Re}}_{{conv}}$ , scaled by $\textit{Ek}^{1/3}$ . (d) Convective power plotted against $\textit{Ek}^{-1/3}{\textit{Re}}_{{conv}}$ .

Figure 9. Long description

Panel A: A scatter plot shows the time-averaged Nusselt and Sherwood numbers on the outer boundary as a function of RaT. The x-axis is labeled RaT and the y-axis is labeled with a logarithmic scale ranging from 10^-4 to 10^1. Different symbols represent different Ekman numbers (Ek = 10^-4, 10^-5, 10^-6). The plot includes lines indicating slopes of 5/4 and 3/4. Panel B: A scatter plot displays the time-averaged thermal, compositional, and total convective power as a function of RaT. The x-axis is labeled RaT and the y-axis is labeled with a logarithmic scale ranging from 10^0 to 10^6. Different symbols represent different Ekman numbers. Vertical dashed lines denote transitions. Panel C: A scatter plot shows Nusselt and Sherwood numbers plotted against the convective Reynolds number, scaled by Ek^1/4. The x-axis is labeled Ek^1/4Re_conv and the y-axis is labeled with a logarithmic scale ranging from 10^-3 to 10^1. Different symbols represent different parameters. A line indicating a slope of 4/3 is included. Panel D: A scatter plot depicts convective power plotted against Ek^-1/5Re_conv. The x-axis is labeled Ek^-1/5Re_conv and the y-axis is labeled with a logarithmic scale ranging from 10^7 to 10^13. Different symbols represent different Ekman numbers. A line indicating a slope of 5/2 is included.

To quantify the turbulent transport in the flow, we show the trend in the Nusselt and Sherwood numbers on the outer boundary, as a function of $\widetilde {{Ra}}_T$ . Figure 9(a) shows that, at all points, $\overline {\textit{Sh}}-1$ is larger than $\overline {{\textit{Nu}}}-1$ (as the diffusive transport is naturally lower for the slower-diffusing scalar). For smaller values of $\widetilde {{Ra}}_T$ , $\overline {{\textit{Nu}}}$ and $\overline {\textit{Sh}}$ are both close to unity, representing very low levels of convection. As $\widetilde {{Ra}}_T$ increases, $\overline {{\textit{Nu}}}$ and $\overline {\textit{Sh}}$ increase, scaling as $(\overline {{\textit{Nu}}}-1),(\overline {\textit{Sh}}-1)\sim \textit{Ek}^{-1/6}\widetilde {{Ra}}_T^{5/4}=\textit{Ek}^{3/2}\textit{Ra}_T^{5/4}$ . The influence of rotation here is relatively small, and the power $-1/6$ is not exact; a relatively good fit is found for $\textit{Ek}^{-0.16\pm 0.1}$ . For comparison, Wood et al. (Reference Wood, Garaud and Stellmach2013) found $\overline {{\textit{Nu}}}-1\sim \textit{Ra}_T^{1/3}$ for non-rotating, layered SC. For rapidly rotating thermal convection, Gastine et al. (Reference Gastine, Wicht and Aubert2016) found the upper limit on the power law to be $\overline {{\textit{Nu}}}\sim \textit{Ek}^2\textit{Ra}_T^{3/2}/{\textit{Pr}}^{1/2}$ , while near to the onset of rotating convection $\textit{Ra}_T^c$ , the scaling was $\overline {{\textit{Nu}}}-1\sim \textit{Ra}_T/\textit{Ra}_T^c-1$ . The $\textit{Ra}_T^{5/4}$ law that we obtain sits between these two limits.

The transition at $\widetilde {{Ra}}_T=50$ is clearly visible in figure 9(a), with a marked jump of a factor of ten in both $\overline {{\textit{Nu}}}-1$ and $\overline {\textit{Sh}}-1$ , and a change in the scaling law to $\textit{Ek}^{1/6}(\overline{\textit{Nu}}-1)\sim \widetilde {{Ra}}_T^{3/4\pm 1/4}$ . In the fingering regime, Monville et al. (Reference Monville, Vidal, Cébron and Schaeffer2019) found that $\overline {{\textit{Nu}}}-1$ always remained small, meaning that the thermal transport was dominated by diffusion, allowing small-scale fingering patterns to develop. Here, by contrast, $\overline {{\textit{Nu}}}-1$ achieves values larger than unity due to the vigorously convecting layer for larger values of $\widetilde {{Ra}}_T$ .

The trend in the convective power $P=P_T+P_C$ against $\widetilde {{Ra}}_T$ , as well as the components $P_T$ and $P_C$ , is shown in figure 9(b). The data collapse to a master curve $\textit{Ek}^{4/3}P\sim \widetilde {{Ra}}_T^2$ at low values of $\widetilde {{Ra}}_T$ , with the slope shallowing to $\widetilde {{Ra}}_T^{3/2}$ above the transition to turbulence at $\widetilde {{Ra}}_T=50$ . Upon closer inspection, it can also be seen that, for $\widetilde {{Ra}}_T\lt 50$ , the total power $\overline {P}=\overline {P_T}+\overline {P_C}$ is very close to $\overline {P_T}$ , with a small contribution from $\overline {P_C}$ . However, for $\widetilde {{Ra}}_T\gt 50$ , $-\overline {P_C}$ becomes closer to $\overline {P_T}$ , leading to a relative reduction in $\overline {P}$ compared with $\overline {P_T}$ . This matches the description of a convection-dominated flow, where both thermal and compositional fields are strongly disrupted by fluid motions, giving a value of $\overline {P_C}$ that is close to $-\overline {P_T}$ . Notably, the transition at $\widetilde {{Ra}}_T=300$ is not visible in either the Nusselt/Sherwood numbers or the convective powers, similarly to what was observed in the flow snapshots.

Figures 9(c) and 9(d) show, respectively, the dependence of the Nusselt and Sherwood numbers, and the convective power, on the convective Reynolds number ${\textit{Re}}_{{conv}}$ . Fairly simple scalings are visible, with $\overline {{\textit{Nu}}}-1$ and $ \overline {\textit{Sh}}-1\sim (\textit{Ek}^{1/4}\overline {{\textit{Re}}}_{{conv}} )^{4/3}$ , and $\overline {P}\sim (\textit{Ek}^{-1/5}\overline {{\textit{Re}}}_{{conv}} )^{5/2}$ . These simple scalings, based on ${\textit{Re}}_{{conv}}$ , rather than $\textit{Re}$ , and in particular showing no significant change at the jet–convection transition, confirm that the majority of turbulent transport is done by convection, with a smaller contribution from the zonal jets.

5.1. Dynamo action: onset and self-sustained regime

Figure 10(a) summarises a number of simulations for a range of $\textit{Ek}$ , $\textit{Ra}_T$ , $\textit{Ra}_C$ and $\textit{Pm}$ , showing the value of $\textit{Rm}$ achieved, and whether or not a dynamo developed. It is clear that in general, larger values of $\overline {\textit{Rm}}$ produce dynamo action; although the onset varies strongly with $\textit{Ek}$ and $\textit{Pm}$ . However, defining the magnetic Reynolds number based on the full domain neglects the fact that the convective flow is rather weak in the outer SSL. We instead define a convective magnetic Reynolds number.

(5.3) \begin{equation} \textit{Rm}_{{conv}}= \textit{Pm}{\textit{Re}}_{{conv}}\left (\frac {R_o^3-R_i^3}{R_2^3-R_1^3}\right )^{1/2}\frac {R_2-R_1}{R_o-R_i}. \end{equation}

where $R_1$ and $R_2$ are the inner and outer boundaries of the convective zone as defined by (4.16)–(4.17), and $R_i$ and $R_o$ are the inner and outer boundaries of the domain (such that $h_{\textit{SSL}}=R_o-R_2$ ). This new parameter is based on an assumption that the convective part of the kinetic energy is contained entirely in the convective region, introducing a factor of $ ((R_o^3-R_i^3)/(R_2^3-R_1^3) )^{1/2}$ into the root -mean-square velocity. The factor of $(R_2-R_1)/ (R_o-R_i)$ is due to rescaling the length in the definition of $\textit{Re}$ .

Figure 10(b) shows the same data as panel (a), plotted with this new convective magnetic Reynolds number $\overline {\textit{Rm}}_{{conv}}$ . Here, there is a more consistent trend in the onset of dynamo action across the values of $\textit{Ek}$ , confirming that the dynamo is only active in the convective zone, and not in the SSL. For the range of $\textit{Ek}$ studied, the critical value of $\textit{Rm}_{{conv}}$ for dynamo action scales approximately with $\textit{Ek}^{1/3}$ . This is a favourable regime for a geophysical setting, suggesting only small $\textit{Rm}$ may be necessary for dynamo when planetary values of $\textit{Ek}$ are considered. However, we note that the trend cannot hold to very low values of $\textit{Ek}$ , otherwise the critical value of $\textit{Rm}$ could become arbitrarily small.

Figure 10.

Summary of a large number of dynamo simulations, for ${\textit{Pr}}=0.3$ , $\textit{Sc}=3$ , $R_\rho =1.2$ and a range of values of $\textit{Ek}$ , $\textit{Ra}_T$ , $\textit{Ra}_C$ and $\textit{Pm}$ , with (a) $\overline {\textit{Rm}}$ and (b) $\overline {\textit{Rm}}_{{conv}}$ plotted against $\textit{Ek}$ . Dots and crosses (coloured by $\textit{Pm}$ ) designate simulations with successful and unsuccessful dynamos, respectively. Black dashed line in (b) represents an approximate scaling for the onset of dynamo.

Our simulations exist primarily in the weak field dynamo regime. Figure 11(a) shows the Rossby number $\textit{Ro}=\textit{Ek}{\textit{Re}}$ against the reduced Rayleigh number $\widetilde {{Ra}}_T$ , coloured by the ratio $E_b/E_u$ (with hydrodynamic simulations shown as black crosses). The maximum value of $E_b/E_u$ is $7.4$ , with all but $4$ strong field simulations having $E_b/E_u\lt 1$ , and $60\,\%$ of the simulations having $E_b/E_u\lt 0.1$ . For reference, we expect a value in the range $10^{-4}\lt E_b/E_u\lt 1$ for Jupiter (Pružina et al. Reference Pružina, Cébron and Schaeffer2025). In all cases, the Rossby number is not significantly affected by the magnetic field, with the MHD simulations broadly collapsing onto the same curve as the hydrodynamics simulations. Figure 11(b) shows $\textit{Ro}$ against the convective Rossby number $\textit{Ro}_C = \textit{Ek}\sqrt {\textit{Ra}_T/{\textit{Pr}}}$ , coloured by $h_{\textit{SSL}}$ , which collapses better than in figure 11(a). The rotationally constrained ( $\textit{Ro}\lt 0.1$ ) regime characteristic of planets is associated with a wide SSL (and small $\widetilde {{Ra}}_T$ and $\textit{Ro}_C$ ), although all simulations except one produced $\textit{Ro}\lt 1$ .

Figure 11.

Rossby number $\textit{Ro}=\textit{Ek}{\textit{Re}}$ plotted against (a) reduced Rayleigh number $\widetilde {{Ra}}_T$ , coloured by the ratio of the magnetic to kinetic energies $E_b/E_u$ , and (b) convective Rossby number $\textit{Ro}_C$ , coloured by the size of the SSL. In (a), black crosses mark hydrodynamic simulations, while black and red dashed lines show the position of the transitions at $\widetilde {{Ra}}_T=50$ and $300$ , respectively. Black dotted lines show $\textit{Ro}=0.1$ ; below this value the flow is considered to be rotationally constrained.

Figure 11. Long description

Two scatter plots illustrate the relationship between Rossby number and other variables in convective fluids. Panel A: A scatter plot shows the Rossby number (Ro) on the vertical axis and the reduced Rayleigh number (RaT) on the horizontal axis. The data points are colored by the ratio of magnetic to kinetic energies (Eb/Eu). Black crosses mark hydrodynamic simulations, while black and red dashed lines indicate transitions at specific values. Black dotted lines show Ro = 0.1, below which the flow is considered rotationally constrained. Panel B: Another scatter plot displays the Rossby number (Ro) on the vertical axis and the convective Rossby number (RoC) on the horizontal axis. The data points are colored by the size of the SSL (h_SSL/ΔR).

The majority of our simulations were carried out with $\Delta R=0.5$ , ${\textit{Pr}}=0.3$ and $\textit{Sc}=3$ . While we conducted a small number of simulations with different values of $\Delta R$ and $\textit{Pr}$ , the Lewis number ${\textit{Le}}=\textit{Sc}/{\textit{Pr}}=10$ was kept constant. While this value is not far from the relevant value for gaseous planets (e.g. ${\textit{Le}} \approx 30$ in Jupiter, see Pružina et al. Reference Pružina, Cébron and Schaeffer2025), we cannot exclude a Lewis number dependence in the curves seen in figure 11. This moderate value of ${\textit{Le}}=10$ also somewhat limits the region of parameters where SC can be studied (see (3.10)); a larger value of ${\textit{Le}}$ would increase the range of $R_\rho$ allowed in simulations from $1 \leqslant R_\rho \leqslant 3.25$ in this paper to $1 \leqslant R_\rho \lesssim 10$ for Jupiter.

5.2. Dipolarity: effect of flow structure

The magnetic fields of giant planets such as Jupiter and Saturn, as well as the Earth, are strongly dipolar, with Saturn’s field also being very axisymmetric. As such, it is of considerable interest to identify trends in the form of the magnetic field generated by the semi-convective dynamo. We quantify the dipolarity and axisymmetry using the dipole and axisymmetry fractions

(5.4) \begin{eqnarray} f_{\textit{dip}} &=& \sqrt {\frac {\sum _{m=0}^1\left (2-\delta _{m0}\right )\left |\boldsymbol{b}\left (l=1,m,t\right )\right |^2}{\sum _{l\geqslant 1,m\geqslant 0}\left (2-\delta _{m0}\right )\left |\boldsymbol{b}\left (l,m,t\right )\right |^2}}, \end{eqnarray}

(5.5) \begin{eqnarray} f_{{\textit{axi}}} &=& \sqrt {\frac {\sum _{l\geqslant 0}\left (2-\delta _{m0}\right )\left |\boldsymbol{b}\left (l,m=0,t\right )\right |^2}{\sum _{l\geqslant 1,m\geqslant 0}\left (2-\delta _{m0}\right )\left |\boldsymbol{b}\left (l,m,t\right )\right |^2}}, \end{eqnarray}

where $l$ and $m$ are spherical harmonic degree and order, $\boldsymbol{b}(l,m,t)$ is the $(l,m)$ component of the magnetic field at the surface of the simulation $r=R_o$ and $\delta _{ij}$ is the Kronecker delta function. Note that these definitions are not totally standardised: it is common to define $f_{\textit{dip}}$ using a spectrum truncated at the maximum $l$ of the current models of the target planet’s magnetic field, while $f_{{\textit{axi}}}$ is not a widely used measure.

Figure 12.

(a) Rescaled squared Reynolds number ratio $(\overline {{\textit{Re}}_{\textit{jet}}}/\overline { {\textit{Re}}_{{conv}}})^2$ , and dipolar and axisymmetric fractions $\overline {f_{\textit{dip}}}$ , $\overline {f_{{\textit{axi}}}}$ , along a transect of dynamo simulations with $\textit{Pm}=1$ , ${\textit{Pr}}=0.3$ , $\textit{Sc}=3$ , $\Delta R=0.5$ , $\textit{Ek}=10^{-4}$ , $R_\rho =1.2$ and $6.2\times 10^6\geqslant \textit{Ra}_T\geqslant 2.6\times 10^8$ . (b–g) Snapshots of the surface radial magnetic field $b_r$ , axisymmetric part of the radial velocity field $u_r$ , axisymmetric part of the azimuthal velocity $u_\phi$ and squared Brunt–Väisälä frequency $N^2$ , for the two simulations indicated in black boxes in (a). (h–k) Magnetic energy spectra in spherical harmonic degree $l$ and order $m$ for each simulation, at a range of radii from the deep interior to the surface.

Figure 12(a) shows how the dipole and axisymmetry fractions vary along a transect of dynamo simulations with $\textit{Pm}=1$ , as well as the ratio of the squared jet/convective Reynolds numbers (representing the ratio of jet/convective kinetic energies). It appears that high values of $\overline {f_{\textit{dip}}}$ and $\overline {f_{{\textit{axi}}}}$ are associated with the jet-dominated regime. Figure 12(b–g) shows snapshots of the magnetic field on the surface, the axisymmetric part of the velocity, and the squared Brunt–Väisälä frequency, for two simulations, one in the jet-dominated regime, and one where convection dominates. In the first simulation, the surface magnetic field is clearly strongly dipolar, with a clear partition between the two hemispheres. The velocity field is dominated by strong jets, with $u_r$ significantly smaller in magnitude than $u_\phi$ . Figure 12(f) shows a wide outer SSL. This simulation closely resembles that considered by Pružina et al. (Reference Pružina, Cébron and Schaeffer2025). By contrast, in the second simulation, the magnetic field (figure 12 c) displays much more fine structure, with large contributions from higher- $l$ modes. The velocity jets are less columnar, and the radial velocity is closer in magnitude to the azimuthal velocity, showing a much higher proportion of energy in convection compared with jets than in figure 12(d). In figure 12(g), the convective zone takes up the entire domain, with much finer-scale density structures visible.

Figure 12(h–k) shows the magnetic energy spectra at a range of radii for each simulation; in both cases, the spectrum is rather flat in the interior, and does not change significantly in the convective region. Further out, in the SSL, the spectrum decays with radius, with the dipole ( $l=1$ ) component decaying at the slowest rate. In the jet-dominated regime, a wide SSL means that the field has decayed significantly and become strongly dipolar by the time it reaches the surface. By contrast, in the convection-dominated regime, the SSL is very thin, so the structure of the field at the surface is rather similar to that deep inside the shell. The effect of radius on the $m$ -spectra is similar, with the $m\gt 0$ modes decaying more rapidly than $m=0$ as radius increases, resulting in an axisymmetric-dominated field at the surface.

5.2.1. A recipe for planetary magnetic fields?

It has been proposed (e.g. Stevenson Reference Stevenson1982; Cao et al. Reference Cao, Russell, Christensen, Dougherty and Burton2011) that the magnetic field of Saturn is strongly axisymmetric and dipolar because of a SSL overlaying a convective zone where the field is generated. Simulations with an imposed SSL above the dynamo region confirm that this is possible (Yadav et al. Reference Yadav, Cao and Bloxham2022). Our simulations of SC show that this structure can be generated spontaneously by SC. With this in mind, we can identify the correct regions of parameter space to produce planet-like magnetic fields. A thick SSL is needed to filter the field, so a small-as-possible value of $\widetilde {{Ra}}_T$ is preferable. However, velocities must also be large enough to generate a dynamo, so $\widetilde {{Ra}}_T$ must also be large enough to provide these turbulent flows. As such, a value $\widetilde {{Ra}}_T$ just greater than $50$ should be ideal. Figure 6 showed that the scaling for $h/\Delta R$ did not depend significantly on domain depth; therefore, a larger SSL can also be obtained by simulating a thicker shell.

Figure 13.

Results from a simulation with $\textit{Ek}=10^{-5}$ , $\textit{Ra}_C=4.2\times 10^9$ , $\textit{Ra}_T=3.5\times 10^8$ , ${\textit{Pr}}=0.3$ , $\textit{Sc}=3$ , $\textit{Pm}=1$ , $\Delta R=0.8$ , showing (a) surface magnetic field $b_r$ , (b) surface field when the fluid in $0.7\leqslant r \leqslant 1$ is insulating and axisymmetric parts of the (c) magnitude of the magnetic field $|\boldsymbol{b}|$ , and radial magnetic field $|b_r|$ (d) square of the Brunt–Väisälä frequency $N^2$ , (e) radial and azimuthal velocity, (f) $l$ -spectrum and (g) $m$ -spectrum of the magnetic energy at a range of radii, with the surface spectra for the insulating outer layer case shown as blue dashed lines. The boundary between convective and stably stratified regions, and the inner boundary of the fluid domain, are marked on the polar plots in white/black.

Figure 13. Long description

Panel A: Two ellipsoid visuals show the surface magnetic field with a color gradient ranging from -40 to 40. Panel B: Another ellipsoid visual displays the surface field when the fluid is insulating, using a similar color gradient. Panel C: A polar plot illustrates the axisymmetric part of the magnetic field magnitude, with a color scale from 10^-1 to 10^3. Panel D: A polar plot shows the square of the Brunt-Väisälä frequency, with a color scale from -2.5 to 7.5 x 10^20. Panel E: A polar plot depicts radial and azimuthal velocity, with a color scale from -500 to 500. Panel F: A line graph shows the l-spectrum of magnetic energy at various radii, with different colored lines representing different radii and blue dashed lines for the surface spectra in the insulating outer layer case. Panel G: Another line graph displays the m-spectrum of magnetic energy at various radii, with similar color coding as Panel F. The boundary between convective and stably stratified regions, and the inner boundary of the fluid domain, are marked on the polar plots in white and black.

Based on these criteria, we present the results of a simulation for $\textit{Ek}=10^{-5}$ , $\textit{Ra}_C=4.2\times 10^9$ , $\textit{Ra}_T=3.5\times 10^8$ , ${\textit{Pr}}=0.3$ , $\textit{Sc}=3$ , $\textit{Pm}=1$ , in a shell with $\Delta R=0.8$ . Note that these parameter values are rather moderate, so the simulation is not particularly computationally challenging, with a resolution of only $(N_R,N_l,N_m)=(256,200,200)$ in the radial direction and spherical harmonic $l$ and $m$ . These values are chosen such that $\widetilde {{Ra}}_T=75$ ; just beyond the transition. Figure 13(a) shows the surface radial magnetic field for this simulation; the field is very strongly dipolar and axisymmetric. This field is also rather strong, compared with the majority of our simulations, with an energy ratio of $E_b/E_u = 1.52$ . For comparison, figure 13(b) shows the surface field if the fluid in the layer $r_L\leqslant r\leqslant 1$ is insulating (calculated by setting the toroidal component of $\boldsymbol{b}$ to zero for $r\gt r_L$ , and the imposing poloidal component $P(r) = P(r_L) \times (r/r_L)^{-l-1}$ , with $r_L=0.7$ ). This ‘truncated’ field is noticeably less dipolar and axisymmetric than the full field, and in general shows larger-scale structures.

Figure 13(c–e) shows the axisymmetric ( $m=0$ ) part of the magnitude of the total and radial magnetic fields, squared Brunt–Väisälä frequency and azimuthal and radial velocity fields. The black line shows the boundary between the convective zone and SSL. Inside the convective zone, the radial and azimuthal velocities are of a similar magnitude, and there is a strong magnetic field; in the SSL, the radial velocity becomes small, and the strength of the magnetic field drops (particularly the radial component).

Figure 13(f–g) shows spectra of the magnetic energy at a range of radii; deep in the interior, the $l$ and $m$ spectra are both rather flat, with no variation with depth. Closer to the surface, the power drops off sharply with radius, resulting in a surface field with very strongly dominant dipole ( $l=1$ ) and axisymmetric ( $m=0$ ) components, with $\overline {f_{\textit{dip}}}=0.95$ and $\overline {f_{{\textit{axi}}}}=0.99$ . Saturn’s magnetic field has dipole and axisymmetric fractions of $f_{\textit{dip}}\approx 0.99$ and $f_{{\textit{axi}}}\approx 0.985$ (Cao et al. Reference Cao, Dougherty, Hunt, Provan, Cowley, Bunce, Kellock and Stevenson2020), so our simulation compares rather favourably. The blue dashed lines in figure 13(f–g) show the surface power spectra of the truncated magnetic field (with an insulator for $r_L\leqslant r\leqslant 1$ ). Compared with the full spectra, these profiles show a slower decay at low values of $l$ and $m$ , with values of $\overline {f_{\textit{dip}}}=0.87$ and $\overline {f_{{\textit{axi}}}}=0.92$ – markedly smaller than for the full magnetic field. This greater dominance of $l=1$ and $m=0$ of the full field demonstrates that the flow in the SSL actively filters the magnetic field, rather than simply acting as an insulating layer where the field decays (as discussed by, e.g. Stevenson Reference Stevenson1982).