# Effects of stratification in spherical shell convection

###### Abstract

We report on simulations of mildly turbulent convection in spherical wedge geometry with varying density stratification. We vary the density contrast within the convection zone by a factor of 20 and study the influence of rotation on the solutions. We demonstrate that the size of convective cells decreases and the anisotropy of turbulence increases as the stratification is increased. Differential rotation is found to change from anti-solar (slow equator) to solar-like (fast equator) at roughly the same Coriolis number for all stratifications. The largest stratification runs, however, are sensitive to changes of the Reynolds number. Evidence for a near-surface shear layer is found in runs with strong stratification and large Reynolds numbers.

4350 \Yearpublication2011 \Yearsubmission2010 \Month9 \Volume332 \Issue1 \DOI10.1002/asna.201012345

2012 Jan 12

## 1 Introduction

Numerical simulations of turbulent convection in spherical geometry have become a standard tool in the study of differential rotation and magnetism in the solar and stellar context (Miesch & Toomre 2009). The current state of the art models can reproduce many aspects of the solar internal rotation (e.g. Miesch et al. 2006), and large-scale oscillatory dynamo action occurs when rotation is rapid enough (Brown et al. 2010, 2011). However, reproducing the solar cycle has turned out to be elusive, even though large-scale oscillatory fields are now seen in some simulations with the solar rotation rate, too (e.g. Ghizaru et al. 2010; Racine et al. 2011). The reason for the remaining discrepancies may lie in the fact that much of the physics have either to be simplified or neglected altogether due to severe numerical constraints (e.g. Käpylä 2011). Furthermore, the simulations that are being carried out are so demanding that often only a single or a few representative cases can be done. Although many results, such as the change from anti-solar (slow equator) to solar-like (fast equator) differential rotation as the rotation rate increases (e.g. Chan 2010; Käpylä et al. 2011a), and the appearance of mostly axisymmetric large-scale magnetic fields (e.g. Gilman 1983; Glatzmaier 1985, Browning et al. 2006; Brown et al. 2010, 2011; Käpylä et al. 2010) appear robust, their exact dependence on different simulation parameters has not been explored in detail.

Here we study the effect of density stratification on rotating spherical shell convection. Our main goal is to study how the transition from anti-solar to solar-like rotation is affected. This is relevant because most convection models in spherical shells still have rather modest density stratification in comparison to the Sun. Although most of the mass within the convection zone is located near the base, fast downflows at the vertices of convection cells originate near the surface. The effect of these downflows on angular momentum transport is yet unclear. We are also interested in the statistical properties, such as anisotropy, of turbulence as the stratification is increased.

## 2 Model

Our model is based on that used by Käpylä et al. (2010, 2011a). We model a segment of a star, i.e. a “wedge”, in spherical polar coordinates, where denote the radius, colatitude, and longitude. The radial, latitudinal, and longitudinal extents of the computational domain are given by , , and , respectively, where is the radius of the star. In all of our runs we take and .

We solve the following equations of compressible hydrodynamics in a frame of reference rotating with angular velocity ,

(1) |

(2) |

(3) |

where is the advective time derivative, is the density, is the velocity, is the specific entropy, is the temperature, and is the pressure. The fluid obeys the ideal gas law with , where is the ratio of specific heats at constant pressure and volume, respectively, and is the internal energy.

Furthermore, is the kinematic viscosity, is the radiative heat conductivity, is the unresolved turbulent heat conductivity, and is the gravitational acceleration given by

(4) |

where is the gravitational constant, is the mass of the star, and is the unit vector in the radial direction. We omit the centrifugal force in our models. The rate of strain tensor is given by

(5) |

where the semicolons denote covariant differentiation; see Mitra et al. (2009) for details. Unlike in our previous studies (Käpylä et al. 2010, 2011a), we omit stably stratified layers below and above the convectively unstable layer.

### 2.1 Initial and boundary conditions

In the initial state the atmosphere is adiabatic and the hydrostatic temperature gradient is given by

(6) |

where is the polytropic index. We use Eq. (6) as the lower boundary condition for the temperature. This gives the logarithmic temperature gradient (not to be confused with the operator ) as

(7) |

Density stratification is obtained by requiring hydrostatic equilibrium. The heat conduction profile is chosen so that radiative diffusion is responsible for supplying the energy flux in the system, and decreases rapidly within the convection zone (see, Fig. 1).

The radial and latitudinal boundaries are taken to be impenetrable and stress free, according to

(8) | |||||

(9) |

On the latitudinal boundaries we assume that the thermodynamic quantities have zero first derivative, thus suppressing heat fluxes through the boundaries.

On the upper boundary we apply a black body condition

(10) |

where is the Stefan–Boltzmann constant. In our runs we use a modified value for that takes into account that our Reynolds and Rayleigh numbers are much smaller than in reality, so is much larger and therefore the flux too high. The black body boundary for the temperature has previously been used in mean-field models of Brandenburg et al. (1992). In our runs is negligibly small near the surface so that the unresolved convective energy flux transports practically all of the energy through the upper boundary. This is similar to what is commonly used in the ASH simulations (e.g. Brun et al. 2004).

### 2.2 Dimensionless parameters

We obtain non-dimensional quantities by choosing

(11) |

where is the density at . The units of length, velocity, density, and entropy are then given by

(12) |

The simulations are governed by the Prandtl, Reynolds, Coriolis, and Rayleigh numbers, defined by

(13) | |||||

(14) |

where is the turbulent thermal conductivity in the middle of the convection zone (i.e. at ), is an estimate of the wavenumber of the energy-carrying eddies, is the thickness of the layer, and is the rms velocity, where the angle brackets denote volume averaging. In our definition of we omit the contribution from the -velocity, because its value is dominated by effects from the differential rotation (Käpylä et al. 2011a). Sometimes we show which is the fluctuating rms velocity as a function of radius and from which we have subtracted the azimuthally averaged velocities. The entropy gradient, measured at , is given by

(15) |

where , and is the pressure scale height at . Due to the fact that the initial stratification is isentropic, we quote values of from the thermally saturated state of the runs.

The energy that is deposited into the domain at the base is controlled by the luminosity parameter

(16) |

where is the constant luminosity, and is the energy flux imposed at the bottom boundary. Furthermore, the stratification is determined by the normalized pressure scale height at the surface

(17) |

where . Similar parameter definitions were used by Dobler et al. (2006). We use three different values of which result in density contrasts of 5, 30, and , respectively (see Fig. 1). Now the convection zones span between roughly 2.5 and 7.5 pressure scale heights.

The simulations were performed using the Pencil Code^{1}^{1}1http://pencil-code.googlecode.com/, which uses
sixth-order explicit finite differences in space and a third-order
accurate time stepping method; see Mitra et al. (2009) for further
information regarding the adaptation of the Pencil Code to
spherical coordinates.

Run | |||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

A0 | 2.5 | 0.025 | 41 | —- | 0.116 | —- | —- | —- | —- | —- | |||

A1 | 2.5 | 0.020 | 33 | 1.37 | 0.003 | 0.944 | 0.78 | ||||||

A2 | 2.5 | 0.016 | 26 | 3.48 | 0.000 | 0.798 | 1.06 | ||||||

A3 | 2.5 | 0.014 | 23 | 5.97 | 0.000 | 0.889 | 1.06 | ||||||

A4 | 2.5 | 0.013 | 21 | 8.78 | 0.000 | 0.914 | 1.05 | ||||||

B0 | 5 | 0.027 | 22 | —- | 0.037 | —- | —- | —- | —- | —- | |||

B1 | 5 | 0.026 | 22 | 1.04 | 0.005 | 0.906 | 0.79 | ||||||

B2 | 5 | 0.025 | 20 | 2.24 | 0.004 | 0.879 | 1.06 | ||||||

B3 | 2.5 | 0.024 | 40 | 4.54 | 0.001 | 0.786 | 1.06 | ||||||

B4 | 2.5 | 0.021 | 36 | 7.61 | 0.001 | 0.581 | 1.05 | ||||||

C0 | 5 | 0.017 | 26 | —- | 0.075 | —- | —- | —- | —- | —- | |||

C1 | 5 | 0.017 | 25 | 0.91 | 0.003 | 0.772 | 0.92 | ||||||

C2 | 5 | 0.016 | 25 | 1.85 | 0.002 | 0.760 | 0.90 | ||||||

C3 | 5 | 0.015 | 22 | 4.12 | 0.003 | 0.333 | 1.00 | ||||||

C4 | 5 | 0.013 | 19 | 7.22 | 0.003 | 0.250 | 1.00 | ||||||

D0 | 2 | 0.019 | 73 | —- | 0.078 | —- | —- | —- | —- | —- | |||

D1 | 2 | 0.019 | 72 | 0.79 | 0.002 | 0.826 | 0.94 | ||||||

D2 | 2 | 0.018 | 66 | 1.72 | 0.002 | 0.956 | 0.81 | ||||||

D3 | 2 | 0.018 | 68 | 3.35 | 0.002 | 0.589 | 0.99 | ||||||

D4 | 2.5 | 0.018 | 54 | 5.11 | 0.001 | 0.718 | 1.02 |

## 3 Results

We have performed four sets of simulations differing by their density stratification (see Table 1) and Reynolds number. We vary the rotation rate within each set so that the Coriolis number changes by roughly an order of magnitude. We increase the gravity by a factor of in Sets C and D in order to limit the Mach number to roughly 0.1 near the surface. The grid resolution in Sets A–C is . In Set C this means that the ratio of the pressure scale height to the radial grid spacing at the surface is , which is on the limit of resolving the structure. We have remeshed snapshots from the saturated states of the runs in Set C to double resolution (Set D) where we are also able to increase the Reynolds number.

### 3.1 Flux balance

In contrast to our earlier studies using a polytropic setup with (Käpylä et al. 2010, 2011a), we now use a setup in which convection transports the majority of the flux. This is achieved by decreasing the heat conductivity within the convection zone and introducing a turbulent heat conductivity which is responsible for unresolved convective transport of heat (e.g. Chan & Sofia 1996; Brun et al. 2004). We apply a constant value of , of the order of the kinematic viscosity , in the bulk of the convection zone () and an order of magnitude larger value above in order to transport the flux through the upper boundary. Below , goes smoothly to zero, see Fig. 1.

To verify that the system is in thermal equilibrium, we consider the radiative, convective, kinetic, viscous, and turbulent energy fluxes, defined as

(18) | |||||

(19) | |||||

(20) | |||||

(21) | |||||

(22) |

where the averages are taken over and . Representative results from Run B3 are shown in Fig. 2. Radiative diffusion transports the total flux through the lower boundary and decreases rapidly as a function of . The radiative flux is less than 10 per cent above . The flux due to resolved convection is responsible for transporting the majority of the luminosity within the convection zone. The flux of kinetic energy is directed downwards and is responsible for roughly 10 per cent of the flux near the surface. Note that the maxima of and are significantly larger in the non-rotating cases. The unresolved turbulent flux is small in the bulk of the convection zone and carries the flux out through the outer boundary. The viscous flux is small in all of our runs. The flux balance is similar to the ASH simulations (e.g. Brun et al. 2004; Miesch et al. 2008) which employ stratification from a 1D solar model and somewhat different profiles of the diffusion coefficients.

### 3.2 Properties of convection

Visualizations of the radial velocity near the surface of the star for all of our models are shown in Fig. 3. We find that as the stratification increases the size of convection cells decreases. This is a consequence of the decreasing pressure scale height near the surface. In the case of the highest stratification, Sets C and D, the granulation pattern is similar to the high-resolution run reported by Miesch et al. (2008) with a comparable stratification. As noted above, the resolution in Set C is close to critical in resolving the stratification properly which is also manifested by numerical artefacts in Fig. 3. The higher resolution runs in Set D, however, are well behaved and show similar convection patterns in the non-rotating and slowly rotating cases.

Furthermore, as the rotation rate is increased, one sees the formation of a cartridge belt-like pattern that is also known as banana cells. However, these structures become less pronounced at larger stratification in Set C. In Set D, on the other hand, strong banana cells are again observed. The main difference between Sets C and D is that in the latter the Reynolds number is modestly increased. The smaller size of convection cells at high stratification leads to a smaller “effective” Reynolds number based on the horizontal eddy size, which is not well reflected by our definition of . It is possible that in Run C4 this effective is below critical to excite the formation of banana cells and strong prograde differential rotation.

Another way to see the difference between weak and strong stratification is to visualize the flows in the meridional plane. In Fig. 4 we show a cut of the radial velocity at from Runs A0, B0, and D0 in the saturated state. Whereas the downflows in Run A0 go easily through the whole convection zone, smaller-scale structures originating from the surface appear already in Run B0. For the strongest stratification the anisotropy of the flow is clear also to the naked eye with smooth large-scale flows in the deep layers and small-scale irregular flows near the surface. However, the strongest downflows are still able to go all the way through the convection zone.

The squares of the fluctuating velocity components are shown in Fig. 5. For the weakest stratification the profiles are almost symmetrical with respect to . The apparent asymmetry of the and velocities is likely due to the large size of the convective cells in comparison to the domain size. For larger stratification the velocities increase near the surface and the anisotropy of the horizontal velocities decreases. An important effect that follows is that the convective turnover time, defined as

(23) |

where is the local pressure scale height, changes substantially between the bottom and the surface. In the Run A0 is almost constant in the whole layer whereas in Runs B0 and C0 it varies by factors of 9.1 and 31, respectively. Provided that mixing length arguments hold, the rotational influence on the flow, measured by the Coriolis number, has the same variation as a function of radius.

We define the vertical anisotropy parameter as

(24) |

where the averages are taken over and , and the primes denote that -averaged mean velocities are subtracted. For weak stratification, is at most 0.4 in the middle of the convection zone. Note that at the boundaries due to the impenetrable boundary conditions. For Runs B0 and C0, the turbulence is more anisotropic and the point at which changes from negative to positive moves closer to the boundaries. The anisotropy measured by also increases and in Runs B0 and C0 it is more than double the value of Run A0. In comparison to forced turbulence simulations of Brandenburg et al. (2011a), our results correspond best to cases with small scale separation, i.e. large-scale forcing. This is consistent with large convective cells that span the whole depth of the convection zone. The importance of is that in rotating turbulence, it acts as a source for the vertical -effect (e.g. Rüdiger 1989; Käpylä & Brandenburg 2008) which drives radial differential rotation.

### 3.3 Differential rotation

The differential rotation profiles from all runs with are shown in Fig. 6. We find that an anti-solar differential rotation pattern with strong meridional circulation forms at the lowest rotation rates. As increases, equatorial acceleration gradually develops. However, in many cases (e.g. Runs B3, B4, C4, D4) there is a minimum of the local angular velocity, , at mid-latitudes and a polar vortex at high latitudes. Similar profiles have been reported by Miesch et al. (2000) and Elliott et al. (2000). Large-scale vortices arise in Cartesian convection simulations at sufficiently high rotation rates (Chan 2007; Käpylä et al. 2011b; Mantere et al. 2011). It is unclear whether the polar vortices in spherical geometry are related to the vortex instability but it is an intriguing possibility.

We quantify the horizontal differential rotation by the parameter

(25) |

where , and . The results for Sets A to D, along with corresponding data from Käpylä et al. (2011a) are shown in Fig. 8. We find that for , is the largest for the smallest stratification. The results for Sets B and D seem to converge for more rapid rotation and produce solar-like rotation () for .

In Set C, however, the transition to a solar-like profile does not yet occur in the parameter range studied here, although the largest Coriolis number is of the order of 7. Another factor that comes into play is the fact that the effective Reynolds number based on the typical scale of convection cells is reduced in the runs with the largest stratification. This seems to be confirmed by the simulations in Set D, which are the higher Reynolds number counterparts of the runs in Set C, although for Run D4, is still negative. This is surprising given the profile seen in Fig. 6 which clearly shows a rapidly rotating equator. The discrepancy is due to a sharp negative radial gradient of near the surface at the equatorial regions in Run D4 (see Fig. 7). This is similar to the near-surface shear layer observed in the Sun (e.g., Benevolenskaya et al. 1999).

The formation of a negative near-surface shear layer only occurs for strong stratification. Furthermore, Runs C4 and D4 (where is twice as large), suggests that negative near-surface shear also requires large Reynolds numbers. It is still unclear whether the current simulations can really capture the physics of the solar near-surface shear layer (e.g. Miesch & Hindman 2011), but the present results might indicate a path that is worth following. Measuring from a little deeper down at gives which is similar to the results from Runs B3 and B4 (see Table 1).

We also note that the results from Käpylä et al. (2011a) roughly fall in line with Sets B and D for . However, in the rapid rotation regime, from Käpylä et al. (2011a) is consistently larger than in the current results with the exception of Set A. This is probably due to the difference in the setups, i.e. we omit here the stably stratified overshoot layer below the convection zone and the isothermal cooling layer near the surface.

It is customary to represent the latitudinal profile of the angular velocity in terms of Gegenbauer polynomials, i.e.,

(26) |

(We use here a definition where all associated Legendre polynomials with odd values are positive, i.e., , for example.) We have determined the coefficients , , and via a fitting procedure using as weighting factor to put more emphasis on the equatorial regions. The resulting coefficients are given in the last 3 columns of Table 1. Note that changes sign for . This is also consistent with Fig. 9, where we plot the latitudinal profiles of together with their fits from Eq. (26). The change of sign of gives a more robust indicator of the transition from antisolar to solar-like rotation than just the parameter.

## 4 Conclusions

We study turbulent convection in spherical shells with varying density stratification. We find that the typical size of convection cells decreases as the stratification increases, which is in accordance with mixing length arguments. At the same time the anisotropy of turbulence increases and the turnover time varies by more than an order of magnitude from the base to the top of the convection zone.

Although convection seemingly changes greatly as the stratification increases, the rotation profiles and their qualitative trend as a function of the Coriolis number change surprisingly little. However, we find that the results for the largest density stratification are sensitive to changes of the Reynolds number and that apparently smaller latitudinal differential rotation for large Coriolis numbers is obtained for large stratification. Measuring the differential rotation simply as a difference between the surface values of at the equator and at high latitudes turns out to give misleading results for our high-resolution runs with the largest stratification. This is due to the self-consistent generation of a sharp radial gradient of near the surface, reminiscent of the near-surface shear layer in the Sun. A similar feature is discernible in the highly stratified simulations of Bessolaz & Brun (2011). However, it remains to be seen whether this is a robust feature.

Another interesting aspect of increasing stratification is related to the generation of large-scale magnetic fields, because some of the contributions to the -effect of mean-field dynamo theory are proportional to density stratification (e.g. Krause & Rädler 1980). Furthermore, the negative magnetic pressure instability (e.g. Brandenburg et al. 2011b; Käpylä et al. 2011c), that can lead to magnetic field concentrations of the form of active regions, becomes stronger when stratification increases. We plan to revisit these issues in forthcoming papers.

###### Acknowledgements.

We thank the referee for the suggestion to add the expansion in terms of Gegenbauer polynomials. Computational resources granted by CSC – IT Center for Science, who are financed by the Ministry of Education, and financial support from the Academy of Finland grants No. 136189, 140970 (PJK), 218159 and 141017 (MJM), and the European Research Council under the AstroDyn Research Project 227952 (AB) are acknowledged.## References

- [1999] Benevolenskaya, E.E., Hoeksema, J.T., Kosovichev, A.G., Scherrer, P.H.: 1999, ApJ 517, L163
- [2011] Bessolaz, N., Brun, A.S.: 2011, ApJ 728, 115
- [1992] Brandenburg, A., Moss, D., Tuominen, I.: 1992, A&A 265, 328
- [2011] Brandenburg, A., Rädler, K.-H., Kemel, K.: 2011a, A&A (submitted), arXiv:1108.2264
- [2011] Brandenburg, A., Kemel, K., Kleeorin, N., Mitra, D., Rogachevskii, I.: 2011b, ApJL 740, L40
- [2010] Brown, B.P., Browning, M.K., Miesch, M.S., Brun, A.S., Toomre, J.: 2010, ApJ 711, 424
- [2011] Brown, B.P., Miesch, M.S., Browning, M.K., Brun, A.S., Toomre, J.: 2011, ApJ 731, 69
- [2006] Browning, M.K., Miesch, M.S., Brun, A.S., Toomre, J.: 2006, \apj 648, L157
- [2004] Brun, A.S., Miesch, M.S., Toomre, J.: 2004, \apj 614, 1073
- [1996] Chan, K.-L., Sofia, S.: 1996, \apj 466, 372
- [2007] Chan, K.-L.: 2007, AN 328, 1059
- [2010] Chan, K.-L.: 2010, IAUS 264, 219
- [2006] Dobler, W., Stix, M., Brandenburg, A.: 2006, \apj 638, 336
- [1999] Elliott, J.R., Miesch, M.S., Toomre, J.: 2000, \apj533, 546
- [2010] Ghizaru, M., Charbonneau, P., Smolarkiewicz, P.K.: 2010, ApJ 715, L133
- [1983] Gilman, P.A.: 1983, ApJS 53, 243
- [1985] Glatzmaier, G.A.: 1985, ApJ 291, 300
- [2008] Käpylä, P.J., Brandenburg, A.: 2008, A&A 488, 9
- [2010] Käpylä, P.J., Korpi, M.J., Brandenburg, A., Mitra, D., Tavakol, R.: 2010, AN 331, 73
- [2011] Käpylä, P.J., Mantere, M.J., Guerrero, G., Brandenburg, A., Chatterjee, P.: 2011a, A&A, 531, A162
- [2011] Käpylä, P.J., Mantere, M.J., Hackman, T.: 2011b, \apj, 742, 34
- [2011] Käpylä, P.J., Brandenburg, A., Kleeorin, N., Mantere, M. J., Rogachevskii, I.: 2011c, MNRAS (submitted), arXiv:1104.4541
- [2011] Käpylä, P.J.: 2011, AN, 332, 43
- [1980] Krause, F., Rädler, K.-H.: 1980, Mean-field Magnetohydrodynamics and Dynamo Theory, Pergamon Press, Oxford
- [2011] Mantere, M.J., Käpylä, P.J., Hackman, T.: 2011, AN (submitted), arXiv:1109.4317
- [2000] Miesch, M.S., Elliott, J.R., Toomre, J., et al.: 2000, \apj 532, 593
- [2006] Miesch, M.S., Brun, A.S., Toomre, J.: 2006, \apj 641, 618
- [2008] Miesch, M.S., Brun, A.S., DeRosa, M.L., Toomre, J.: 2008, \apj 673, 557
- [2009] Miesch, M.S., Toomre, J.: 2009, AnRFM, 41, 317
- [2011] Miesch, M.S., Hindman, B.W.: 2011, ApJ (submitted), arXiv:1106.4107
- [2009] Mitra, D., Tavakol, R., Brandenburg, A., Moss, D.: 2009, ApJ 697, 923
- [2011] Racine, ., Charbonneau, P., Ghizaru, M., Bouchat, A., Smolarkiewicz, P. K.: 2011, \apj735, 46
- [1989] Rüdiger, G.: 1989, Differential Rotation and Stellar Convection: Sun and Solar-type Stars (Akademie Verlag, Berlin)