last term on the right-hand side of this equation is the viscous heating term.

The domain used for the simulations is a spherical shell with an inner

radius of 5 — 105 km and an outer radius of 6.9 — 105 km. The density

contrast from the bottom of the domain to the top is about 30, which,

although substantially less than the actual contrast across the convection

zone, allows for signi¬cant strati¬cation e¬ects to be felt. Thermal forcing

is modelled by assuming that the mean strati¬cation in the bulk of the

Elliott

320

Fig. 21.2. The divergence µ of the radiative ¬‚ux (dotted line: actual, dashed line:

that assumed in the model), relative to its value µ0 at the bottom of the simulation

domain, and the density (solid line); z is the fractional height from the bottom.

convection zone is always nearly adiabatic, allowing the divergence of the

radiative ¬‚ux, µ, to be assumed to be ¬xed ” in the simulations described

in this paper, the following form is used for µ:

ζ1 e’ζ1 z ζ2 e’ζ2 (1’z)

’

µ=± , (21.5)

1 ’ e’ζ1 1 ’ e’ζ2

where ± is an appropriate constant of proportionality, z is the non-

dimensional height within the domain, and ζ1 and ζ2 are non-dimensional

constants, which take the values 5.5 and 13.0 respectively. Figure 21.2 com-

pares the actual radiative ¬‚ux divergence in the convection zone with the

formula used, and also shows the density pro¬le assumed. Near the surface,

the assumption of a near-adiabatic strati¬cation breaks down, and there is

an extremely narrow region of strong radiative cooling ” this is arti¬cially

broadened (by making the scale of the region comparable to the local density

scale height) to enable it to be resolved on the computational grid used.

Consistent with the DNS framework, the Navier-Stokes equations are

solved with entropy di¬usion and viscous terms having coe¬cients that de-

pend only on radius. Tests have been carried out with several di¬erent

dependencies of these coe¬cients on radius”the e¬ects on the results are

not too substantial. Su¬ciently large values of entropy di¬usion and vis-

Numerical simulations of the solar convection zone 321

Fig. 21.3. The vertical velocity on three horizontal surfaces (from left to right) near

the top, middle, and bottom of the shell, in a DNS simulation of solar convection.

Dark shades denote down¬‚ows and light shades up¬‚ows. The dashed line marks

the equator.

cosity are required to produce stable, well-behaved solutions; typical values

used lead to Reynolds numbers in the region of 50.

In a linear analysis, the most unstable modes of the system are convective

rolls aligned with the rotation axis. These are also seen in the non-linear

solutions, albeit with considerable distortion from non-linear e¬ects. Figure

21.3 shows typical snapshots of the vertical velocity on three horizontal

levels of a solution from the DNS code. The convective rolls, or “banana

cells”, can clearly be seen in the left and middle panels of this ¬gure (which

correspond to levels near the surface, and near the centre of the convection

zone, respectively), with rotationally aligned lanes of down¬‚owing material

interspersed by up¬‚owing material.

This simulation used 64 radial grid points and spherical harmonics up

to degree 85. The Reynolds number for the largest eddies was about 30.

Higher resolution would enable higher Reynolds numbers to be attained;

however, the computational e¬ort to increase the resolution is very high,

since the number of ¬‚oating point operations scales with the ¬fth power of

resolution (since the spherical harmonic transforms scale with the cube of

the horizontal resolution, and the number of time steps also scales linearly

with the resolution).

Turning to other aspects of the results, the time-averaged di¬erential ro-

tation for the same simulation is shown in Fig. 21.4. Comparing with

Fig. 21.1, which shows the solar di¬erential rotation as deduced by helio-

seismology, some similarities are clear. First, the simulation gives a fast

equator and slow poles. Secondly, the contrast in angular velocity is similar

to that seen in the sun. Unfortunately, there are signi¬cant di¬erences in

Elliott

322

Fig. 21.4. The time-averaged angular velocity from the same simulation as shown

in Fig. 21.3. The scale on the right is in nHz.

other respects. Firstly, most of the angular velocity contrast in the sim-

ulation is seen near the equator, whereas in the sun there is a continuous

change in angular velocity from the equator to the poles. Also, the sun shows

considerable tilting of the lines of constant angular velocity relative to the

rotation axis, whereas the simulation shows very little. Other simulations

(Gilman 1977; Glatzmaier 1984) have also shown little of such tilting, with

angular velocity contours nearly parallel to the rotation axis. The reasons

for these discrepancies are currently unclear, although what does seem to be

the case is that Reynolds stresses alone (in this case correlations of the radial

and latitudinal velocity components) are not strong enough, and typically

of the wrong sign, to produce the observed tilting”buoyancy forces almost

certainly play a part, implying that the poles are warm with respect to the

equator.

21.3 VLES results

The simulations described so far rely on large, ¬xed viscosities and thermal

di¬usivities to obtain a stable solution. An alternative option is to use an

SGS model within an LES framework. This avoids the use of a large amount

of dissipation where it is unnecessary, the only drawback being the di¬culty

Numerical simulations of the solar convection zone 323

of ¬nding an appropriate SGS model (especially in the presence of rota-

tion and strati¬cation). Another way of avoiding unnecessary dissipation,

without this drawback, is to use an appropriate numerical advection scheme

that is itself responsible for the dissipation. Such advection schemes belong

to the general class of NFT methods for which ensuring linear stability (in

particular the familiar Courant-Friedrichs-Lewy condition) also assures non-

linear stability. We refer to numerical simulations of ¬‚uid ¬‚ow employing

such advection schemes as VLES, to distinguish them from LES.

This section describes VLES simulations of solar convection carried out us-

ing a code described in Smolarkiewicz et al. (2001). The code is an anelastic,

grid-based code, with options for carrying out advection either by means of

a semi-Lagrangian scheme or an Eulerian scheme (MPDATA, Smolarkiewicz

& Margolin 1998). Because NFT methods are inherently two time-level (i.e.

the current and the next time step), accurate (e.g. second-order) time in-

tegration couples the non-linear terms into the implicit solve, necessitating

the solution of a large non-symmetric linear system representing the complex

non self-adjoint elliptic pressure equation. For deep ¬‚uids, such as the sun™s

convection zone, such a problem can be solved using standard Krylov sub-

space methods (although note that in the case of the Earth™s atmosphere,

the resulting elliptic equations are extremely sti¬, and preconditioning is

necessary) ” that approach is adopted in this code.

In terms of the equation set, potential temperature is used instead of en-

tropy, leading to the following momentum and potential temperature equa-

tions:

‚v p θ 1

+ (v · ∇) v = ’∇ + 2v — „¦ ’ ∇ · (νρ0 ∇v) ,

+g (21.6)

‚t ρ0 θ0 ρ0

‚θ θ0 1 1

’ ∇ · kρ0 ∇θ + 2νρ0 eij eij ’ (∇ · v)2 ,

+ (v · ∇) θ =

‚t ρ0 KT0 ρ0 3

(21.7)

where θ is the ¬‚uctuation of potential temperature away from the reference

state, K is the thermal capacity of the ¬‚uid, and k is the thermal di¬usivity.

Note the somewhat di¬erent form of the viscous stress in the momentum

equation. This term is calculated by means of scalar Laplacians on each

component, rather than taking the full vector Laplacian ” this means that

certain parts of the stress (in particular those arising from the curvature of

the coordinate system) are neglected.

Coordinate singularities, such as that present at the poles in the spherical