Fig. 21.1. Solar angular velocity in a meridional cut as determined by helioseis-

mology. (Reprinted with permission from Thompson et al., Science, 296, 1300.

Copyright 1996 American Association for the Advancement of Science.)

a broad region in radius near the base of the convection zone. The time scale

of the convection is about 30 days, with characteristic velocities of a few

hundred m s’1 and characteristic temperature ¬‚uctuations of a few Kelvin.

Since the convection and rotation time scales are very similar (leading to

a convective Rossby number of order unity), the convection is signi¬cantly

in¬‚uenced by rotation”this leads to the observed di¬erential rotation, with

the poles rotating signi¬cantly more slowly than the equator (see Fig. 21.1).

In addition, since the density in the convection zone varies by over ten orders

of magnitude, the convection is strongly in¬‚uenced by strati¬cation e¬ects.

Since the solar plasma is a very good electrical conductor (leading to small

magnetic di¬usivities and high magnetic Reynolds numbers), magnetic ¬elds

are yet another factor which in¬‚uences the dynamics of the solar interior;

the interaction of ¬‚uid ¬‚ow and magnetic ¬elds leads to the 22-year cycle

of magnetic ¬elds and sunspot activity, although the exact mechanism is as

yet unclear.

The combined e¬ects of rotation, strati¬cation and magnetic ¬elds in a

high Reynolds number (and therefore highly turbulent) ¬‚ow leads to a sys-

tem with a rich spectrum of behaviour that we can only hope to understand

through detailed numerical simulations. Since they are computationally in-

tensive, such simulations have only become possible in the last few decades

with advances in high-performance computing. In spite of these advances,

Numerical simulations of the solar convection zone 317

certain approximations are necessary in carrying out simulations of solar

convection”the most important of these is to ¬lter out fast sound waves

that play a small part in the overall dynamics of the convection zone but

could unreasonably limit the computational time step. In both the codes

described here, this is done using the anelastic approximation (Gough 1969).

Another fundamental simpli¬cation is to neglect magnetic ¬elds”the valid-

ity of this approximation depends on magnetic ¬elds not being too strong.

A ¬nal simpli¬cation made in the simulations described here is to limit the

range of densities modelled, e¬ectively cutting o¬ the simulations near the

solar surface, where the density scale height becomes small; this removes the

need to use multiple grids to model the di¬erent scales of convection and

simpli¬es the problem.

The choice of numerical integration scheme in modelling solar convection

is paramount. Interestingly, thermodynamic considerations provide some

insight. Radiative heating near the base of the convection zone and com-

pensating cooling near the solar surface act as a source of negative entropy,

which in equilibrium must be balanced by the action of irreversible processes.

Two irreversible processes dominate”¬rstly, the transfer of heat from one

¬‚uid element to another by radiative di¬usion, and secondly, the dissipa-

tion of kinetic energy by the action of viscosity. The relative importance

of the two is unclear, although since the Prandtl number is low (the ra-

diative di¬usivity being much higher than the kinematic viscosity) the ¬rst

is probably dominant. Unfortunately, both operate at scales far too small

to be resolved computationally. For this reason, irreversible entropy pro-

duction must be accounted for by arti¬cial dissipation, either in the form of

enhanced values of the radiative di¬usivity and viscosity, or a dissipative nu-

merical scheme. Studies employing non-dissipative numerical schemes (such

as pseudo-spectral and CTS”Centred in Time and Space”methods) have

often justi¬ed the arti¬cial dissipation added as in some way modelling the

e¬ect of small eddies on the large-scale ¬‚ow (e.g. Gilman 1977; Glatzmaier

1984). Since the di¬usion and viscosity took the form of Laplacian operators

with constant coe¬cients, these studies may be classi¬ed as DNS simula-

tions. The next section summarizes some of the results from Elliott at al.

(2000) and Miesch et al. (2000), which fall into this category.

The problem with DNS techniques is that large thermal di¬usivities must

be introduced everywhere, even where the ¬‚ow is laminar. A better approach

is to adopt the LES technique, where the ¬‚ow ¬eld is split into resolved

and unresolved scales with the latter acting back on the former through

the divergence of the SGS (Sub-Grid Scale) stress tensor. In general, the

components of this stress tensor are expressed in terms of the resolved ¬‚ow

Elliott

318

¬eld”this is known as the SGS model. The most well known SGS model

is that postulated by Smagorinsky (1963), in which the SGS stress tensor is

proportional to the local rate of strain of the resolved ¬‚ow.

Historically, SGS models have tended to be judged more on their ability

to suppress false computational oscillations than on their representation of

what is actually happening on sub-grid scales. Once one gives up the notion

that SGS models accurately represent the behaviour of unresolved eddies,

there are simpler, more e¬ective means of achieving the suppression of false

computational oscillations. In particular, there is a class of ¬nite di¬erence

methods “ Non-oscillatory Forward in Time (NFT) “ that have the remark-

able property of representing LES/VLES without recourse to any explicit

SGS model. Such methods are inherently dissipative and provide an irre-

versible entropy source, as they must given the constraints mentioned above.

The use of this class of numerical methods in geophysical ¬‚uid dynamics has

been justi¬ed by simulations of the Earth™s planetary boundary layer (Mar-

golin et al. 1999) that compare well with observations. This paper presents

¬rst results for solar convection obtained from a code based on NFT meth-

ods, in particular the MPDATA (Multidimensional Positive-De¬nite Advec-

tion and Transport Algorithm) scheme (Smolarkiewicz & Margolin 1998).

This is not the ¬rst time that comparisons between di¬erent methods

of simulating solar convection have been made. For example, Cattaneo et

al. (1991) compared three-dimensional Cartesian convection simulated with

a pseudo-spectral code and with a code implementing the PPM method

(Piecewise Parabolic Method, Collela & Woodward 1984); the latter is an-

other method from the NFT class that can be run stably with no explicit

di¬usion or viscosity. Although detailed comparisons of results were not

presented in this study, a description of the philosophy of NFT methods

was included. Much confusion has been introduced into the subject by the

inappropriate use of the word “inviscid” to describe NFT methods.

This paper is divided into three further sections. The next covers re-

sults from a pseudo-spectral DNS model that uses constant di¬usivities and

viscosities. Next, results from a VLES simulation of solar convection are

described. The ¬nal section highlights the conclusions from this work.

21.2 DNS results

This section describes numerical simulations carried out within the DNS

framework using a code (see Clune et al. 1999) that expands the prognostic

variables (entropy, pressure and velocity) in spherical harmonics in the hor-

izontal, and in Chebyshev polynomials in the vertical. Spherical harmonic

Numerical simulations of the solar convection zone 319

expansions have been widely used in meteorological applications and have

the advantage of providing uniform resolution over the sphere (at least when

triangular truncation is used) ” this avoids the so-called pole problem. Non-

linear terms are calculated in physical space (the transform method), using

a grid with su¬cient resolution to avoid quadratic aliasing. Apart from the

uniform resolution, spherical harmonic basis functions also o¬er the advan-

tage of yielding a very easy solution of the elliptic pressure equation, since

the linear terms decouple.

We now describe the equation set used. The prognostic variables are

velocity and entropy. The continuity equation in the anelastic approximation

(which applies to both models described in this paper) is

∇ · (ρ0 v) = 0 , (21.1)

where ρ0 is the density in the reference state (T0 and θ0 are the correspond-

ing temperature and potential temperature respectively) and v is the ¬‚uid

velocity. The momentum equation is

‚v 1 ρ 1

+ (v · ∇) v = ’ ∇p + g + 2v — „¦ ’ ∇ · D , (21.2)

‚t ρ0 ρ0 ρ0

where p and ρ are the departures of the pressure and density from the

reference state values, g is the acceleration due to gravity, „¦ is the rotation

vector, and D is the viscous stress tensor, given in terms of the viscosity ν

and rate-of-strain tensor e by

1

D ij = ’2ρ0 ν eij ’ (∇ · v) δij . (21.3)

3

The entropy equation is given by

‚S 1

(∇ · v)2 ,

+ (v · ∇) S = µ + ∇ · κρ0 T0 ∇S + 2νρ0 eij eij ’

ρ0 T0

‚t 3

(21.4)

where S is the ¬‚uctuation of entropy away from the reference state value, κ