Deep model simulation of polar vortices in gas giant atmospheres

A homogeneous fluid with density $\rho$ and constant physical properties –thermal diffusivity $\kappa$, thermal expansion coefficient $\alpha$, and dynamic viscosity $\mu$– is considered, to study thermal convection in a rotating spherical shell. The latter is defined by the inner and outer radius $r_{i}$ and $r_{o}$, and a constant angular velocity ${\bf\Omega}=\Omega{\mathbf{k}}$ about the vertical axis. Gravity acts in the radial direction ${\bf g}=-\gamma\mathbf{r}$ ($\gamma$ is constant and $\mathbf{r}$ the position vector) and the relation $\rho=\rho_{0}(1-\alpha(T-T_{0}))$ is assumed (the Boussinesq approximation), but just in the gravitational term. For the other terms appearing in the equations a reference state $(\rho_{0},T_{0})$ is considered (see for instance Pedlosky 1979; Chandrasekhar 1981).

At the boundaries, which are perfectly conducting, a temperature difference is imposed $\Delta T=T_{i}-T_{o}$, $T(r_{i})=T_{i}$ and $T(r_{o})=T_{o}$. Stress-free boundary conditions for the velocity field are appropriate for planetary atmospheres (Christensen, 2002; Heimpel et al., 2005). The mass, momentum and energy equations are written in the rotating frame of reference as in Simitev & Busse (2003) and expressed in terms of velocity (${\mathbf{v}}$) and temperature ($\Theta=T-T_{c}$) perturbations of the basic conductive state ${\mathbf{v}}=0$ and $T_{c}(r)=T_{0}+\eta d\Delta T(1-\eta)^{-2}r^{-1}$, $\eta=r_{i}/r_{o}$ being the aspect ratio, $d=r_{o}-r_{i}$ being the gap width, and $T_{0}=T_{i}-\Delta T(1-\eta)^{-1}$ being a reference temperature. With units $d$ for the distance, $\nu^{2}/\gamma\alpha d^{4}$ for the temperature, and $d^{2}/\nu$ for the time, the governing equations read

being $p^{*}$ is a dimensionless scalar containing all the potential forces. Centrifugal effects are neglected by assuming $\Omega^{2}/\gamma\ll 1$. This is usual for the case of spherical shell convection (Schubert & Zhang, 2000). Four non-dimensional parameters -the aspect ratio $\eta$, the Rayleigh ${\rm Ra}$, Prandtl $\Pr$, and Taylor ${\rm Ta}$ numbers- summarise the physics of the problem. The Prandtl number characterises the relative importance of viscous (momentum) diffusivity to thermal diffusivity, the Taylor number is the ratio between rotational and viscous forces, and the Rayleigh number ${\rm Ra}$ is associated with buoyancy forces and thermal forcing. The parameters are defined by

To solve the model equations 1-3 with the prescribed boundary conditions a pseudo-spectral method (Tilgner, 1999) is used (see Garcia et al. 2010 for a detailed description). In the radial direction a collocation method on a Gauss–Lobatto mesh (Sánchez et al., 2016) is considered and spherical harmonics are employed in the angular coordinates. In the Boussinesq approximation the toroidal/poloidal decomposition for the velocity field can be used (Chandrasekhar, 1981). The code is parallelized in the spectral and in the physical space using OpenMP directives. Optimized libraries (FFTW3, Frigo & Johnson, 2005) for the fast Fourier transforms FFTs in longitude and matrix-matrix products (dgemm GOTO, Goto & van de Geijn, 2008) for the Legendre transforms in latitude are implemented for the computation of the nonlinear terms. Aliasing is properly removed in the pseudo-spectral transform method (Orszag, 1970).

High order implicit-explicit backward differentiation formulas (IMEX–BDF, Garcia et al., 2010) are used for the time integration of the discretized equations. The nonlinear terms are considered explicitly, to avoid solving nonlinear equations at each time step. However, the Coriolis term is treated fully implicitly and this allows a time integration with larger time steps (Garcia et al., 2010).

Two different models of Boussinesq thermal convection are considered to explore the transition from Saturn-like (Model 1) to Jupiter-like (Model 2) polar dynamics. In these two models the parameters $\eta$, ${\rm Ta}$ and $\Pr$ are the same. They are similar to previous numerical models of gas giant atmospheres (Heimpel et al., 2005; Heimpel & Aurnou, 2007; Heimpel et al., 2015) focusing on the understanding of jet dynamics away from the polar regions. We select an aspect ratio $\eta=r_{i}/r_{o}=0.9$ which roughly falls between the values estimated for Jupiter and Saturn (Liu et al., 2008). We set a sufficiently large Taylor number, ${\rm Ta}=10^{11}$ - as used in Heimpel et al. (2005); Heimpel & Aurnou (2007); Heimpel et al. (2015) - and a Prandtl number ${\rm Pr}=10^{-2}$, which is a reasonable regime for gas giant atmospheres (French et al., 2012) due the uncertainty in determining $\Pr$. A key difference between the current polar models and those of Heimpel et al. (2005); Heimpel & Aurnou (2007) and Heimpel et al. (2015) devoted to lower latitudes, is the selection of a small Prandtl number, which strongly influences the type of convection (Simitev & Busse, 2003; Kaplan et al., 2017; Garcia et al., 2019). In our case convection directly onsets at high latitudes (Garcia et al., 2018) whereas previous deep models require a strong nonlinear regime to develop polar convection.

With an aspect ratio of $\eta=0.9$, ${\rm Ta}=10^{11}$ and ${\rm Pr}=10^{-2}$ convective onset takes the form of nonaxisymmetric polar modes confined to large latitudes (Garcia et al., 2008, 2018) at a critical Rayleigh number ${\rm Ra_{c}}=4.65\times 10^{4}$. The main difference between our two models comes from the Rayleigh number. Model 1 (M1) corresponds to ${\rm Ra}=5\times 10^{5}$ (${\rm Ra}/{\rm Ra_{c}}=10.8$), whereas Model 2 (M2) corresponds to a slightly larger ${\rm Ra}=8\times 10^{5}$ (${\rm Ra}/{\rm Ra_{c}}=17.2$). For low Prandtl and large Taylor numbers even a small degree of supercriticallity provides strong turbulent flows (Kaplan et al., 2017), characterised by a predominance of advection, rather than diffusion, in the temperature transport.

The DNS presented here have spatial resolutions very similar to those used in previous deep convection models for the same range of parameters (Eg. Heimpel et al. 2005, 2015). Specifically, we consider spherical harmonics up to order and degree $L_{\max}=500$ and a radial mesh with $n_{r}=60$ points. In contrast to previous studies, however, our numerical simulations do not rely on the use of azimuthal symmetry constraints, nor on imposed hyperviscosity, allowing us to capture low wave number coherent dynamics in the flow. These are otherwise filtered if $m=m_{d}>6$-fold symmetry constraints are imposed. However, the system of equations has in return dimension $n=(3L_{\text{max}}^{2}+6L_{\text{max}}+1)(N_{r}-1)\approx 4.4\times 10^{7}$, making the numerical integration of Models 1 and 2 very challenging.

Both models are evolved for around one viscous time unit ($1.5$ for M1 and $1.3$ for M2), including a large initial transient, with a time step of $\Delta t=2\times 10^{-7}$. The initial condition for integrating Model 1 corresponds to a transient flow at the same $\Pr=10^{-2}$ and ${\rm Ta}=10^{11}$ but with ${\rm Ra}=2.32\times 10^{5}$. The latter simulation has been initialised from a model at ${\rm Pr}=3\times 10^{-3}$, ${\rm Ta}=10^{7}$, and ${\rm Ra}=2\times 10^{3}$  (${\rm Ra}/{\rm Ra_{c}}=6.6$) obtained from the chaotic polar flows described in Garcia et al. (2019). Model M2, with larger ${\rm Ra}$, is obtained from M1. This is a common strategy in deep convection models (Christensen, 2002; Heimpel et al., 2005; Garcia et al., 2014). Because of the large dimension of the system (which has been integrated without any symmetry assumption) and the small time step required to obtain the numerical simulations, only two simulations are presented. These are however some of the highest resolution convective models, in rotating thin shells, performed to date.

The two models exhibit features of turbulent flows that occur in gas giant atmospheres. The values shown in Table 1 of the time and volume averaged Rossby $\overline{{\rm Ro}}$ (measuring the balance between Coriolis and inertial forces) or Reynolds $\overline{{\rm Re}}$ (a measure of turbulent flows) numbers are in agreement with those given in previous studies on gas giant atmospheres, for instance those based on the peak zonal flow velocities, given in Table II of Heimpel & Aurnou (2007). The modified Rayleigh number ${\rm Ra}^{*}={\rm Ra}{\rm Ta}^{-1}{\rm Pr}^{-1}\geq 0.005$, measuring the ratio of buoyancy to Coriolis forces, matches reasonably as well. In addition, the Péclet number (measuring the ratio of temperature transport by advection and diffusion) is larger than $10$ (see Kaplan et al. 2017) a value which marked the transition between laminar and turbulent flows at low ${\rm Pr}$ and large ${\rm Ta}$.

Figure 1 displays snapshots of the contour plots of the azimuthal velocity $v_{\varphi}$ at the outer sphere and the radial velocity $v_{r}$ at the sphere with radius $r=(r_{i}+r_{o})/2$ (i.e. in the middle of the shell) for both models. The spheres are viewed from $45^{\circ}$ latitude, and the $70^{\circ}$ and $80^{\circ}$ parallels surrounding the north pole are marked with a dashed line. The flow is strongly axisymmetric for both models (see first row of Fig. 1) as is typical in stress-free convection models (Aurnou & Olson, 2001; Christensen, 2002), with a noticeable wide positive (prograde) equatorial band which is more evident for M2. This is typical for flows with ${\rm Ra}^{*}\ll 1$ (Eg. Aurnou et al. 2007) in which the Coriolis forces dominate over buoyancy effects. In this case Reynolds stresses associated with columnar convection are feeding the zonal flow at mid and low latitudes (Christensen, 2002). In contrast to the latter studies, strong prograde zonal flows are produced at high latitudes in models M1 and M2.

The main difference between the models is in the latitudinal extent of these high latitude zonal circulations and their topology, which for M1 has a noticeable $m=1$ azimuthal wave number component. The existence of small scale vortices associated with turbulent flows is best seen in the contour plots of the radial velocity (second row of Fig. 1). In this case very small scale flow structures at high latitudes coexist with large scale cells elongated in the axial direction near the equator. As evidenced by the figure, the radial vortices within the $70^{\circ}$ parallel are stronger in the case of M2. The geostrophic character of the flows, i.e. the tendency of the flow to align with the rotation axis, is reflected in Figure 2 on the meridional sections of $v_{\varphi}$ and $v_{r}$. The zonal circulation of polar jets extends down to the inner boundary and is almost independent of the $z$ coordinate. This is fulfilled for the small as well as the large scale radial vortices. Notice how these large scale cells are only allowed to grow outside the tangent cylinder.

Figure 3 shows the radial vorticity patterns on the outer sphere for M1 and M2. This figure displays a transition, in agreement with shallow layer models O’Neill et al. (2015, 2016); Brueshaber et al. (2019), in which a single polar vortex (M1) is divided into smaller vortices (M2) as the thermal forcing (${\rm Ra}$) is increased. The flow structure of M1, seen in the first row of Fig. 3, reveals key features similar to those that appear in Saturn’s atmosphere (Cassini mission, Sayanagi et al., 2017). There is a single strong cyclonic vortex almost centered at both poles and with spiralling morphology. The polar vortices are surrounded by circular belts containing small vortices far away from $\pm 70^{\circ}$ latitude. These two strong cyclones are present throughout the simulation (at least one viscous time unit) in agreement with the observations (Fletcher et al., 2015).

In contrast to M1, the DNS corresponding to M2 (2nd row of Fig. 3) has long-lived polar structures resembling the polar features observed in Jupiter’s atmosphere (Adriani et al., 2018). It reveals the prevalence of several strong cyclonic coherent vortices close to the north as well as south pole, some of them with an oval shape. Azimuthally elongated structures exist, in agreement with the observations of the Juno mission (Orton et al., 2017; Adriani et al., 2018). In addition, far away from the polar regions, circular vortices are arranged in parallel lines (east-west) and the strongest lie very close to $25.8^{\circ}$ latitude, in either the northern or southern hemisphere as previously simulated in Heimpel et al. (2015) with stratified convection models.

We note that Figure 3 evidences a multimodal structure as found in the rotating convection experiments of Aurnou et al. (2018) at similar $\Pr$ number. Modes with $|m|\leq 30$ are responsible for long lived polar vortices, which extend down to the inner boundary, as well as low latitude large-scale cells (meridional sections of Fig. 2). In contrast, high wave numbers $|m|>30$ only contribute far away from the latitude circles $\pm 70^{\circ}$. For this component of the flow two different sizes of vortex are observed: those elongated in the latitudinal direction near equatorial latitudes and those with oval shape appearing around $\pm 24^{\circ}$ latitudes.

The latitude profiles at the outer surface shown in Fig. 4 mark the latitude position of the key features for both models seen in the previous figures. The time averaged Rossby number ${\rm Ro}=v_{\varphi}/\Omega r_{o}$ (ratio of inertial to Coriolis forces) and radial vorticity $w_{r}$ are plotted versus latitude for both models. The time series span a interval of around 1300 planetary rotations for Model 1 and 1500 for Model 2, with a time step of around 20 planetary rotations in both models. The time averaged zonal Rossby number ${\rm Ro}_{z}=\langle v_{\varphi}\rangle/\Omega r_{o}$ (where $\langle v_{\varphi}\rangle$ is the azimuthal average of the azimuthal velocity) is not displayed since it is almost equal to the time averaged ${\rm Ro}$, i.e. nonaxisymmetric $v_{\varphi}$ time fluctuations tend to cancel out. For both models the Rossby number ${\rm Ro}$ is clearly large close to the poles (see Fig. 4, panels (a,c)) and remains smaller at equatorial latitudes. This means that our modelling is not reasonable for the study of dynamical behaviour of gas giant planets at low latitudes (Heimpel & Aurnou, 2007; Heimpel et al., 2005, 2015), although it reproduces qualitatively several key features of the equatorial wind such as its latitudinal extent, the negative peaks close to $\pm 25^{\circ}$ latitudes and even a noticeable decrease of magnitude near the equator (see Fig. 4(c)). The polar dynamics exhibited by the present models is a new feature not previously observed in rotating thermal convection in spherical shells. There are two peaks of ${\rm Ro}$ around $76^{\circ}$ and $-78^{\circ}$ for Model 1 and they are located closer to the poles $81^{\circ}$ and $-81^{\circ}$ for Model 2. Indeed, the time dependence away from $\pm 20^{\circ}$ latitudes is significantly more vigorous for Model 1. For this model the maximum and minimum curves of ${\rm Ro}$ clearly display two latitudes $\pm 63^{\circ}$ and $\pm 81^{\circ}$ for which ${\rm Ro}$ remains roughly constant. Besides a transition between different polar vortex characteristics between Model 1 and 2, there is also a transition between equatorial behaviour: the time dependence of ${\rm Ro}$ for Model 2 is enhanced near the equator ($\pm 20^{\circ}$ latitudes) providing maximum values of ${\rm Ro}$ around $\pm 5^{\circ}$ latitude, comparable to the peaks near the poles (see Fig. 4(c)).

The radial vorticity latitude profiles shown in Fig. 4(b,d) clearly display a strong time dependence of vorticity. The main characteristic featured by both models is the strong cyclonic vorticity near both poles. For Model 1 the polar jet extends down to $\pm 70^{\circ}$ latitudes (where $w_{r}$ becomes zero) and a similar extent, down to $\pm 73^{\circ}$, is reached for Model 2. For the latter, the relative maximum/minimum near $\pm 96^{\circ}$ marks the position of the several circumpolar vortices. The radial vorticity of Model 2 exhibits stronger oscillations at lower latitudes, especially around $\pm 24^{\circ}$, where the time average amplitude is noticeable. The relative extrema at intermediate latitudes $\theta\in[25^{\circ},70^{\circ}]$ are larger for Model 1, indicating that vortex activity is more vigorous (compare also the contour plots seen in Fig. 3). The bottom row of Fig. 4 evidences the different mode contribution to vortex activity in Model 2 by considering only the spherical harmonics with $|m|\leq 30$ or with $|m|>30$. The time averaged $w_{r}$ of the latter component is nearly zero, meaning that the main contribution to mean vorticity comes from low order modes. Indeed, the time dependent $w_{r}$ is zero near the poles for the $|m|>30$ component of the flow, and thus the latter does not contribute to the origin of the polar dynamics.

Quantitatively, our models underestimate the zonal wind amplitudes measured in gas giant planets. For instance, in the case of Saturn (Sayanagi et al., 2018) the peak zonal winds close to the poles (around $\pm 85^{\circ}$ latitude, Sánchez-Lavega et al. 2006) are around $150$ m/s giving rise to ${\rm Ro}\sim 0.01$ which is roughly three times that measured for Model 1. As commented previously, values for the parameters are not so far from those obtained from the DNS of Heimpel & Aurnou (2007); Heimpel et al. (2005, 2015) (at larger $\Pr$) modelling wind jets at moderate and low latitudes. However, our models exhibit larger ${\rm Ro}$ and circular vortices at high latitudes, in contrast to previous modelling. In addition, preliminary results for models with larger ${\rm Ra}$ point to a decrease of zonal flows and a disappearance of coherent structures in the polar regions. Thus there exists a flow transition, in which zonal flow and vortices at high latitudes progressively weaken, in a parameter regime relevant for giant planetary atmospheres. This transition also occurs at lower rotation rates when polar modes are preferred at the onset (Garcia et al., 2019). Understanding this transition may provide a simple explanation for the existence of polar convective coherent vortices in planetary atmospheres.

An analysis of the kinetic energy distribution among the spherical harmonic modes with different degree $l$ is performed, to elucidate the inverse cascade mechanism for these three-dimensional rapidly rotating flows (Rubio et al., 2014). In addition, kinetic energy spectra are a tool which helps to validate the spherical harmonic resolution of a given DNS. A rule of thumb (Christensen et al., 1999; Kutzner & Christensen, 2002) indicates that a model is well resolved if the kinetic energy decreases by more than two orders of magnitude from the maximum to the smallest wave-length scale.

Figure 5(a,b) displays the volume averaged kinetic energy $K$ (dashed line) and the surface averaged kinetic energy $K_{o}$ at $r=r_{o}$ (solid line) of a particular mode with degree $l$ for both models M1 and M2. Notice that the energies of few harmonics just before the maximum spectral degree $L_{\text{max}}=500$ increase a little bit. This behaviour is present in some three-dimensional DNS in rotating spherical geometry, but it has been shown (see Christensen et al. 1999 for the geodynamo case) that a reasonable energy increase at the tail does not significantly alter the large-scale structure of the flow.

The spectra exhibit slopes of $-5/3$ and $-3$ which are in reasonable agreement with those obtained from 2D maps of velocity measurements extracted from Cassini images (Choi & Showman, 2011). Furthermore, the distribution of energy among the different degree $l$ displayed in Fig. 5(a,c) resembles the distribution exhibited on Fig.5(a) of Choi & Showman (2011), constructed from Jovian velocity measurements. Because the flow has a strong equatorially symmetric component, odd $l$ have larger energy than even $l$, especially for low values, giving rise to the spiking behaviour of the spectra. Note that Choi & Showman (2011) computes the power spectra of kinetic energy, rather than the energy content of each mode, so that even $l$ have larger energy. The length scale marking the change of slope from $-5/3$ to $-3$ of the Jovian spectra of Choi & Showman (2011), assumed to be the scale at which energy is injected to the system, is in the interval $70-140$ depending whether global or eddy kinetic energy are considered. In our case the situation is similar and the $-5/3$ slope is lost from $l\approx 60$ in the case of M1 or from $l\approx 120$ in the case of M2.

In the context of three-dimensional Rayleigh-Bénard rotating convection an upscale energy transfer mechanism, known as the condensation process, has recently been identified in the regime of GT (Rubio et al., 2014). Large scale barotropic vortices, fulfilling the Taylor-Proudman constraint, are generated through an inverse cascade mechanism. Small scale convective motions act as an energy source for these barotropic vortices, which in turn organise these small convective structures. This is reflected in the energy spectra of the flow when considering the barotropic and baroclinic components separately. Flow motions in regions close to the poles in the case of three dimensional rotating spherical shells can be approximated to flow motions in a three-dimensional rotating plane layer if the shell is sufficiently thin. Because this is the case of M1 and M2, we may compare our kinetic spectra with those obtained on a plane geometry. Following Rubio et al. (2014) we compute the power spectral energy of the barotropic and baroclinic components of the flow. Close to the poles we may assume ${\bf\hat{e}_{r}}\approx{\bf\hat{e}_{z}}$ (the radial direction is almost parallel to the axis of rotation) and the plane spanned by ${\bf\hat{e}_{\theta}}$ and ${\bf\hat{e}_{\varphi}}$ is parallel to the $x-y$ plane. We define the barotropic kinetic energy $K_{\text{bt}}=\frac{1}{2}(\langle v_{\theta}\rangle^{2}+\langle v_{\varphi}\rangle^{2})$ and the baroclinic kinetic energy $K_{\text{bc}}=\frac{1}{2}(v_{r}^{{}^{\prime}2}+v_{\theta}^{{}^{\prime}2}+v_{\varphi}^{{}^{\prime}2})$, $\langle f\rangle$ being the radial average (barotropic component) and $f^{\prime}=f-\langle f\rangle$ the baroclinic component.

Figure 5(c) displays the power spectra for the barotropic and baroclinic kinetic energies, as in Fig. (3) of Rubio et al. (2014). The qualitative agreement between both figures is remarkable. For the smaller scales the barotropic kinetic energy $K_{\text{bt}}$ displays an $l^{-3}$ scaling as expected for a downscale enstrophy cascade (Rubio et al., 2014). In contrast, at larger scales the expected $l^{-5/3}$ scaling is replaced by $l^{-3}$ due to the appearance of a large dipole structure. The scaling $l^{-5/3}$ observed for the baroclinic kinetic energy $K_{\text{bc}}$ is in agreement with a three-dimensional turbulence cascade from the convective scale. The different scalings observed for the barotropic and baroclinic components provide an explanation for the different scalings observed for the volume averaged total kinetic energy $K$ of Fig. 5(a,b). At larger scales $K_{\text{bt}}$ dominates and thus $K$ exhibits the same scaling. In contrast, for $l>30$ the baroclinic flow is dominant and thus there is an interval for which $K$ follows the baroclinic $-5/3$ scaling.

In this section we provide numerical evidence of the formation of a large scale dipole close to the polar regions, as in Rubio et al. (2014). In addition, we interpret the transition between a single polar vortex (as model M1) and multiple polar vortices (as model M2) following ideas developed by O’Neill et al. (2015, 2016) in the context of shallow layer modelling.

Coherent large scale structures are ubiquitous in the atmospheres of gas giant planets (Fletcher et al., 2018; Adriani et al., 2018), and have been simulated within the 2D and 3D GT context (see Rubio et al. 2014 and references therein). The latter study demonstrated the existence of a dipole condensate in GT flows and provided an explanation of the physical mechanism behind for its formation using a 3D reduced model in the limit of small ${\rm Ro}$. As noted in the previous section, the advection of 3D small scale convective baroclinic vortex produces Taylor-Proudman (almost $z$-independent, 2D) barotropic large scale motions, which in turn organise the small scale convective motions in a self sustained process.

The relevant regime for the formation of the dipole condensate found in Rubio et al. (2014) is ${\rm Ra}_{\text{layer}}=O({\rm E}^{-4/3})$ (see also Julien et al. 2012b), corresponding to a regime that is overall dominated by rotation, i. e. GT . For $\Pr=1$ a coherent large scale dipole was found in coexistence with small scale convective vortices. In this regime the volume renderings of axial vorticity of Rubio et al. (2014) help to visualise the organization of large columnar structures extending from the bottom to the top boundary together with very small convective vortices. Taking into account the differences in the definition of ${\rm Ra}={\rm Ra}_{\text{layer}}(1-\eta)$ the condition ${\rm Ra}_{\text{layer}}{\rm E}^{4/3}=O(1)$ becomes ${\rm Ra}{\rm E}^{4/3}=O(0.1)$. In our case ${\rm Ra}{\rm E}^{4/3}=0.023$ for M1 and ${\rm Ra}{\rm E}^{4/3}=0.037$ for M2, values not so far from $O(0.1)$. As noticed in Julien et al. (2012b) for small $\Pr$ the GT regime is attained for smaller values of ${\rm Ra}_{\text{layer}}{\rm E}^{4/3}$ which is consistent with our results, as $\Pr=0.01$.

The three-dimensional structure of the dipole condensate for both models can be visualised in Fig.6 which displays the surfaces of constant radial vorticity $\omega_{r}$ for two different values (roughly half and one tenth of the maximum value). Two different views are provided, containing the $x-y$ and the $y-z$ planes, respectively, which help us to understand the horizontal and vertical structure of $\omega_{r}$. The surfaces of constant $\omega_{r}$ have been restricted to the volume within the coaxial cylinder of radius $3d$ and on the northern hemisphere, for comparison with the results of the rotating plane geometry of Rubio et al. (2014). This comparison makes sense as the rotating spherical shell is very thin and thus well approximated by a plane layer (see Fig.6, views containing the $z$ axis).

The structure of the main polar cyclonic vortex in the case of M1 is strongly anisotrophic, with slightly smaller horizontal scales (compared to the vertical scale) spanning the whole shell. Small scale cyclonic convective vortices of comparable amplitude surround the main cyclone, some of them being very small (see the red surfaces of Fig.6, left column, Model 1). For small vorticity amplitudes the anticyclonic vortices, having relatively small horizontal scales, are less common than cyclonic vortices, evidencing a clear asymmetry between cyclonic/anticyclonic motions. This asymmetry was also found in the experiments of Vorobieff & Ecke (2002) investigating plane Rayleigh-Bénard rotating convection for ${\rm Ro}\lesssim 1$.

The surfaces of constant $\omega_{r}$ for M2 display similar structures to those of M1, the main difference being the disruption of the main cyclone of M1 into the several vortices with smaller horizontal scale of M2, as shown previously in the contour plots of Fig. 3 in Sec. 3.1. The vertical scale of these main cyclonic vortices is of the order of the width of the shell (see the red surfaces of Fig.6, left column, Model 1). Again a clear asymmetry between cyclonic/anticyclonic vortices is present in M2. The asymmetry between cyclonic and anticyclonic vorticity is mainly due to the axisymmetric component, which is not present for plane layer modelling (Rubio et al., 2014) as the latitudinal variation of the Coriolis force is not considered in their models. The contribution to the axisymmetric component to cyclonic vorticity is evidenced in Fig. 7 when considering only $\omega_{r}$ for the nonaxisymmetric component of the flow. In this case the number of cyclonic/anticyclonic vortices is similar, but the cyclonic motions have larger vertical and horizontal scales.

The analysis by O’Neill et al. (2015, 2016) (and also Brueshaber et al. 2019) of shallow atmospheric models has provided evidence for a transition between different dynamical regimes which reproduces important features seen in giant planet atmospheres. The transition between Jupiter-like and Saturn-like polar dynamic regimes is described in terms of the Burger number. To compare our fully convective models with previous shallow models, the Burger number for M1 and M2 must be estimated. For the present global simulations the spherical domain as well as its rotation is fixed, but the characteristic horizontal length scale $L$ of the polar convective region varies between models, giving rise to different Burger numbers. In our models the variation of $L$ is achieved by modifying the Rayleigh number, i.e. the thermal forcing, and keeping all of the other parameters (Eq. 4) fixed. In this sense the Burger number represents an output parameter from the simulations which is computed after the flow is saturated, as suggested in Read (2011). This approach differs from the simulations of O’Neill et al. (2015, 2016); Brueshaber et al. (2019) in which the Burger number represents an input parameter of the models.

According to Pedlosky (1979) the Burger number is defined as ${\rm Bu}=\frac{g\Delta\rho D}{\rho 4\Omega^{2}L^{2}}$, $L$ and $D$ being the horizontal and vertical length scales and $\Delta\rho$ the density difference between the surface and the base of the fluid. From the Boussinesq approximation the latter can be expressed as $\Delta\rho=-\rho_{\text{outer}}\alpha\Delta T$, and assuming $\rho\approx\rho_{\text{outer}}$ the Burger number is

Considering the full spherical shell the only quantity in this latter expression that changes between M1 and M2 is $\Delta T$. However, as we are interested in studying polar dynamics we should restrict the definition of ${\rm Bu}$ to the polar region. This will also permit a comparison with the local shallow atmospheric models of O’Neill et al. (2015, 2016); Brueshaber et al. (2019). For these reasons we take as characteristic vertical and horizontal length scales $D=d$ and $L=\theta_{\text{polar}}r_{o}$, respectively, $\theta_{\text{polar}}$ being the latitudinal angle of the extent of polar motions on the outer surface, which clearly changes between the models, as shown in previous sections. In terms of the nondimensional input parameters the Burger number is then approximated as

Figure 6 provides a way to infer the latitudinal extent of polar motions in order to estimate ${\rm Bu}$ for each model. We focus on the surfaces of constant vorticity (at a value roughly half the maximum) viewed from the north pole that appear in Fig. 6. For M1, the main cyclone and its surrounding small vortices are contained in a rectangle (on the $x-y$ plane, see Fig. 6) with a maximum dimension of $3d$ (in the $x$ direction). Projecting this rectangle onto the outer sphere gives $\theta_{\text{polar}}\approx 17^{\circ}$. In contrast, for M2, vortices having half of the maximum vorticity amplitude spread within and over the whole coaxial cylinder of diameter $6d$. This gives rise to a latitudinal extent $\theta_{\text{polar}}\approx 35^{\circ}$, two times larger than for M1. The corresponding Burger numbers (from Eq. 5) are ${\rm Bu}_{1}=1.4\times 10^{-4}$ and ${\rm Bu}_{2}=5.4\times 10^{-5}$ for M1 and M2, respectively. These values qualitatively agree with and are not so far from those obtained in Brueshaber et al. (2019), which found ${\rm Bu}>10^{-3}$ and ${\rm Bu}<5\times 10^{-4}$ for Saturn and Jupiter like polar dynamics, respectively. The order of magnitude of difference in the computation of ${\rm Bu}$ can be attributed to the uncertainty in the determination of the Rossby deformation radius $L_{\text{d}0}$ of Jupiter and Saturn from which ${\rm Bu}=(L_{\text{d}0}/a)^{2}$ ($a$ being the osculating radius at the poles,) was estimated in Brueshaber et al. (2019). According to Brueshaber et al. (2019) the smallest estimates of $L_{\text{d}0}$ for Jupiter may give rise to ${\rm Bu}_{\text{J}}\lesssim 10^{-4}$, and in the case of Saturn a value of $L_{\text{d}0}=1000$ km estimated in Choi et al. (2009) gives rise to ${\rm Bu}_{\text{J}}\approx 2\times 10^{-4}$ (assuming $a=66810$ km). These values are certainly in reasonable agreement with our estimations.

The physical mechanisms giving rise to a strong vortex centered at the pole for shallow models (O’Neill et al., 2015, 2016) and our thermal convective simulations share two key characteristics. Small scale baroclinic instabilities cascade, creating a vertically aligned strong barotropic vortex. This is demonstrated in the first part of this section as well as in O’Neill et al. (2015, 2016). In addition, the same parameter (see also Brueshaber et al. 2019) is controlling the transition between Jupiter-like and Saturn like polar patterns for both types of modelling. The authors of O’Neill et al. (2016) raise the possibility that the weather layer of Saturn may be coupled with the convective interior below. As discussed in Sec. 8 of O’Neill et al. (2016) convective structures may allow local stability, which could affect the coupling between the molecular zone below the weather layer. In this respect, the present study is providing more support for the idea that polar vortices are generated deep in the atmosphere, as recently discovered for the well-known east-west low and mid latitude zonal jets (see Kaspi et al. 2018).

The analysis of the force balances involved in the generation of geophysical flows helps us to identify and characterise different flow regimes and to shed light on the transitions occurring among them (Oruba & Dormy, 2014; Gastine et al., 2016; Garcia et al., 2017; Julien et al., 2012a). In the context of spherical shell rotating convection the generation of zonal flows strongly depends on the large scale force balance, and more specifically on the balance between Coriolis and buoyancy forces. The nondimensional number ${\rm Ra}^{*}$ was used in Aurnou et al. (2007) to distinguish between rotation dominated regimes ${\rm Ra}^{*}\ll 1$, with strong generation of prograde zonal flows at mid-low latitudes fed by Reynolds stresses associated with columnar convection in spherical geometry (Christensen, 2002), regimes in which ${\rm Ra}^{*}\sim 1$ buoyancy starts to play a role and three dimensional convection acts to homogenise the fluid within the shell, generating retrograde zonal flows in the equatorial region and strong prograde flows at high latitudes. As commented in Glatzmaier (2018) there are studies (for instance Kaspi et al. (2018)) that have attributed the zonal flow generation to a thermal wind balance, in which the curl of Coriolis force balances the curl of buoyancy. In convective models of rotating spherical shells (Glatzmaier, 2018), and in our models M1 and M2, this balance does not hold as the curl of Coriolis force dominates over the other terms.

Because ${\rm Ra}^{*}\ll 1$ for the numerical models M1 and M2 (see Sec. 2.2), Coriolis forces are governing the global dynamics, and prograde zonal flows are generated at low latitudes. The time and volume averages of the forces $\mathcal{F}=\frac{1}{V}\int_{V}(\mathcal{F}_{r}^{2}+\mathcal{F}_{\theta}^{2}+\mathcal{F}_{\varphi}^{2})^{1/2}\text{d}V$ are shown in Table 2 and verify the predominance of the Coriolis forces. The mean values of the Coriolis $\mathcal{F}_{\text{C}}$, inertial $\mathcal{F}_{\text{I}}$, and viscous $\mathcal{F}_{\text{V}}$ forces verify $\mathcal{F}_{\text{Coriolis}}>\mathcal{F}_{\text{inertial}}\gg\mathcal{F}_{\text{viscous}}$ which is the force balance believed to operate in Jupiter’s convective atmosphere (Schneider & Liu, 2009).

However, and in contrast to previous numerical studies (Eg. Aurnou et al. 2007; Gastine et al. 2013), strong prograde zonal flows are generated at high latitudes as well, which are preferred for ${\rm Ra}^{*}\sim 1$ due to an angular momentum mixing produced by three-dimensional convection (Aurnou et al., 2007). As the angular momentum depends only on the cylindrical coordinate the zonal flow generated at polar latitudes is still geostrophic, as in models M1 and M2. For these models the ${\rm Ra}^{*}\ll 1$ condition seems to be relaxed when considering the force balance only in the azimuthal direction, allowing the development of high latitude prograde zonal flows. Table. 3 summarises the balance for each component $(\mathcal{F}_{r},\mathcal{F}_{\theta},\mathcal{F}_{\varphi})$ of the forces. The maximum absolute value within the shell is listed. We note that the values are picked up at the same time instant as the snapshots of the contour plots and surfaces of constant radial vorticity shown previously. Whereas the Coriolis force is clearly dominant for the radial and colatitudinal components, in the azimuthal direction the inertial forces are comparable.

The instantaneous spatial distribution of the azimuthal component of the Coriolis and inertial forces at the outer surface is illustrated in Fig. 8 for both models at the same time instant as the previous figures of contour plots. The figure shows that clear correlation of both forces in the polar regions, indicating the importance of three-dimensional convection. Figure 8 also suggests that the Coriolis predominance in the azimuthal direction is relaxed close to the equatorial regions as well, particularly for M1. This may be the reason for the relatively small amplitude of the prograde zonal flow of M1 and M2 (see Fig. 4) at equatorial latitudes with respect to previous models of Jupiter and Saturn (see for instance Heimpel & Aurnou 2007).

The strong dipolar character of the force balance in the case of M1 is lost for M2. The ratio $(\mathcal{F}_{\text{I}})_{\varphi}/(\mathcal{F}_{\text{C}})_{\varphi}$ between the maximum values over the shell (see Table. 3) is $0.23$ for M1 and $2.67$ for M2. This suggests that the value $(\mathcal{F}_{\text{I}})_{\varphi}/(\mathcal{F}_{\text{C}})_{\varphi}\approx 1$ seems to control the transition from flows with a single polar cyclone like M1, to flows with multiple polar cyclones like M2. At the outer surface the situation is similar, the ratio being $0.13$ for M1 and $0.6$ for M2 (see colour palette of Fig. 8). Although the ratio for M1 is still smaller than unity, we note that in the framework of Rayleigh-Bénard rotating convection with parallel horizontal boundaries (Julien et al., 2012a), which approximates our model close to the polar regions, the GT regime is attained for values of ${\rm Ra}^{*}$ which can be smaller than $O(1)$.

Previous studies in rapidly rotating thin spherical shells (Aurnou et al., 2008) have analysed the heat transfer mechanisms in the context of planetary atmospheres. In the equatorial regions the heat flux is inhibited in the cylindrical direction as a consequence of the strong generation of zonal flows. In contrast, in the polar regions three-dimensional thermal plumes develop which are not sensitive to the geostrophic zonal flow (Sreenivasan & Jones, 2006), and thus at large latitudes the heat transfer resembles that predicted for plane layer convection. Figure 9 displays the instantaneous contour plots of the temperature perturbation $\Theta$ on a meridional section for both models M1 and M2. Thermal perturbations in the equatorial regions are affected by the flow spiralling in the azimuthal direction and extend parallel to the rotation axis in a wide region outside the tangent cylinder. The structure of $\Theta$ close to the poles is different, $\Theta$ being positive close to the outer boundary and negative close to the inner, in agreement previous studies in rotating spherical geometry of gas giant atmospheres, see for instance Fig. 4 of Aurnou et al. (2008). However, in contrast to the latter study the Nusselt number of our models is significantly smaller, as in the polar regions the regime of GT is already attained due to the low Prandtl number of our models (Julien et al., 2012b).

For Rayleigh-Bénard rotating convection within parallel horizontal boundaries a recent study Julien et al. (2012b) provided an analysis of the heat transfer in terms of a reduced model in the regime of GT (Julien et al., 2012a; Rubio et al., 2014). As described prevously, our models M1 and M2 exhibit features, such as the dipole condensate (Sec. 3.4) and the kinetic energy spectra (Sec. 3.3) characteristic of GT flows, thus the thermal properties of M1 and M2 are in concordance with Julien et al. (2012b). Figure 10 displays the surfaces of constant nonaxisymmetric temperature perturbation. The latter quantity can be compared to the temperature fluctuations with respect to the time averaged temperature (strongly axisymmetric) of the modelling of Julien et al. (2012b), yielding qualitative agreement of the tridimensional structure of the thermal deviations (see Fig. 1 of Julien et al. 2012b).

A scaling law for the heat transport (measured by the Nusselt number ${\rm Nu}$) was derived in Julien et al. (2012b) for GT flows. In this regime turbulent 3D convection damps the heat transport in the bulk of the fluid rather than in the boundary layers, as a consequence of the weak $z$-dependence of the flow. The results of Julien et al. (2012b) point to the independence of heat transport with respect to microscopic diffusion coefficients, leading to the scaling

with ${\rm E}_{\text{layer}}=1/2{\rm Ta}^{1/2}$ for sufficiently small ${\rm Pr}\leq 1$ and large ${\rm Ra}_{\text{layer}}$. In terms of the definition of ${\rm Ra}={\rm Ra}_{\text{layer}}(1-\eta)$ used in the present study the Nusselt scaling becomes

The time averaged Nusselt number ${\rm Nu}$, defined as the ratio of the average of the total radial heat flux to the conductive heat flux (both through the outer surface), has been computed for M1 and M2 in the same time interval as the mean physical properties shown in Sec. 2.2, Table. 1. For M1 ${\rm Nu}_{1}-1=0.00827$ whereas for M2 ${\rm Nu}_{2}-1=0.0174$, giving rise to almost the same constant $C_{1}^{1}=0.0296$ and $C_{1}^{2}=0.0308$, meaning that our models verify the GT heat transfer scaling derived in Julien et al. (2012b). Indeed, the value of the constant is quite similar to the value of $C_{1}=0.04$ obtained in Julien et al. (2012b) by means of a detailed exploration of the parameter space, and thus our results fit reasonably well with the theory of heat transfer developed for plane rotating Rayleigh-Bénard convection.

A significant achievement of the full-globe general circulation models of Liu & Schneider (2010) was to reproduce a meandering jet with hexagonal-like structure, as observed in Saturn (Godfrey, 1988). Our convective Model 1 generates a similar structure: belts of weak vortices surrounding the north pole (Fig. 3, first row, left section) appear to have polygonal shape. The relatively strong $m=6$ component is reflected in the plot of volume averaged kinetic energy contained in each wave number $m$ shown in Fig. 11. For M1 the $m=4$ and $m=6$ components are dominant for $m>2$ and they give rise to polygonal structures surrounding the poles. A peak at around $m=90$ is noticeable in the kinetic energy spectra. This corresponds to a higher harmonic multiple of $m=6$, which is responsible for the edges and straight lines forming the hexagonal pattern. A north pole hexagonal boundary can clearly be identified if we consider only the component of the flow containing the $m=6k$, $k\in\mathbb{Z}$, azimuthal wave numbers in its spherical harmonic expansion. This flow component is strong, providing $82.7\%$ of the kinetic energy of the flow. The contour plots of the kinetic energy density at the outer surface ($r=r_{o}$) are displayed on Fig. 12, the 1st and 2nd left sections corresponding to the north and south poles, respectively. As described for Saturn (E.g. the review (Sayanagi et al., 2018)), the north hexagonal boundary is very close to $75^{\circ}$ latitude. Our results account for a noticeable asymmetry between north and south dynamics, which is a characteristic feature of Saturn (see discussion in Sec. 12.4 of Sayanagi et al. 2018). The contour plots of the temperature perturbation shown in the two rightmost sections of Fig. 12 enable us to make a qualitative comparison bewteen our results and real measurements such as the brightness temperature maps (from Cassini/CIRS spectroscopy) obtained in Fletcher et al. (2015), with reasonable agreement.

More than 6 years of Cassini observations have revealed that the hexagonal pattern has rotated around 30^∘ in the azimuthal direction (Fletcher et al., 2015). Thus the pattern remains nearly stationary in the planetary reference system (Sánchez-Lavega et al., 2014). Our preliminary investigations indicate that azimuthal drift (in the planetary rotating frame) of the hexagonal pattern seems to be slow. This is displayed in Figure 13 showing the variation of the isosurfaces of kinetic energy density of the flow after 2 and 4 planetary rotations, from a polar viewpoint. As observed for Saturn (E.g. Figure 1(a) of Sánchez-Lavega et al. (2014)) the hexagonal pattern is bounding small horizontal scale vortices at around 75^∘ latitude (marked with segments on the axis) and the bottom vertex (that on the $x$ axis) remains stationary. The vortices extend down to the deep interior in a spiralling fashion, so its vertical scale is $d$.

The analysis of Fletcher et al. (2015) showed that the main polar vortex and its surrounding hexagon have been persistent features on the troposphere for over a decade, despite the seasonal evolution of the temperature and composition. The stability of the hexagonal pattern despite seasonal variations of insolation induced the authors of Sánchez-Lavega et al. (2014) (see also Fletcher et al. 2018) to propose its origin as consequence of a Rossby wave extending deep into the planetary atmosphere. The simulated deep structures for M1 persist over more than $10^{4}$ planetary rotations, so the hexagonal pattern can survive in a deep atmospheric model. This reinforces the idea of a deep origin of the hexagonal pattern, as suggested by Sánchez-Lavega et al. (2014); Fletcher et al. (2018).