A&A 459, 935-944 (2006)
DOI: 10.1051/0004-6361:20065519

Reynolds stresses and meridional circulation from rotating cylinder simulations

C. Hupfer1 - P. J. Käpylä1,2,3 - M. Stix1


1 - Kiepenheuer-Institut für Sonnenphysik, Schöneckstr. 6, 79104 Freiburg, Germany
2 - Observatory, University of Helsinki, PO Box 14, 00014 University of Helsinki, Finland
3 - Astronomy Division, Department of Physical Sciences, PO Box 3000, 90014 University of Oulu, Finland

Received 28 April 2006 / Accepted 4 August 2006

Abstract
Aims. The latitude and depth dependences of Reynolds stresses are obtained from numerical simulations of a solar-type convection zone where the star is assumed to rotate with a uniform angular velocity.
Methods. A two-dimensional model, using a cylindrical annulus with axis perpendicular to the axis of rotation, is introduced as an approximation of a spherical section along the meridional plane. We solve the fully compressible Navier-Stokes equations numerically and use Cartesian and cylindrical geometries to simulate convection under the influence of rotation.
Results. For moderate Coriolis numbers both models yield strong extrema of the Reynolds stress component $Q_{\theta \varphi }$ at low latitudes near the equator, and a meridional cell pattern is found in the cylindrical model. For Coriolis numbers larger than about 10 the flow becomes aligned with the direction of the rotation axis.

Key words: hydrodynamics - convection - Sun: rotation

1 Introduction

The differential rotation of the Sun and solar-type stars is still under examination both in observation and theory. There are three contributions capable of redistributing the angular momentum in a convection zone: meridional circulation, turbulent transport, and magnetic fields. In the present work the treatment is restricted to the non-magnetic case, mainly based on the assumption that the magnetic field in the present Sun has no great effect on the internal rotation - this may not be true for young stars, and therefore for the young Sun as well, e.g. Pinsonneault et al. (1989).

In spherical coordinates $(r,\theta, \varphi)$ the equation of angular momentum transport reads, see e.g. Rüdiger (1989),

 \begin{displaymath}{\frac{\partial }{\partial t} \left({ \rho \Omega s^2 }\right... ...}} - \rho s^2 \nu \vec{\nabla}_{{\rm m}} \Omega}\right) = 0}, \end{displaymath} (1)

where $\Omega$, ${\nu}$, ${s = r \sin \theta}$, and ${\vec{u}_{{\rm m}} = ( u_r, u_\theta,0)^T}$ denote the angular velocity, kinematic viscosity, lever arm, and meridional velocity, respectively, and axisymmetry is assumed; ${\vec{\nabla}_{{\rm m}} = \vec{e}_r \frac{\partial} {\partial r} + \vec{e}_{\theta} \frac{1}{r}\frac{\partial}{\partial \theta}}$ is a meridional gradient.

Since turbulent diffusion is much greater than molecular diffusion, a statistical approach via mean-field theory is appropriate, where the time evolution of, say, the mean velocity ${\langle \vec{u} \rangle}$ is calculated using the averaged Navier-Stokes equation, which is often called the Reynolds equation. Averaging of the equation of angular momentum transport introduces the Reynolds stress tensor ${\underline{Q}}$ as a possible generator of differential rotation (Rüdiger 1989). The Reynolds stresses are defined as

\begin{displaymath}{Q_{ij} = \langle {u_i}' {u_j}'\rangle}, \end{displaymath} (2)

where the brackets denote averaging and ${{u_i}' = u_i - \langle{u_i} \rangle}$ is the fluctuation of the ith velocity component. For the mean angular momentum transport equation one obtains (neglecting the stress due to molecular diffusion)

 \begin{displaymath}{\frac{\partial }{\partial t}\left({\rho \Omega s^2}\right) +... ...angle {u'_{\varphi} \vec{u}'_{{\rm m}}\rangle }}\right) = 0} , \end{displaymath} (3)

which contains the Reynolds tensor components ${Q_{r\varphi}}$, responsible for radial, and $Q_{\theta \varphi }$, responsible for horizontal transport of angular momentum. If this transport is from higher to lower latitudes towards the equator, then ${Q_{\theta \varphi} > 0}$ on the northern hemisphere of the Sun. Similarly, ${Q_{r\varphi} > 0}$ means outward transport of angular momentum.

The quantity $Q_{\theta \varphi }$ has been measured, for example, by Ward (1965) from the motion of sunspot groups, by Schröter & Wöhl (1976) from the Ca+ network, and more recently by Vrsnak et al. (2003) also for higher latitudes from coronal bright points. Typical values for $Q_{\theta \varphi }$ are of the order ${ 1~000\; {\rm m^2~s^{-2}}}$. The experimental determination of the horizontal correlation is difficult due to the large rotation velocity ${u_{\varphi} \approx 2000\; {\rm m~s^{-1}}}$, compared to the meridional velocity ${\vert u_{\theta}\vert \approx 10 \; {\rm m~s^{-1}}}$. There is no observational information on ${Q_{r\varphi}}$, nor on the depth dependence of both correlations. Reynolds stress values for other stars are unknown.

Meridional circulation appears in both forms of the angular momentum transport equation. It seems to be confirmed by helioseismology that there is a meridional circulation pattern down to a depth of at least ${15~000\;{\rm km}}$(Zhao & Kosovichev 2004). Direct observation of the surface flows yields results indicating a northward flow on the northern hemisphere of the order of ${-10\;{\rm m~s^{-1}}}$. From recent observations of K-type giants (Weber et al. 2005) there are indications of strong meridional surface velocities, of order ${1000\;{\rm m~s^{-1}}}$, and anti-solar differential rotation.

1.1 The closure model of Boussinesq and the ${\Lambda }$ effect

From a theoretical point of view it is difficult to calculate Reynolds stresses due to the fact that the derived expressions generally involve higher order correlations. This gives rise to the closure problem which requires that the higher order terms need to be related to lower order correlations or mean quantities in order to make the equations tractable. An early ansatz comes from Boussinesq (1897):

 \begin{displaymath} Q_{ij} = -\nu_{\rm t}\left({\frac{\partial\langle u_i\rangl... ... +\frac{\partial\langle u_j\rangle}{\partial x_i}}\right)\cdot \end{displaymath} (4)

Using this in spherical coordinates one obtains
  
$\displaystyle {{Q}_{r\varphi}}={-\nu_{\rm t}r\sin\theta\frac{\partial\Omega}{\partial r}} ,$     (5)
$\displaystyle { {Q}_{\theta \varphi}} = {- \nu_{\rm t} \sin \theta \frac{ \partial \Omega}{\partial \theta}}\cdot$     (6)

These expressions resemble the terms in Eq. (1) if the molecular diffusion ${\nu}$ is replaced by the turbulent diffusion $\nu_{\rm t}$; but the Boussinesq terms neither provide generation of differential rotation if the star is rotating uniformly initially, nor can they maintain the angular velocity gradient against the strong diffusion implied by the turbulent viscosity in a star which already rotates differentially. Moreover, the observed signs of  $Q_{\theta \varphi }$ and $\partial\Omega/\partial\theta$ are incompatible with Eq. (6). There must be an additional term that is non-diffusive and acting even in a uniformly rotating star. The non-diffusive term is usually called the ${\Lambda }$-effect. It can be written (Rüdiger 1989) as
  
$\displaystyle {Q_{r\varphi}} = { - \nu_{\rm V} r \sin \theta \frac{ \partial \Omega}{\partial r} + \Lambda_{\rm V} \Omega \sin \theta }~,$     (7)
$\displaystyle {Q_{\theta\varphi}} = { - \nu_{\rm H} \sin \theta \frac{ \partial \Omega}{\partial \theta} + \Lambda_{\rm H} \Omega \cos \theta } ,$     (8)

where a distinction between a vertical (radial) and horizontal ${\Lambda }$ effect with correspondingly different turbulent viscosities ${\nu_{\rm V}, \nu_{\rm H}}$ has been made.

1.2 Numerical approaches

Reynolds stresses have been determined already by Hathaway & Somerville (1983) with box simulations in the moderate rotation regime at a near-equator latitude; they calculated both ${Q_{r\varphi}}$ and $Q_{\theta \varphi }$. More recently similar box calculations were performed by Pulkkinen et al. (1993), Chan (2001), Käpylä et al. (2004), Hupfer et al. (2005), and Rüdiger et al. (2005) to explore the Reynolds stresses and the ${\Lambda }$ effect, and the possibility of a latitude dependency of the turbulent heat transport ("pole-equator temperature difference''). Robinson & Chan (2001) made simulations in a section of a spherical shell, and global spherical simulations were run by Miesch et al. (2000) and Elliott et al. (2000) with a spherical harmonics code.

In the subsequent two sections we present our models and the equations that have been solved; in Sect. 4 we briefly describe the numerical method. Sections 5 and 6 contain the results and a discussion.

2 The models

In this section we briefly describe our two models. Both have an initially piecewise polytropic stratification in hydrostatic equilibrium, which yields constant thermal conduction coefficients and linear temperature profiles in each layer. The total vertical extent of the box models and the cylinder models CA is 2.65 density scale heights initially; for the cylinder models CB and CC it is 2.9 density scale heights.

  
2.1 The box model

The calculations of this model are carried out within a rectangular box that consists of three layers, as shown in Fig. 1. The middle layer, of thickness d, is convectively unstable; the lower layer, of thickness 0.68d, is stably stratified, but the flow is allowed to overshoot into it from above. The upper layer, of thickness 0.32d, is also stable; in Hupfer et al. (2005) we presented results for Reynolds stresses from models in which we cooled this layer effectively, so that it remains nearly isothermal; we considered this upper layer mainly as a technical means to remove the heat that enters the box from below. The present calculations, however, have shown that, at least for the parameter regime under consideration, such cooling is not necessary; it was therefore omitted in the present simulations. In order to explore the influence of rotation we vary the strength of the Coriolis force and place the box at different latitudes. The box is assumed to be small so that the co-moving unit vectors ${\hat{\vec{e}}_x, \hat{\vec{e}}_y, \hat{\vec{e}}_r}$ are identified with the local spherical coordinate system ${\hat{\vec{e}}_\theta, \hat{\vec{e}}_\varphi, \hat{\vec{e}}_r}$, respectively. In contrast to some of the above-mentioned publications the vertical axis  ${\hat{\vec{e}}_r}$ has an outward orientation. The layer boundaries are located at ${r_0 = r_{\rm b},\ r_1 = r_0 + 0.68d,\ r_2 = r_1 + d,\ r_3 = r_{\rm t}}$, where  ${r_{\rm b},r_{\rm t}}$ denote the radial positions of the bottom and the top of the box.


  \begin{figure} \par\includegraphics[width=7.5cm,clip]{5519fg01.eps}\end{figure} Figure 1: The rectangular domain of computation, with three layers, here located at an intermediate latitude on the southern hemisphere.
Open with DEXTER

Boundary conditions for the box model

The horizontal extent of the box is 4d. The boundary conditions in the horizontal directions (x,y) are periodic for all variables. The energy flux and thereby the temperature gradient is prescribed at the bottom of the box. Contrary to the simulations presented in Käpylä et al. (2004) or Hupfer et al. (2005) we apply this condition at the top of the box as well, with the same heat flux that enters the box at the bottom leaving it at the top. The heat conduction coefficients in the lower and upper stable layer are the same, resulting in the same temperature gradient at the vertical boundaries. The box is assumed to be impenetrable in the vertical direction, and for the horizontal velocity we use stress-free boundary conditions:

\begin{displaymath}{\left({u_{r}}\right)_{\vert r=r_{\rm b},r_{\rm t}}} = {0}. \end{displaymath} (9)


\begin{displaymath}{\left({\frac{\partial u_{x}}{\partial r}}\right)_{\vert r=r_... ...{y}}{\partial r}}\right)_{\vert r=r_{\rm b},r_{\rm t}}} = {0}. \end{displaymath} (10)

The box model has two notable deficiencies:


  \begin{figure} \par\includegraphics[width=7.5cm,clip]{5519fg02CMJN.eps} \end{figure} Figure 2: Rotating cylinder model with three layers. The latitude ${\beta } = 90\hbox {$^\circ $ }-\theta $ is measured from the equator, and rb, rt are the radii of the bottom and the top of the vertical computation domain.
Open with DEXTER

  
2.2 The cylinder model

In order to compare whether geometry and connected flows influence the results obtained by box calculations we have developed an approach that may be considered as an intermediate step between local and global models: a cylindrical annulus is assumed to represent a spherical section in the ${r, \theta}$-plane; the cylinder axis is perpendicular to the axis of rotation (Fig. 2). We restrict ourselves to a two-dimensional model, depending on the coordinates r and $\beta$, with 3 velocity components ${u_r, u_{\beta}, u_{\phi}}$. In this setup the azimuthal coordinate of the cylindrical geometry corresponds to the meridional direction ${\theta}$, and the direction along the cylinder axis serves as an approximation of the spherical ${\varphi}$ coordinate, however, with sign changes due to the different orientations of the coordinate systems. We use the latitude $\beta$ rather than the co-latitude and the symbol ${\phi}$ as the axial direction of the cylinder model (${\phi}$ has the dimension of a length). The orientation is such that $\hat{\vec{e}}_{\theta} = - \hat{\vec{e}}_{\beta}$ and $\hat{\vec{e}}_{\varphi} = - \hat{\vec{e}}_{\phi}$ (eastern hemisphere of the corresponding embedded sphere). Thus, there are sign changes in the velocities as well, namely $u_\varphi = u_y = -u_\phi$ for the azimuthal direction, and $u_\theta = u_x = -u_\beta$ for the meridional direction. The separation into three layers, the intermediate layer boundaries, and the stratification are the same as in the box models. Efficient calculations could be restricted to one quadrant or one hemisphere. However, our simulations cover the full $360\hbox{$^\circ$ }$ range, with 721 grid points in the meridional direction and 64 grid points for the radial direction. The advantage of such "global'' calculations is the possibility to check for symmetry or anti-symmetry of flows on different quadrants or hemispheres.

Since the thickness of the convective layer is our length scale, of order 105 km in the deep convection zone, we choose in our cylinder model series CA the radii $r_{\rm b} = 5$ and $r_{\rm t} = 7$, which would correspond to ${500~000\; {\rm km}}$ and ${700~000 \;{\rm km}}$ on the Sun. The radius ratio is therefore $\eta = \frac{r_{\rm b}}{r_{\rm t}} = \frac{5}{7}$. The series CA consists of 26 calculations from no rotation to very strong rotation, with Coriolis numbers - defined by (25) below - ranging from 0 to $\approx$40. In addition, we performed simulations with either different radius ratios ${\eta}$ (models CAS), or with a convective layer of vertical extent 2d (models CB, CC) in order to compare with the CA results. In these models the upper and lower layers remain unchanged, so that the total depth is 3d, cf. Fig. 7 below. An overview of the cylinder models is given in Table 2.

A deficiency of the cylinder model is that there is no true polar region due to the Cartesian geometry along the cylinder axis. In an axisymmetric sphere there would be a converging or diverging meridional flow at the poles; this is not necessarily the case in the model of a cylinder rotating around an axis perpendicular to its own axis.

Boundary conditions for the cylinder model
Since the full circle is calculated, the mathematical periodicity assumed for the meridional coordinate x in the box models turns into a physical periodicity - there are no obstacles for a meridional flow, so we are using periodic boundary conditions also in this case. The vertical boundary conditions are the same as for the box model.

3 Equations and dimensionless parameters

We use ${\frac{D}{Dt} = \frac{\partial }{\partial t} + \vec{u}\cdot \vec{\nabla}}$ for the substantial derivative, and the Einstein summation in the case of index notation. The symbols $\Delta$, $\vec{\nabla} \cdot \vec{A }$, and $\nabla_k \Psi$ denote the Laplacian operator, the divergence of a vector $\vec{A}$, and the kth component of the gradient of a scalar function ${\Psi}$. We give the hydrodynamical equations for Cartesian and cylindrical geometry separately and use the usual symbols ${\rho, e, \vec{u}}$ for the mass density, the energy density and the velocity, respectively. Since the density varies by more than two scale heights over the radial extent of the computational domain it is preferable to use ${\ln \rho}$ rather than the density itself; this also simplifies the equations and enhances the speed of the computation.

3.1 Cartesian geometry

For the box model we use Cartesian coordinates:

   
    $\displaystyle {\frac{{\rm D} \ln \rho}{{\rm D} t}} = { - \vec{\nabla} \cdot \vec{u}},$ (11)
    $\displaystyle {\frac{{\rm D} e}{{\rm D} t}} = { -\frac{P}{\rho}\vec{\nabla} \cd... ...vec{\nabla} \cdot \left({\kappa \vec{\nabla} e}\right) + 2\nu \underline{S}^2},$ (12)
    $\displaystyle { \frac{{\rm D} u_i}{{\rm D} t}} = { - \frac{1}{\rho}\frac{\parti... ...S_{ij}}{\partial x_j} + S_{ik}\frac{\partial \ln \rho}{\partial x_k}}\right)} ,$ (13)

where $\underline{S}$ denotes the symmetric strain tensor

 \begin{displaymath}{S_{ij}} = { \frac{1}{2}\left({\frac{\partial u_i}{\partial x... ...ight) - \frac{1}{3}\delta_{ij} \vec{\nabla} \cdot \vec{u}} , \end{displaymath} (14)

with ${\underline{S}^2 = S_{ik}S_{ik}}$; ${\nu}$ is the kinematic viscosity, and ${\kappa}$ the heat conduction coefficient (divided by $c_{\rm V}$). The influence of rotation is given by the Coriolis force ${\vec{F}_{{\rm Cor}} = - 2~\rho~ \vec{\Omega} \times \vec{u}}$, which introduces a dependence on the latitude $\beta$, sometimes called "$\beta$ - effect'' in Geophysics:


 
$\displaystyle {\vec{F}_{{\rm Cor}} = 2~\rho~\Omega \left( \begin{array}{c} u_{y... ...~ \beta + u_{x}\sin~ \beta }\right) \\ u_{y} \cos~ \beta \end{array}\right) .}$     (15)

The angle $\beta$ is measured from the equatorial plane to the North pole with positive values, i.e. from ${0\hbox{$^\circ$ }}$ to ${90\hbox{$^\circ$ }}$. The gravity is assumed to have only a radial component, i.e. ${\vec{g} = (0, 0, g_r)}$.

The conservation equations are supplemented by the equation of state for a perfect gas,

\begin{displaymath}{P = \left({\gamma - 1}\right) \rho e ,} \end{displaymath} (16)

where ${\gamma = c_{\rm P}/c_{\rm V} = 5/3}$.

3.2 Cylindrical geometry

For cylindrical geometry both the mass and energy conservation Eqs. (11) and (12) remain unchanged if one uses vector notation. The three components of the Navier-Stokes equation read:

                            $\displaystyle {\frac{{\rm D} u_r}{{\rm D} t}}$ = $\displaystyle { \frac{u^2_{\beta}}{r} - \frac{1}{\rho}\frac{\partial P}{\partia... ...ho}\left(\vec{F}_{{\rm Cor}}\right)_r } +2\nu\left(S_{rk}\nabla_k\ln\rho\right)$  
       
    $\displaystyle {+ \nu \left[{\Delta u_r - \frac{u_r}{r^2} - \frac{2}{r^2}\frac{\... ...ac{\partial }{\partial r} \left({\vec{\nabla} \cdot \vec{u}}\right)}\right] } ,$ (17)


                            $\displaystyle {\frac{{\rm D} u_\beta}{{\rm D} t}}$ = $\displaystyle { - \frac{u_{\beta}u_r}{r} - \frac{1}{\rho r}\frac{\partial P}{\p... ...c{F}_{{\rm Cor}}\right)_{\beta} } +2\nu\left(S_{\beta k}\nabla_k \ln\rho\right)$  
    $\displaystyle {+ \nu \left[{\Delta u_{\beta} - \frac{u_\beta}{r^2} + \frac{2}{r... ...\partial }{\partial \beta} \left({\vec{\nabla} \cdot \vec{u}}\right)}\right]} ,$ (18)


                            $\displaystyle {\frac{{\rm D} u_{\phi}}{{\rm D} t}}$ = $\displaystyle {- \frac{1}{\rho }\frac{\partial P}{\partial \phi} + \frac{1}{\rh... ...\vec{F}_{{\rm Cor}}\right)_{\phi} } +2\nu\left(S_{\phi k}\nabla_k\ln\rho\right)$  
    $\displaystyle {+ \nu \left[{\Delta u_{\phi} + \frac{1}{3}\frac{\partial }{\partial \phi} \left({\vec{\nabla} \cdot \vec{u}}\right)}\right]} ,$ (19)

where
$\displaystyle {\underline{S} = \left( \begin{array}{ccc} S_{rr} & S_{r\beta} & ... ... \phi} \\ S _{\phi r} & S_{\phi \beta} & S_{\phi \phi} \\ \end{array}\right)}$     (20)

(note that ${\phi}$ has the dimension of length). The precise form of the cylindrical stress tensor for fully compressible flow can be found in Appendix A. The Coriolis force in this system reads
 
$\displaystyle {\vec{F}_{{\rm Cor}}} = {2~\rho~\Omega \left( \begin{array}{c} - ... ...\phi} \sin \beta \\ u_{r} \cos\beta - u_{\beta} \sin\beta\end{array}\right)} ,$     (21)

in the order of the unit vectors ${\vec{\hat{e}}_{r}, \vec{\hat{e}}_{\beta}, \vec{\hat{e}}_{\phi}}$.

3.3 Parameters

3.3.1 Scaling
Since only hydrodynamics is considered, basic mechanical units are sufficient for the scaling of the equations. We choose d as unit of length, and $g={\vert\vec{g}\vert = \vert g_r\vert}$ as unit of gravity; it follows that ${t_0= \sqrt{d/g}}$ and ${u_0 = \sqrt{d g}}$ are the units of time and velocity. In order to compare our results with solar values we use ${g = 370\;{\rm m~s^{-2}}}$ as a typical value for the gravitational acceleration in the middle of the solar convection zone, and ${d \approx 100~000\;{\rm km}}$. This yields ${t_0 \approx 10 \;{\rm min}}$ and ${u_0 \approx 200\; {\rm km~s^{-1}}}$, which is of the order of the sound speed at the bottom of the convection zone. The unit of the viscosity is given by $\nu_0 = u_0 d$, the unit of density is ${\rho_1 = \rho(r_1)}$, the density at the interface between the lower and middle layers.

  
3.3.2 Dimensionless parameters

A simulation is controlled by the following dimensionless parameters: The Rayleigh number

\begin{displaymath}{{\rm Ra} = \frac{L^4~g ~\delta} {\chi \nu H_{\rm P}} ,} \end{displaymath} (22)

the Prandtl number

\begin{displaymath}{{\rm Pr} = \frac{\nu}{\chi} ,} \end{displaymath} (23)

and the Taylor number

\begin{displaymath}{{\rm Ta} = \left(\frac{2~\Omega~ L^2}{\nu}\right)^2\cdot} \end{displaymath} (24)

The vertical extent L of the unstable layer, which is used in the definitions of Ra and Ta, is equal to the length scale d in all models except CB and CC. The parameter ${\delta = 1/(m_2+1) - (\gamma - 1)/\gamma = 0.1}$is the superadiabaticity of the middle (unstable) layer in the initial state, m2 is the polytropic index in the unstable layer, ${H_{\rm P}}$ is the pressure scale height in the middle of that same layer, and $\chi =\kappa/\rho_1$ is the thermal diffusivity. For our choice of parameters, we have ${H_{\rm P} = 0.45d}$initially, and this remains approximately so during the course of the calculation. The heat conductivity ${\kappa}$ is also from the middle layer. As far as the rotational influence is concerned, the inverse Rossby or Coriolis number,

 \begin{displaymath}{\rm Co} = \frac{2~\Omega~ d}{u_{{\rm rms}}}, \end{displaymath} (25)

is of more interest than the Taylor number. In our model Co is a derived quantity. For the box models the chosen parameters yield ${{\rm Co} \approx 4}$, which is representative for the lower part of the solar convection zone; for the cylinder models we explored a wider range, from ${{\rm Co} = 0 }$ to ${\rm Co} \approx 40$. From ${{\rm Ra, Pr, Ta}, d, g,\gamma }$ and the polytropic indices m1, m2, m3 of the three-layer models (see Sects. 2.1 and 2.2) the equation parameters ${\nu, \kappa}$ and $\Omega$are calculated. We fixed ${d \equiv 1, g_r \equiv -1, \gamma = 5/3, m_1 = 3, m_2 = 1, m_3 = 3}$ for all models. The condition of hydrostatic equilibrium yields the stratification of the initial state. For details of the initial stratification we refer to Brandenburg et al. (1996), Ossendrijver et al. (2001), or Käpylä et al. (2004), where the same three-layer models and parameter systems are used.

4 Numerics

The hydrodynamical equations were solved with a semi-discrete approach, using 6th-order-accurate finite-difference approximations for spatial derivatives, and either the Adams-Bashforth-Moulton 3rd-order predictor-corrector scheme described in Caunt & Korpi (2001), or the 2N-Runge-Kutta 3rd-order method mentioned in Brandenburg (2003) for the time evolution; the latter method allows larger time steps. In all simulations the time step was determined as the smaller of the diffusive and sound travel time over one spatial step, multiplied by a safety factor 0.2 ("Courant parameter''). For the box models we used the finite-difference code already described in Hupfer et al. (2005); for the cylinder models we used cylindrical coordinates and implemented the corresponding geometry into our code. The module containing the hydrodynamic source terms was completely rewritten for this new geometry.

The simulations where run on the 16-processor Linux Cluster KABUL of the Kiepenheuer-Institut and on private personal computers. For the box models the code was validated by determining critical Rayleigh numbers for the onset of convection, and comparison with results of Gough et al. (1976); moreover, comparison simulations were made by M. Ossendrijver using a code derived from one originally written by A. Brandenburg, and by W. Dobler using the Pencil code. No such tests were possible for the cylinder models: we know of no other such models (with cylinder axis perpendicular to the axis of rotation). A measure of the statistical error in box-model results is the behavior of  $Q_{\theta \varphi }$ at the pole ( $\beta = 90^\circ$); we verified that integrating for a longer time decreased this error. But again this is not true for the cylinder models where $Q_{\theta \varphi }$ generally has a finite value at the pole.

Table 1: Box models. Common parameters for all series are ${\rm Ra} = 250~000$, ${\rm Pr} = 0.4$, and $\nu = 5.963\times 10^{-4} \; \nu_0$. Due to the CFL condition and the larger grid the time step for BD is smaller than for BC, otherwise these models are identical. The Coriolis number ${{\rm Co}}$ is determined according to (25) after the data analysis has been done.

5 Results

For both the box and the cylinder models we performed simulations, varying the input parameters. Table 1 summarizes the box models. Each series consists of calculations on different latitudes either on the southern (BA, BB) or on the northern hemisphere. We were interested especially in the behavior of the Reynolds stresses near the equator, where we did additional simulations in small steps of  $\Delta\beta=2.5\hbox{$^\circ$ }$.

The parameters of the cylinder models are summarized in Table 2. Table 3 lists the variation of the Taylor number for the series CA. The values of  $u_{\rm rms}$as well as the values of  $Q_{\theta \varphi }$ listed in that table clearly demonstrate the increasing rotational effect at moderate Coriolis number, and the quenching at large Coriolis number. The models of the series CAS are identical to the CA models, except for different radii $r_{\rm b}$ and $r_{\rm t}$, as listed in Table 4. Tables 5 and 6 contain models where the thickness of the convectively unstable layer is 2d, for two values of the Rayleigh number.

5.1 Averaging method

As a number of previous studies have shown, the Reynolds stresses are extremely fluctuating quantities in space and time. The momentary change in time does not reveal anything but turbulence. Since differential rotation is a long term phenomenon, averaging both in space (horizontal coordinates) and time is done. We used up to 1500 individual snapshots to calculate time averages, equally spaced by one time unit for the cylinder models, and by up to four time units for the box models. Since the cylinder simulations use only two coordinates, no horizontal average is performed in order to keep the information on the latitudinal dependence - the results presented in the corresponding figures are time averages only. Sliding averages over a latitude range do not yield significantly different results and are not presented here.

Table 2: Cylinder models. The Prandtl number is always ${{\rm Pr} = 1}$, the polytropic indices are identical to those used in the box models BC and BD.The quantity L is the thickness of the convective layer.

Table 3: Cylinder model series CA, with ${\rm Ra} = 10^6$, ${\rm Pr} = 1.0$. The value for $Q_{\theta \varphi }$ is calculated from volume and time averaging over the unstable region, the latitude range 0-35 $\hbox {$^\circ $ }$, and 1500 time units.

Table 4: Models CAS - cylinder runs with different radius ratios. The models CAS03, CAS04, CAS05 are identical to CA12, CA14, CA16 and are included here for completeness. The angular velocities can be found in Table 3 for the respective Taylor numbers.

Table 5: Models CB - cylinder runs with a convective layer of thickness 2d and Rayleigh number $16\times 10^6$.

Table 6: Models CC - cylinder runs with a convective layer of thickness 2d and Rayleigh number $8\times 10^6$.

Table 7: Position of the maximum of $Q_{\theta \varphi }$ from box calculations, and comparison with other authors. All results are adapted to the northern hemisphere. The depth position is indicated by CL = cooling layer, UCL = upper convective layer, BCL = bottom of convective layer.


  \begin{figure} \par\includegraphics[angle=-90,width=8cm,clip]{5519fg03.eps} \end{figure} Figure 3: Reynolds stress component $Q_{\theta \varphi }$ for different northern latitudes as a function of depth from the box model BD, with ${{\rm Co} \approx 4}$. The vertical lines give the boundaries of the unstable layer. (This and the following figures are available in color in electronic form.)
Open with DEXTER


  \begin{figure} \par\includegraphics[angle=-90,width=8cm,clip]{5519fg04.eps} \end{figure} Figure 4: As Fig. 3, but for the cylinder model CA16, with ${{\rm Co}\approx 3}$. The largest peak appears at ${\beta = 12.5\hbox {$^\circ $ }}$.
Open with DEXTER

5.2 Results for Q ${_{\theta \varphi }}$

In Table 7 we compare the box models BB, BC, BD with results published in other papers. The common result is that there is a peak of  ${Q_{\theta\varphi}(r)}$ near the surface of the convective layer and at low latitude, near the equator. The latitude of this peak seems to be ${\beta_{{\rm max}} = 10\hbox{$^\circ$ }\pm 2.5\hbox{$^\circ$ }}$. In Fig. 3 we show the depth dependence of $Q_{\theta \varphi }$ for boxes placed at different latitudes from model BD. The latitudes from  $2.5\hbox{$^\circ$ }$ to  $20\hbox{$^\circ$ }$ dominate obviously over higher latitudes and the equator itself. The latitude and depth dependence from model BD is - apart from the different hemisphere - similar to that presented in Fig. 3 in Hupfer et al. (2005), which was produced with data from model BB. The main differences between BB and BD are slightly smaller amplitudes and the fact that all peaks are shifted into the convective layer for the latter model. These differences must be due to a slightly different stratification and a modified heat transport in the upper stable layer, which was not cooled in model BD.

In order to compare the results from box and cylinder simulations it is important to use models which have about the same Coriolis number. We have done box simulations up to ${{\rm Co} \approx 4}$. These are compared with results from the moderate rotation regime, ${{\rm Co}\approx 1.8 {\rm\ to\ Co} \approx 4.0}$, of our cylinder models (models CA12 to CA18 from Table 3). A typical behavior of ${Q_{\theta\varphi}(r)}$ for different latitudes on the northern hemisphere (first quadrant) is shown in Fig. 4 for model CA16, with a maximum peak at ${\beta = 12.5\hbox {$^\circ $ }}$.

The amplitude of ${Q_{\theta\varphi}(r)}$ obtained from cylinder models is always greater than for results from box models. The most probable reason is that the cylinder model has a two-dimensional geometry - there is no third dimension to calculate horizontal averages. Averages over a range of latitudes reduce the amplitudes only slightly. Comparing with Fig. 11 in Käpylä et al. (2004), obtained from box calculations with ${{\rm Co} \approx 10}$, we found no such peak behavior for low latitudes in our cylinder model CA22, which has a similar Coriolis number. In the range ${{\rm Co} = 4}$ to 10 the dominant low-latitude peak shifts both to higher latitude, around ${\beta \approx 40\hbox{$^\circ$ }}$, and more into the middle of the unstable layer. When ${{\rm Co}}$ exceeds 20 a sign reversal occurs at all latitudes, and the peaks are now minima, located at the top of the convective layer and at ${\beta = 45\hbox {$^\circ $ }}$ (Fig. 5). However, the absolute values are smaller by a factor of about five than in model CA16, which indicates quenching of the ${\Lambda }$ effect, as described by Kichatinov & Rüdiger (1993). A third dimension used in horizontal averages could also reduce the root mean square velocity  $u_{\rm rms}$, shifting the Coriolis numbers from the cylinder models closer to the values reported by Käpylä et al. (2004).

Apart from smaller amplitudes there is no significant change for $Q_{\theta \varphi }$ if the ratio of the radii is altered. A result for ${{\rm Co} = 2.9}$ and ${\eta = \frac{2}{4}}$ is shown in Fig. 6. There is also no large deviation if the thickness of the unstable layer is doubled as has been done in the series CB and CC; the main change consists of broader peaks near the transition region between the convective and the upper stable layer, see Fig. 7. Since both the Rayleigh and the Taylor numbers are proportional to d4 it is necessary to adapt both numbers in order to get the same mechanical and thermal viscosities as in the CA cases, otherwise the results would be hard to compare to each other. A common result for all cylinder models is a non-vanishing Reynolds stress component $Q_{\theta \varphi }$ at the poles, ${\beta = \pm 90\hbox{$^\circ$ }}$. This is because in the cylinder models neither  ${u_{\phi}}$ nor  ${u_{\beta}}$ necessarily vanishes at the poles.


  \begin{figure} \par\includegraphics[angle=-90,width=8cm,clip]{5519fg05.eps} \end{figure} Figure 5: As Fig. 3, but for the cylinder model CA24 for fast rotation with ${{\rm Co} = 24.1}$. The minimum peak is at ${\beta = 45\hbox {$^\circ $ }}$. Note the different scaling as compared to Fig. 4.
Open with DEXTER


  \begin{figure} \par\includegraphics[angle=-90,width=8cm,clip]{5519fg06.eps} \end{figure} Figure 6: As Fig. 3, but for the cylinder model CAS02, with ${{\rm Co} = 2.9}$. Note the different radii of the lower and upper boundary, as compared to Fig. 4. The largest peak appears at ${\beta = 10\hbox {$^\circ $ }}$.
Open with DEXTER


  \begin{figure} \par\includegraphics[angle=-90,width=8cm,clip]{5519fg07.eps} \end{figure} Figure 7: As Fig. 3, but for the cylinder model CC02, with ${{\rm Co} = 2.96}$. The thickness of the convective layer has been doubled compared to the CA series. See Fig. 4 for differences of this model with CA16 which has virtually the same Coriolis number. The largest peak appears at ${\beta = 12.5\hbox {$^\circ $ }}$.
Open with DEXTER

  
Comparison with solar observations

Ward (1965) determined the correlation $Q_{\theta \varphi }$ from the motion of sunspot groups on the northern hemisphere, with a maximum at ${\beta \approx 25\hbox{$^\circ$ }}$ and ${Q_{\theta\varphi} \approx 2000\;{\rm m^2~s^{-2}}}$. Gilman & Howard (1984) found similar values from the motion of individual sunspots, with a maximum at ${\beta = 30\hbox{$^\circ$ }}$. In a more recent paper Vrsnak et al. (2003) published results derived from the motion of coronal bright points, whereby the latitudinal range was extended up to ${\beta \approx 60\hbox{$^\circ$ }}$. They confirmed that ${Q_{\theta \varphi} > 0}$ on the northern hemisphere and the order of magnitude of this correlation, however with large standard deviations (see Fig. 6a and Table 3 in Vrsnak et al. (2003); note that this component of the Reynolds stress changes sign when latitude instead of co-latitude is used).

Since ${u_0 \approx 200\; {\rm km~s^{-1}}}$ the Reynolds stress components determined by the present numerical models are too large by at least three orders of magnitude if the observed solar surface values are used for comparison. A similar discrepancy occurs if we evaluate the turbulent viscosity according to $\nu_{\rm t} = u_{{\rm rms}}d/3$, where  $u_{\rm rms}$ is obtained from a volume and time average over the unstable layer; values are listed in Tables 3 to 6. These discrepancies are due to the unrealistic parameter regime which is enforced by the constraints of solving the Navier-Stokes equations numerically (essentially the too large input energy flux). Apart from the Coriolis number all other dimensionless parameters given in Sect. 3.3.2 are far away from solar values; e.g., Ossendrijver (2003) lists ${\rm Ra} = 10^{16}\dots 10^{20}$, ${\rm Pr} = 10^{-7}$, and ${\rm Ta} = 10^{14}\dots 10^{27}$ for the Sun. It is therefore only possible to compare the results qualitatively. Both observation and numerical experiment with various codes and geometric setups show that the Reynolds stress $Q_{\theta \varphi }$ is on the whole positive on the northern hemisphere and has a significant peak in the latitude band from  ${10\hbox{$^\circ$ }}$ to  $30\hbox{$^\circ$ }$, with analogous results for the southern hemisphere. However, the numerical results suggest the peak at the transition region between the unstable and the upper stable layer; this is not to be identified with the surface of the Sun if the calculations are considered ro represent the lower half of the convection zone, which is the interesting case because at that depth the rotational influence is significant and the overshooting below the unstable layer is included. Unfortunately, there is no observational information about the dependence of  $Q_{\theta \varphi }$ on the depth in the solar convection zone.

The indications for the angular momentum transport are obvious: a positive correlation (on the northern hemisphere) means an equatorward transport, from higher to lower latitudes. The occurrence of the strong peaks in the equator zone should lead also to an enhanced transport of angular momentum and to an increase of the angular velocity, leaving latitude bands in which the rotation is faster than average, and as well bands in which it is smaller than average. Indeed, Vrsnak et al. (2003) presented results for such bands, especially the stronger rotation region between  ${10\hbox{$^\circ$ }}$ and  ${25\hbox{$^\circ$ }}$, cf. Fig. 1b in Vrsnak et al. (2003), which remarkably well coincides with the results from numerical simulations in this paper or in Pulkkinen et al. (1993), Chan (2001), and Käpylä et al. (2004).


  \begin{figure} \par\includegraphics[angle=-90,width=8cm,clip]{5519fg08.eps} \end{figure} Figure 8: Reynolds stress components ${Q_{\theta \varphi }, Q_{r\theta }, Q_{r\varphi }}$, averaged over the range 0-35 $\hbox {$^\circ $ }$ of latitude, for the cylinder models CA as functions of the Coriolis number.
Open with DEXTER

Table 8: Maximum position for $Q_{\theta \varphi }$, as function of Coriolis number ${{\rm Co}}$ for moderate rotation (models CA and CAS). The subscripts 1 and 2 indicate maxima near the surface and near the bottom of the convective layer, respectively.

5.3 Averages over latitude

In order to show more clearly the qualitative behavior of the Reynolds stresses derived from the cylinder model, we display in Fig. 8 averages over the latitude range 0-35 $\hbox {$^\circ $ }$. This range is probably the most relevant to the Sun; and the cylindrical geometry does not deviate too much from the spherical geometry.

Except for very slow rotation, the horizontal transport is positive in the range of Coriolis numbers relevant to the Sun ( ${{\rm Co} \leq 20}$). As far as the sign is concerned, this is consistent with the result of Kichatinov & Rüdiger (1993) who used a given anisotropic turbulence model; the magnitude of our result is, however, much larger, for the reasons given above (Sect. 5.2).

The radial transport is negative, directed to greater depths, in about the same range of Coriolis numbers. This is consistent with other recent numerical simulations, e.g. Rüdiger et al. (2005), but in contrast to the earlier analytical result of Kichatinov & Rüdiger (1993). The meridional correlation  $Q_{r\theta}$ is as well negative, as indicated also by the orientation of the circulation pattern.

At very large Coriolis number all three correlations decrease in magnitude, a result that is consistent with Kichatinov & Rüdiger (1993): very fast rotation quenches the ${\Lambda }$ effect.

5.4 Meridional flow


  \begin{figure} \par\includegraphics[width=8cm,clip]{5519fg09.eps} \end{figure} Figure 9: Contour plot of the mean meridional velocity ${\langle u_{\theta}\rangle }$ from the cylinder model CA16. The labels indicate the angles in degrees, counted from the eastern equator. The flow of the dominant cell is directed polewards in the upper, and equatorwards in the lower part of the model.
Open with DEXTER

Information on meridional flow and a possible circulation can be obtained from a global model like the cylinder approach. We found meridional circulation already for the case of no rotation, in the model CA00. In this model the circulation pattern is partitioned into multiple small cells along the meridional direction, having an extent of about  ${15\hbox{$^\circ$ }}$. For small Coriolis numbers there is only one cell in the radial direction, concentrated in the convective layer and - due to strong overshooting - in the lower stable layer. As long as ${{\rm Co} < 0.6}$ the cell near the equator has always a counter-clockwise orientation, i.e. a poleward oriented flow in the upper part of the model. In the range ${{\rm Co} = 0.6{-}1.2}$ there is clockwise circulation near the equator, combined with merging of cells, i.e. the number of cells is reduced. For ${{\rm Co} > 1.2}$ the cells near the poles vanish and two cells in radial direction appear; the outer cell is confined to the upper stable layer and is very thin, its orientation is clockwise, and the flow of the inner cell is contained almost completely within the convective layer, the circulation being counter-clockwise, see Fig. 9. The cells have a latitudinal extent up to ${45\hbox{$^\circ$ }}$; this angle decreases with increasing Coriolis number to about $30\hbox{$^\circ$ }$. In the strong-rotation regime ( ${{\rm Co} > 4 }$) the meridional cell pattern vanishes and we find Taylor columns (Fig. 10) as a consequence of the Taylor-Proudman theorem, cf. Chandrasekhar (1961).

The box models cover only a narrow range of latitude; nevertheless we may use these models to extract some information on meridional flow, simply by calculation of $\langle u_x \rangle$ for each box. In Figs. 11 and 12 the meridional velocity is shown for the northern hemisphere as obtained from the box models BD and from CA16, a cylinder model with comparable Coriolis number. Positive values indicate an equatorward flow. In Fig. 12 the poleward orientation of the flow in the bulk of the convective layer is clearly discernible, as well as the equatorward orientation at the bottom of the unstable layer and at the top of the computational domain. For the box model a partitioning in (at least) two cells of circulation occurs: the low-latitude cell has a clockwise orientation, whereas the cell at higher latitudes has a counter-clockwise orientation and a smaller magnitude - the lower cell dominates. This is one of the main discrepancies between plane-layer and cylinder calculations. The difference in the form the circulation is probably due to the periodic boundary conditions applied in the meridional direction of the box models.

  \begin{figure} \par\includegraphics[width=8cm,clip]{5519fg10.eps} \end{figure} Figure 10: Contour plot of the mean meridional velocity ${\langle u_{\theta}\rangle }$ from the cylinder model CA26. The flow is aligned in layers of Taylor columns.
Open with DEXTER

The circulation pattern vanishes completely if the Coriolis number exceeds a value of $\approx$4 which also seems to be the limit for the existence of the peak in the Reynolds stress $Q_{\theta \varphi }$. This could be a particular feature of the cylinder model; unfortunately we do not know whether in box calculations the circulation also vanishes for Coriolis numbers beyond 10. From box models with the same stratification specifications as have been applied in this paper, Käpylä et al. (2004) in their Fig. 10 for ${{\rm Co} \approx 10}$ show the mean meridional velocity ux, which is very weak in the convective layer. If there is any circulation, it is rather suppressed in the rapid-rotation regime - a result which fits our findings from the cylinder simulations.


  \begin{figure} \par\includegraphics[angle=-90,width=8cm,clip]{5519fg11.eps} \end{figure} Figure 11: Average over horizontal coordinates and time of the meridional velocity ux at different northern latitudes from the box model BD.
Open with DEXTER


  \begin{figure} \par\includegraphics[angle=-90,width=8cm,clip]{5519fg12.eps} \end{figure} Figure 12: Time average of the meridional velocity ${u_{\theta }}$ from model CA16 for different northern latitudes.
Open with DEXTER

6 Conclusion and summary

We have performed numerical simulations of convection zones of solar-type stars with two different geometries in order to explore the influence of rotation and the resulting Reynolds stresses. We found that the horizontal transport coefficient $Q_{\theta \varphi }$ is positive on the northern hemisphere at almost all depths and latitudes, and a range of Coriolis numbers that is relevant to the solar convection zone. This agrees with results determined at the solar surface by various observers, e.g., Ward (1965), Gilman & Howard (1984), Vrsnak et al. (2003). However, in the regime of very rapid rotation the sign of $Q_{\theta \varphi }$ is reversed, and the horizontal transport is polewards. Comparing both geometries we confirmed the existence of a strong peak of ${Q_{\theta\varphi}(r)}$ at a low latitude and near the surface of the convection layer, but this peak disappears in the cylinder models when the Coriolis number exceeds a value of about 4.

Regarding the meridional flow we found that there is a circulation in the weak and moderate rotation regime, ${{\rm Co} = 0 \dots 4}$; for more rapid rotation one obtains Taylor columns in the meridional plane. The circulation is partitioned into multiple latitudinal cells for weak rotation and forms one meridional cell for moderate rotation, with counter-clockwise orientation in the convective layer.

Although the two-dimensional cylinder model has (geometric) problems in the polar regions it proves to be an efficient way to get information on connected flows with a high latitudinal resolution.

Acknowledgements
We gratefully acknowledge discussions with M. Ossendrijver and W. Dobler, and suggestions for improving the text by Reiner Hammer and an unknown referee. CH was supported by the Deutsche Forschungsgemeinschaft, and PJK acknowledges support from the Finnish graduate school for astronomy and space physics.

   
Appendix A: The hydrodynamical stress tensor in cylindrical coordinates

The stress tensor is symmetric; the six independent components  ${S_{rr}, S_{r\beta}, S_{r\phi}, S_{\beta\beta}, S_{\beta\phi}, S_{\phi\phi}}$ are given by:

                     Srr = $\displaystyle {\frac{\partial u_r}{\partial r} - \frac{1}{3} \vec{\nabla} \cdot \vec{u}} ,$ (A.1)
$\displaystyle {S_{r\beta}}$ = $\displaystyle {\frac{1}{2} \left[{\frac{1}{r}\frac{\partial u_r}{\partial \beta}+ r \frac{\partial}{\partial r}\left({\frac{u_{\beta}}{r}}\right)}\right]} ,$ (A.2)
$\displaystyle {S_{r\phi}}$ = $\displaystyle {\frac{1}{2} \left({\frac{\partial u_r}{\partial \phi} + \frac{\partial u_\phi}{\partial r}}\right)} ,$ (A.3)
$\displaystyle {S_{\beta\beta}}$ = $\displaystyle { \frac{1}{r}\frac{\partial u_{\beta}}{\partial \beta} + \frac{u_r}{r} - \frac{1}{3} \vec{\nabla} \cdot \vec{u}} ,$ (A.4)
$\displaystyle {S_{\beta\phi}}$ = $\displaystyle {\frac{1}{2} \left({\frac{1}{r}\frac{\partial u_\phi}{\partial \beta} + \frac{\partial u_\beta}{\partial \phi}}\right)} ,$ (A.5)
$\displaystyle {S_{\phi\phi}}$ = $\displaystyle { \frac{\partial u_\phi}{\partial \phi} - \frac{1}{3} \vec{\nabla} \cdot \vec{u}} .$ (A.6)

References

 
Copyright ESO 2006