A&A 488, 9-23 (2008)
DOI: 10.1051/0004-6361:20079098

## Lambda effect from forced turbulence simulations

P. J. Käpylä1,2 - A. Brandenburg2

1 - Observatory, PO Box 14, 00014 University of Helsinki, Finland
2 - NORDITA, AlbaNova University Center, Roslagstullsbacken 23, 10691 Stockholm, Sweden

Received 19 November 2007 / Accepted 22 May 2008

Abstract
Aims. We determine the components of the -effect tensor that quantifies the contributions to the turbulent momentum transport even for uniform rotation.
Methods. Three-dimensional numerical simulations are used to study turbulent transport in triply periodic cubes under the influence of rotation and anisotropic forcing. Comparison is made with analytical results obtained via the so-called minimal tau-approximation.
Results. In the case where the turbulence intensity in the vertical direction dominates, the vertical stress is always negative. This situation is expected to occur in stellar convection zones. The horizontal component of the stress is weaker and exhibits a maximum at latitude   - regardless of how rapid the rotation is. The minimal tau-approximation captures many of the qualitative features of the numerical results, provided the relaxation time tau is close to the turnover time, i.e. the Strouhal number is of order unity.

Key words: hydrodynamics - turbulence - Sun: rotation - stars: rotation

### 1 Introduction

Differential rotation plays a crucial role in dynamo processes that sustain large-scale magnetic activity in stars like the Sun (e.g. Moffatt 1978; Krause & Rädler 1980). The internal rotation of the Sun is familiar from helioseismology (e.g. Thompson et al. 2003), but the processes sustaining the observed rotation profile are not understood well. The angular momentum balance of a rotating star is determined by the conservation equation

 (1)

where is the meridional flow, s the cylindrical radius, the density (neglecting however its fluctuations), the local angular velocity, and the zonal component of the Reynolds stress. Overbars denote averages over the azimuthal direction.

The meridional flow can also be directly affected by the Reynolds stresses (e.g. Rüdiger 1989), but it is more strongly determined by the baroclinic term that arises if the isocontours of density and pressure do not coincide. Such a configuration can appear because of latitude-dependent turbulent heat fluxes that arise naturally in rotating convection (e.g. Rüdiger 1982; Pulkkinen et al. 1993; Käpylä et al. 2004; Rüdiger et al. 2005a) or from a subadiabatic tachocline (Rempel 2005) which is likely to occur in the Sun (Rempel 2004; Käpylä et al. 2006b). The flows due to thermodynamic effects are likely to be needed to avoid the Taylor-Proudman balance in the Sun (e.g. Durney 1989; Brandenburg et al. 1992; Kitchatinov & Rüdiger 1995; Rempel 2005). The overall importance of the meridional flow in the angular momentum balance of the Sun is, however, still unclear since no definite observational information about it is available below roughly 20 Mm depth (e.g. Zhao & Kosovichev 2004).

Although not much more is known about the Reynolds stresses from observations, already this limited knowledge can be used to gain insight into the theory of turbulent transport. Solar surface observations indicate that there is an equatorward flux of angular momentum, as described by the Reynolds stress component , of several 103 m2 s-2 in the latitude range where sunspots are observable (e.g. Ward 1965; Pulkkinen & Tuominen 1998; Stix 2002). In mean-field theory the simplest approximation that can be made concerning the Reynolds stresses is to assume them proportional to the gradient of mean velocity (the Boussinesq ansatz):

 (2)

In the Sun this ansatz turns out to be insufficient because the observed and solar surface differential rotation profile indicate that the turbulent viscosity is negative. Thus, in analogy to mean-field dynamo theory, additional contributions to the Reynolds stress were conjectured to appear (e.g. Wasiuty ski 1946; Krause & Rüdiger 1974), leading to the present description

 (3)

where is a third-rank tensor describing the -effect that contributes to the Reynolds stress even in the case of rigid rotation. These terms are often referred to as non-diffusive'' contributions to the Reynolds stress. The zonal components of the stress can be written in the form (e.g. Stix 2002)

 (4) (5)

where and describe the non-diffusive transport and where is the turbulent viscosity. The factor  in the latter equation indicates that the turbulent viscosity can be anisotropic. Furthermore, the coefficients  , , and can vary as functions of the spatial coordinates.

Much effort has been devoted to determining Reynolds stresses from convection simulations (Hathaway & Somerville 1983; Pulkkinen et al. 1993; Rieutord et al. 1994; Brummell et al. 1998; Chan 2001; Käpylä et al. 2004; Rüdiger et al. 2005b; Hupfer et al. 2005, 2006; Giesecke 2007). These studies have confirmed the existence of the -effect and revealed some surprising results that are at odds with theoretical considerations (e.g. Kitchatinov & Rüdiger 1993, 2005) derived under the second-order correlation approximation (SOCA). The discrepancies are most prominent in the rapid rotation regime, , where

 (6)

is the Coriolis (or the inverse Rossby) number. Here, is the angular momentum-averaged rotation rate and the convective turnover time. In the solar convection zone, varies between 10-3 near the surface to ten or more in the deep layers. Convection simulations in the latter regime show that the horizontal angular momentum flux is directed toward the equator, corresponding to positive in the northern hemisphere, and that it peaks very sharply near the equator (Chan 2001; Käpylä et al. 2004; Hupfer et al. 2005). On the other hand, the vertical stress can be directed outward (Käpylä et al. 2004), contradicting the theory for vertically dominated turbulence (e.g. Biermann 1951; Rüdiger 1980, 1989). So far, these results remain without proper explanation.

Often the Reynolds stress realized in the simulation is taken to solely represent the -effect. This approach seems like a reasonable starting point but in an inhomogeneous system large scale mean flows are generated when rotation becomes important. These flows affect the Reynolds stresses via the turbulent viscosity. Furthermore, in the presence of stratification, heat fluxes can also significantly affect the stresses (Kleeorin & Rogachevskii 2006). In the present study we simplify the situation as much as possible in order to disentangle the effect of the turbulent velocity field from other effects. Thus we neglect stratification and heat fluxes by adopting a periodic isothermal setup. Turbulence is driven by external forcing, which provides clear scale separation between the turbulent eddies and the system size, which is typically not achieved in convection simulations. Further insight is sought from comparison of simple analytical closure models with numerical data.

Preliminary results on the -effect are presented in Käpylä & Brandenburg (2007). In the present paper, numerical datasets covering a significantly larger part of the parameter space are analyzed, and a more thorough study of the validity and results of the minimal tau-approximation are presented.

The remainder of the paper is organized as follows. Section 2 summarizes the model and the methods of the study, and in Sects. 3 and 4 the results and the conclusions are given.

### 2 The model and methods

#### 2.1 Basic equations

We model compressible hydrodynamic turbulence in a triply periodic cube of size . The gas obeys an isothermal equation of state characterized by a constant speed of sound, . The continuity and Navier-Stokes equations read

 (7)

 (8)

where denotes the advective derivative, is the velocity, the density, the viscous force, and the forcing function. Due to the periodic boundaries, mass is conserved and the average density retains its initial value at all times. Compressibility is retained but we consider low Mach number flows, . The viscous force is given by

 (9)

where is the kinematic viscosity and

 (10)

the traceless rate of strain tensor. The forcing function is given by
 (11)

where is the position vector, is a normalization factor, and describes the direction-dependent amplitude of the forcing (see below), , is the length of the time step, and a random delta-correlated phase. The vector is given by
 (12)

where is an arbitrary unit vector. Thus, describes nonhelical transversal waves with , where k is chosen randomly from a predefined range in the vicinity of the average wavenumber at each time step, where k1 is the wavenumber corresponding to the domain size, and the wavenumber of the energy-carrying scale.

To make the resulting turbulence anisotropic the forcing amplitude depends on direction

 (13)

where f0 and f1 are the amplitudes of the isotropic and anisotropic parts of the forcing, the unit vector in the vertical direction, and the angle between the vertical direction and the wave vector k. In the following we assume , but even then the amount of anisotropy depends on the Reynolds number and is in any case only quite modest (see, e.g. Table B.1).

At this point it is important to emphasize that our forcing function is designed to capture the effects that lead to a finite -effect. The implementation of anisotropy should therefore be a simple one. In stars, anisotropy is produced by stratification and convection. Our goal is clearly not to simulate properties of convection other than its tendency to make the turbulence anisotropic.

The numerical computations were made with the P ENCIL-C ODE, which uses sixth-order accurate finite differences in space, and a third-order accurate time-stepping scheme (Brandenburg & Dobler 2002; Brandenburg 2003). Resolutions up to 2563 grid points were used in the simulations.

#### 2.2 Nondimensional units

In the following we use non-dimensional variables by setting

 (14)

This means that the units of length, time, and density are
 (15)

However, in most of the plots we present the results in explicitly non-dimensional form using the quantities above.

#### 2.3 Coordinate system, averaging, and error estimates

The simulated domain is thought to represent a small rectangular portion of a spherical body of gas. We choose (x,y,z) to correspond to of spherical coordinates. With this choice the rotation vector can be written as

 (16)

where is the angle between the rotation axis and the local vertical direction, i.e. the colatitude.

Since the turbulence is homogeneous, volume averages are employed and denoted by overbars. An additional time average over the statistically saturated state of the calculation is also taken. We define the Coriolis number as

 (17)

where is the rms-velocity. In comparison to the commonly used definition (6), definition (17) is smaller by a factor of , i.e. . The other dimensionless number relevant in this study is the Reynolds number based on the forcing scale
 (18)

See Fig. 1 for a snapshot from a typical run with .

Errors are estimated by dividing the time series into three equally long parts and computing mean values for each part individually. The largest departure from the mean value computed for the whole time series is taken to represent the error.

 Figure 1: Ux at the periphery of the simulation domain from a slowly rotating run with , , and , resolution 2563. Open with DEXTER

#### 2.4 The -effect from the minimal tau-approximation

To have some understanding of the numerical results, we compare with the simplistic tau-approximation (hereafter MTA) in real space (e.g. Blackman & Field 2002; Brandenburg et al. 2004). Unlike the usual first-order smoothing approximation where nonlinearities in the fluctuations are neglected, they are retained in an approximate manner in MTA.

In order to develop a theory for the Reynolds stress, , we derive an equation for its time derivative,

 (19)

In the absence of large-scale flows, i.e. , the equation for the fluctuating part can be written as

 (20)

where is a nonlinear term. Multiplying Eq. (20) by uj gives

 (21)

Inserting this into Eq. (19) yields

 (22)

where are the triple correlations. Under the assumption of periodic boundary conditions, this term can be written in the form

 (23)

In MTA the higher than second-order terms, i.e. Tij, in the equations of turbulent correlations are retained in a collective manner by parametrizing them with a term that is equal to the original correlation divided by a relaxation time, i.e.

 (24)

Now the equations can be solved for the stresses in terms of the forcing. A simple ansatz for parameterizing the forcing is given in terms of the non-rotating equilibrium solution,

 (25)

Inserting the parameterizations (24) and (25) into Eq. (22) yields

 (26)

see Appendix A for more details on the equations used in the MTA-model. The model equations are very similar to those of Ogilvie (2003; see also Garaud & Ogilvie 2005) who studied hydrodynamic and magnetohydrodynamic turbulence and angular momentum transport due to shear flows and the magnetorotational instability. Among other things, they also introduced an isotropization term that causes decaying turbulence to become isotropic. Some of our models give explicit evidence of such a term.

The rotational influence in the MTA-model is measured by the Coriolis number

 (27)

where

 (28)

is the Strouhal number, which is the main free parameter in the model.

On account of previous results on similar systems of forced turbulence (e.g. Brandenburg & Subramanian 2005) is used as our reference model. To reproduce the numerical results, we introduce an empirical isotropization term in the MTA-model; see Sect. 3.1 and in particular Eqs. (29) and (30). The use of this term is regulated by another free parameter , which can obtain the values zero or unity.

The MTA-model results were obtained by advancing the time-dependent Eqs. (22) with the parameterizations (24) and (25) until a stationary solution was reached.

### 3 Results

#### 3.1 Diagonal components of the stress

Consider turbulence where the intensity in the vertical direction is stronger than the intensities in the horizontal directions. This situation is encountered if the turbulence is caused by the convective instability. First we consider the case with the maximum anisotropy that was achieved with the present model. Table B.1 summarizes the results for seven values of the Coriolis numbers ranging from 0.06 to roughly 5.4. Calculations at seven latitudes from the pole ( ) to the equator ( ) with equidistant intervals of were made with each Coriolis number, in all runs. In this set of runs, f0 = 10-6 and f1 = 0.2 were used in Eq. (13).

 Figure 2: Volume-averaged Reynolds stress components Qxx ( top), Qyy ( middle), and Qzz ( bottom), normalized by the square of the rms-velocity, as functions of latitude and rotation from the turbulence simulations listed in Table B.1. Coriolis number, as defined in Eq. (17), varies as indicated in the legend in the middle panel. Open with DEXTER

 Figure 3: Diagonal Reynolds stress components from the MTA-model. Here the Coriolis number is defined via Eq. (27) with and . Open with DEXTER

 Figure 4: Same as Fig. 3 but with . Compare with Figs. 2 and 3. Open with DEXTER

Figure 2 shows the diagonal components of the Reynolds tensor as functions of latitude and Coriolis number from the numerical turbulence simulations. As rotation is increased, the magnitudes of the horizontal components Qxx and Qyy increase monotonically while Qzz decreases, which illustrates an isotropizing effect of rotation on the turbulence. In the following, we refer to this effect as the rotational isotropization of turbulence. This effect is a purely empirical and refers to the observation that increased rotation leads to stronger mixing which washes out anisotropies that are caused by other effects. Of course, rotation itself can cause the turbulence to become anisotropic, but this seems to play a role only at much higher rotation rates. This is indeed seen in the most rapidly rotating case where the behavior is more complex. The MTA-model, Fig. 3, on the other hand, fails to reproduce isotropization of turbulence when rotation is included. The behavior is most obvious at the pole where rotation contributes no net effect to linear order, see Appendix A, Eqs. (A.1) to (A.6). The same is true for the Qxx component at the equator.

Some persistent trends arise as a function of latitude in the numerical simulations; for instance, Qxx peaks at the pole and decreases toward the equator except again for the fastest rotation. An approximately opposite trend is seen for Qyy for slow and rapid rotation, whereas in the intermediate range a minimum appears at mid-latitudes. For slow rotation, Qzz behaves approximately in the opposite way to Qyy, having a maximum at the pole and a minimum at the equator, whereas for rapid rotation the trend is reversed. There is also a persistent maximum at mid-latitudes. The MTA-model, however, fails to reproduce most of the characteristics of the latitude distribution of the diagonal stresses. This indicates that nonlinear terms, which are manifestly not described adequately well by the MTA-assumption, are the deciding factor in determining the behavior of the diagonal stresses.

To capture the rotational isotropization with the MTA-model at least qualitatively, we experimented by adding a term

 (29)

on the rhs of Eq. (26). Here Q is the trace of Qij and

 (30)

The functional form of is chosen purely empirically so that the magnitudes of the off-diagonal stresses in comparison to the simulations are fairly accurately reproduced with the MTA-model (see Fig. 7).

Figure 4 shows the results for the diagonal components with . Now the magnitudes of the turbulence intensities are more in line with the full numerical simulations, although the latitude distribution is still manifestly wrong. Although the off-diagonal components are much better represented by the linear terms appearing in the equation of the Reynolds stress (see Sect. 3.3), the rotational isotropization term helps for reducing their magnitudes closer to the levels seen in the direct simulations also in that case.

The functional form of in Eq. (29) indicates that for rapid rotation. This behaviour cannot be justified based on the present numerical data (see Fig. 2).

#### 3.2 Anisotropy of the turbulence

Turbulence anisotropies can be characterized by the quantities

 (31) (32)

The importance of these quantities is that, for slow rotation, they can be considered as proxies of the -effect according to and (e.g. Rüdiger 1980, 1989), where is the correlation time of the turbulence.

Table B.1 shows that for slow rotation increases monotonically from the pole to the equator. For , exhibits a negative minimum at mid-latitudes and reaches a positive maximum at the equator. The maximum at the equator can be explained as resulting from the surviving (y-)component of the Coriolis force at low latitudes. The relation between the horizontal -effect and the anisotropy parameter  is at best poorly confirmed by the numerical results. It must also be noted that is smaller by at least an order of magnitude in comparison to , which also enters the equation for Qxy but with a higher order in (see Eq. (A.7)).

The vertical anisotropy, , on the other hand, retains its sign in all models. For slow rotation the absolute value of decreases monotonically from the pole toward the equator. The latitude distribution for rapid rotation is approximately opposite, although an additional minimum ( ) or a maximum ( ) can occur at mid-latitudes. Here the correspondence between and holds at least for the sign for all calculations, bar the few cases in the intermediate rotation regime where Qyz is positive (see below).

#### 3.3 Off-diagonal components of the stress

In stellar convection zones, the off-diagonal Reynolds stresses contribute to the angular momentum transport and work to generate (-effect) or to smooth out (turbulent viscosity) differential rotation. In the present case only the former effect is in operation. Figure 5 summarizes the results for the runs listed in Table B.1.

The Qxy component of the stress corresponds to in spherical coordinates and is responsible for latitudinal angular momentum transport. In the simulations this component is always positive and the latitude distribution peaks at latitude  (see the uppermost panel of Fig. 5). The sign is in accordance with solar observations (Ward 1965; Pulkkinen & Tuominen 1998) and analytical turbulence models (Kitchatinov & Rüdiger 1993, 2005). The latitude distribution in the rapid rotation regime is significantly different from convection simulations where Qxy is sharply concentrated near the equator (Chan 2001; Käpylä et al. 2004; Hupfer et al. 2005; Rüdiger et al. 2005a). The reason for this difference is still unclear but is evidently related to the physics that have been omitted in the present study. The MTA-model captures the latitude dependence rather well, but the magnitude of the stress is clearly too great (see Fig. 6). If the rotational isotropization term, Eq. (29), is taken into account (see the uppermost panel of Fig. 7), the agreement is also better for the magnitude. Keeping the forcing and rotation rate fixed, the best agreement with the 3D models is found if is used, see Fig. 8.

 Figure 5: Same as Fig. 2 but for the off-diagonal stress components Qxy ( top), Qxz ( middle), and Qyz ( bottom). Linestyles as in Fig. 2. Open with DEXTER

 Figure 6: Same as Fig. 3 but for the off-diagonal stress components Qxy ( top), Qxz ( middle), and Qyz ( bottom). , . Linestyles as in Fig. 3. Open with DEXTER

 Figure 7: Same as Fig. 6 but with . Linestyles as in Fig. 3. Open with DEXTER

 Figure 8: Same as Fig. 6 but with , . Linestyles as in Fig. 3. Open with DEXTER

While the analogue of the Qxz component does not play a direct role in the angular momentum balance in stars, it can still contribute via generating meridional flows (e.g. Rüdiger 1989). In the numerical simulations we find that, for all cases except the most rapidly rotating one, Qxz is negative and peaks at (middle panel of Fig. 5). For , however, the sign changes near the equator, with positive values toward the equator and negative ones toward the pole. This is at odds with the analytical result of Rüdiger et al. (2005b), but does agree with the results of convection simulations (Pulkkinen et al. 1993; Käpylä et al. 2004).

The MTA-model gives qualitatively similar results, although the sign change occurs at significantly slower rotation; see the middle panel of Fig. 6. Rotational isotropization helps to correct the magnitude, but not the earlier occurrence of the sign change (Figs. 7 and 8).

Since the turbulent intensity of the vertical (z-)motions is greater than the horizontal ones, the expectation is that the stress component Qyz is negative, i.e. that (Biermann 1951). This is indeed seen in the simulations quite consistently (lowermost panel of Fig. 5), although at intermediate rotation low positive values can occur at high latitudes. The highest values of Qyz occur for as opposed to for Qxy. Rotational quenching of Qyz seems to be stronger and occur for lower  than for Qxy(see also Fig. 10). A similar trend was seen in convection simulations by Käpylä et al. (2004). The consistently positive values of Qyz for rapid rotation seen in convection simulations (Käpylä et al. 2004; Chan 2007, private communication) do not occur in the present calculations.

The value of Qyz always reaches a maximum at the equator. This contradicts with analytical results derived under first-order smoothing (Kitchatinov & Rüdiger 1993, 2005) and the MTA-model that predicts a maximum around for intermediate and rapid rotation. The rotational isotropization term is again needed to reduce the magnitude of the stress. It seems that to reproduce Qyz correctly, one would need to apply somewhat stronger isotropization than what is presently used (see Figs. 7 and 8).

One conclusion that can be drawn from the MTA-model results is that, although the diagonal Reynolds stresses are quite poorly reproduced in comparison to the 3D simulations, the off-diagonals have most of the qualitative features correct. This seems to imply that, to model the off-diagonals, the exact parameterization of the nonlinearities in the equation of the Reynolds stress is not crucial. The empirical rotational isotropization term helps to capture some of the missing features for the diagonal components and reduces the magnitudes of the off-diagonals to the level that is also seen in the simulations.

##### 3.3.1 Comparison with inhomogeneous simulations

The homogeneous setup used so far prevents any mean flows from being generated. This is good for the purpose of testing the sole effect of turbulent velocity field on the Reynolds stresses, but can be argued to be unphysical because astronomical objects where turbulence is important; e.g., the solar convection zone, have boundaries and cannot be considered homogeneous.

To test how much the assumption of homogeneity affects the results, we made a set of simulations with a setup where the z-boundaries are impenetrable. For the horizontal velocity components, stress-free boundary conditions are used, i.e.

 Ux,z = Uy,z = Uz = 0. (33)

To compare with the homogeneous simulations and at the same time minimize the effects of the boundaries, we average the results over and the full (x,y)-extents.

Mean flows are generated in the runs that, however, are small in comparison to the rms-velocity of the turbulence. The exception to this trend is the equator, i.e. , where a large mean shear flow develops, for which

 (34)

Similar flows have been seen in convection simulations (Chan 2001; Käpylä et al. 2004; Brandenburg 2007). The volume average of this flow is zero. For the present simulations the mean velocities near the z-boundaries can approach the speed of sound or even exceed it. This prevents us from computing models with until a statistically saturated state.

The results for the off-diagonal stresses are compared in Fig. 9. For slow rotation, , the differences are minute at all other latitudes except at the equator. The tendency for Qyz to have a minimum at the equator is reminiscent of results from convection simulations (e.g. Chan 2001; Käpylä et al. 2004; Rüdiger et al. 2005b). For more rapid rotation, the trends for different components seem to diverge: Qxy is somewhat reduced whereas Qxz seems to increase somewhat and there is hardly any change for Qyz apart from the equatorial case. For the stresses at the equator are not saturated because the large scale flow is not fully developed.

##### 3.3.2 Comparison to SOCA results

To connect to earlier studies, the 3D simulation data and MTA-model results are compared to the analytical SOCA results of Kitchatinov & Rüdiger (2005, hereafter KR05; see also Kitchatinov & Rüdiger 1993). Since analytical results are only available for the components relevant to the -effect, only the Qxy and Qyz components are thus considered. Furthermore, we restrict the comparison to a subset of the models with .

 Figure 9: Same as Fig. 5 but for homogeneous ( left panel) and inhomogeneous ( right) simulations. Open with DEXTER

 Figure 10: Comparison of the numerical simulations (solid lines), two MTA-models with , (dotted), and , (dashed), with the SOCA results of KR05 (dash-dotted). Left ( right) hand panels show ( ) for Coriolis numbers 0.15 (uppermost panels), 1.5 ( middle), and 5.4 ( lower panels). Open with DEXTER

The conventional way of writing the -effect is (e.g. Rüdiger 1989)

 (35) (36)

where is the turbulent viscosity, and the dimensionless quantities H and V are given by

 (37) (38)

where H(1), V(0), and V(1)=-H(1) depend on the Coriolis number.

The normalized stresses are known from the simulations and the MTA-model, whereas H and V can be computed analytically for the turbulence model of KR05 (see Appendix B). The results are shown in Fig. 10 for three Coriolis numbers. Our Coriolis number is smaller than that of KR05 by a factor of .

For the horizontal stress, the SOCA results and the rotationally quenched MTA-model seem to fare similarly well. The former underestimates the magnitude for small and overestimates it for large , whereas for the latter the trend is exactly the opposite. If rotational isotropization is not taken into account, the agreement is poor for all Coriolis numbers considered here.

For Qyz the standard MTA-model is almost spot on for , but overestimates the magnitudes by at least a factor of two for the other cases. As discussed in the previous section, the latitude distribution shows a mid-latitude maximum that is not present in the numerical simulations. When rotational isotropization is taken into account, at least the magnitude can be reconciled with the numerical results. The SOCA result does not fare very well in this case, predicting, in general, values that are too low and an incompatible latitude distribution with zero stress at the equator.

 Figure 11: Volume-averaged Reynolds stress components Qxy ( top), Qxz ( middle), and Qyz ( bottom) as functions of latitude and vertical turbulence anisotropy , see Eq. (32). Line styles as indicated by the legend in the lower-most panel. Coriolis number in each case is roughly 0.3. Open with DEXTER

#### 3.4 Dependence on the amount of anisotropy

Table B.2 and Fig. 11 show the results for four sets of calculations in which the Coriolis number is kept approximately constant at 0.3 whereas the turbulence anisotropy  is varied between -0.07 and -0.47. This is done by choosing suitable values for f0 and f1 so that the rms-velocity stays approximately constant.

From the figure it is seen that the ratio of the stress to the amount of anisotropy decreases monotonically as a function of . Although part of the difference can be explained by the somewhat smaller Co (0.31 as opposed to 0.34 in the other cases) in the calculation with the largest , the trend still persists.

This trend can be understood as follows: from the approximate relation we obtain . The decreasing trend of seen in the results suggests that the correlation time changes when the turbulence anisotropy is varied, i.e. decreases when is increased. This is plausible since to change different values of the forcing amplitudes, f0 and f1 (see Eq. (13)), need to be used resulting in differences in the turbulence.

 Figure 12: Off-diagonal Reynolds stresses Qxy ( top), Qxz ( middle), and Qyz ( bottom) as functions of latitude and Reynolds number, see the legend in the lower-most panel. Open with DEXTER

#### 3.5 Dependence on Reynolds number

Figure 12 shows the off-diagonal stresses as functions of latitude and Reynolds number for a constant (see also Table B.3). There is no a priori reason to expect that the stresses should depend on the Reynolds number if and if the turbulence anisotropy is kept constant. This is essentially what is borne out of the simulations, although there seems to be a weak decreasing trend as a function of Re for for Qxy and Qxz, although the results are within error bars. For Qyz, however, the largest Reynolds number case shows a distinct drop in comparison to the less turbulent cases.

The decrease in the stresses as a function of the Reynolds number is likely to have the same origin as the decrease seen when the turbulence anisotropy is increased (see the previous section). In order to obtain the same in the simulations different values, f0 and f1 are needed for different values of , the kinematic viscosity. More precisely, the less the viscosity, the more difficult it is to obtain large anisotropy. Thus, for smaller , higher ratio f1/f0 is required to achieve a given . Thus, if the interpretation that the higher the ratio f1/f0, the smaller  becomes is correct, the decreasing trend seen as a function of might be an artefact due to the small differences in the forcing between the different runs.

#### 3.6 Strouhal number from simulations

The MTA-model reproduces the simulation results for the off-diagonal stresses reasonably well when is used. Now we turn to the numerical simulations to determine independently. Although we have not been able to provide any direct support of the basic MTA-assumption, , the MTA-model is still able to reproduce many of the features of the numerical simulations adequately. This implies that a value for could be extracted from the MTA-relations.

We consider two methods to determine the Strouhal number from the simulations: (i) MTA-relations derived for the Reynolds stresses; and (ii) similar relations for passive scalar transport under the influence of rotation.

##### 3.6.1 Off-diagonal Reynolds stresses versus MTA-relations

Using the minimal tau-approximation and assuming a stationary state where , the off-diagonal Reynolds stresses can be derived from Eq. (8), yielding

 (39) (40) (41)

where are a set of relaxation times that allow for the possibility that there can be different values of for different components of the Reynolds stress. Here Qij are the components of the Reynolds stress tensor in the simulations. These equations can be solved for the via

 (42) (43) (44)

The first and the last of these equations, i.e. and , yield reasonable values for most cases, whereas is less well-behaved, producing sign changes and showing no clear trend in magnitude as a function of latitude or rotation. This is because the two terms in the denominator of Eq. (43) tend to nearly cancel each other; see the upper panels of Fig. 5. If, however, the values of in Eqs. (39) to (41) are considered the same, which is a basic MTA-assumption, it is possible to solve for Qxz in terms of the diagonal stresses

 (45)

from which it follows that

 (46)

Although the assumption that all are equal in Eqs. (39) to (41) is not exactly realized in the numerical simulations, relation (46) gives results that are compatible with those from Eqs. (42) and (44). The Strouhal number can now be computed from

 (47)

where are given by Eqs. (42), (46), and (44), respectively. Representative results are given in Fig. 13.

 Figure 13: Strouhal numbers computed from Eq. (47) at colatitude as functions of the Coriolis number for the runs listed in Table B.1. A power law proportional to shown for reference. Open with DEXTER

##### 3.6.2 Passive scalar transport with rotation

As an independent check of the dependence of the Strouhal number on rotation, we expand the passive scalar transport case, which was studied by Brandenburg et al. (2004) to cases where rotation is included. The numerical model is the same as in the runs presented so far, except that isotropic forcing is used with f0=0.01-0.03 and f1 = 0 in Eq. (13).

The turbulent passive scalar flux is denoted by , where c is the fluctuation of passive scalar density, i.e. the passive scalar concentration per unit volume. Following the MTA-approach, we solve first for the time derivative

 (48)

where the fluctuation of the passive scalar field, neglecting diffusive terms, is given by

 (49)

Here and are the mean passive scalar concentration and the mean velocity, respectively. Following Brandenburg et al. (2004) we impose a large-scale gradient of passive scalar concentration according to . In what follows G=0.1 is used. Now, assuming incompressibility and using the fact that , we arrive at

 (50)

where the last three terms denote the triple correlations. For the z-component of Eq. (50) the triple correlations are given by
 (51)

where is the reduced pressure (or enthalpy). In the non-rotating case, , and the contributions from the momentum equation also cancel on average, i.e. . The former relation follows from the periodic boundary conditions used and remains valid when rotation is added. The latter relation, however, is no longer true and is now balanced by the Coriolis term. Thus, the MTA-assumption should be applied to

 (52)

Assuming a stationary state in Eq. (50), the passive scalar fluxes can now be written as

 (53) (54) (55)

where we have retained the possibility that the values of from different equations are unequal.

In the passive scalar cases we consider isotropically forced turbulence for which even when rotation is included. Equations (53)-(55) yield

 (56) (57) (58)

which can be used to compute the Strouhal numbers

 (59)

where are given by Eqs. (56) to (58).

In the passive scalar case, the second and third order terms are indeed correlated (Brandenburg et al. 2004) according to the basic MTA-assumption, and a Strouhal number can be thus computed using

 (60)

where

 (61)

See Fig. 14 for representative results.

 Figure 14: Uppermost panel: St from triple correlation , middle panel: , corresponding to the tau in Eq. (58), lowermost panel: Strouhal numbers from Eq. (59) as functions of the Coriolis number for a constant Reynolds number of roughly nine. Open with DEXTER

##### 3.6.3 Discussion

For slow rotation, the Strouhal number is consistently between one and three when computed from the Reynolds stresses (Fig. 13). Similar values are obtained from the passive scalar transport (Fig. 14) when the Reynolds number is ten or larger. The Strouhal number computed from the triple correlations, , is more strongly dependent on the Reynolds number, but it seems to converge slowly towards a constant value near unity. These results are in line with the values required in the MTA-model and earlier studies in different contexts employing similar turbulence calculations (Brandenburg et al. 2004; Brandenburg & Subramanian 2005, 2007).

However, when the Coriolis number approaches or exceeds unity, the Strouhal numbers computed from the equations of the Reynolds stresses, Eqs. (42)-(44), decrease rapidly so that for it has dropped at least by an order of magnitude (see Fig. 13). Similar results are obtained for the passive scalar transport under the influence of rotation, see Eqs. (56)-(58) and the lower panels of Fig. 14. The trend is clearer for and , whereas for the decreasing trend is seen only for rapid rotation, i.e. when . For slow rotation, however, is almost constant and increases somewhat when the Coriolis number approaches unity. These results seem to confirm the trend seen earlier in convection simulations (Käpylä et al. 2005, 2006a).

The Strouhal number from the triple correlations follows a trend similar to , with increasing values up to after which there is a rapid decrease, see the uppermost panel of Fig. 14. For low Reynolds number St can become negative in the range , hence the gaps in the corresponding data in Fig. 14.

### 4 Conclusions

Turbulent momentum fluxes, which are described by the Reynolds stresses, were determined from numerical simulations of homogeneous rotating anisotropic turbulence. Since no large-scale shear is present, the generated Reynolds stresses correspond to contributions that are already present for uniform rotation. The resulting term is known as the -effect (Krause & Rüdiger 1974). The component responsible for the horizontal transport, Qxy, is positive and peaks around latitude  regardless of the Coriolis number. The vertical component is predominantly negative and it always peaks at the equator.

Although the numerical results for the -effect broadly agree with analytical SOCA calculations (Kitchatinov & Rüdiger 1993, 2005), the MTA-model seems to reproduce certain features of the numerical results somewhat more closely. The present numerical results do not show the enigmatic results, such as the extreme latitude distribution of Qxy or a positive Qyz for rapid rotation, which have been reported from convection simulations (e.g. Chan 2001; Käpylä et al. 2004). The difference lies most likely in our neglecting stratification and heat fluxes. The exact manner in which they affect the Reynolds stresses is not within the scope of the present paper, but should be investigated more closely in the future.

By applying the minimal tau approximation closure relation to the Reynolds stress equation, qualitatively similar results are obtained, but the magnitude of the stresses is in general too large. The vertical flux in the MTA-model, however, has a maximum at mid-latitudes for intermediate and rapid rotation. Adding an empirical rotational isotropization term (motivated in Sect. 3.1) also brings the magnitude in line with the 3D simulations. Although adding this term with this particular form has no rigorous theoretical basis, we can see that phenomenological effects of isotropization of turbulence due to rotation are indeed seen in the simulations and that the term is thus justified.

Another drawback of the MTA-model is that the diagonal components of the Reynolds tensor are rather badly reproduced since the nonlinear effects of rotation manifest in the numerical simulations are not explicitly taken into account. The empirically added rotational isotropization term augments the magnitudes, but not the latitude distribution. Furthermore, no direct evidence of the validity of the MTA-assumption was found in the numerical simulations. Contrasting the behavior of the diagonal components to the fairly good correspondence between the numerical simulations and the MTA-model for the off-diagonal components leads us to conclude that, whereas the behavior of the diagonal components is dominated by the inadequately modeled nonlinear effects, the off-diagonals are fairly well presented by the linear terms.

A Strouhal number of order unity in the MTA-model gives best fits to the numerical results. Fitting the numerical results to expressions derived under the MTA, similar values of  are found for slow rotation. For Coriolis numbers of order unity or larger, however, the Strouhal number obtained in this manner decreases rapidly. In the passive scalar case, the situation is somewhat more complex, although a similar decreasing trend of the Strouhal number is recovered for rapid rotation, see Fig. 14. These results are in accordance with earlier results from convection simulations (Käpylä et al. 2005, 2006a) using Reynolds stresses or correlation analysis of the velocity field.

A related aspect in turbulent transport that requires closer study is the turbulent viscosity (see preliminary results in Käpylä & Brandenburg 2007) and the possibility of a -effect due to the anisotropy induced by a large-scale shear flow (Leprovost & Kim 2007). These matters will be considered in more detail in a future publication.

Acknowledgements
The computations were performed on the facilities hosted by the Center of Scientific Computing in Espoo, Finland, who are financed by the Finnish ministry of education. P.J.K. acknowledges the financial support from Helsingin Sanomat foundation and the Academy of Finland grant No. 121431. P.J.K. acknowledges the hospitality of Nordita during the program Turbulence and Dynamos'' during which this work was finalized. The anonymous referee is acknowledged for the critical reading and helpful comments on the manuscript.

### Appendix A: Stationary solutions for the stresses from MTA

Setting and using the MTA-closure , and parameterizing the contributions of the forcing as in Eq. (22) gives six equations for the six unknowns of the symmetric tensor Qij, i.e.

 (A.1) (A.2) (A.3) (A.4) (A.5) (A.6)

It is possible to solve for Qij in terms of , , , and Qij(0). From Eqs. (A.1) to (A.6) it is clear that only three components of Qij are independent. Thus it is sufficient to solve for the off-diagonal components. After some algebra we arrive at

 (A.7) (A.8) (A.9)

where

 (A.10) (A.11)

and is given by Eq. (27).

In the present study the forcing is such that Qyy(0) - Qxx(0) = 0, so the equations reduce to

 (A.12) (A.13) (A.14)

### Appendix B: Coefficients of the -effect from SOCA

The fluxes of angular momentum have commonly been parameterized by Eqs. (35) and (36), and the normalized fluxes by Eqs. (37) and (38). Kitchatinov & Rüdiger (2005) computed these coefficients using SOCA

 (B.1) (B.2)

where is their definition of the Coriolis number, and the turnover time. Note that there is a difference of in comparison to our definition, Eq. (17). For simplicity, we retain in the expressions that follow. The coefficients H(1) and V(0) are given by

 (B.3) (B.4)

where is the correlation length, the density scale height, and a = 2 an anisotropy parameter'' that reduces the amount of anisotropy for slow rotation.

The functions Ii and Ji are given by

 (B.5) (B.6) (B.7) (B.8)

In our model there is no density stratification, so is taken to be a free parameter that we choose to be equal to unity.

## References

### Online Material

Table B.1: Summary of the turbulence anisotropies and normalized Reynolds stresses, , for the calculations. , , and grid resolution 643 was used in all runs.

Table B.2: Summary of the turbulence anisotropies and normalized Reynolds stresses, , for a set of runs with varying turbulence anisotropy. , , , and grid resolution 643 was used in all runs.

Table B.3: Turbulence anisotropies and Reynolds stresses for varying Reynolds numbers. and were used in all runs. .