A&A 386, 331-346 (2002)
DOI: 10.1051/0004-6361:20020126

Box simulations of rotating magnetoconvection

Spatiotemporal evolution

U. Ziegler

Astrophysikalisches Institut Potsdam, An der Sternwarte 16, 14482 Potsdam, Germany

Received 31 July 2001 / Accepted 7 January 2002

Abstract
This paper reports 3D numerical simulations of compressible, thermal convection in a small rectilinear domain placed tangentially on a rotating sphere at latitude $45^{\circ}$. The spatiotemporal evolution of convection is studied in an initially non-homogeneous toroidal magnetic field located at the interface of a 2-layer unstable/stable polytropic stratification. The effects from a variation of the rotation rate and magnetic field strength are explored. In weak field convection the solutions bear close resemblance to the field-free situation with the magnetic field treated as a passive ingredient. In this case a significant amount of magnetic energy is transported downwards into the stable layer by penetrative motions where the field is concentrated in small-scale tube-like features. In the case of a dominating magnetic field, the overall structure of convection changes dramatically towards a two-dimensional, more laminar flow with convective motions occurring in columnar cells aligned with the mean magnetic field. The latter quickly becomes flat in the convection zone due to magnetic buoyancy effects. Magnetic quenching of the flow heavily influences the mixing properties of convection, thereby, reducing the fluxes of kinetic energy and enthalpy and suppressing the downward transport of magnetic energy. Adding rotation in strong field convection generates streaming motions parallel to the mean field. These motions contain a large fraction of the total kinetic energy. The horizontally-averaged (mean) flows, on the other hand, are less energetic. In the presence of strong rotation the horizontal mean flows are triggered by inertial oscillations. Averaging the mean flows over time gives an estimate for steady flows persistent in a convection zone. Whereas the vertical flow component obtained in this way shows a systematic dependence on the rotation rate and field strength, the flows in the meridional and zonal directions depend in a more complex manner on these parameters. This includes reversals in their orientation when going from moderate to strong rotation and when increasing the magnetic field strength.

Key words: convection - magnetohydrodynamics - turbulence - magnetic fields

1 Introduction

Observations of the solar atmosphere reveal the existence of a rich set of magnetic flux elements of various scales making up active regions on the surface. Diameters of these elements range from about 100 km for intense flux tubes up to several thousand km for sunspots. It is widely accepted that the appearance of the surface magnetic field is largely controlled by the dynamical processes in the layers beneath the surface where thermal convection plays a dominant role. Any serious attempt to explain the field topology requires thus an understanding of the phenomenon of convection, its interaction with the magnetic field and, of course, the mechanisms of magnetic field generation in the interior. The current picture postulates an ordered toroidal magnetic field which derives from a dynamo process taking place in a thin shear layer located at the base of the convection zone. The existence of such a shear layer has indeed been proposed on the basis of helioseismological measurements. These yield an angular velocity profile which is characterized by a nearly rigid rotation in the deeper layers and which is almost constant on radial lines throughout the convection zone. The toroidal field is continuously regenerated out of poloidal field by differential rotation ($\Omega$-effect) whereas the poloidal field is produced out of expelled toroidal flux subject to cyclonic convective motions ($\alpha$-effect). Of course, the whole process relies on efficient transport mechanisms in order to allow both expulsion of toroidal field from the shear layer into the convection zone and, vice versa, downward transport of poloidal field from the convection zone into the shear layer. The strong toroidal field may occasionally disrupt due to the action of magnetic buoyancy instabilities producing smaller field fragments. These fragments, ideally in the form of discrete flux tubes, eventually rise through the convection zone to the surface where it finally emerge as the features one can observe.

The linear theory of magnetoconvection has lead to considerable insight under which circumstances convection sets in and what the growth rates are. Exploration of the non-linear regime, however, exclusively depends upon numerical simulations. The computational modeling of the convective outer envelope of a star suffer from the enormous complexity of the problem. Many different physical processes are involved and flows are expected to be in a highly turbulent state. Reynolds numbers, ${\rm Re}$, of the order of 10¹², like in the solar convection zone or even higher are the rule. Such high Reynolds number flows are, of course, intractable numerically as it would require a grid resolution far beyond what can be achieved even on supercomputers. This lack of attainable resolution might be a point of criticism, especially in global models that consider convection in a full spherical shell (see e.g.Glatzmaier & Toomre 1995). Expected small Prandtl numbers, Pr(ratio of viscous dissipation to thermal dissipation), and small magnetic Prandtl numbers, ${Pr}_{\rm m}$(ratio of viscous dissipation to Joule dissipation), are another severe problem. ${Pr}\ll 1$ means that viscous dissipation occurs on much smaller spatial scales than thermal dissipation. Therefore, in computer simulations the Prandtl number has to be restricted to values ${Pr}\,\raise 0.3ex\hbox{$\space > $ }\kern -0.8em \lower 0.7ex\hbox{$\sim$ }\,0.1$ in order to resolve both scales accurately on a numerical mesh with typically 100 grid points per coordinate direction. In total, with present available computer resources it is impossible to take care of all aspects of convection in an accurate manner. Instead, idealized models must stand to racket trying to be as realistic as possible for the problem under study. To cover at least the moderately turbulent regime of compressible convection, a local model is best. Even then, however, Reynolds numbers much in excess of order 10³ are hard to realize.

In this series of papers, numerical results based on such an idealized, local model of thermal convection are presented, including the effects of stratification, rotation and magnetic field. The applied model in many respects takes guidance from previous work and tries to combine the advantages of all those approaches. For example, Cattaneo et al.(1991) (hereafter CBTMH) studied the transition from laminar to turbulent compressible convection without rotation and without an overshoot region in a polytropic stratification identical to our unstable layer. Brummell et al.(1996) (hereafter BHT) extended the work of CBTMH by including the effects of the Coriolis force. Penetrative, non-rotating convection in two space dimensions has been investigated by Hurlburt et al.(1994) using a similar 2-layer stable/unstable configuration as in this study. The current model does not, however, account for a background shear flow which is present at the base of a real convection zone. Therefore, the self-consistent generation of magnetic field by an $\Omega\alpha$-dynamo is ab initio prohibited. Although such an investigation would be of considerable interest, it is not the primary aim here but must await future studies. Rather, the present paper focuses on the spatiotemporal evolution of convection, the development of mean flows and the properties of energy transport. The effects from systematically varying the strength of rotation and strength of the applied magnetic field are explored. In forthcoming papers the dynamics of convective penetration and dynamo properties will be examined.

2 Basic equations and numerical solver

The evolution equations describing thermal magnetoconvection in a stratified, visco-resistive medium including the effect of the Coriolis force are given by:

$$\partial_{\rm t}\rho = -{\nabla\cdot}(\rho \vec{v})$$ (1)

$$\partial_{\rm t}\left(\rho \vec{v}\right) -{\nabla\cdot}(\rho\vec{vv})-\nabla p +\nabla\cdot\tau -2\rho\vec{\Omega}\times\vec{v}$$

$$+\frac{1}{\mu}(\nabla\times\vec{B})\times\vec{B} +\rho \vec{g}$$ (2)

$$\partial_{\rm t}\vec{B} \nabla\times(\vec{v}\times\vec{B} -\eta\nabla\times\vec{B})$$ (3)

$$\partial_{\rm t} e -{\nabla\cdot}(e\vec{v})-p\nabla\cdot\vec{v} +Q_{\rm vis} +{\nabla\cdot}(\kappa\nabla T)$$

$$+\frac{\eta}{\mu}\vert\nabla\times\vec{B}\vert^2.$$ (4)

The dependent variables of the problem are the mass density $\rho$, pressure p, thermal energy density e, temperature T, velocity $\vec{v}$ and magnetic field $\vec{B}$. The remaining quantities in the above equations are the viscous stress tensor $\tau$ given by

$${\tau}_{ij}=\nu\rho\left(\nabla\vec{v}_{ij} +\nabla\vec{v}_{ji} -\frac{2}{3}{\nabla\cdot}\vec{v}\delta_{ij}\right)$$ (5)

where $\nu$ denotes the kinematic viscosity coefficient, the viscous heating term

$$Q_{\rm vis}=\sum_{ij}\tau_{ij}\nabla\vec{v}_{ij}\,,$$ (6)

the gravitational field $\vec{g}$, the magnetic permeability $\mu$, the thermal conductivity coefficient $\kappa$, and the magnetic diffusivity $\eta =1/\sigma \mu$ where $\sigma$ is the electrical conductivity.

Ionization of matter and radiation transfer are neglected. Hence, the mathematical description is suited to model the deeper layers of a star's convection zone where both processes are believed to be relatively unimportant. It is assumed that the thermodynamic variables are related through the ideal gas equation of state

$$p=\frac{k}{\overline{\mu}m_{\rm u}}\rho T=(\gamma -1)e$$ (7)

where k is Boltzmann's constant, $m_{\rm u}$ the atomic mass unit, $\overline{\mu}$ the mean molecular weight and $\gamma =c_{\rm p}/c_{\rm V}=5/3$ the ratio of specific heats, $c_{\rm p}$ and $c_{\rm V}$. Any chemical composition of the gas is ignored ( $\overline{\mu}=1$). Furthermore, as usual, the fluid is assumed to be non-magnetizable with $\mu=4\pi\times 10^{-7}$ Vs A^-1 m^-1.

The set of Eqs. (1)-(4) is solved with the finite-difference, fractional-step code NIRVANA (Ziegler 1998; Ziegler 1999). For the purpose of simulating thermal convection at low Prandtl number the original code has been extended by a time-implicit solver for the heat conduction term in Eq. (4) in order to avoid an otherwise overly restrictive time-step. An ADI method is used which differs somewhat from the standard ADI procedure in the way the dimensional splitting is performed (Ames 1977).

NIRVANA, in addition, has been improved in the numerical treatment of the Coriolis force term. Inertial oscillations described by

$$\partial_{\rm t}\vec{v}=-2\vec{\Omega}\times\vec{v}$$ (8)

are handled in a way taking advantage of the analytical solution to (8). Denoting $\vec{v}_0=\vec{v}(t_0)$ the velocity at a time instant t₀the solution of (8) for later times t>t₀ is given by

$$\vec{v}(t) (\hat{\vec{\Omega}}\cdot\vec{v}_0)\hat{\vec{\Omega}} +(\vec{v}_0\times\hat{\vec{\Omega}}) \sin(2\Omega(t-t_0))$$ =

$$+(\hat{\vec{\Omega}}\times\vec{v}_0) \times\hat{\vec{\Omega}} \cos(2\Omega(t-t_0))$$ (9)

where $\hat{\vec{\Omega}}$ is the angular velocity unit vector. The corresponding partial update in the operator-split scheme of NIRVANA from a time level tⁿ to a time level $t^{n+1}=t^n+\delta t$then reads

$$\vec{v}^{n+1} (\hat{\vec{\Omega}}\cdot\vec{v}^n)\hat{\vec{\Omega}} +(\vec{v}^n\times\hat{\vec{\Omega}}) \sin(2\Omega\delta t)$$ =

$$+(\hat{\vec{\Omega}}\times\vec{v}^n) \times\hat{\vec{\Omega}} \cos(2\Omega\delta t)$$ (10)

where $\delta t$ denotes the actual time-step. Former use of simple forward Euler time-stepping for (8) is seen to fail retaining initial hydrostatic equilibrium in the case of fast rotation. The failure relies upon a kind of oscillating instability arising when $\delta t$ does not explicitly take into account the time-scale for inertial oscillations which is given by $\pi/\Omega$. On the contrary, Eq. (10) provides the correct velocity update for arbitrary $\delta t$. Note, however, that for very large $\Omega$ the Coriolis force term becomes stiff and this may nevertheless lead to inaccurate results when combined with less stiff terms within the fractional-step scheme. Rotation rates are not that large here, hence, such problems never arise.

3 The convection model

	Figure 1: Sketch of the 2-layer local model representing a portion of a spherical shell at latitude $\theta =45^{\circ }$.
Open with DEXTER

A rectangular domain placed tangentially on a sphere at north latitude $\theta =45^{\circ }$ is considered. The spatial extent of this box is assumed to be small compared to the sphere's radius so that it is justified to describe the undisturbed atmosphere by a planar stratification in a constant gravitational field $\vec{g}=-g\hat{\vec{z}}$, g>0. Cartesian coordinates with unit vectors $\hat{\vec{x}}$ (pointing equatorward), $\hat{\vec{y}}$ (pointing in azimuthal direction) and $\hat{\vec{z}}$ (pointing in radial direction) are adopted. A 2-layer configuration is constructed which consists of an upper, convectively unstable region of depth d placed on top of a stable region with the same vertical extent. At the bottom of the stable layer a constant heat flux is driven. The model setup is illustrated schematically in Fig. 1. The influence of the sphere's rotation on the ensuing convection will be characterized by the Taylor number, ${Ta}=4\Omega^2d^4/\nu^2$, which defines the ratio of the viscous timescale to the rotation period squared. In the box coordinate system the angular velocity vector is given by $\vec{\Omega} =-\Omega\sin\theta \hat{\vec{x}} +\Omega\cos\theta \hat{\vec{z}}$ and, thus, has two components. The computational domain spans a volume $(x,y,z)\in [0,4d]\times[0,4d]\times[-2d,0]$discretized by $96\times 96\times 128$ grid points. The grid points are uniformly distributed in each coordinate direction. The aspect ratio of grid spacings is $1:1:\frac{3}{8}$. The smaller grid spacing in z accounts for the presence of vertical stratification.

The thermal conductivity coefficient $\kappa$ which controls the initial temperature gradient and, hence, the stability of the layer against convection is assumed to be a piecewise constant function:

$$\kappa(z)=\left\{\begin{array}{lll} \kappa_{\rm s} & \quad -... ...uad -d<z<0 & \quad \mbox{(unstable region)} \end{array}\right.$$ (11)

where the subscript "s'' ("c'') denotes the stable (convectively unstable) region. With this assumption the initial static state, which must be a solution of the equations

$$0 = -\frac{{\rm d}p}{{\rm d}z}-\rho g$$ (12)

$$-\kappa\frac{{\rm d}T}{{\rm d}z} = {\rm const.} = F_0$$ (13)

$$p = \frac{k}{m_{\rm u}}\rho T$$ (14)

is a piecewise polytropic stratification expressed by

$$\left\{\begin{array}{ll} \frac{F_0}{\kappa_{\rm s}}\left(z_1-z\ri... ...ac{F_0}{\kappa_{\rm c}}\left(z_2-z\right)+T_2 & \quad -d<z<0 \end{array}\right.$$ (15)

$$\rho \left\{\begin{array}{ll} \rho_1\left(\frac{T}{T_1}\right)^{m_{\rm... ... \rho_2\left(\frac{T}{T_2}\right)^{m_{\rm c}} & \quad -d<z<0 \end{array}\right.$$ (16)

$$\left\{\begin{array}{ll} \frac{k}{m_{\rm u}}\rho_1 T_1\left(\frac... ... T_2 \left(\frac{T}{T_2}\right)^{m_{\rm c}+1} & \quad -d<z<0 \end{array}\right.$$ (17)

where the subscript 0 (1, 2) refers to values at vertical position z₀=-2d (z₁=-d, z₂=0). Instead of prescribing the gravity field strength g the polytropic indices $m_{\rm s}$ and $m_{\rm c}$ which determine the degree of stratification in both layers are input parameters of the problem. These parameters are related to the quantity g by

$$m_{\rm s} \frac{m_{\rm u}}{k}\frac{\kappa_{\rm s}}{F_0}g-1$$ = (18)

$$m_{\rm c} \frac{m_{\rm u}}{k}\frac{\kappa_{\rm c}}{F_0}g-1$$ = (19)

from which the relation

$$\frac{\kappa_{\rm s}}{\kappa_{\rm c}}=\frac{m_{\rm s}+1}{m_{\rm c}+1}$$ (20)

can be derived. It is therefore sufficient to specify only one thermal diffusion coefficient, say $\kappa_{\rm c}$. F₀ denotes the heat flux injected at the bottom z-boundary. With $m_{\rm s}$, $m_{\rm c}$, F₀ and $\kappa_{\rm c}$ given g is fixed. Continuity of the solution (15) and (16) at z₁ requires

$$\frac{F_0}{\kappa_{\rm c}}d+T_2$$ T₁ = (21)

$$\rho_1 \rho_2\left(\frac{T_1}{T_2}\right)^{m_{\rm c}}.$$ = (22)

The level of forcing of the convective instability is determined by the Rayleigh number, Ra, which for a viscous, thermal conducting fluid is given by (ignoring modifications due to rotational effects and magnetic fields; Spiegel 1965)

$${Ra}=\frac{\rho}{T}\frac{gc_{\rm p}d^4}{\kappa\nu}\left(\frac{{\rm d}T}{{\rm d}z} -\frac{g}{c_{\rm p}}\right)\cdot$$ (23)

For convection to set in at all, Ra must exceed a critical value which, very roughly, is of the order of 10³ (Gough et al.1976). A polytropic layer with index greater than the adiabatic one, $1/(\gamma -1)=1.5$, is convectively stable (negative Ra). For the lower layer we choose $m_{\rm s}=9$ which is such high that overshooting convective motions do not reach the lower z-boundary of the domain. A polytropic index of $m_{\rm c}=1$ is adopted for the convection zone. The convection zone spans about 5 pressure scale heights and its density contrast is $\approx$11 for our choice of Ra and Pr(see below). The density contrast over the full domain is $\approx$50.

To complete the model description the quantities F₀, $\rho_2$ and T₂ must be specified and the problem parameters must be declared. First are asserted as boundary conditions. Parameters of the problem are the polytropic indices (see above), the conductivity coefficient $\kappa_{\rm c}$ (recall that $\kappa_{\rm s}$ is then fixed by (20)), the viscosity coefficient $\nu$ and the rotation rate $\Omega$. However, instead of prescribing $\kappa_{\rm c}$, $\nu$ and $\Omega$ the characteristic numbers Ra, Pr and Ta serve to control the simulations. Furthermore, in order to facilitate comparison with earlier work reference values are introduced: coordinates are measured in units of d, time in units of the sound crossing time $\tau_{\rm S2}=d/c_{\rm S2}$ where $c_{\rm S2}=(p_2/\rho_2)^{1/2}$ is the isothermal sound speed at z₂, velocity in units of  $c_{\rm S2}$, density in units of $\rho_2$, pressure in units of p₂ and temperature in units of T₂. Periodic boundary conditions are assumed in the horizontal directions. At the lower and upper z-boundary stress-free, impermeable conditions are imposed for the velocity i.e. $\partial_zv_x=\partial_zv_y=v_z=0$at z=z₀,z₂. At z=z₂ the temperature is fixed, T₂=1. Also, the heat flux F₀ is fixed through the temperature gradient $\partial_zT=-10\kappa_{\rm c}/\kappa_{\rm s}=-2$ at z=z₀.

The magnetoconvection simulations start with initial conditions taken from the evolved states of non-magnetic convection superimposed by a magnetic field. Most of magnetoconvection models to date assume an initially homogeneous magnetic field for the sake of convenience (Brandenburg et al.1990; Nordlund et al.1992; Hurlburt et al.1996; Weiss et al.1996). Here the choice is guided from the idea of a dynamo-generated toroidal magnetic field located at the base of the convection zone. The field is restricted in vertical direction to the region -3d/2<z<-d/2 and has the profile of a shifted Gauß-function:

$$\vec{B}(z)=\frac{B_1}{9}\left(10\exp\left\{-4\frac{(z-z_1)^2}{d^2} \ln10\right\}-1\right)\hat{\vec{y}}.$$ (24)

It is symmetric with respect to the stable/unstable interface with maximum B₁ on the interface (z=z₁) and vanishes at the ends of the given interval. The strength of the field will be expressed by $\beta=2\mu p_1/B_1^2$ which is a further control parameter of the simulations. The magnetic field is introduced in a way to keep the total (thermal + magnetic) pressure continuous. For that the thermal pressure is locally reduced to compensate magnetic pressure. Also, the temperature structure is let unchanged by decreasing the density. In dimensionless notation the magnetic field is measured in units of B₁.

Of course, periodic boundary conditions are adopted for the magnetic field in the horizontal directions. Some kind of "open'' boundary conditions are used at the z-boundaries z=z₀,z₂. These can be formulated as follows. First,

$$\partial_zB_x=\partial_zB_y=0$$ (25)

which, on the staggered grid representation, means that B_x(x,y,z₊)=B_x(x,y,z_-) and B_y(x,y,z₊)=B_y(x,y,z_-) where z₊ (z_-) denotes the z-coordinate of the corresponding variable in the first ghost zone (first active zone) at the boundary. To calculate B_z(x,y,z₊) one can make use of the discretized form of ${\nabla\cdot}{\vec{B}}=0$ which gives

$$B_z(x,y,z_+)=B_z(x,y,z_-)-\left(\partial_xB_x+\partial_yB_y\right)_- (z_+-z_-).$$

(26)

Note that, because of the above condition for B_x and B_y, $(\partial_xB_x+\partial_yB_y)_-= (\partial_xB_x+\partial_yB_y)_+$ and B_z(x,y,z₊) is therefore uniquely determined. The chosen boundary condition is neither a perfect conductor boundary condition (which requires B_z=0) nor of "pseudo-vacuum'' type. Rather, it is a condition which suggests an approximate extrapolation of a magnetic field line across the boundary. The drawback is that, in contrast to the other mentioned boundary conditions, it gives a non-vanishing Poynting flux in vertical direction which then has to be considered in the magnetic energy budget. Although a clear physical interpretation of this boundary condition is missing, it seems quite convenient to the author to use it for the current problem relaxing the restricting condition of a perfectly conducting boundary.

Ultimately, interpretation of the numerical data is aided by defining mean values of a space- and time-dependent variable f. This is done in different ways: horizontally-averaged quantities are denoted by $\langle f\rangle(z,t)$, volume-averaged quantities by $\langle\langle f\rangle\rangle(t)$ and time-averaged quantities by $\overline{f}(\vec{x})$. In case of volume averages a subscript "s'' ("c'') indicates that the averaging process has been restricted to the stable layer z<-1 (convectively unstable layer z>-1). Finally, averages over space and time can be combined e.g. $\overline{\langle f\rangle}(z)$.

4 Numerical results

4.1 Preparing initial conditions: non-magnetic convection

	Figure 2: RMS velocity in the convection zone as a function of time for non-rotating convection (solid line), ${Ta}=6\times 10^4$ (dash-dotted line) and ${Ta}=6\times 10^5$ (dashed line).
Open with DEXTER

The primary aim of the following non-magnetic models is to provide initial conditions for the more complex magnetoconvection problem. However, since this is the first attempt to model convection with NIRVANA, it deserves a detailed discussion in its own. A series of runs has been performed for Taylor numbers ${Ta}=\{0,6\times 10^4,6\times 10^5\}$, Rayleigh number ${Ra}=3\times 10^5$ and Prandtl number Pr=0.1. Note that, since $\nu\rho$ and $\kappa$ are assumed to be constant (layer-wise in case of $\kappa$), ${Pr}\propto \nu\rho/\kappa$ does not depend on z but Ra and Ta in general do so because their definitions involve the density and/or temperature which are both z-dependent functions. The above quoted values correspond in particular to a height z=-1/2i.e.the parameters are evaluated in the middle of the unperturbed convection zone. Ra lies well above the critical value which places convection in the supercritical regime.

	Figure 3: Mean profiles of entropy, density, turbulent Mach number and Reynolds number at $t_{\rm e}$ for Ta=0 (solid line), ${Ta}=6\times 10^4$ (dash-dotted line) and ${Ta}=6\times 10^5$ (dashed line). In the entropy- and density plots the initial distribution is also shown (thin solid line).
Open with DEXTER

Convection is initiated by perturbations in the thermal energy density and the system is evolved to a time $t_{\rm e}=27.3$ corresponding to roughly 13 convective turnover times. A statistically steady state solution quickly develops, corroborated by the fact that the resulting RMS velocity $\langle\langle(\delta\vec{v})^2\rangle\rangle^{1/2}_{\rm c}$becomes almost constant in time (Fig. 2). The initial growth phase starts somewhat later for ${Ta}=6\times 10^5$ because of the stabilizing effect of rotation for the onset of convection (the critical Rayleigh number increases with rotation; Chandrasekhar 1961). The strength of convection as measured by $\langle\langle(\delta\vec{v})^2\rangle\rangle^{1/2}_{\rm c}$is nearly independent of rotation for the considered range of Tawith a slight trend towards larger values for increasing Ta. There is, however, a significant bias in vertical kinetic energy in non-rotating convection reflecting the buoyancy-driven generation of vertical momentum whereas rotating convection promotes equipartition of kinetic energy. As pointed out by BHT, this is related to the fact that the Coriolis force provides an effective mechanism to convert vertical momentum into horizontal momentum through linear coupling. Such a mechanism is missing in non-rotating (or non-tilted) convection where momentum conversion can occur only through non-linear effects.

	Figure 4: Perspective view of the computational box displaying the velocity structure at $t_{\rm e}$ for Ta=0 (left panel) and for ${Ta}=6\times 10^5$ (right panel). The grey-scale image shows vzat or near the domain faces and in the horizontal plane z=-0.9 (displaced vertically), and the arrows represent the 3D velocity field.
Open with DEXTER

Figure 3 displays the z-profiles of some selected mean quantities at $t_{\rm e}$. The entropy $s=\ln(p \rho^{-\gamma})$, density, turbulent Mach number $\langle (\delta\vec{v})^2\rangle^{1/2}/c_{\rm S}$and Reynolds number (definition follows later) is shown. Different regions can be distinguished with depth. Near the top boundary a thermal layer of thickness $\approx$0.2 has developed followed by a nearly isentropic zone which extends up to the stable/unstable interface (consult the entropy plot). The build-up of such an isentropic region is peculiar to convection which serves to equalize entropy. The deviations from ideal isentropy in the interior are somewhat larger for rotating convection which is likely due to an enhanced horizontal mixing of thermodynamic quantities with the consequence of a decrease of the vertical enthalpy flux (see below), hence, delaying the homogenization process. A similar effect has been found for rotating convection in Boussinesq approximation (Julien et al.1996). The relaxation towards an isentropic state is associated with a smart restructuring of the density stratification (consult the density plot). The changes are of significance in the interior of the convection zone where, for instance, $\langle \rho\rangle-\langle \rho \rangle_0 \approx 0.15\langle \rho\rangle_0$at z=-0.3 but almost negligible elsewhere.

	Figure 5: Scatter plots of the kinetic energy flux $F_{\rm kin}$ (left), enthalpy flux $F_{\rm e}$ (middle) and convected energy flux $F_{\rm c}$ (right) at depth z=-0.5 and time te for Ta=0. The fluxes are scaled to the initial heat flux F0.
Open with DEXTER

The Mach number remains below unity everywhere which characterizes the flow as subsonic. This does not exclude the possibility that it occasionally becomes supersonic. A more detailed data analysis shows, however, that such events are very rare and, if present, are transient phenomena preferentially occurring near the top boundary. There is no evidence for the existence of permanent shocks. Supersonic convection has been studied by Malagoli et al. (1990) using a higher-order Godunov scheme (PPM, Woodward & Colella 1984). PPM is especially suited in that context because its numerical dissipation is very small, hence, allowing simulations with lower effective Prandtl numbers compared to our models. Indeed, Malagoli et al.report on the existence of short-lived shock structures which form spontaneously in the vicinity of strong downflows near the top layer where the sound speed is smallest. Such features are not observed here, probably because of insufficient resolution.

A measure of the degree of non-linearity of the flow is the Reynolds number. It seems convenient to define

$${Re}(z)=\frac{\langle(\delta\vec{v})^2\rangle^{1/2}d}{\nu}$$ (27)

where the depth d of the convection zone is used as the typical length scale and the RMS velocity as the typical velocity. In accordance with previous work we consider a value of 10³ or greater as indication for a moderately turbulent flow. This is the case in the spatial region from $z\approx -1.45$(-1.35, -1.2) to $z\approx -0.5$ for Ta=0 ( $6\times 10^4$, $6\times 10^5$). The turbulent zone thus extends into the stable layer because of the penetration effect of convective motions. The depth of penetration depends on the rotation rate and decreases for increasing Ta (Ziegler & Rüdiger 2001, hereafter ZR). This explains the rightward shift of the Re(z)-profile for larger rotation rates. To be complete, the dynamics below the penetration zone is dominated by gravity waves which are excited by the overshooting motions.

To get an impression of the three-dimensional structure of convection Fig. 4 (left panel) shows an instantaneous view of the velocity pattern for Ta=0. The grey-scale image represents the vertical velocity component where negative values are shown by darker tones and positive ones by lighter tones. Many characteristic features observed in previous simulations of compressible convection are also found here: near the top of the domain motions appear in a cellular-like network which is made up of broader regions where warmer material ascends surrounded by narrower downflows which carry cooler gas. The downflows are connected and form a polygonal structure. This quasi-laminar network only exists in the thermal boundary layer consistent with lower Reynolds numbers there. With increasing depth the network becomes more and more disintegrated. The downflow lanes disconnect into sheet-like structures and fast moving plumes embedded in a much more disorganized flow which makes up the turbulent region below the thermal boundary layer. The fast-moving plumes have a coherence length of the order of d and are comparatively long-lived objects.

The picture emanating from non-rotating convection largely persists for rotating convection (Fig. 4, right panel). That is a cellular structure at the surface, a turbulent interior region and the existence of vertically coherent downflows. Main differences are found in the temporal behavior and topology of the surface network. The eddies are smaller in size and the downflow lanes appear more curvaceous in shape and less connected with some of it ending abruptly in upflow islands. The surface network is enduringly influenced by inertial oscillations of larger-scale flows existing in the solution which, in effect, lead to a distortion and reorganization of the overall pattern on a time scale given by the inertial time scale.

An another important aspect of convection is its mixing properties. Because of the asymmetry between downflow regions and upflow regions typical for convection in a stratified medium it makes sense to work out the contributions of the downflows and upflows to the convective transport separately. It is convenient to use a phase-space representation where the various energy fluxes are plotted as a function of vertical velocity v_z. The kinetic energy flux $F_{\rm kin}$ and enthalpy flux $F_{\rm e}$ are

$$F_{\rm kin} \frac{1}{2}\rho\vec{v}^2v_z$$ = (28)

$$F_{\rm e} c_V\rho \delta Tv_z$$ = (29)

where $\delta T=T-\langle T\rangle$ denotes temperature fluctuations relative to the mean temperature $\langle T\rangle$. The sum of both of these fluxes define the convected energy flux, $F_{\rm c}$. Figure 5 shows scatter plots of all three kinds of energy fluxes evaluated in the horizontal plane z=-0.5 and at time $t_{\rm e}$ for Ta=0. The $(F_{\rm kin},v_z)$-diagram exhibits a cubic distribution which is expected because of the v³-dependence of  $F_{\rm kin}$. Note that each dot represents the value at a certain grid point in the plane and their sum gives the horizontal average $\langle F_{\rm kin}\rangle(z=-0.5)$. Apparently, a significant contribution to $\langle F_{\rm kin}\rangle(z=-0.5)$ comes from the faster downflows. These flow features thus provide channels through which kinetic energy can efficiently be transported downwards. Both upflows and downflows give net positive contributions to the enthalpy flux but with a clear preponderance on the part of downflows. It are the fast downflows again which enable an efficient transport of thermal energy, however, upwards. One remarkable property of non-rotating convection is seen when looking at the phase-space representation of the convected energy flux. The $(F_{\rm c},v_z)$-diagram reveal that the downflows transport almost as many kinetic energy downwards as thermal energy upwards expressed by the nearly symmetric appearance of the dot distribution on the downflow side. This balance infers that the net contribution to $\langle F_{\rm c}\rangle(z=-0.5)$ must derive from the broader upflows. This surprising effect has already been detected by CBTMH using a mixed pseudo-spectral/finite-difference code. We basically confirm their results using NIRVANA.

	Figure 6: Mean convected energy flux $\langle F_{\rm c}\rangle$ normalized to F0 at $t_{\rm e}$ for Ta=0 (solid line), ${Ta}=6\times 10^4$ (dash-dotted line) and ${Ta}=6\times 10^5$ (dashed line).
Open with DEXTER

Energy mixing turns out not much sensitive to slow rotation. Larger differences exist, however, for ${Ta}=6\times 10^5$. One of the more prominent changes is outlined in Fig. 6 which presents the mean convected energy flux as a function of z. For ${Ta}=6\times 10^5$ the convected energy flux is enhanced in the interior of the convection zone which is due to a larger contribution of the upflows. The sharp decline at $z\approx -1$ and subsequent sign reversal is a consequence of buoyancy braking in the stable layer which inhibits the upward transport of enthalpy. $F_{\rm c}$ is then dominated by the negative contribution of the kinetic energy flux. The detailed behavior of $\langle F_{\rm c}\rangle (z)$ below the stable/unstable interface therefore reflects to some degree the dependence of convective penetration on the rotation rate suggesting a smaller penetration depth for increasing rotation rate (see ZR).

4.2 Magnetoconvection

Figure 7: Transition from non-magnetic convection to magnetoconvection illustrated by the time evolutions of the horizontally-averaged RMS velocity (grey-scale image; darker tones mean larger velocities) and its volume-averaged counterpart (line plot; right axis). The cases ( ${Ta}=6\times 10^5,\beta =50$) (top) and ( ${Ta}=6\times 10^5,\beta =5$) (bottom) are shown. The magnetic field is introduced at $t_{\rm e}$.

Open with DEXTER

The simulations of magnetoconvection start with initial conditions given by the final states of non-magnetic convection just presented with Pr=0.1, ${Ra}=3\times 10^5$ and ${Ta}=\{0,6\times 10^4,6\times 10^5\}$. These states are imposed by a magnetic field of type (24). In the absence of convection the unperturbed magnetic stratification is stable against buoyancy effects but becomes unstable if perturbations are attached inducing local density gradient inversions. To integrate the role of such (non-convective) instabilities for the later convection process special runs have been performed where heat conduction has been switched off but the finite perturbations have been retained. The following behavior is noted without going much into detail. First, while the top layer (former convection zone) still has a superadiabatic temperature profile, no convective state develops in these non-conductive simulations because superadiabaticity turns out not sufficient to overcome diffusive effects. Second, for the strongest field case ($\beta =5$) a weak magnetic instability of Rayleigh-Taylor type is observed. The developing magnetic field structure is 2D (y-independent) with the tendency of flux to be concentrated in a few tube-like features aligned along the y-direction. There is no evidence for the formation of arched magnetic flux tubes or vortex tubes as found in 3D numerical simulations of a thin magnetic slab (Matthews et al.1995; Brandenburg & Schmitt 1998). For our smoothly varying magnetic field the resulting instability seems to be too weak to generate secondary Kelvin-Helmholtz instabilities due to the shear between descending and ascending fluid elements which could drive the generation of such vortices.

	Figure 8: Series of snapshots illustrating the time evolution of vz at z=-0.2 for $\beta =500$ (top) and $\beta =5$ (bottom). Time instances are t=28.8,30.2,31.5,33.0,34.4 ( 27.4,28.7,29.9,31.1,32.4) for $\beta =500$ (5) from left to right.
Open with DEXTER

In the following fully convective simulations such magnetic buoyancy effects are superimposed and will be part of the system's dynamics. Four distinct cases with different magnetic field strengths, $\beta=\{5000,500,50,5\}$, are considered. The magnetic Prandtl number is fixed to ${Pr}_{\rm m}=1$. Note that the Taylor number is varied for each $\beta$-case which sums to a total of 12 distinct runs. The dynamical importance of the various $\beta$-cases becomes obvious when one considers the initial ratio of magnetic energy to kinetic energy of the pre-existing flow. One obtains $E_{\rm mag}/E_{\rm kin}\approx 45/\beta$ for ${Ta}=6\times 10^5$ (and similar ratios for other Ta). Ratios larger than unity let expect a strong influence on convection whereas ratios much smaller than unity a weak influence. In this sense, the models with $\beta =5$ may be regarded as strong field cases, the models with $\beta =5000$ as weak field cases and the simulations with $\beta=\{50,500\}$ lie somewhere in between. All simulations are evolved from time $t_{\rm e}=27.3$ to time $t_{\rm e}^*=54.6$.

4.2.1 Spatiotemporal evolution

Initially the MHD system is not in statistically steady state convection and a transition phase occurs. Figure 7 illustrates that transition for two the models with parameter set ( ${a}=6\times 10^5, \beta=50$) respective ( ${Ta}=6\times 10^5,\beta =5$). The time evolution of the RMS velocity is shown in different ways: i) as grey-scale (t,z)-image of $\langle(\delta\vec{v})^2\rangle^{1/2}$and ii) as line plot of $\langle\langle(\delta\vec{v})^2\rangle\rangle^{1/2}_{\rm c}$. For $\beta =50$ the transition is relatively smooth. In this case convection settles into a new statistically steady state within a time span of 10-15 time units (corresponding to 4-6 convective turnover times in the new state). During the phase of readjustment the RMS velocity decreases everywhere. The volume-averaged value falls by $\approx$30% from $\approx$0.55 to $\approx$0.38. A somewhat more strange behavior is found for the strong field case $\beta =5$. Because of rapid redistribution of magnetic flux due to magnetic buoyancy effects, the flow likewise reorganizes on a short time scale. The toroidal magnetic field with initial Gaussian z-profile develops into a rather flat configuration in the unstable region (see below) affecting overall convection. This violent dynamical relaxation, perceptible in Fig. 7 by a peak in $\langle\langle(\delta\vec{v})^2\rangle\rangle^{1/2}_{\rm c}$ at $t\approx 30$, lasts a few time units. In the later course of evolution, the RMS velocity continuously decreases. It may appear that statistical equilibrium is still not fully reached at $t_{\rm e}^*$ but should be close to it with a terminal RMS velocity of $\approx$0.2. The drastic decrease of the RMS velocity is due to strong quenching by the magnetic field aligned perpendicular to the buoyancy force driving convection. The resulting change in the spatial structure of convection is expected to depend crucially on $\beta$. This is illustrated in Fig. 8 showing an early-stage sequence of snapshots of the vertical velocity in the surface-near plane z=-0.2 for a weak ( $\beta =500$) and strong ($\beta =5$) magnetic field. Whereas for $\beta =500$ the typical network structure found for non-magnetic convection persists, it is quickly destroyed for $\beta =5$. In the latter case the curvaceous downflow lanes dissolve and are finally replaced by sheet-like filaments oriented along the y-axis.

	Figure 9: Velocity pattern at t=54.6 for the models ( ${Ta}=6\times 10^5,\beta =5000$) (upper left), ( ${Ta}=6\times 10^5,\beta =50$) (upper right), ( ${Ta}=0,\beta =5$) (lower left) and ( ${Ta}=6\times 10^5,\beta =5$) (lower right). The grey-scale image shows vzand the arrows represent the 3D velocity field.
Open with DEXTER

Figure 9 presents in more detail the resulting velocity field structure at $t_{\rm e}^*$ for various models. The solution for ( ${Ta}=6\times 10^5,\beta =5000$) obviously bears close resemblance to the corresponding non-magnetic model which underlines the passive role of the magnetic field. For intermediate field strength ($\beta =50$) the granular-like cells have become significantly elongated as a consequence of the existing mean horizontal magnetic field which is primarily oriented in y-direction. Due to the Lorentz force, diverging motions in upflows are squeezed in the x-direction giving rise to their ellipsoidal-like shape. As anticipated in the corresponding plot and in accord with a reduced Reynolds number, the velocity structure turns out smoother both in the more laminar surface layer, recognizable by the broadening of the downflow lanes, and in the more turbulent interior region. The strong field solutions ($\beta =5$) are best described as quasi-laminar throughout the domain. This is consistent with a decrease in convective forcing associated with the increase of the critical Rayleigh number in the presence of the magnetic field. The final configuration shows the trend towards two-dimensionality with the direction of the mean magnetic field (y-direction) as the direction of invariance. The maximum value of Reynolds number has decreased to $\approx$500; a value which can no longer be considered representing turbulence (recall that ${Re}_{\max}\approx 1400$ for the turbulent $\beta =5000$ models). Convective motions now occur roughly in the form of columnar cells aligned along the mean magnetic field. The columnar structures are highly vertically coherent and there is little or no evidence for their corrugation with depth. In fact, no punctuated plume-like features emerge through progressive disintegration of the downflow lanes as have been found for weak-field convection. Such a disintegration with plume formation is suppressed by the strong magnetic field.

Even in strong-field convection the flow is time-dependent but shows much more regularity compared to the more turbulent cases. If rotation is absent a high level of temporal coherence exists expressed by the relative invariability of the flow pattern over many turnover times. The influence of rotation in strong field convection manifests in the details of the velocity field structure. In that situation the Coriolis force favors streaming motions parallel to the direction of the mean magnetic field caused by momentum transfer from the x-component and z-component to the y-component via linear coupling. The bulk of streaming motions occur in a couple of pipes with flow directions in +y and -y (Fig. 10). A significant amount of the total kinetic energy is stored in such motions. For example, one yields $E_{\rm kin,{\it x}}:E_{\rm kin,{\it y}}:E_{\rm kin,{\it z}}=6:12:4$ for ${Ta}=6\times 10^5$ but $E_{\rm kin,{\it x}}:E_{\rm kin,{\it y}}:E_{\rm kin,{\it z}}=23:5:26$ for Ta=0 where the numbers are scaled in such a way to allow a direct comparison. Hence, most of the kinetic energy in the rotating case resides in the y-component whereas y-motions are relatively weak in the non-rotating case. These streaming motions are the dominant flow structures in strong field convection if rotation rates are sufficiently large. They are much less pronounced or absent in weak-field convection. A remarkable anti-correlation exists between the strength of downflows and strength of streaming motions. Whereas for Ta=0 one can still distinguish between narrow, strong downflows and broader, weak upflows such a separation is much less obvious for ${Ta}=6\times 10^5$ in which both flow components have comparable strength (see Fig. 9). The diminishing asymmetry with increasing rate of rotation is because a growing fraction of kinetic energy is extracted from the faster downflows and fed into horizontal streaming motions through Coriolis force coupling. This process has strong influence on the convective transport properties (see Sect. 4.2.3).

	Figure 10: y-averaged vy-profile at $t_{\rm e}^*$ for model $({Ta}=6\times 10^5, \beta =5)$. Solid contours denote positive values whereas dashed contours negative values.
Open with DEXTER

Figure 11: Distribution of magnetic energy at $t_{\rm e}^*$ in case of rotating convection ( ${Ta}=6\times 10^5$) for initial field strength $\beta =5000$ (top) and $\beta =5$ (bottom). Isosurfaces belong to values $\vec{B}^2=\vec{B}^2_{\rm max}/6$ for $\beta =5000$ and $\vec{B}^2=\vec{B}^2_{\rm max}/2$ for $\beta =5$. The sideplots show the mean field components $\langle B_x\rangle$ (dashed), $\langle B_y\rangle$(solid) and $\langle B_z\rangle$ (dash-dotted) at $t_{\rm e}^*$. The thin solid line shows the initial magnetic field.

Open with DEXTER

The distribution of magnetic energy at $t_{\rm e}^*$ is presented in Fig. 11 for the models $({Ta}=6\times 10^5, \beta =5000)$ and $({Ta}=6\times 10^5, \beta =5)$. In the weak field case magnetic energy is concentrated in small-scale elements irregularly distributed in space. Such spatial appearance is in agreement with the picture of turbulent convection. The turbulent flow interacts with the magnetic field, wrapping it up by gradients in the velocity and stretching it when advected along with the strong downflows. Since the magnetic field thereby is amplified this mechanism may infer local dynamo action (Brandenburg et al.1996). We dispense here with a more detailed discussion of dynamo action which will be addressed in a separate paper. A significant part of magnetic energy is transported downwards deep into the stable layer which is due to the penetration effect of convective motions. As indicated by the sideplots in Fig. 11 representing the initial- and resulting mean magnetic fields a significant redistribution of magnetic flux has been taken place during the coarse of evolution. Such mean field transport is in part due to the large-scale vertical flow present in the solution and in part caused by turbulent pumping investigated in more detail in ZR. For $\beta =5000$the developed mean field has two components, $\langle B_y\rangle$ and $\langle B_x\rangle$, where latter shows multiple sign reversals with depth. The vertical component $\langle B_z\rangle =0$ because of the choice of boundary conditions. Both relevant mean field components are time-dependent because of the time-dependence of the vertical mean flow (see Sect. 4.2.2). The larger-scale structures found for $\beta =5$, on the other hand, underline the more laminar behavior of strong field convection. Magnetic energy accumulates in a few ropes which show a high degree of spatial coherence in accord with the columnar convective eddies. The mean field is primarily oriented in y-direction ( $\langle B_x\rangle \ll \langle B_y\rangle$). It has become nearly constant within the convection zone and, to a fairly good approximation, has become time-independent. The downward transport of mean magnetic field into the stable layer is dominated by diffusion rather than advective transport or turbulent pumping both of which play an important role in weak field convection. Latter processes, however, are substantially suppressed in the strong field limit (ZR).

4.2.2 Mean flow generation

There are several routes to generate mean flows in the system under study. To see this more clearly we derive the evolution equation for the mean momentum $\langle \rho \vec{v}\rangle (z,t)$ by averaging the Navier-Stokes Eq. (2) over the horizontal (periodic) directions. It follows

$$\partial_{\rm t}\langle \rho v_x\rangle -\partial_z\left\langle \rho v_xv_z -\frac{1}{\mu}B_xB_z\right\rangle +\nu_{\rm dyn}\partial_z^2\langle v_x\rangle$$ =

$$+2\Omega\cos\theta \langle \rho v_y\rangle$$ (30)

$$\partial_{\rm t}\langle \rho v_y\rangle -\partial_z\left\langle \rho v_yv_z -\frac{1}{\mu}B_yB_z\right\rangle +\nu_{\rm dyn}\partial_z^2\langle v_y\rangle$$ =

$$-2\Omega\cos\theta \langle \rho v_x\rangle -2\Omega\sin\theta \langle \rho v_z\rangle$$ (31)

$$\partial_{\rm t}\langle \rho v_z\rangle -\partial_z\left\langle \rho v_zv_z +\frac{1}{2\mu}\left(B_x^2+B_y^2-B_z^2\right)\right\rangle$$ =

$$+\frac{4}{3}\nu_{\rm dyn}\partial_z^2\langle v_z\rangle +2\Omega\sin\theta \langle \rho v_y\rangle$$ (32)

where $\nu_{\rm dyn}=\nu\rho =\rm const.$ is the dynamic viscosity coefficient. In the absence of rotation mean flows can therefore be produced through nonlinear coupling processes via the Reynolds stress terms and Maxwell stress terms. These terms are the sources of mean flow generation. Their structure depend on the flow- respective magnetic field topology as well as on externally imprinted conditions such as the spatial extent of the computational box or the choice of boundary conditions. Indeed, as shown by Rucklidge & Matthews (1993) and Brummell & Julien (1997) latter effects can lead to symmetry breaking of the flow and, consequently, mean flow production.

Mean flows are observed in the simulations of non-rotating, non-magnetic convection. This is illustrated in Fig. 12 (first row of plots) which shows the spatiotemporal evolution of the mean horizontal flows $\langle v_x\rangle (z,t)$ and $\langle v_y\rangle (z,t)$ in form of grey-scale coded (t,z)-diagrams. The mean flows become noticeable after convection sets in. Both flow components exhibit an irregular behavior in time. Their amplitudes are typically of the order of a few percent of the reference sound speed $c_{\rm S2}$ and, therefore, contain only a small fraction of the total kinetic energy. A relative strong negative $\langle v_y\rangle$ develops after $t\approx 15$ at the top boundary which persists up to the end of simulation. Note that the existence of such boundaries flows are consistent with the imposed stress-free boundary condition. Brummell et al.(1998) found no clear evidence for mean flows in their non-rotating models which is likely due to the different setup of the problem. Note in particular that the effect of penetrative convection, ignored in Brummell et al., is to induce a net vertical flow in the vicinity of the overshoot region. This is because penetrative plumes give rise to a non-vanishing $\langle \rho v_zv_z\rangle$.

	Figure 12: Space-time dependence of the mean horizontal flows $\langle v_x\rangle$ (left) and $\langle v_y\rangle$ (right). Various models with parameters $({Ta},\beta )$ as indicated in the plot are shown. The inherent time irregularity of the mean flows is demonstrated by an overlaid 1D time plot taken at z=-0.5 (right axis in the plots).
Open with DEXTER

	Figure 13: Representation of $\langle v_y\rangle (z,t)$ obtained from a long-term simulation with parameters ${Ta}=6\times 10^5$ and $\beta =5000$ accompanied by the corresponding frequency spectrum for that time series at height z=-0.5.
Open with DEXTER

	Figure 14: Persistent flows for models with ${Ta}=6\times 10^5$ and varying $\beta$. The dashed line refers to the long-term run with $\beta =5000$. The vertical bars attached at some height in the convection zone indicate the peak variation of the corresponding mean flow at that location within the interval of time averaging.
Open with DEXTER

A natural mechanism for flow symmetry breaking is tilted rotation or the presence of a large-scale magnetic field. In addition of being a potential source of mean flow generation, rotation serves to exchange momentum among the different flow components through the Coriolis force. In order not to overload the presentation by too many figures it is focused on a few typical cases to analyze the effects from rotation. Figure 12 (last three rows of plots) shows (t,z)-diagrams for the runs with parameters $({Ta},\beta)=(6\times 10^5, 5000),(6\times 10^5, 5),(6\times 10^4, 50)$. In the weak field case with ${Ta}=6\times 10^5$ the presence of rotation introduces inertial oscillations of the horizontal flows with a characteristic period of $P=\pi/\Omega\cos\theta \approx 6$ corresponding to a frequency $P^{-1}\approx 0.17$. The effect is more clearly seen in Fig. 13 (left plot) which presents the (t,z)-diagram of $\langle v_y\rangle$ for a long-term run with $({Ta},\beta )=(6\times 10^5,5000)$stopped at $t_{\rm e}^{**}=268$ which corresponds to $\approx$130 convective turnover times. The frequency spectrum of $\langle v_y\rangle$ is also shown in Fig. 13 (right plot). The distribution has indeed a strong peak at $P^{-1}\approx 0.17$ in agreement with the predicted frequency for inertial oscillations. In all cases studied with ${Ta}=6\times 10^5$ such oscillations are present in the mean horizontal flows. If rotation becomes moderate, however, (third row of plots in Fig. 12) inertial oscillations turn out no longer significant. The influence of rotation on larger-scale motions in the convection can be measured in terms of the Rossby number ${Ro}=({Ra}/{Ta\cdot Pr})^{1/2}$. One yields ${Ro}\approx 7$ for ${Ta}=6\times 10^4$ and ${Ro}\approx 2$ for ${Ta}=6\times 10^5$. The above result can therefore also be understood in terms of Rosince smaller values imply a larger influence of rotation on the mean flows.

The resulting mean horizontal flows in case of rotating convection with ${Ta}=6\times 10^5$ in a strong toroidal magnetic field ($\beta =5$) is presented in the fourth row of Fig. 12. Again one finds clear evidence for inertial oscillations. In the convection zone $\langle v_y\rangle$ also exhibits oscillations with higher frequency. These high-frequency oscillations with a period of $\approx$1 time unit are the result of the Coriolis force interaction between the y- and z-component of mean momentum. Namely, in the vertical mean flow $\langle v_z\rangle$ such oscillations are the dominant dynamical feature. They represent Alfvénicmotions in the presence of the strong mean magnetic field. The amplitude of the Alfvénicoscillations is initially large compared to the typical time-averaged value of $\langle v_z\rangle$(see below) but it are damped significantly as time progresses. Since there is no linear coupling between the x- and z-component those high-frequency oscillations are much less pronounced in $\langle v_x\rangle$. Except in the initial phase of evolution the mean flow patterns for $\langle v_x\rangle$ and $\langle v_y\rangle$, dominated by inertial motions, exhibit a rather regular behavior in time. There is a phase shift of $\approx$180$^{\circ}$ roughly constant in time between the top and bottom of the convection region. This is reminiscent of a propagating wave in vertical direction having a wavelength of $\lambda \approx 2$ and phase speed of $\approx$ $4\times 10^{-3}$ measured in units of $c_{\rm S0}$. Note, in addition, that there is a phase shift of approximately $90^{\circ}$ between $\langle v_x\rangle$ and $\langle v_y\rangle$. Consequently, the mean velocity field vector, at each height in the convection zone, rotates in horizontal planes. The sense of rotation is anticylonic i.e.opposite to the sense of $\Omega_z\hat{\vec{z}}$.

One can remove the time-dependence of the mean flows by defining "persistent'' flows through time-averaging. These persistent flows may then stand for the typical steady velocity field in a convection zone. Recall in this context that the computational box is positioned at north latitude $45^{\circ}$ so that the calculated flows represent characteristic values at that location in spherical shell convection. The persistent flow components will be denoted by $\overline{\langle v_x\rangle}$ (meridional flow; north-south direction), $\overline{\langle v_y\rangle}$ (zonal flow; west-east direction) and $\overline{\langle v_z\rangle}$ (vertical flow; radial direction) where the overbar in most cases means time-averaging over the interval $t\in[t_{\rm e}+(t_{\rm e}^*-t_{\rm e})/2,t_{\rm e}^*]$. We dispense here with showing results for the pure hydrodynamic convection simulations and lay focus on the MHD cases with ${Ta}=6\times 10^5$. Note first that, except for the special long-term run, the averaging interval covers only two inertial periods. The persistent flows defined in that way are therefore not truly persistent but are expected to be sensitive to the interval length. This effect can be seen in Fig. 14 by comparing the dashed curves (long-term run) with the curves label by "5000'' (short-term run). The persistent flows predicted from the long-term run are the more reliable ones and are somewhat smaller in magnitude than their short-term counterparts. The vertical bars at certain heights indicate the range within the mean flows vary for the duration of time averaging. These departures reflect the presence of the inertial oscillations in case of the horizontal flows and Alfvénicoscillations in case of the vertical flow. $\overline{\langle v_z\rangle}$ is the most weak (note the different scaling compared to the plots of the horizontal persistent flows). It reverses sign at the stable/unstable interface from negative in the stable layer to positive in the unstable layer. This feature is independent of the magnetic field strength. There is a second sign reversal near the top boundary with the point of reflection shifted downwards with increasing field strength. In the strongest field case no such reversal near the top occurs, however. The amplitude of $\overline{\langle v_z\rangle}$ clearly decreases with increasing field strength. Sign reversals are also observed for the horizontal persistent flows. For $\beta =\{5000,500,50\}$ the meridional flow $\overline{\langle v_x\rangle}$ is negative (retrograde i.e.in opposite sense to rotation) in the bulk of the convection zone and positive (prograde i.e.like rotation) below. Except for $\beta =5$ the zonal flow $\overline{\langle v_y\rangle}$ turns out largest in the thermal boundary layer where it is prograde but is relatively weak elsewhere. It should not be concealed that the situation becomes somewhat more complicated in case of moderate rotation with ${Ta}=6\times 10^4$. For example, in weak-field convection ( $\beta =5000,500$) the surface-near zonal flow becomes retrograde instead of prograde, as has been observed for ${Ta}=6\times 10^5$. For stronger magnetic fields, however, it becomes prograde again. Altogether, we find a rather complex dependence of the spatial structure of persistent flows on the strength of rotation and magnetic field strength.

	Figure 15: Scatter plots of the kinetic energy flux $F_{\rm kin}$ (left), enthalpy flux $F_{\rm e}$ (middle) and convected energy flux $F_{\rm c}$ (right) at depth z=-0.5 and time $t_{\rm e}^*$. Upper plots are for Ta=0 and lower plots for ${Ta}=6\times 10^5$. The same scaling as in Fig. 5 is used.
Open with DEXTER

4.2.3 Magnetoconvective energy transport

One remarkable feature observed in the models of non-magnetic, turbulent convection is that downflows carry roughly as much kinetic energy downwards as thermal energy upwards. In the presence of a strong magnetic field this property does not persist. To work out the differences the following discussion focuses on the simulations with $\beta =5$. A survey of the results for the remaining $\beta$-parameters indicate a behavior which lies somewhere in between the extreme cases of non-magnetic convection examined in Sect. 4.1 and strong field convection. Figure 15 presents scatter plots of the various fluxes. The same axes range and scaling as in Fig. 5 is used to enable a direct comparison with those results. The fluxes are calculated in the horizontal plane z=-0.5 at time $t_{\rm e}^*$. Note first that the sample of dots in phase-space show a high degree of correlation which is an expression of the laminar character of convection enforced by the strong magnetic field. Both the kinetic energy flux and enthalpy flux has decreased significantly although individual velocities on the downflow side still reach magnitudes of the order of unity in the non-rotating case. For ${Ta}=6\times 10^5$, in contrast, peak velocities are much smaller and are roughly equal in downflows and upflows. In both cases the contribution of the downflows to the convected energy flux is non-negligible and definitely positive in contrast to the behavior observed in turbulent convection. A possible explanation for that difference stems from the idea that the enthalpy flux is likely less affected by magnetic quenching than is the kinetic energy flux being proportional to the third power of velocity. This would imply that the oppositely directed fluxes of kinetic energy and enthalpy no longer compensate on the downflow side in the diagram as observed. In non-magnetic convection the effect of approximate flux cancellation in downflows has been attributed to the intermittent structure of those flow features. In strong field convection intermittency is significantly reduced in the rotating case but not for Ta=0 where fast, downward-directed sheet-like flows are still present. In both cases, however, flux cancellation does not occur. One is therefore willing to argue that the differences do not arise because of intermittency but result from fundamentally distinct transport properties inherent to laminar convection on the one hand and turbulent convection on the other hand.

5 Conclusions

The spatiotemporal evolution, mean flow generation and mixing properties of rotating, compressible convection have been studied in a local Cartesian subset of a spherical shell. To account for overshooting effects a 2-layer polytropic stratification has been constructed consisting of a convectively unstable region on top of a stable region. Starting with a fully developed state of non-magnetic, turbulent convection the system's reorganization after imposing a non-homogeneous, toroidal magnetic field localized at the stable/unstable interface has been examined.

The initiated transition to magnetoconvection turns out to be a complicated process. It has been demonstrated that the resulting spatial structure, generated mean flows and mixing properties of convection crucially depend on the problem parameters $({Ta},\beta )$. In weak field convection ( $\beta =5000$) the structure is essentially like in non-magnetic convection with the magnetic field treated as a passive ingredient advected with the flow. Also, the properties of convective energy transport are comparable. The solution is characterized by a quasi-laminar surface layer disguising the more turbulent interior of the convection zone. The flow is made up of narrow downflows surrounded by broader upflows; a typical behavior for compressible convection in a stratified medium. Consistent with the turbulent nature of convection the magnetic energy is stored in small-scale structures to some part located in the stable layer beneath the convection zone where it has been carried by overshooting motions. The mean magnetic field is found to be time-dependent even at the latest stages of evolution and has both a y-component, $\langle B_y\rangle$, and x-component, $\langle B_x\rangle$, where latter involves multiple sign reversals when plotted as a function of z. If the field strength is increased the changes in the topological structure of the flow are towards two-dimensionality and laminar-type. More specifically, for the strongest field case considered ($\beta =5$) the resulting solution is best described by a laminar flow throughout the domain with convective motions occurring roughly in the form of cylindrical rolls aligned with the mean magnetic field which is primarily oriented in y-direction in that case. The mean magnetic field has become spatially flat and quasi-stationary in the convection zone which is due to magnetic buoyancy. Magnetic quenching effects strongly affect the mixing properties of convection involving a significant reduction of the kinetic energy flux and enthalpy flux as well as a suppression of the downward transport of magnetic energy by overshooting plumes. The laminar character of strong field convection infers a significant net positive contribution of the downflows to the convected (kinetic and enthalpy) energy flux unlike the situation of turbulent convection where it has been found that the downflows carry as much kinetic energy downwards as enthalpy upwards and, hence, cancel approximately.

One major finding of the study is the existence of mean flows. Their characteristic properties depend in a complex fashion on the system parameters Ta and $\beta$. Even for $({Ta},\beta )=(0,\infty)$ mean flows are present in all coordinate directions although mechanisms for flow symmetry breaking are not apparent. We believe that, at least for the vertical mean flow, this is in part a result of permitting overshooting which is clearly different from the situation of an impermeable boundary. The finite horizontal extent of the computational box may also lead to symmetry breaking and, hence, to a (artificial) contribution to the Reynolds stresses being the source of mean flow production in the mean-field equations. In the presence of strong rotation the mean horizontal flows show inertial oscillations with period $P=\pi/\Omega_z$. This feature turns out independent from $\beta$ although the stronger the magnetic field is the more clearly defined are the temporal variations. One another obvious influence of a strong toroidal field is to generate Alfvénicoscillations in the vertical mean flow. These oscillations have a higher frequency than the inertial motions by a factor of $\approx$6. The Alfvénicoscillations are also seen in the horizontal flows (in particular in $\langle v_y\rangle$) but with much weaker intensity where they are excited by Coriolis force coupling. The time-averaged mean flows ("persistent'' flows) are an estimate for local steady motions in a star's convection zone. The vertical flow shows a systematic behavior with a variation of Ta and $\beta$. It is negative in the stable layer and positive in the convection zone except near the top boundary where it becomes negative again. Its strength is almost one order of magnitude smaller than the horizontal flows and decreases with increasing magnetic field. The dependence of the meridional flow (north-south direction, x) and zonal flow (west-east direction, y) on these parameter, however, is more delicate. For example, in the case of strong rotation the meridional flow typically changes its orientation from prograde in the vicinity of the stable/unstable interface to retrograde in the bulk of convection zone. However, a strong magnetic field forces the meridional flow prograde in the top half of the convection zone. In the regime of moderate rotation no systematic behavior is observed. One further aspect which may need more discussion is the possible influence of the magnetic field boundary condition at the top. In general, we expect the detailed structure of the mean flows to depend on the type of boundary condition. It is neither a perfect conductor condition nor is it of pseudo-vacuum type. Yet, it is the belief of the author that the adopted boundary condition is most suitable among the local boundary conditions for the current problem. To further strengthen confidence in the obtained results it is planned in a future work to relax the top boundary condition by adding a fiducial stable layer.

Acknowledgements

The computations were performed on workstations at the Astrophysikalisches Institut Potsdam.

References

 
• Ames, W. F. 1977, in Numerical methods for partial differential equations (Academic Press, Inc.) In the text  
• Brandenburg, A., Nordlund, A., Pulkkinen, P., Stein, R. F., & Tuominen, I. 1990, A&A, 232, 277 In the text NASA ADS  
• Brandenburg, A., Jennings, R. L., Nordlund, A., et al. 1996, J. Fluid Mech., 306, 325 In the text  
• Brandenburg, A., & Schmitt, D. 1998, A&A, 338, L55 In the text NASA ADS  
• Brummell, N. H., Hurlburt, N. E., & Toomre, J. 1996, ApJ, 473, 494 In the text NASA ADS  
• Brummell, N. H., Hurlburt, N. E., & Toomre, J. 1998, ApJ, 493, 955 In the text NASA ADS  
• Cattaneo, F., Brummell, N. H., Toomre, J., Malagoli, A., & Hurlburt, N .E. 1991, ApJ, 370, 282 In the text NASA ADS  
• Chandrasekhar, S. 1961, in Hydrodynamic and Hydromagnetic Stability (Oxford: Clarendon) In the text  
• Glatzmaier, G. A., & Toomre, J. 1995, in GONG 94: Helio- and Asteroseismology from the Earth and Space, ed. R. K. Ulrich, E. J. Rhodes Jr., & W. Däppen (San Francisco: ASP), ASP Conf. Ser., 76, 200 In the text  
• Gough, D. O., Moore, D. R., Spiegel, E. A., & Weiss, N. O. 1976, ApJ, 206, 536 In the text NASA ADS  
• Hurlburt, N. E., Toomre, J., Massaguer, J. M., & Zahn, J.-P. 1994, ApJ, 421, 245 In the text NASA ADS  
• Hurlburt, N. E., Matthews, P. C., & Proctor, M. R. E. 1996, ApJ, 457, 933 In the text NASA ADS  
• Julien, K., Legg, S., McWilliams, J., & Werne, J. 1996, J. Fluid Mech., 322, 243 In the text  
• Malagoli, A., Cattaneo, F., & Brummell, N. H. 1990, ApJ, 361, L33 In the text NASA ADS  
• Matthews, P. C., Hughes, D. W., & Proctor, M. R. E. 1995, ApJ, 448, 938 In the text NASA ADS  
• Nordlund, A., Brandenburg, A., Jennings, R. L., et al. 1992, ApJ, 392, 647 In the text NASA ADS  
• Spiegel, E. A. 1965, ApJ, 141, 1068 In the text NASA ADS  
• Weiss, N. O., Brownjohn, D. P., Matthews, P. C., & Proctor, M. R. E. 1996, MNRAS, 283, 1153 In the text NASA ADS  
• Woodward, P. R., & Colella, P. 1984, J. Comp. Phys., 54, 115 In the text  
• Ziegler, U. 1998, Comp. Phys. Comm., 109, 111 In the text  
• Ziegler, U. 1999, Comp. Phys, Comm., 116, 65 In the text  
• Ziegler, U., & Rüdiger, G. 2001, A&A, in prep. In the text