Penetrative Convection in Partly Stratified Rapidly Rotating Spherical Shells

Celestial objects host interfaces between convective and stable stratified interior regions. The interaction between both, e.g., the transfer of heat, mass, or angular momentum depends on whether and how flows penetrate into the stable layer. Powered from the unstable, convective regions, radial flows can pierce into the stable region depending on their inertia (overshooting). In rapidly rotating systems, the dynamics are strongly influenced by the Coriolis force and radial flows penetrate in stratified regions due to the geostrophic invariance of columnar convection even in the limit of vanishing inertia. Within this study, we numerically investigate both mechanisms and hence explore the nature of penetrative convection in rapidly rotating spherical shells. The study covers a broad range of system parameters, such as the strength of the stratification relative to the Coriolis force or the inertia. Guided by the application to Saturn, we model a sandwiched stable stratified layer (SSL) surrounded by two convective zones. A comprehensive analysis of the damping behavior of convective flows at the edges of the SSL showed that the mean penetration depth is controlled by the ratio of stratified and unstratified buoyancy gradients and is hence independent of rotation. A scaling law is derived and suggests that the penetration depth decreases with the square root of the ratio of unstabilizing and stabilizing entropy gradients. The influence of the Coriolis force, however, is evident by a modulation of the penetration depth along latitude, since convective columns are elongated vertically and hence pierce predominantly into the SSL around mid-latitudes and outside the tangent cylinder. Our result also show that the penetration depth decreases linearly with the flow length scale (low pass filter), confirming predictions from the linear theory of rotating partially stratified convection.


Introduction
When the local temperature gradient is steeper than the one associated with an adiabat, small perturbations from the hydrostatic equilibrium amplify to the well-known Rayleigh-Taylor convective instability. This leads to vigorous convection that very rapidly reestablishes bulk adiabatic gradients of density and temperature due to the inherent mixing and heat transport efficiency.
However, stable stratified regions, in which the heat flux is conductive or radiative, are widespread phenomena in stars and planets. Those regions are caused by either subadiabatic temperature or positive heavy element gradient. In stars, the large efficiency of radiative heat transport induces stratified regions, e.g. in the radiative part of the solar interior (e.g. Zahn, 1991;Brun et al., 2011). Stable stratification due to temperature inversion or subadiabatic temperature lapse rates are known for the Gas Giants (Jupiter Galileo probe, (Atkinson et al., 1996)) and outermost atmosphere of terrestrial planets, where radiation is most efficient. In the Earth's liquid outer core, the outermost layer seems stratified caused by either a subadiabatic temperature gradient or a local enrichment of light elements (Fearn and Loper, 1981;Lister and Buffett, 1998). Such a layer would have profound consequences for core convection and dynamo action in the core and be traceable in the geomagnetic field (Nakagawa, 2011;Buffett, 2014). A potential effect of electrically conducting, stable stratified layers surrounding a dynamo region are weak amplitude, strongly axisymmetric surface magnetic fields and are often proposed for Mercury's or Saturn's dynamo regions (Stevenson, 1982;Schubert et al., 1988;Christensen, 2006).
For Saturn a particular interesting scenario, consisting of a stratified layer between two convective ones is favored by the H/He demixing behavior (Stevenson and Salpeter, 1977;Püstow et al., 2016;Schöttler and Redmer, 2018) and in agreement with several observational constraints, such as the peculiarly axisymmetric magnetic field (e.g. Cao et al., 2012), the apparently impeded thermal evolution (Leconte and Chabrier, 2013;Nettelmann et al., 2013) and the detection of gravity wave induced pulsations in Saturn's rings (Hedman and Nicholson, 2013;Fuller, 2014). Saturn's interior might have been homogeneously mixed in the early, hot stages of the thermal evolution. Once the adiabatic temperature gradient crosses the demixing curve, He droplets form, sink downwards and are remixed at greater depth. This process builds up a compositional gradient roughly around mid-depth bridging from He-depleted outer to He-enriched inner convective zone (Stevenson and Salpeter, 1977;Schöttler and Redmer, 2018).
Radial, convective flows originating from the unstable regions might penetrate into the stratified region due to inertia or the geostrophy of the columnar convective structures.
It is useful to distinguish qualitatively the different end-members scenarios in terms of the dominant forces. In rotating, partially stratified convection, the three main forces, each associated with a timescale, are the Coriolis force (τ rot ), the buoyancy of the rising fluid parcel (τ buo ) and the inverse buoyancy due to the stratification (τ strat ). The time scales are given by the rotation rate Ω, τ rot = (2Ω) −1 , the Brunt-Väisälä frequency N , τ strat = √ N −2 = (−g/ρ dρ/dr) −1/2 and the time scale of buoyancy to accelerate fluid parcels τ buo = (g/ρ ρ /D) −1/2 , respectively. Here ρ is the density fluctuation of a buoyant fluid parcel, ρ the ambient density, D a typical length scale and g the local gravity. For our setup, τ buo characterizes the convective and τ strat the stratified regions.
The ratio of the time scales can be expressed as the convective Rossby number Ro c and the stratification parameter I s : In slowly rotating systems, rotational forces are less important (I s < 1, Ro c > 1) and a hot fluid parcel is accelerated by the buoyancy in the convective region, rises radially upwards and eventually pierces into the SSL, where the depth of penetration depends on the previously gained inertia. This represents the classical case of non-rotating penetrative convection, typically called overshooting or inertia penetration. For this end-member, the penetration depth should be independent of colatitude, and may be the dominant form of penetration close to the rotation axis for Ro c < 1 and in the equatorial regions in our models.
In rapidly rotating systems, convective flows are strongly modified by the rotation (Ro c < 1). The dominant Coriolis force constrains the flow to be invariant along the rotation axis leading to the well-known convective columns and strongly geostrophic zonal flows. This implies that the flows extend into the stratified layer independent of inertia, but with a characteristic colatitudinal modulation. When the Coriolis force is stronger than the gravity force associated with the stratification (i.e. I s > 1, Ro c < 1), rotational penetration is strong and acts predominantly at mid-latitudes and outside the tangent cylinder (TC), where convective columns can vertically extend through the whole spherical shell. Since in the equatorial regions and close to the rotation axis no convective columns can be extended into the SSL, rotational penetration yields a characteristic colatitude dependence. If, on the other hand, I s < 1 and Ro c < 1 the SSL should be devoid of radial flows since the stratification is strong enough to efficiently break the vertical stiffness of geostrophic flows.
A fundamental understanding of penetrative convection in rapidly rotating spherical shells is of ample importance in geo-and astrophysics. The efficiency or vigor of penetrative convection controls the transport of heat, mass, angular momentum and magnetic fields across the interface between unstable and stable stratified regions. If the penetration is strong and the associated heat transport efficient, adiabatic regions are extended into the stratified layers (e.g. Browning et al., 2004). This means that stratified layers can be eroded by the permanent entrainment of convective flows and subsequent efficient mixing (e.g. Ellison and Turner, 1959;Levy and Fernando, 2002). Another complication arises from the combination of spherical geometry and rapid rotation. If rotational forces dominate over the inverse buoyancy associated with the stratification, piercing radial flows are non-uniform along colatitude leading to latitudinal entropy (density) gradients and hence drive baroclinic instabilities. Those in turn will act as a boundary condition for the differential rotation in the convective regions and potentially alter the magnetic field (Stevenson, 1982).
Classic studies of penetrative convection dealt with non-rotating cartesian setups of Rayleigh-Bénard convection (Veronis, 1963). The effect of rotation was firstly considered in the framework of oceanic dynamics (e.g. Julien et al., 1996). The investigation of rotating convection in spherical shells below or above a stable stratified layer was mainly driven by studies on the interaction and efficiency of radiative and/or convective heat transport in solar or stellar interiors (e.g. Zahn, 1991;Brummell et al., 2002;Rogers et al., 2006). For rapidly rotating A-type stars, where a convective core is surrounded by a radiative outer envelope, Browning et al. (2004) reported that rotational and overshooting penetrative convection generates adiabatic regions in the radiative zone preferentially at higher latitudes (prolate adiabatic core).
Several numerical investigations tuned to the solar setup with a convective envelope enclosing a deeper radiation dominated zone, investigate the parameter dependence of penetrative convection (Zahn, 1991;Hurlburt et al., 1994;Rogers et al., 2006). Since the stratification originates from enhanced radiative heat transport, it is typically assumed that penetration depth is largely characterized by two parameters, one being the Péclet number given by P e = U D/κ, where U is a typical flow speed, D a length scale and κ the thermal diffusivity measuring the radiative heat transport (Zahn, 1991;Brummell et al., unstable region (Hurlburt et al., 1994). However, different scaling relations between the penetration depth and S have been reported (Hurlburt et al., 1994;Brummell et al., 2002;Rogers et al., 2006). Interestingly, in models that take spherical geometry and rotation into account, the associated Coriolis forces appear unimportant for the magnitude of the penetration depth.
In the planetary context, the Galileo Probe revealed subadiabatic temperature gradients and non-ceasing zonal flows in Jupiter's outermost atmosphere down to roughly 20 bar of atmospheric pressure. This suggests that zonal flows are not damped, but rather maintained in stratified regions yet disconnected from their potential sources, such as convection or irradiation gradients. Assuming that the axisymmetric winds are driven from correlations of the deep convection (i.e. Reynolds stresses), Schubert (1996, 1997) showed that zonal mean flows are vertically extended by the Coriolis forces in accordance with the Taylor-Proudman theorem, but non-axisymmetric convective flows are damped by the inverse buoyancy gradients due to the stratification. From a more theoretical point of view, the linearized model of inviscid, inertia-less, rapidly rotating penetrative convection as studied by Takehiro and Lister (2001) shows a linear length scale dependence of penetration depth such that larger scale flows (like differential rotation) are less damped than short length scale flows typically associated with convective motions. The overall damping amplitude depends on the gravity force associated with the stable stratification relative to the Coriolis force (I s ). Subsequent studies including nonlinear inertia showed a more complex damping behavior (Takehiro and Lister, 2002), but the results offered an attractive explanation how zonal flows are extended through stratified regions.
The major aim of the present study is to understand how overshoot and rotational penetration act across the interface between stratified and unstratified regions in rapidly rotating, spherical shell models, which are most suitable for interiors of planets. Even though the effects of stratification on differential rotation and magnetic fields as far as the emerging waves inside the SSL are left aside in this study, the investigation of fundamental properties of rapidly rotating, penetrative convection serves as a basis of upcoming research. Our models cover a comprehensive part of the parameter space, such that the different regimes appropriate for planets are reached in terms of I s and Ro c .
The general setup featuring a thin, sandwiched stratified layer centered at mid-depth and surrounded by two thick convective zones is motivated from the application to Saturn.

Model
The non-dimensional governing equations for conservation of mass, momentum and thermal energy for an ideal gas in the anelastic approximation are given by (Jones, 2011;Gastine and Wicht, 2012;Verhoeven et al., 2015;Wicht et al., 2018): where p is the reduced pressure, F ν is the viscous force, and Q ν represents viscous heating given by where σ ij is the stress tensor. The non-dimensional parameters are given by These equations were formulated for an adiabatic background state. In our model, a prescribed non-adiabaticity of the background state powers or suppresses convection.
The amplitude of these destabilizing and stabilizing entropy gradients must remain small relative to the adiabat to make sure the equations still hold. The non-adiabaticity of the new background state is thus scaled with s 1, where ds/dr is the dimensionless analytically prescribed stratification profile and Γ the Grüneisenparameter which is in an ideal gas related to the polytropic index n and specific heats by the dissipation number Di, which sets the background density variation where |ds/dr| r i is the dimensional reference entropy gradient at the inner boundary. The background entropy gradient is prescribed analytically and sets convective (unstable stratified, ds/dr < 0) and stable (ds/dr > 0) regions. The ds/dr-profile is given by where A SSL is the amplitude of stable stratification (basically setting the BV frequency, see where a polynomial of second order is used for the former: with Γ characterizes the nature of the ideal gas and is found by fitting p-T -curves from Nettelmann et al. (2013) with polytropic laws such that yielding a best fitting Γ = 0.513 for the deep interior. The mean density contrast is set by Di = 3, corresponding to a top-to-bottom ratio of ρ i /ρ o ≈ 38. Then the temperature and density gradient are set by a chosen entropy gradient profile and eq. 8. Fig. 1, panel b) shows the so derived background state, indicating sub-and superadiabatic temperature gradients (middle panel, upper row). The black, dashed line corresponds to an isentropic (ds/dr = 0) model for reference. Note, that the non-adiabaticity in the figure is strongly exaggerated using s = 0.5 and A SSL = 1.0. For production runs, the relative deviation from an adiabatic background must remain small ( s A SSL 1). To ensure that the a)̃� 1 � � isentropic stable layer g(r) from N13 2 nd order fit p(T) from N13 power law fit (metallic/molecular) newly defined background state is steady, the heat flux must be constant on every radius.
Hence a thermal conductivity profile (k =ρc pκ ) is set such that where Q o is the heat flux at the outer boundary. The profile ofk is shown in fig. 1 convective Rossby (Ro c ) and the stratification relative to rotation rate (I s ).
The colors indicate simulations from different groups. The green triangles mark the cases used in Takehiro and Lister (2002).
frequency and the non-dimensional stratification parameter are 3 Numerical results Eqs. 3 and 4 are solved numerically using MagIC 5.6 (Wicht, 2002;Gastine and Wicht, 2012;Schaeffer, 2013) which is modified such that the non-adiabatic background state relations (eq. 8) are integrated. The mechanical confinement yields impenetrable and free-slip conditions at both walls. The entropy boundary conditions are fixed-flux in accordance with the stratification profile (eq. 11), which sets unstable and stable regions.
For a more detailed analysis of the apparent colatitudinal variation and the length scales dependencies, a FFT transforms from v r (r, ϑ) to v r (r, ), where is the spherical harmonic degree. For a few radii, this is shown as a function of either colatitude or spectral degree in fig. 6. The radial flow intensity drops slightly before reaching the lower SSL edge at r lb /r o = 0.57 (top left), though with weak modulation along colatitude ϑ.
At the edge (blue line), the damping is strongest close to the rotation axis and around the equatorial region. In between rotational penetration somewhat reduces the damping of radial flows (profiles are closer to each other). This is the expected behavior, as the

penetration depth
For a qualitative assessment we assume that the penetration depth of radial flows is a linear function of the radial flow intensity itself (Takehiro and Lister, 2001)  varying along colatitude. This implies, that v r (r, ϑ) ∝ exp −r/δ , where the penetration depth δ is defined by the e-fold decay scale height Fig. 3, panel b) show that the radial profiles v r drops exponentially in the vicinity of the neutral buoyancy radii located at r lb /r o = 0.57 and r ub /r o = 0.67, respectively. The penetration depth is calculated as a function of radius using eq. 19 while averaging over the colatitude ϑ 20) or as function of colatitude taken at the lower edge of the stratified region Finally, we also investigate the length scale dependence in terms of the spherical harmonic degree as used in the linear theory study by Takehiro and Lister (2001): where v r (r, ) = FFT {v r (r, ϑ)} is the Fourier-transformed. is reduced. All models seem to share significant higher penetration depths at midlatitudes. Inside the TC and in equatorial regions the penetration depth is typically smaller. However, the penetration depth is bounded by the viscosity (horizontal green lines, δ ν = √ E) and apparently cannot be reduced far below this value. We conclude that the penetration depth at intermediate colatitudes of 30 Takehiro and Lister (2001) concluded that the stable layer acted as a spectral filter with low pass nature, more precisely they predict that δ should decay like 1/| | . Fig. 8 shows δ for the nine cases from group 4 (see tab. 1 and fig. 7). The blue horizontal lines denote the mean value listed in tab. 1, the green one marks the thickness of the Ekman layer.
It can be seen, that penetration depth clearly decays for smaller length scales even for cases with rather large inertia and weak stratification (top right). As shown for the case 4.10 in fig. 6, the largest (small ) and smallest length scales (large ) contain little kinetic energy and are excluded from the analysis. Therefore, power laws of the form Finally, the dependence of the mean penetration depth on the model parameters (Ra, E, P r and A SSL ) is studied. The penetration depth at the lower and upper boundary of the SSL are indicated by δ + r and δ − r . Those are derived by taking the minimal value of δ r in the vicinity (±0.01r o ) of r lb and the maximal value of −δ r near r ub (see also fig. 9). Fig. 10 shows δ + r and δ − r for all models from tab. 1. It can be clearly seen, that δ r decreases with increasing stratification and settles somewhat below δ r = 0.01. This lower boundary is due to viscosity and hence depending on the Ekman number. For all studied Ekman numbers δ r can be decreased to values approximatively the Ekman layer thickness (δ E = E 1/2 ) suggesting that our numerical setup is only limited by the viscosity. Models with δ r larger than the width of SSL (0.1 r o ) are robustly constrained, as the mean radial flow intensity drops noticeable across the SSL, yet less than e-fold.
As shown by the dark-blue profile in the left panel of fig. 3, even models with neutral buoyancy gradient (A SSL = 0) in the SSL create a strong dimple in the poloidal energy.
Apparently δ + r and δ − r are nearly identical, suggesting that there is no major dynamic difference between a scenario where the stable region is on top of the convective one and the opposite situation. As the Coriolis force clearly dominates the flows in convective regions (Ro c < 1), it appears reasonable to expect that I s = τ strat /τ rot is the appropriate force ratio to estimate the penetration depth ( fig. 10, left). However, this does not seem Previous studies targeted to the radiative-convective boundary in the Sun, such as Hurlburt et al. (1994); Brummell et al. (2002); Rogers et al. (2006) controlled the degree of stratification by reducing or increasing the local polytropic index n i below or above the adiabatic one n a . Using the polytropic relations (e.g. ρ = ρ o ξ n , where ξ is a function of r), the deviation from the adiabat can be expressed by logarithmic density gradients and polytropic indices where index a refers to the adiabat. Consequently, the stratification measure S defines the ratio of subadiabatic (positive) to superadiabatic (negative) entropy gradients and is compatible with our stratification amplitude S ≈ A SSL . Whereas the numerical models in the previous study cover stratifications between S = 1 and S = 30, our parameter space with respect to A SSL spans from 0 to 750 (see tab. 1). Further the scaling above suggested that δ r ∝ I s Ro c = 2 g(r) ds dr Hence, our results show that the penetration depth increases with the square root of the ratio of positive and negative entropy gradients. Firstly reported by (Hurlburt et al., 1994), comparable models found different scaling exponents, which changed as a function of S. Typically for small S the penetration depth scales with δ r ∝ S −1 , but for larger S with a significantly shallower S −1/4 (Hurlburt et al., 1994). Other studies only found the weak scaling throughout the studied S-range (Brummell et al., 2002;Rogers et al., 2006).
For guidance, both scaling laws are drawn into fig. 10, right panel, thin grey dashed lines. Indeed those predictions seem to approximate the data in local sub-regions of the parameter space. However, the general trend appears much better reproduced by a single scaling law with an exponent of roughly S −1/2 .

summary and discussion
We have performed an extensive numerical modeling campaign of rapidly rotating convection in a spherical shell that exhibits a sandwiched stable stratified layer between two convective zones. Such a system is most suitable for Saturn, where the H/He demixing generates a compositional gradient around mid-depth and though suppresses thermal convection locally. Our numerical models cover a large fraction of the appropriate parameter space, e.g. in the terms of the stratification relative to rotational forces 2Ω/N , where N is the Brunt-Väisälä frequency or relative to the inertia characterized by Ro c (see also tab. 1).
Stratification acts directly to suppress radial flows, but radial flows can pierce across the interface by different mechanisms. The classic overshooting penetration depends on the inertia a fluid parcel has gained from buoyancy instabilities. Furthermore, rapid rotation and the spherical geometry allows convective columns to be vertically extended into the stratified regions in accordance with the Taylor-Proudman theorem. We investigate and distinguish the potential effects of both penetration mechanisms by altering the leading order force balances. Our results indicate that the mean penetration depth depends linearly on the square root of the ratio of destabilizing and stabilizing entropy gradients (δ r ≈ A −1/2 SSL ). This builds on the rather diverse results derived in the solar context (Zahn, 1991;Hurlburt et al., 1994;Brummell et al., 2002;Rogers et al., 2006), where different scaling exponents with stronger scaling for weaker stratification were reported. The exponents in those studies range from −1/4 to −1, but the robust −1/2-scaling suggested in the present study reproduces the results convincingly. In addition, and in contrast to the expectations from Takehiro and Lister (2001), the magnitude of the penetration depth appears largely unaffected by the rotation.
The Coriolis force, however, adds a characteristic latitudinal modulation of the penetration depth amplifying to ±(40 − 50)% of the mean value, where the strongest penetration resides at mid-latitudes (outside TC). This clearly reflects the action of columnar structures vertically extended by the Coriolis force and hence represents rotational penetration. Therefore, latitudinal density gradients are expected to build up at the edges of stratified regions. Such gradients are known to drive baroclinic instabilities, e.g. thermal winds.
Alongside original studies using linear waves, such as Takehiro and Lister (2001), a spectral (length scale) dependence of the penetration depth is expected and this is indeed found in our strongly non-linear simulations. The decay exponent γ approaches unity, what confirms the expected value from linear theory (Takehiro and Lister, 2001). The SSL though clearly shows the quality of a low pass filter. Consequently, radial flow in the SSL, even though weak in amplitude, must be expected to be of larger length scale than the convective flows pummeling into the SSL. However, the striking agreement with the theoretical predictions is limited to an intermediate length scale range.
As a peculiar result, we have identified wave-like radial flows that seem to somewhat resist the inverse buoyancy forces in the stratified layer. Those are non-geostrophic, of large length scale, weak in amplitude and predominantly driven in the equatorial regions of the stratified layer. We showed that the required additional source of radial flow is likely not due to the length scale filtering effect of the rotational penetration or due to inertia penetration, but is potentially powered from the convergence of horizontal flows in the shape of gravito-inertial waves.
Finally, we have shown that viscous effects only play a role for very shallow penetration depth and at the smallest length scales. The overall results should not depend on the fact that we have only explored an Ekman number rage limited by the numerical constraints.
More important is the spherical geometry. We speculate that a thin stable stratified layer will be significantly penetrated by Coriolis force induced columnar flows vertically extending into the SSL and hence more closely obey the linear theory by Takehiro and Lister (2001).
As an outlook for future work, a thorough investigation of the gravito-inertial waves generated in the SSL, the emergence of thermal winds due to latitude-dependent penetrative convection and consequently how differential rotation is paused through stratified zones. In addition, the effect on a dynamo process originating from underneath will help to understand the dynamics and magnetic fields of partly stratified rapidly rotating convection, which is most suitable for the planetary atmospheres and interiors.