text:
# Solar Multi-Scale Convection and Rotation Gradients Studied in Shallow Spherical Shells

## 1 Introduction

Helioseismology has revealed that the differential rotation profile observed at the solar photosphere roughly extends throughout the bulk of the convection zone (Thompson et al., 1996; Schou et al., 1998). From about 0.75 $`R`$ to 0.95 $`R`$ (where $`R`$ is the solar radius), the angular velocity $`\mathrm{\Omega }`$ has a small radial gradient, particularly at mid latitudes, as seen in Figure 1. In contrast, regions of strong radial shear are evident near both the bottom and top of the convection zone (shown shaded in Fig. 1), and these shear layers are believed to play important dynamical roles within the solar convection zone. While the tachocline region at the base of the convection zone has commanded much recent attention (as it is likely the seat of the global solar dynamo and the associated 22-year magnetic activity cycle), the dynamics within the near-surface shear layer, extending from 0.95 $`R`$ to 1.00 $`R`$, are also likely to have additional dynamical consequences that affect the appearance and evolution of flows and magnetic structures visible at the surface. Such dynamics are presently not well understood, but are now becoming accessible to study through direct numerical simulations that capture the important effects of sphericity, compressibility, and rotation.

Several questions arise about the near-surface shear layer. First, what dynamical mechanisms within the coupling of turbulent convection with rotation leads to such a boundary layer involving negative radial gradients of $`\mathrm{\Omega }`$? In contrast, in the bulk of the convection zone the gradients are much smaller and positive. Second, why does this boundary layer have a depth of about 5% in solar radius? Third, does the presence of such a strong radial shear zone play a significant role in the complex large-scale meandering flows and reversing meridional circulations that have been shown by helioseismology to coexist with the intense smaller-scale convection in the upper reaches of the convection zone? These three issues motivate our studies as we seek to resolve both supergranulation and global-scale responses in our simulations of convection in rotating spherical shells. Computational constraints have encouraged us to begin by studying thin shells of such turbulent convection, encompassing at this stage only the upper portion of the solar convective envelope. A preliminary account of such modeling is presented in DeRosa & Toomre (2001). We shall here show that the resulting multi-scale convection is able to redistribute angular momentum so as to yield radial gradients in $`\mathrm{\Omega }`$ that are largely in accord with the helioseismic findings. These first steps are important in defining more complex simulations to be undertaken within deep convection shells that capture much of the depth range of the solar convection zone.

The stratification within the near-surface shear layer serves to drive vigorous motions possessing a wide range of spatial and temporal scales, visible at the surface as the convective patterns of supergranulation, mesogranulation, and granulation (Spruit et al., 1990). Many aspects of such small-scale but intensely turbulent convection, influenced by radiative transfer effects and complex equations of state and opacities, have been studied through three-dimensional simulations within localized planar domains positioned near the solar surface (e.g. Stein & Nordlund 1998, 2000, 2001). The driving in such convection is enhanced by the latent heat released within the ionization zones of helium and hydrogen that are present in the near-surface shear layer (e.g. Rast et al. 1993). These small-scale turbulent convective motions are likely to facilitate the transport of angular momentum along both radial and latitudinal velocity gradients within the shear layer, and thus may be able to affect the dynamics within the convection zone on a more global scale. In particular, the horizontal extent and overturning time of supergranular flows suggest that such convection will be at least weakly influenced by rotational effects, which can yield Reynolds stresses of significance in transporting angular momentum within the layer. The coupling of turbulent compressible convection with rotation has also been studied extensively in localized $`f`$-plane domains (e.g. Brummell et al. 1995, 1996, 1998, 2002; Brandenburg et al. 1996; Chan 2001) using perfect gases, revealing that the presence of coherent structures associated with strong downflow plumes and networks play a crucial role in the redistribution of angular momentum. Such studies are complemented by a broad range of other simulations of compressible convection that exhibit intrinsic asymmetries between upflows and downflows, of complex vorticity structures that influence the transport of heat, momentum and magnetic fields, and of rich time-dependence involving a broad range of time scales (e.g. Cattaneo et al. 1991; Porter & Woodward 1994, 2000; Saikia et al. 2000; Robinson & Chan 2001; Tobias et al. 2001).

Velocity features larger than the spatial scale of solar supergranulation are also in evidence in the near-surface shear layer. Bands of slightly faster rotation, or torsional oscillations, that gradually propagate toward the equator as the magnetic activity cycle advances are detected in Doppler measurements of the surface (LaBonte & Howard, 1982; Hathaway et al., 1996; Ulrich, 1998). They are also seen in global $`f`$\- and $`p`$-mode helioseismic studies (e.g. Kosovichev & Schou 1997; Schou 1999; Howe et al. 2000b; Vorontsov et al. 2002), and are present over at least the outer 8% in radius. Even more complex flows within the near-surface shear layer, now called solar subsurface weather (Toomre, 2002), are revealed by local-domain helioseismic techniques such as ring-diagram analyses (e.g. Hill 1988; Haber et al. 1998, 2000; Basu et al. 1999) and time-distance methods (e.g. Duvall et al. 1993; Giles et al. 1997; Duvall & Gizon 2000; Chou & Dai 2001). Mappings of subsurface flow fields over a range of depths reveal evolving large-scale horizontal flows that are somewhat reminiscent of jet streams, meridional circulations that may possess multi-celled structures in one hemisphere and not in the other, and distinctive flow deflection in the vicinity of active complexes (Haber et al., 2002). Although the flow speeds in the meridional circulations are only of order 20 m s<sup>-1</sup>, they may be quite effective in redistributing angular momentum in latitude, thereby coupling widely separated regions within the near-surface shear layer and possibly having a role in its existence.

Photospheric magnetic field observations reveal structured concentrations that also possess a wide range of spatial scales, including active regions, sunspots, pores, and emergent flux elements. On the smallest observable scales, concentrations of filamentary magnetic flux elements are found to be laterally advected by larger-scale surface flow patterns. Outflows associated with the convective patterns of supergranulation and granulation in particular are observed to readily advect such small-scale flux toward intercellular lanes and concentrate these fields on scales small enough for dissipation to occur (Schrijver et al., 1997; Berger et al., 1998).

Figure 1 indicates that the radial gradient of angular velocity $`\mathrm{\Omega }`$ is largely negative at low and mid latitudes within the near-surface shear layer, such that the rotation rate decreases by about 2–4% as one moves outward across the layer. Such a radial gradient in $`\mathrm{\Omega }`$ may be interpreted as a tendency for fluid parcels in the convection zone to partially conserve their angular momentum as they move toward or away from the axis of rotation. This idea was originally suggested by Foukal & Jokipii (1975), and may explain why larger-scale magnetic tracers at the surface have a faster rotation rate relative to the photospheric plasma, assuming that these magnetic features are anchored at a radius below the photosphere where the rotation rate is faster. Numerical simulations of Boussinesq fluids confined to thin shells (Gilman & Foukal, 1979) showed that angular momentum is roughly conserved along radial lines, and small-scale convective motions are able to transport angular momentum inward, thereby maintaining the negative radial gradient of rotation rate with radius for such an incompressible fluid. Whether the same is true for a compressible fluid is one of the topics addressed in this paper.

The maintenance of the relatively small radial gradients of $`\mathrm{\Omega }`$ throughout the bulk of the convection zone must be a direct consequence of the interaction of rotation with the turbulent fluid motions that exist within the solar convection zone. Recent three-dimensional numerical simulations of such deep convection within rotating spherical shells (e.g. Miesch et al. 2000; Elliott et al. 2000; Brun & Toomre 2001, 2002) indicate that for a range of parameter values, solar-like differential rotation profiles can be established, even with viscous and thermal diffusivities in a regime far removed from their solar values. Many of these simulations possess about a 30% contrast in angular velocity $`\mathrm{\Omega }`$ between the equator and high latitudes and have small radial gradients of $`\mathrm{\Omega }`$ in the mid-latitude regions, features that roughly match the helioseismic determinations of the interior rotation profile within the bulk of the convection zone. Analyses of the angular momentum transport within these simulations indicate that the fast equatorial rotation relative to the higher latitudes is primarily maintained by a complex interplay between global meridional circulation and Reynolds stresses achieved within the domains, both of which contribute to the equatorward transport of angular momentum with latitude.

The radial velocity planforms within the more laminar convection zone simulations take the form of rotationally aligned banana cell structures, with the downflowing fluid lanes extending throughout most of the radial extent of the domain. As the level of turbulence is increased, these organized patterns become less prominent, giving way to a network of narrower downflow lanes that form plume-like structures at the interstices in the network. Such plumes tend to possess significant vortical motion and span the entire domain in radius. The influence of rotation on these vertical plumes preferentially tilts these structures such that they are partially aligned with the axis of rotation, which in turn creates the Reynolds stresses that facilitate angular momentum transport within the domain. These deep-shell simulations of the convection zone typically place the upper boundary at about 0.96 $`R`$, and thus do not capture the smaller scales of convection that exist closer to the surface. Consequently, the convection in even the most turbulent of these simulations involves overall pattern scales of order $`20^{}`$$`30^{}`$, or several hundred Mm, although the sheets and plumes associated with the downflow network are individually narrower and more concentrated.

To understand more clearly some of the physical processes occurring within the near-surface shear layer of the solar convection zone, we have constructed numerical simulations of compressible fluids within thin spherical shells that extend up to 0.98 $`R`$, encompassing solely the near-surface shear layer region. Continual advances in supercomputing technology now permit three-dimensional compressible fluid simulations that explicitly resolve spatial and temporal scales spanning several orders of magnitude. As a result, we are able for the first time to employ direct numerical simulations to investigate the effects of supergranular-sized convection on the more global dynamics within the near-surface shear layer of the sun.

In formulating our simulations, we have adopted the viewpoint that the latitudinal variation of the angular velocity in the sun, with the equatorial regions rotating more rapidly than the poles as in Figure 1, is established and maintained within the bulk of the convection zone somewhere below the lower boundary of our thin shell models. We have thus imposed a solar-like differential rotation profile as a no-slip lower boundary in three of the four simulations presented here, in order to capture some of the dynamical effects related to such an angular velocity structure. In so doing, we are implicitly assuming that the global differential rotation profile is not substantially affected by the convection within our thin shells, even though in the sun such shearing layers could have subtle effects on these dynamics.

Our primary focus in this paper will be to investigate the angular momentum transport achieved by multi-scale convection, involving both global and supergranular scales, within shearing layers analogous to the near-surface shear layer of the sun. We shall consider radial stratifications that resemble ones deduced from stellar structure models over the depth range being studied, though the physics of the gases is highly simplified. After briefly discussing in §2 the governing equations and numerical approach used in solving them, we review the parameters used to initialize our thin-shell simulations in §3. We next examine in §4 the multi-scale convective velocity patterns of the mature solutions, and discuss the meridional circulation, time-evolution and angular momentum balance achieved within the thin shell domains. Lastly, we discuss the connection between these simulations and the near-surface shear layer of the solar convection zone, and present possible directions for future research in §5.

## 2 Governing Equations and Numerical Approach

The convection simulations described here are carried out using the anelastic spherical harmonic (ASH) computer code (Clune et al., 1999), which solves the anelastic equations of hydrodynamics describing a compressible fluid confined to a spherical shell heated from below. The fluid motions are calculated with respect to a rotating frame. The complex structures and intricate behavior of the resulting convection require high spatial resolution, and the flows must be studied over extended periods of time for statistical equilibrium to be achieved. As a result, the ASH code is designed to run efficiently on massively parallel architectures such as the Cray T3E, SGI Origin2000, and IBM SP-3 machines (e.g. Miesch et al. 2000; Elliott et al. 2000). The current multi-processor version follows the numerical approach first implemented by Glatzmaier (1984).

The ASH code employs a pseudo-spectral approach, where all fluid velocities and state variables are projected onto orthogonal basis functions in each of the three spatial dimensions. The radial structure of the solution variables is represented by an expansion based on Chebyshev polynomials characterized by radial order $`n`$, while functions of latitude and longitude are expanded over spherical harmonic basis functions $`Y_{\mathrm{}}^m`$ characterized by angular degree $`\mathrm{}`$ and azimuthal order $`m`$. This discretization scheme ensures that the horizontal resolution is approximately uniform everywhere on the spherical domain when all $`(\mathrm{},m)`$-pairs for a given maximum degree $`\mathrm{}_{\mathrm{max}}`$ are retained in the modal expansion (such as the triangular truncation used here). Conversely, the simplest finite-difference scheme, where computational gridpoints are uniformly distributed along lines of latitude and longitude, suffers from the problem that the spatial resolution varies with latitude such that the gridpoints are more closely spaced near the poles compared to equatorial regions (colloquially known as the pole problem). In addition, spectral methods also provide increased computational accuracy for a given grid size when compared to finite-difference representations.

The ASH code solves an approximate form of the Navier-Stokes equations known as the anelastic equations. The anelastic approximation (Gough, 1969) allows us to handle the effects of compressibility while filtering out acoustic perturbations that would otherwise severely limit the computational time step. This approximation is valid when the convective fluid velocities are subsonic, which in turn requires that the stratification of the fluid be only slightly superadiabatic. Such acoustic filtering is achieved by insisting that the time derivative of density vanishes in the continuity equation, or equivalently that the divergence of the momentum be zero, or that the momentum vector be solenoidal. This approximation is effectively equivalent to allowing pressure disturbances to equilibrate instantaneously, forcing the system to evolve on convective rather than sound-speed time scales. It is therefore implicitly assumed that sound waves do not play a significant role in the dynamical evolution of the system, which is in agreement with the expectation that the coupling of convection, stratification, and rotation are the major dynamical influences within the convection zone.

The convective motions are computed relative to a spherically symmetric mean state having a density $`\overline{\rho }`$, temperature $`\overline{T}`$, pressure $`\overline{p}`$, and specific entropy $`\overline{s}`$, much as discussed in Miesch et al. (2000); Brun & Toomre (2002). These mean quantities initially satisfy the equation of hydrostatic equilibrium,

$$\frac{d\overline{p}}{dr}=\overline{\rho }g,$$
(1)
where $`g`$ is the acceleration due to gravity, and a perfect-gas equation of state

$$\overline{p}=\frac{\gamma 1}{\gamma }c_p\overline{\rho }\overline{T},$$
(2)
from which the mean specific entropy $`\overline{s}`$ is defined (to within an arbitrary constant) by

$$\frac{d\overline{s}}{dr}=c_p\left(\frac{1}{\gamma \overline{p}}\frac{d\overline{p}}{dr}\frac{1}{\overline{\rho }}\frac{d\overline{\rho }}{dr}\right),$$
(3)
where the specific heats at constant pressure and volume are represented by $`c_p`$ and $`c_v`$ respectively, with their ratio $`\gamma =c_p/c_v`$. The mean thermodynamic state is allowed to evolve in time.

The anelastic approximation involves neglecting the time-derivative of density in the mass continuity equation (as described above), such that

$$\mathbf{}\mathbf{}(\overline{\rho }𝒖)=0,$$
(4)
and expanding the Navier-Stokes equations around the spherically symmetric mean state. Perturbations to the state variables are denoted by $`\rho `$, $`T`$, $`p`$, and $`s`$. The ASH code solves the equations describing the evolution of the fluid velocity,

$$\overline{\rho }\frac{𝒖}{t}=\overline{\rho }(𝒖\mathbf{}\mathbf{})𝒖+2\overline{\rho }(𝒖\mathbf{\times }𝛀_0)\mathbf{}p\rho g\widehat{𝒓}+\mathbf{}\mathbf{}𝓓,$$
(5)
and of the specific entropy,

$$\overline{\rho }\overline{T}\frac{s}{t}=\overline{\rho }\overline{T}(𝒖\mathbf{}\mathbf{})(\overline{s}+s)+\mathbf{}\mathbf{}𝒒+\mathrm{\Phi },$$
(6)
where the viscous stress tensor $`𝓓`$ and viscous heating term $`\mathrm{\Phi }`$ are defined

$$𝒟_{ij}=2\overline{\rho }\nu \left[e_{ij}\frac{1}{3}(\mathbf{}\mathbf{}𝒖)\delta _{ij}\right]$$
(7)
and

$$\mathrm{\Phi }=2\overline{\rho }\nu \left[e_{ij}e_{ij}\frac{1}{3}(\mathbf{}\mathbf{}𝒖)^2\right],$$
(8)
with $`e_{ij}`$ denoting the strain rate tensor. The diffusive heat flux $`𝒒`$ is defined

$$𝒒=\kappa _s\overline{\rho }\overline{T}\mathbf{}(\overline{s}+s)+\kappa _rc_p\overline{\rho }\mathbf{}(\overline{T}+T).$$
(9)
Furthermore, these equations are subject to the linearized equation of state,

$$\frac{\rho }{\overline{\rho }}=\frac{p}{\overline{p}}\frac{T}{\overline{T}}=\frac{p}{\gamma \overline{p}}\frac{s}{c_p}.$$
(10)
As with the temporal scales of motion, fully resolving all spatial scales of motion in a numerical simulation of the solar convection zone is infeasible at this time, as the dynamically active scales range from $`10^2`$ Mm (of order the depth of the zone) to $`10^4`$ Mm (typical dissipation scale), thereby encompassing a factor of at least $`10^6`$ in scale. Because current simulations can cope with a range of only about $`10^3`$ in each of the three physical dimensions, the ASH code adopts the common approach of parameterizing the transport properties of sub-grid scale (SGS) turbulent eddies and resolving only the largest scales of convection, thus becoming a large eddy simulation (LES).

All LES-SGS simulations require a prescription for representing the effects of SGS convective motions not explicitly resolved in the model. Such a scheme may incorporate characteristics of the resolved flows into their functional forms (see the reviews by Canuto 1996; Lesieur 1997; Canuto & Christensen-Dalsgaard 1998), or may simply enhance the viscous and thermal diffusivities relative to their molecular values. We have adopted the latter approach for simplicity, yet recognize that this aspect requires considerable attention in the future. The main drawback of this scheme is that the enhanced diffusion draws energy from larger resolved scales of motion which should be unaffected by such dissipative effects. In one alternative approach, known as hyperviscosity, one allows the enhanced eddy diffusivities to act on fourth- (or higher-) order derivatives of the velocity field, thereby confining the diffusive effects more toward the smaller end of the energy spectrum. Another class of SGS models involves adding extra stress terms to the equations of motion. Evolution equations for these additional contributions can then be constructed once functional forms for the correlations between second-order variables are specified using some kind of a closure hypothesis. As is true of all LES-SGS studies, one hopes that the specific form by which the SGS motions are parameterized has a relatively small effect on the global dynamics of the system, but this is a property that cannot be readily verified at this stage.

Time-stepping in the ASH code is performed using an implicit second-order Crank-Nicholson procedure for the linear terms and a fully explicit second-order Adams-Bashforth procedure for the nonlinear terms, with the exception of the linear Coriolis forcing terms in the momentum equation which are updated explicitly. Because the explicit time-stepping procedure cannot be performed in the spectral domain, such a pseudo-spectral scheme necessitates conversions between the physical and spectral representations during each time step when switching between solving the implicit and explicit terms in the evolution equations. However, the benefits gained by solving the equations in the spectral domain outweigh the added computational time spent in performing the transformations between the physical and spectral domains.

## 3 Model Formulation

### 3.1 Initializing the Spherically Symmetric Mean State

We have constructed two shallow-shell simulations (Cases S1 and S2) that span a radial extent of 0.94–0.98 $`R`$, equivalent to a shell thickness of 28 Mm. Two additional simulations (Cases D2 and D3) span 0.90–0.98 $`R`$ or 56 Mm, and are twice as thick as the two shallow-shell simulations. In Cases S2, D2, and D3, a solar-like differential rotation profile is applied as a no-slip lower boundary. For comparison purposes, the lower boundary of Case S1 is maintained at a uniform rotation rate equal to the angular velocity $`\mathrm{\Omega }_0`$ of the computational frame, with all other attributes equivalent to Case S2. In each of the four cases, the upper boundary is stress-free and both the lower and upper boundaries are impenetrable. The thermal driving in each case is accomplished by setting the heat flux at the lower boundary equal to the solar value, while the upper boundary is held at a uniform and constant entropy. The parameters of the four thin-shell simulations are summarized in Table 1.

As discussed in the previous section, the anelastic equations of motion are advanced in time by the ASH code using a pseudo-spectral approach, wherein functions of colatitude $`\theta `$ and longitude $`\varphi `$ are expanded over spherical harmonic functions characterized by angular degree $`\mathrm{}`$ and azimuthal order $`m`$. Functions of radius are projected onto Chebyshev polynomials characterized by radial order $`n`$. The four simulations presented here are calculated using spherical harmonic functions with a maximum angular degree $`\mathrm{}_{\mathrm{max}}=340`$, so that horizontal scales as small as about 10 Mm can be explicitly resolved. The highest order Chebyshev polynomial used is $`n=64`$ in Cases S1 and S2, and $`n=128`$ in Cases D2 and D3. Since we expect the resulting convection to have a limited longitudinal scale and seek computational economy, we impose a four-fold azimuthal symmetry by keeping every fourth $`m`$ value in the spherical harmonic expansions for Cases S1, S2, and D2. Such an imposed symmetry is primarily noticeable only in the high-latitude regions, where the convergence of meridional lines near the poles limits the longitudinal scale of the convective structures to sizes smaller than are present at lower latitudes. In addition, flows across the pole are not permitted. For comparison, Case D3 was computed without such an angular periodicity, but is otherwise equivalent to Case D2.

During initialization, the radial profiles of the mean density $`\overline{\rho }(r)`$, temperature $`\overline{T}(r)`$, and pressure $`\overline{p}(r)`$ are determined by jointly solving the equations (1)–(3), given radial profiles of the gravitational acceleration $`g(r)`$ and specific entropy gradient $`d\overline{s}/dr`$ throughout the domain. We specify the initial entropy gradient to have a slightly superadiabatic value (e.g. $`d\overline{s}/dr=10^7`$), while the function $`g(r)`$ is taken from the one-dimensional solar model of Brun et al. (1999). The initial values of $`\overline{\rho }(r)`$ and $`\overline{T}(r)`$ obtained in this way are shown in Figure 2 to compare favorably with the one-dimensional structure model, with the slight discrepancy in density resulting in a greater density scale height $`\lambda _\rho `$ in our simulations than in the structure model.

With our current angular resolution ($`\mathrm{}_{\mathrm{max}}=340`$), we are able to accommodate convection possessing horizontal size scales as small as a few angular degrees. Consequently, we place the upper boundary of each simulation at 0.98 $`R`$ where the imposed solar-like stratification should drive convective structures similar in size to solar supergranulation. The even greater degree of stratification present above this radius is likely to drive modes of convection that are too small in physical size to be resolved with our current grid size. In addition, as the convection becomes less efficient closer to the surface, typical convective velocities become a significant fraction of the speed of sound thereby making the anelastic approximation inappropriate. Furthermore, the ideal gas equation of state and the diffusive treatment of radiative transfer currently used in these calculations would not capture the effects resulting from the ionization of hydrogen. Consequently, the upper boundary of each simulation is located at 0.98 $`R`$.

The ASH code is a LES-SGS simulation and thus requires a prescription to account for the transport of energy and momentum by turbulent motions not explicitly resolved in the simulations. We have adopted the simplest approach of enhancing the molecular values of the viscous and thermal diffusivities, while recognizing that such an SGS treatment is unlikely to capture all dynamical effects of small-scale turbulence. The eddy diffusivities $`\nu (r)`$ and $`\kappa _s(r)`$ used in these simulations are chosen to vary inversely as the square root of the mean density profile, so that $`\beta =\frac{1}{2}`$ in

$$\nu =\nu _{\text{top}}\left(\frac{\overline{\rho }}{\overline{\rho }_{\text{top}}}\right)^\beta \text{and}\kappa _s=\kappa _{\text{top}}\left(\frac{\overline{\rho }}{\overline{\rho }_{\text{top}}}\right)^\beta .$$
(11)
This particular value of $`\beta `$ was chosen to allow some variation of the dissipation scale with the density scale height $`\lambda _\rho `$, but at the same time to prevent numerical instabilities near the bottom of the domain where the dissipation length scale is smallest. The free parameters $`\nu _{\text{top}}`$ and $`\kappa _{\text{top}}`$ in equation (11) are the viscous and thermal diffusivities at the upper boundary, which we have set to $`\nu _{\text{top}}=\kappa _{\text{top}}=1\times 10^{12}`$ cm<sup>2</sup> s<sup>-1</sup> in each case. Consequently, the Prandtl number $`P_r=\nu /\kappa _s=\nu _{\text{top}}/\kappa _{\text{top}}`$ throughout each domain is unity. We note that as computing technology becomes more advanced, the parameterization scheme used to account for SGS transport effects in simulations of highly turbulent fluids should become less of an issue because such global simulations will be able to explicitly resolve the convective motions and associated energy transport at smaller dynamical scales.

In order to prevent the formation of thin diffusive thermal layers having a radial thickness below the current radial resolution, we have introduced an unresolved enthalpy flux to the mean state, much as in Miesch et al. (2000). Without such enhancement, steep entropy gradients would otherwise be required to carry the imposed heat flux, owing to the small radial fluid velocities (and thus the convective transport of heat) in the vicinity of the impenetrable boundaries.

The characteristic size scales of the convective structures that appear near the top of these thin-shell simulations are somewhat sensitive to the degree of driving they experience, which is in part determined by the functional form of the unresolved enthalpy flux. The form of this quantity essentially determines what fraction of the total energy flux must be transported via convection near the top of the domain, which in turn feeds back on the entropy gradient and affects the convection. We find that adjusting either the diffusivities or the superadiabaticity within the domain (that is, altering the Rayleigh number of the convection) will alter the appearance of the convection, especially near the upper boundary. We have chosen a somewhat gentle functional form for this enthalpy flux, having an $`e`$-folding depth away from each boundary of about 0.01 $`R`$ for both the shallow and deep shell cases.

Finally, the radiative diffusivity $`\kappa _r(r)`$ throughout each domain is taken from the structure model used above; however, we note that the radiative heat flux is several orders of magnitude smaller in our simulations than convective heat flux throughout the bulk of the domain (as it is in the sun).

### 3.2 Approach to Thermal Equilibrium

Once the spherically symmetric mean state has been arranged, small seed perturbations are introduced into the fluctuating entropy field $`s`$, and the simulations are advanced in time using the evolution equations (5) and (6). The seed entropy perturbations are soon reflected in the fluctuating density field $`\rho `$, whose variations quickly provide the unstable density imbalance that buoyantly accelerates the fluid from rest. After an initial period of adjustment during which the convective kinetic energy ramps up (typically about 30 days, as shown in Fig. 3 for Case D2), an approximate thermal equilibrium is reached. The mean thermodynamic state will evolve in time, as reflected in the $`\mathrm{}=0`$ component of the thermodynamic variables.

In total thermal equilibrium, the outward energy transport in these simulations must be achieved by a balance of radiative, kinetic, enthalpy, and eddy diffusive fluxes, as represented by the following:

$`F_k`$ $`=u_r\left({\displaystyle \frac{\overline{\rho }u^2}{2}}\right),`$ (12)
$`F_e`$ $`=u_r\overline{\rho }c_p(TT),`$ (13)
$`F_r`$ $`=\kappa _r\overline{\rho }c_p{\displaystyle \frac{(\overline{T}+T)}{r}},`$ (14)
$`F_u`$ $`=\kappa _s\overline{\rho }\overline{T}{\displaystyle \frac{(\overline{s}+s)}{r}},`$ (15)
where the kinetic, enthalpy, radiative, and unresolved eddy fluxes are denoted by $`F_k`$$`F_e`$, $`F_r`$, and $`F_u`$ respectively. The quantity $`TT`$ appearing in the definition of $`F_e`$ is the temperature excess relative to the mean (spherically symmetric) value of the temperature field at each radial level. In a true steady state, the total energy transport at each radius within the domain must equal the imposed energy influx at the lower boundary:

$$4\pi r^2(F_k+F_e+F_r+F_u)=L_{},$$
(16)
where $`L_{}`$ is the solar luminosity.

Figure 4 shows the time-averaged energy transport within Cases S2 and D2 as a function of radius at a late stage in the simulations. In both cases, the enhanced energy transport by the unresolved eddies near the boundaries reflects the additional enthalpy flux applied to the mean state. Figure 4 also shows that Case D2 is not yet in total thermal equilibrium, as the energy output at the upper boundary is approximately 80% of the input applied to the lower boundary. The radial energy transport within Case S1 is qualitatively similar to Case S2, and likewise for Cases D2 and D3.

## 4 Characteristics of Resulting Flows

### 4.1 Multi-Scale Convection

The convective flow patterns in these simulations are intricate, containing complex evolving structures occurring on multiple size scales. We illustrate typical velocity patterns in Figures 5 and 6 by showing the radial velocity for Cases S2 and D2 at several depths within each domain. The horizontal structure of the radial velocity fields realized in Cases S1 and D3 are qualitatively similar to Cases S2 and D2, respectively.

Figures 5$`a`$ and 6$`a`$ show that the largest scale of convection visible near the top of both Cases S2 and D2 is associated with a connected network of downflow lanes (green-blue colors) having a spatial scale of about 200 Mm. The large areas enclosed by the downflow lanes each contain several smaller-scale upflows (orange-red tones) measuring about 15–30 Mm across. Although the upflow cells in the shallow shell (Case S2) tend to be slightly larger than those in the deeper shell (Case D2), this general surface pattern of a network of connected downflow lanes enclosing several distinct smaller upflows appears to be a robust property of convection within our thin shells. The size of the smaller upflows in these simulations is approximately 20 Mm, similar to the horizontal size scale of solar supergranulation.

Figures 5$`b`$,$`c`$ and 6$`b`$,$`c`$ illustrate how the horizontal planforms change with depth within Cases S2 and D2. As the downwelling fluid reaches deeper layers, the degree of connectivity of the downflow network appears to decrease, since the vertical velocities along each downflow lane becomes much less uniform. Near the bottom of each domain, the downflows have fragmented into more isolated and compact plume-like structures while retaining a vestige of the connected network of downflow lanes seen near the surface. This narrowing of scale with depth appears to be related to the larger densities found near the bottom of each domain.

The strongest downflow lanes in the equatorial region of Case D2 possess a noticeable north-south orientation that is somewhat reminiscent of the banana-cell modes evident in more laminar spherical shell convection simulations (e.g. simulation LAM in Miesch et al. 2000). However, the columnar structures seen in Case D2 are much thinner, have more spatial variation, and evolve more quickly than the banana cells in the more laminar cases, as a result of the strong small-scale convection driven throughout the domain.

The smaller-scale upflows visible near the upper boundary of these simulations also change their character with depth. The distinct upflows enclosed by downflow lanes seen in the upper layers gradually become more uniform, forming broad regions of upwelling fluid surrounded by an incomplete network of downflow lanes. Near the lower boundary, these broad upflow regions have largely disappeared, with the fastest upward velocities occurring in shrouds surrounding each of the strong downflows in response to the impenetrable lower boundary that exists at the bottom of the domains. Figure 7 contains vertical cuts of the radial velocity field showing the variation with latitude and radius for Cases S2 and D2. The radial structure of the upflow regions is most evident in Case D2 (Fig. 7$`b`$), wherein the smaller-scale yet faster upflows are only evident near the upper boundary of the domain.

Figure 8 shows the time evolution of the radial velocity field of Case S2. Such time series illustrate the tendency for features in the flow field to be systematically advected by velocity patterns possessing a larger scale. One small upflow (indicated by the arrow) as well as the downflow lane immediately to its left are both advected laterally by the horizontal outflow motions associated with the broader cell. The center of this broad outflow cell is indicated by a cross. Such lateral transport of velocity features by larger scales of convection is most apparent in movie sequences showing the time evolution of the radial velocity field.

In addition to the network of downflow lanes possessing larger spatial scales than are observed at the solar photosphere, we also find that these structures persist for longer time intervals than their counterparts on the sun. For example, the time series shown in Figure 8 illustrates that the prominent downflow network visible near the upper surface of Case S2 evolves on time scales longer than about 10 days. Similar evolutionary time scales are evident in Case D2 as well. Such discrepancies between the numerical simulations and the observations is not surprising, since the level of turbulence is much lower (by several orders of magnitude) in our simulations than for the sun. Furthermore, at 0.98 $`R`$ the upper boundaries of the simulations are located below the radius of the photosphere, and consequently may not accurately depict the convective cell patterns that would exist above this radius.

In Figure 9 we present space-time diagrams of radial velocity sampled in time for four specific latitudes (and all longitudes) at a radius of 0.978 $`R`$ for Case S2, plotted with respect to the rotation rate of the computational frame $`\mathrm{\Omega }_0`$. The advection of radial velocity structures appear as slanted features in each panel, with the higher latitudes rotating more slowly (retrograde) than those near the equator. The mean advection rate of these radial velocity patterns evident in the figure approximately equals the differential rotation rate of the fluid for the same depth; consequently, the latitudinal variation of angular velocity shown in the figure results largely from the imposed differential rotation at the lower boundary. However, it is interesting to note from Figure 9 that the pattern speed of the advected structures does vary with longitude, as structures having the same latitude but separated in longitude may exhibit slightly different rates or propagation.

At the mid-point in radius of their respective domains, the rms radial velocities are measured to be 105 m s<sup>-1</sup> for Case S2 and 140 m s<sup>-1</sup> for Case D2, with overturning times on the order of 6–10 days, depending on the shell depth. These overturning times suggest that the large-scale convective pattern may be weakly sensitive to rotational effects, as both Cases S2 and D2 are rotating at the solar-like mean rate of one rotation per 28 days. We will discuss the rotational influence of these convective overturning motions in more detail in §4.3.

Figure 10 shows the fluctuating temperature field for Case D2, after removing the spherically symmetric component of the temperature at each radius. We find that the locations of the large regions of warm and cool temperatures correlate well with the locations of the broad upflows and narrow downflow lanes visible in the radial velocity field shown in Fig. 6, with the primary difference being that the cooler regions tend to be much broader than their associated radial velocity counterparts. This feature is somewhat surprising given the equal thermal and viscous diffusivities (i.e. $`P_r=1`$) within Case D2, and is likely indicative of the degree to which the cool fluid associated with the strong downflow lanes undergoes mixing.

Near the equator, the temperature field is dominated by columns of alternating warm and cold fluid associated with the weak banana-cell-like structures visible in the radial velocity images. In both temperature and radial velocity, these structures are sheared slightly by the differential rotation within the domain, and extend up to about $`\pm 30^{}`$ of latitude. (For comparison, the cylinder tangent to the inner radius intersects the outer radius at about $`\pm 20^{}`$ of latitude.) These large-scale columnar temperature structures are broken up by smaller-scale variations on the temperature field which also tend to correlate well with some of the small-scale radial velocity features. For example, the localized hot and cold spots particularly evident in the close-up views of Figure 10$`a`$ are coincident with some of the fastest fluid motions visible in the radial velocity image of Figure 6.

Near the upper boundary of Case D2, there exists a significant latitudinal temperature contrast between the equator and the poles, as shown in Figure 11. The temperatures in the near-polar regions are about 10–15 K warmer than near the equator, although over half of the total equator-to-pole difference occurs within $`10^{}`$ of the pole. As stated earlier, we believe that many characteristics of the fluid in the near-polar regions are most likely artifacts of the four-fold azimuthal periodicity imposed in these simulations, and should be interpreted with care.

### 4.2 Axisymmetric Flow Patterns

The time-averaged axisymmetric (longitudinally averaged) profile of angular velocity $`\mathrm{\Omega }/2\pi `$ achieved within each simulation is shown in Figure 12, with the quantity $`\mathrm{\Omega }`$ now corresponding to the angular velocity in the rotating frame after the angular velocity of the computational frame $`\mathrm{\Omega }_0`$ is subtracted out. The no-slip rotating boundary is imposed such that equatorial rate corresponds to the rotation rate of the computational frame, such that $`\mathrm{\Omega }/2\pi =0`$ at these points. The differential rotation of the lower boundary imposed in Cases S2, D2, and D3 decreases from $`\mathrm{\Omega }/2\pi =0`$ nHz at the equator to about $`\mathrm{\Omega }/2\pi =120`$ nHz at a latitude of $`75^{}`$, and is similar in contrast and functional form to the latitudinal variation of the photospheric plasma rate measured by Snodgrass (1984). For comparison purposes, the no-slip lower boundary imposed in Case S1 is uniformly rotating at $`\mathrm{\Omega }/2\pi =0`$ nHz.

With the exception of the polar regions (not shown in the figure), the angular velocity profile within each simulation is retrograde with respect to the rotating coordinate system, with the fastest rotation rates occurring near the bottom of each shell at any given latitude. Figure 12 shows that Cases S1 and S2 possesses a largely constant negative radial gradient of angular velocity with radius at each point (see also Fig. 14), with the overall magnitude of $`\mathrm{\Omega }`$ determined by the rotation rate imposed at the corresponding latitude on the lower boundary. Within Cases D2 and D3, the negative radial gradients in angular velocity throughout the bulk of each shell are smaller, except within the thin viscous boundary layers located near the lower boundaries of both simulations.

The regions poleward of $`\pm 75^{}`$ of latitude are not shown in Figure 12 due to the high angular velocities present there. Even though these regions exhibit reasonable (linear) zonal velocities, the short moment arm at such high latitudes produces higher values of $`\mathrm{\Omega }`$ than are seen elsewhere within the domain at lower latitudes. Furthermore, the polar regions of Cases S1, S2, and D2 are dominated by effects related to the four-fold azimuthal periodicity imposed in these simulations, which disallows flows across the poles and limits the physical size of convective structures. Such effects are likely to produce dynamics which may not be otherwise present in analogous simulations without the four-fold azimuthal symmetry. Although all longitudes are computed in Case D3, it was started from an evolved solution of Case D2 and still shows some effects of the angular periodicity imposed in that simulation.

The time-averaged meridional circulation achieved in our simulations is illustrated in Figure 13, where the profile of the meridional velocity streamfunction $`\mathrm{\Psi }`$ for each simulation is shown. Most of the kinetic energy associated with the meridional flows is contained in $`15^{}`$-wide latitudinal rolls distributed across the mid-latitude and equatorial regions. These rolls span the entire depth of the domain both in the shallow- and deep-shell simulations and possess typical flow speeds of 50–75 m s<sup>-1</sup>. Rolls having poleward velocities in the surface layers appear to be preferred over rolls having the opposite rotational sense, especially in the mid- and low-latitude regions, as they are generally more extended horizontally and possess faster fluid velocities.

Alternatively, one could view these meridional flow patterns as an array of symmetric clockwise and counterclockwise rolls superimposed on a single meridional cell in each hemisphere. We shall see that the magnitude of the single hemispheric meridional cell is larger relative to the smaller-scale rolls when the lower boundary is differentially rotating, as in Cases S2, D2, and D3. We note, however, that because these simulations have impenetrable boundaries, any meridional flow achieved in our simulations is forced to close within the (relatively shallow) domain. In this regard, the profiles of meridional circulation obtained in these four simulations are likely to be unrealistic, as the sun obviously does not contain impenetrable boundaries above and below the near-surface shear layer.

### 4.3 Angular Momentum Balance

In order to understand the maintenance of the differential rotation profiles shown in Figure 12, we will now examine the angular momentum balance within each of the four simulations in more detail. We will consider only axisymmetric quantities, defined by

$$\widehat{A}(r,\theta )=\frac{1}{2\pi }_0^{2\pi }𝑑\varphi A(r,\theta ,\varphi ),$$
(17)
where the hat $`\widehat{}`$ signifies that the quantity $`A(r,\theta ,\varphi )`$ has been averaged over longitude $`\varphi `$ and further over time. Consequently, this operation allows the decomposition

$$A(r,\theta ,\varphi )=\widehat{A}(r,\theta )+A^{}(r,\theta ,\varphi )\text{such that}\widehat{A^{}}=0,$$
(18)
where the prime on $`A^{}`$ denotes the non-axisymmetric component of $`A`$.

By multiplying the zonal component of the momentum evolution equation (5) by $`r\mathrm{sin}\theta `$, we can derive an equation describing the evolution of the angular momentum $`\widehat{L}=\overline{\rho }r\mathrm{sin}\theta \widehat{u_\varphi }`$, written in conservative form as

$$\frac{\widehat{L}}{t}=\mathbf{}\mathbf{}𝑭,$$
(19)
where $`𝑭`$ is the angular momentum flux vector. In a statistically steady state, we must have $`\widehat{L}/t=0`$, such that $`\mathbf{}\mathbf{}𝑭=0`$ throughout the domain. Symbolically, the components of the total flux in the radial and latitudinal directions can be written as

$`F_r`$ $`=F_r^{\text{DIF}}+F_r^{\text{RS}}+F_r^{\text{MC}}`$ (20)
and
$`F_\theta `$ $`=F_\theta ^{\text{DIF}}+F_\theta ^{\text{RS}}+F_\theta ^{\text{MC}},`$ (21)
where we have used the abbreviations in capital letters to signify the contributions to the total angular momentum flux resulting from viscous diffusion (DIF), non-axisymmetric Reynolds stresses (RS), and axisymmetric meridional circulation (MC):

$`F_r^{\text{DIF}}`$ $`=\overline{\rho }r\mathrm{sin}\theta \left[\nu r{\displaystyle \frac{}{r}}\left({\displaystyle \frac{\widehat{u_\varphi }}{r}}\right)\right]`$ (22)
$`F_r^{\text{RS}}`$ $`=\overline{\rho }r\mathrm{sin}\theta \widehat{u_r^{}u_\varphi ^{}}`$ (23)
$`F_r^{\text{MC}}`$ $`=\overline{\rho }r\mathrm{sin}\theta \widehat{u_r}(\widehat{u_\varphi }+\mathrm{\Omega }_0r\mathrm{sin}\theta )`$ (24)
and
$`F_\theta ^{\text{DIF}}`$ $`=\overline{\rho }r\mathrm{sin}\theta \left[\nu {\displaystyle \frac{\mathrm{sin}\theta }{r}}{\displaystyle \frac{}{\theta }}\left({\displaystyle \frac{\widehat{u_\varphi }}{\mathrm{sin}\theta }}\right)\right]`$ (25)
$`F_\theta ^{\text{RS}}`$ $`=\overline{\rho }r\mathrm{sin}\theta \widehat{u_\theta ^{}u_\varphi ^{}}`$ (26)
$`F_\theta ^{\text{MC}}`$ $`=\overline{\rho }r\mathrm{sin}\theta \widehat{u_\theta }(\widehat{u_\varphi }+\mathrm{\Omega }_0r\mathrm{sin}\theta ).`$ (27)
By examining the contributions to the total angular momentum flux from each of these components, we will show that the radial gradients of angular velocity realized in each of the four simulations are supported against diffusion by Reynolds stresses associated with non-axisymmetric convective motions. This behavior can be thought of as the tendency of convective fluid elements to partially conserve their angular momentum per unit mass $`\lambda =\mathrm{\Omega }r^2\mathrm{sin}^2\theta `$ as they move toward or away from the axis of rotation. As suggested by Foukal & Jokipii (1975), constancy of $`\lambda `$ along radial lines may also explain why surface magnetic tracers on the sun have a faster rotation rate relative to the surface fluid, if one assumes that the magnetic features are anchored at a radius slightly below the photosphere where the rotation rate is faster. Gilman & Foukal (1979) tested this notion by numerically modeling Boussinesq convection confined to a thin shell, and found that angular momentum was roughly conserved along local radii for the case of an incompressible fluid. They demonstrated that the convective motions were able to transport angular momentum inward, thereby maintaining the negative radial gradient of rotation rate with radius. The models presented here indicate that compressible convection behaves similarly.

We begin by examining the angular momentum balance within the two shallow-shell simulations, Cases S1 and S2. These two simulations are identical, except that Case S2 has a differentially rotating lower boundary, whereas in Case S1 it is forced to rotate uniformly. The upper boundary in both cases is stress-free. From Figure 12 it is apparent that the radial angular velocity gradient $`\mathrm{\Omega }/r`$ is negative over a broad range of latitudes. This effect is further illustrated in Figure 14, where the angular velocity profiles shown in Figure 12 for Cases S1 and S2 are plotted as a function of radius for selected latitudes. The dash-dot-dot-dot lines in each panel indicate the angular velocity of an isolated fluid parcel, assuming it were moving in a purely radial direction and conserving $`\lambda `$ (making $`\mathrm{\Omega }r^2`$).

Figure 15 illustrates the contributions to the total radial and latitudinal angular momentum fluxes from each of the component fluxes listed in equations (20) and (21). The radial fluxes of Cases S1 and S2 are shown in Figures 15$`a`$,$`b`$ after integrating over latitude. In both cases, the average angular momentum fluxes associated with viscous diffusion $`F_r^{\text{DIF}}`$ and meridional circulation $`F_r^{\text{MC}}`$ are positive at all radii, signifying a transport of angular momentum toward the top of the domain. We find that this upward transport is balanced almost entirely by the downward transport achieved by non-axisymmetric Reynolds stresses associated with the nonaxisymmetric convective motions. For this effect to occur in this manner, the velocity correlation $`\widehat{u_r^{}u_\varphi ^{}}`$ must be negative on average by equation (23) for $`F_r^{\text{RS}}<0`$, indicating that small-scale non-axisymmetric radial motions, such as the fast downflows, tend to possess a retrograde tilt in the $`r\varphi `$-plane. Such a tendency is exactly what one would expect if non-axisymmetric convective motions were to partially conserve their angular momentum per unit mass $`\lambda `$ in radial motion.

In the sun, while $`\mathrm{\Omega }/r`$ is negative between 0.94–0.98 $`R`$, its magnitude is much smaller than in Cases S1 and S2, as shown in Figure 14. In the figure, the dash-dot lines denote the run of $`\mathrm{\Omega }`$ with radius as determined from helioseismology. Since the viscosity of the sun is much lower, and since the solar meridional circulation is of the same magnitude as in our shallow-shell cases, the tendency for non-axisymmetric convective motions to conserve $`\lambda `$ is likely to be less effective in the sun between 0.94–0.98 $`R`$ than in the numerical simulations presented here. The helioseismic results surfaceward of 0.98 $`R`$, also plotted in the figure, suggest that the smaller-scale convective motions located closer to the photosphere may tend to conserve their angular momentum more fully, thereby maintaining a steeper angular velocity gradient $`\mathrm{\Omega }/r`$ in this region.

Figure 14 shows that both Cases S1 and S2 possess viscous boundary layers near the bottom of each domain, especially at higher latitudes. These Ekman-type layers are formed in response to the zonal velocity imposed at the lower boundary in each case and have no physical analog in the sun. Within this boundary, the influence of viscous dissipation is large enough to flatten out the angular velocity gradient at low latitudes and even produce a positive $`\mathrm{\Omega }/r`$ at the higher latitudes (e.g. at $`60^{}`$ in Fig. 14). We shall see that this effect is even more pronounced in Cases D2 and D3.

The latitudinal flux of angular momentum in Cases S1 and S2 is shown in Figures 15$`c`$,$`d`$ after integrating in radius. We find that Case S2, having a differentially rotating lower boundary, has a much stronger poleward transport of angular momentum relative to Case S1 in which the lower boundary is uniformly rotating. This effect is an expected consequence due the additional angular momentum being removed from the system at the lower boundary in the near-polar regions of Case S2, since the lower boundary in this case is differentially rotating rather than uniformly rotating as in Case S1. Such an enhanced poleward transport of angular momentum within Case S2 is achieved primarily by the convective motions, and to a lesser degree by the meridional circulation, as illustrated by the magnitudes of the fluxes $`F_\theta ^{\text{RS}}`$ and $`F_\theta ^{\text{MC}}`$ in Figure 15$`d`$. A recent analysis of photospheric velocity images suggests that there may be observational evidence of a similar effect occurring on the sun, where the Reynolds stresses associated with solar supergranulation yield a poleward transport of angular momentum within the near-surface layers of the sun (Hathaway et al., 2001).

Superimposed on the broad poleward transport, the angular momentum fluxes $`F_\theta ^{\text{RS}}`$ and $`F_\theta ^{\text{MC}}`$ possess significant latitudinal fluctuations that are directly correlated with the axisymmetric rolls of Figure 13. We find that the contributions to $`F_\theta ^{\text{MC}}`$ from the latitudinal rolls is partially offset by matching contributions to $`F_\theta ^{\text{RS}}`$, such that the net angular momentum transport by the rolls is somewhat limited.

Turning to the deep-shell simulations of Cases D2 and D3, we find that the angular velocity $`\mathrm{\Omega }`$ possesses radial gradients similar to those achieved within Cases S1 and S2. Radial cuts through each domain, such as those presented in Figure 16, show that the radial gradient of angular velocity is negative above 0.91 $`R`$ for each of the latitudes shown. The primary difference between the shallow and deeper shells is the magnitude of the viscous boundary layer adjacent to the lower boundary. The radial gradients of $`\mathrm{\Omega }`$ are larger in magnitude across the viscous boundary layer for the deeper shells, causing the diffusive transport of angular momentum down the radial gradient of $`\mathrm{\Omega }`$ to be correspondingly greater in Cases D2 and D3 than for Cases S1 and S2. This enhanced diffusive transport is illustrated in the curves of $`F_r^{\text{DIF}}`$ of Figure 17$`a`$,$`b`$, where the greater torques applied to the fluid via the no-slip lower boundary are evident in the enhanced diffusive fluxes $`F_r^{\text{DIF}}`$ within Cases D2 and D3. Elsewhere in the domain, viscous effects play a reduced role compared to the radial angular momentum transport achieved by the meridional circulation and Reynolds stresses. As in Cases S1 and S2, the radial Reynolds stress flux $`F_r^{\text{RS}}`$ is negative and the radial meridional circulation flux $`F_r^{\text{MC}}`$ is positive on average throughout the entire domain.

Figures 17$`c`$,$`d`$ show the latitudinal fluxes of angular momentum within Cases D2 and D3, after integrating over radius. These two cases indicate that angular momentum is primarily transported poleward throughout the domain by the meridional circulation (rather than by the Reynolds stresses as in Cases S1 and S2). However, the $`15^{}`$ rolls now contain more kinetic energy than their counterparts in the shallow-shell simulations, and consequently the variations superimposed on $`F_\theta ^{\text{MC}}`$ are larger. As with the shallower cases, the Reynolds stress flux $`F_\theta ^{\text{RS}}`$ due to these rolls tends to partially offset the latitudinal transport of angular momentum by meridional flows. By comparing each of the four simulations, we suspect that the strength of the meridional rolls in our thin-shell domains is affected by the shell depth, which effectively limits the radial extent of any meridional circulation within the domain.

In summary, the time scales of the largest overturning motions in our simulations suggest that they are at least weakly influenced by rotational effects, which in turn may enable Reynolds stresses to facilitate transport angular momentum inward. This inward angular momentum transport balances the outward transport by diffusion and meridional flows on average, thereby maintaining a negative angular velocity gradient throughout much of the layer at low and mid latitudes. This effect may contribute to the observed decrease of $`\mathrm{\Omega }`$ with radius in the near-surface shear layer of the sun as deduced from helioseismic analyses. Behavior at high latitudes is somewhat more complex due to the presence of a viscous boundary layer near the lower boundary, and likewise there is some uncertainty in the helioseismic inferences about the radial gradient in $`\mathrm{\Omega }`$ achieved at latitudes of 60 or greater.

## 5 Conclusions

We have presented results of three-dimensional numerical simulations of turbulence confined to thin rotating spherical shells, seeking to understand some of the dynamical effects that supergranular scales of motion within thin shearing layers might have within the analogous layer located near the top of the solar convection zone. We have focused our analysis on the physical processes that enable the transport of angular momentum within the thin shell domains, in order to investigate analogous angular momentum transport on the sun, and to determine the cause of the negative radial angular velocity gradients shown to exist within the near-surface shear layer of the sun.

The high-resolution simulations presented here allow horizontal structures of order 10 Mm to be explicitly resolved, thereby allowing us for the first time to incorporate dynamical scales on the order of solar supergranulation within global simulations of solar convection. We find that the broad spectrum of scales of motion, while typically much smaller in size than the largest characteristic length scales of the convection zone, are able to influence the large-scale dynamics of the system through their ability to transport angular momentum within the shells on a global scale.

The vigorous convection realized in each of the four simulations presented here is driven by imposing the solar heat flux at the lower boundary of each domain. Shell thicknesses of 4% and of 8% of the solar radius $`R`$ are considered, with the upper boundary in each case located at a radius of 0.98 $`R`$. We have imposed a differentially rotating no-slip lower boundary in three of the four simulations. Both the lower and upper boundaries are impenetrable, and the upper boundary is forced to remain at a constant entropy.

We find in all cases that near the middle of the domain the convection takes the form of a connected network of fast yet narrow downflow lanes that enclose broad regions of warmer, more slowly rising fluid. The cells enclosed by the downflow network typically measure 100–200 Mm across, with the lanes themselves about 20 Mm wide. In the deeper layers where the density is greater, this network loses much of its horizontal connectivity, instead forming more plume-like structures that approach the bottom of the domain, at which point the impenetrability of the lower boundary forces the fluid to be diverted horizontally.

Closer to the surface, the broad cells of upwelling fluid are found to segment into several smaller upflows having a horizontal scale comparable to that of solar supergranulation, of order 20–40 Mm. These small-scale upflow cells appear in both the shallow- and deep-shell simulations, suggesting that the more superadiabatic stratification present near the top of each domain, rather than the depth of the shell, is the primary factor that determines the morphology of the convection near the surface. Time series of the evolving near-surface flow field show that both the smaller upflow cells as well as the narrow downflow lanes are horizontally advected in a sustained fashion as they respond to larger-scale sweeping flows that develop nearby.

Averaging the flow fields in longitude reveal that the angular velocity decreases with radius in the low- and mid-latitude regions of each domain, with the exception of a thin viscous boundary layer that forms near the no-slip lower boundary in each case. An analysis of the angular momentum balance shows that such negative radial gradients of angular velocity are maintained by an inward transport of angular momentum, achieved by Reynolds stresses associated with the convective motions that balance the outward transport of angular momentum from viscous diffusion and the global meridional circulation. Such an inward transport is achieved if radially moving fluid motions, such as the broad upflows and strong downflows seen here, have the tendency to conserve their angular momentum per unit mass while moving radially throughout the shell.

The longitudinally averaged meridional velocity patterns take the form of a series of $`15^{}`$ latitudinal rolls that span the full radial thickness of each shell, with cells having a poleward surface flow tending to have a broader latitudinal extent than cells with the opposite sense. The net effect of these rolls is to transport angular momentum poleward, as required by latitudinally varying angular momentum flux imposed by the differentially rotating lower boundary. However, we reemphasize that the profiles of meridional circulation within these thin-shell simulations are significantly influenced by the impenetrable radial boundaries of our simulations, effectively forcing any circulation in the meridional plane to be completely enclosed within the domain. As a result, the meridional flow profiles realized here are not expected to resemble the solar case.

The continual advance of supercomputing technology will allow simulations of convection within thin spherical shells to be extended to deeper layers in the future. Such models would preclude having to add angular momentum to the system via no-slip boundary conditions, as the differential rotation and near-surface shear layer could then be computed in a self-consistent manner. In addition, large-scale flows driven by the small-scale convective patterns that would not normally be confined to the near-surface layers (such as the meridional circulation) would then be allowed to feed back on the deeper layers below. Preliminary attempts to construct such global models, encompassing the bulk of the convection zone as well as a more highly resolved layer where convection on supergranular scales can exist, are currently underway.

We also believe the treatment of sub-grid scale (SGS) convective motions not explicitly resolved in our simulations deserves considerable attention in the future. The current prescription, whereby the diffusivities are enhanced over their thermal and molecular values, is adopted only for simplicity, and likely does not capture all of the relevant effects of the unresolved scales on the global dynamics. Other treatments, such as those discussed in §2 whose functional forms depend on the shearing properties of the resolved flows, may be more appropriate.

Nevertheless, the thin-shell simulations presented here contain highly evolving, multi-scale convective motions that are able to efficiently redistribute angular momentum in both radius and latitude. Such motions are found to maintain the radial shear within the domains, even in this idealized environment that only approximates the near-surface shear layer of sun by decoupling it from the bulk of the convection zone. We also speculate that convection in the near-surface layers of the sun may behave in a similar fashion, maintaining the negative radial gradients in the near-surface shear layer of the sun as deduced from helioseismology. While we are admittedly still far removed from directly modeling a convective layer with realistic solar parameters, it is encouraging that the flow patterns realized in these simulations exhibit the multiple scales of supergranulation and the more global convection cells being deduced from local helioseismic probing of the near-surface shear layer.

We thank Mark Miesch, Sacha Brun, and Nicholas Brummell for several rounds of useful discussions and constructive feedback during the analysis and writing phases of this paper. This work was partly supported by NASA through grants NAG 5-7996, NAG 5-8133, and NAG 5-10483, and by NSF through grant ATM 9731676. Various phases of the simulations with ASH were carried out through NSF PACI support of the NRAC allocation MCA93S005S at the National Center for Supercomputing Applications (NCSA), the San Diego Supercomputing Center (SDSC), and the Pittsburgh Supercomputing Center (PSC). Much of the analysis of the extensive data sets was carried out in the Laboratory for Computational Dynamics (LCD) within JILA.


format: text

text:
# Properties of Circum(sub)stellar Accretion Disks

## 1. A Survey for Brown Dwarf Disks

The existence of circumstellar disks and their role in planet formation are well established for young solar-type stars. However, little is known about disks around young substellar objects. Such circum(sub)stellar disks might provide a laboratory for studying the physical processes of disks over a wide range of (central object and disk) mass. In addition, the presence of disks around young substellar objects may be an important clue to the origin of brown dwarfs.

We have recently completed a large thermal IR ($`L^{}`$-band; 3.8 $`\mu `$m) survey to study the frequency and properties of disks around young brown dwarfs and very low mass stars (Liu et al. 2002; see also Liu 2002). Our sample comprises young ($``$1–3 Myr) objects in nearby star-forming regions which have been spectroscopically classified to be very cool, corresponding to masses of $``$15 to $``$100 $`M_{\mathrm{Jup}}`$ based on current models. As described in our paper, the objects constitute a well-defined sample and are largely free of selection biases.

A priori, brown dwarf disks are expected to be harder to detect than disks around stars because of lower contrast. Substellar objects are less luminous and have shallower gravitational potentials. Hence, their disks should be relatively cool and may have negligible excesses in the commonly used $`JHK`$ (1.1$``$2.4 $`\mu `$m) bands. In fact, we find that thermal IR data are required to detect most disks around young brown dwarfs. Our survey is also sensitive enough to detect brown dwarf photospheres — hence the absence of a disk can be discerned, and the disk frequency of young brown dwarfs can be measured for the first time.

## 2. Disk Frequency and Properties

We find $``$75% of our sample have intrinsic IR excesses: disks around young brown dwarfs and very low mass stars appear to be common. This disk fraction is similar to that for T Tauri stars in the same star-forming regions. Figure 1 presents histograms of the IR excesses for our sample as a function of spectral type. The three earliest spectral type bins (M5.7$``$M6.4, M6.5$``$M7.4, and M7.5$``$M8.4) are all very similar, based on the K-S test. Since the evolutionary tracks at fixed mass are roughly constant in $`T_{eff}`$ for young ages, spectral types provide a relative mass scale. Hence, we find that the excesses are largely independent of central mass. The exceptions are the coolest (lowest mass) objects, types M8.5$``$M9.4, where the excesses are small, consistent with non-existent.

Figure 2 examines the age dependence of the disks. If we separate objects by spectral subclass, the absolute $`K_S`$-band magnitude can be used as a surrogate for age, due to the fact that model isochrones are roughly horizontal in the HR diagram. We find no statistically significant correlation between the IR excesses and $`M_{K_S}`$ — this suggests that the inner regions of young brown dwarf disks do not evolve substantially over the first $``$3 Myr. This timescale is in accord with studies of disks around T Tauri stars (e.g. Strom et al. 1989).

The $`L^{}`$-band flux arises from the inner disk regions ($`<`$0.1 AU) and hence is sensitive to the presence of inner holes. The disk emission from a specific object depends on the viewing angle to the observer, which is of course unknown. This limitation can be overcome with a large unbiased sample of objects — meaningful constraints on the inner holes can then be found from (1) the maximum observed IR excess and (2) the observed distribution of IR excesses. An example of the latter is shown in Figure 3, which compares the observed IR excess distribution with passive (reprocessing) disk models having different inner hole sizes. Inner holes of $`2R_{}`$ are favored by the observations; disks without inner holes would produce many more objects with large IR excesses than are actually observed.

What is the origin of these holes? For T Tauri stars, inner holes are believed to be created through the truncation of the disk by strong, closed stellar magnetic fields; disk matter then reaches the star by accretion along the field lines. A similar situation may be relevant to young brown dwarfs. Rough estimates based on magnetospheric accretion models say that fields of a few to several hundred Gauss are needed to produce the inner holes (see Liu et al. 2002). These values are modest compared to those measured for pre-main-sequence stars.

## 3. Connections to Disks around Young Solar-Type Stars

It appears that much of the observational paradigm developed for T Tauri disks can be extended to disks around young brown dwarfs. Based on data from the literature, much of our sample shows H$`\alpha `$ emission, suggesting that most young brown dwarfs are accreting, although the inferred rates are much lower than for T Tauri stars (e.g. Muzerolle et al. 2000). Most of our objects show both IR excesses and strong H$`\alpha `$ emission, like the higher mass classical T Tauri stars. Others have little/no IR excess and little H$`\alpha `$ emission, analogous to weak-line T Tauri stars.

We also find that brown dwarf disks, like those of T Tauri stars, possess inner holes. While the hole sizes are very different for T Tauri stars and brown dwarfs, the estimated disk temperatures at the inner radii are similar, of order 1200 K. Moreover, the inner disks are hot enough for significant thermal ionization (of Na and K), which is needed to permit coupling of the magnetic field to the disk. Thus, our inferred hole sizes suggest that the magnetic accretion paradigm developed for T Tauri stars may extend to substellar masses.

Finally, the high disk fraction of young brown dwarfs raises the possibility of forming planets around these objects. Such planetary systems would represent a fascinating alternative to the numerous systems found around solar-type stars.

## References

Liu, M. 2002, in Proc. of IAU 211: Brown Dwarfs, in press (astro-ph/0207477)

Liu, M., Najita, J. & Tokunaga, A. 2002, ApJ, submitted

Muzerolle, J. et al. 2000, ApJ, 545, L141

Strom, K., Strom, S., Edwards, S., Cabrit, S. & Struksie, M. 1989, AJ, 97, 1451