On the Penetration of Meridional Circulation below the Solar Convection Zone II: Models with Convection Zone, the Taylor-Proudman constraint and Applications to Other Stars
aa r X i v : . [ a s t r o - ph . S R ] J un On the Penetration of Meridional Circulation below the SolarConvection Zone II: Models with Convection Zone, theTaylor-Proudman constraint and Applications to Other Stars.
P. Garaud & L. Acevedo-Arreguin
Department of Applied Mathematics and Statistics, Baskin School of Engineering,University of California Santa Cruz, 1156 High Street, CA 95064 Santa Cruz, USA
Abstract
The solar convection zone exhibits a strong level of differential rotation, whereby therotation period of the polar regions is about 25-30% longer than the equatorial regions. TheCoriolis force associated with these zonal flows perpetually “pumps” the convection zonefluid, and maintains a quasi-steady circulation, poleward near the surface. What is theinfluence of this meridional circulation on the underlying radiative zone, and in particular,does it provide a significant source of mixing between the two regions? In Paper I, we beganto study this question by assuming a fixed meridional flow pattern in the convection zone andcalculating its penetration depth into the radiative zone. We found that the amount of mixingcaused depends very sensitively on the assumed flow structure near the radiative–convectiveinterface. We continue this study here by including a simple model for the convection zone“pump”, and calculating in a self-consistent manner the meridional flows generated in thewhole Sun. We find that the global circulation timescale depends in a crucial way on twofactors: the overall stratification of the radiative zone as measured by the Rossby numbertimes the square root of the Prandtl number, and, for weakly stratified systems, the presenceor absence of stresses within the radiative zone capable of breaking the Taylor-Proudmanconstraint. We conclude by discussing the consequences of our findings for the solar interiorand argue that a potentially important mechanism for mixing in Main Sequence stars hasso far been neglected.
1. Introduction
Various related mechanisms are thought to contribute to the generation and maintenanceof large-scale meridional flows in the solar convection zone. The effect of rotation on turbulentconvection induces a relatively strong anisotropy in the Reynolds stresses (Kippenhahn 1963),in particular near the base of the convection zone where the convective turnover time is of the 2 –order of the solar rotation period. The divergence of these anisotropic stresses can directlydrive large-scale meridional flows (see R¨udiger, 1989, for a discussion of this effect). It alsodrives large-scale zonal flows (more commonly referred to as differential rotation) which theninduce meridional forcing through the biais of the Coriolis force, a mechanism referred to as“gyroscopic pumping” (McIntyre 2007). Indeed, the polar regions of the convection zone areobserved to be more slowly rotating than the bulk of the Sun (Schou et al. towards the polar axis. Meanwhile,equatorial regions are rotating more rapidly than the average, and are therefore subjectto a Coriolis force pushing fluid away from the polar axis. The most likely flow patternresulting from the combination of these forces is one with an equatorial upwelling, a surfacepoleward flow and a deep return flow. This pattern is indeed observed near the solar surface:poleward surface and sub-surface flows with velocities up to a few tens of meters per secondhave been observed by measurements of photospheric line-shifts (Labonte & Howard 1982)and by time-distance helioseismology (Giles et al. 1997) respectively.The amplitude and spatial distribution of these meridional flows deeper in the convectionzone remains essentially unknown, as the sensitivity of helioseismic methods rapidly dropsbelow the surface. As a result, the question of whether some of the pumped mass flux actuallypenetrates into the underlying radiative zone is still open, despite its obvious importance formixing of chemical species (Pinsonneault, 1997; Elliott & Gough, 1999), and its presumedrole in the dynamical balance of the solar interior (Gough & McIntyre 1998, McIntyre,2007, Garaud, 2007, Garaud & Garaud, 2008) and in some models of the solar dynamo (seeCharbonneau, 2005 for a review).In Paper I (Garaud & Brummell, 2008), we began a systematic study of the penetrationof meridional flows from the convection zone into the radiative zone by considering a relatedbut easier question: assuming that the amplitude and geometry of meridional flows in theconvection zone are both known , what is their influence on the underlying radiative zone?This simpler question enabled us to study the dynamics of the radiative zone only by assum-ing a flow profile at the radiative–convective interface (instead of having to include the morecomplex convection zone in the calculation). The overwhelming conclusion of that first studywas that the degree to which flows penetrate into the stratified interior (in the model) is very sensitive to the interfacial conditions selected. Hence, great care must be taken when usinga “radiative-zone-only” model to make definite predictions about interior flow amplitudes.In addition, that approach makes the implicit assumption that the dynamics of the radiativezone do not in return influence those of the convection zone, but the only way to verify thisis to construct a model which includes both regions. This was the original purpose of thepresent study; as we shall see, the combined radiative–convective model we construct hereprovides insight into a much broader class of problems. 3 –We therefore propose a simplified model of the Sun which includes both a “convective”region and a “radiative” region, where the convective region is forced in such a way as topromote gyroscopic pumping of meridional flows. We calculate the flow solution everywhereand characterize how it scales in terms of governing parameters (e.g. stellar rotation rate,stratification, diffusivities, etc..), focusing in particular on the flows which are entering theradiative zone. We begin with a simple Boussinesq Cartesian model (Section 2), first in theunstratified limit (Section 2.3) and then in the more realistic case of a radiative–convectivestratification (Section 2.4). Although the Cartesian results essentially illustrate most of therelevant physical phenomena, we confirm our analysis with numerical solutions of the full setof equations in a spherical geometry in Section 3. We then use this information in Section4 to discuss the effects of mixing by meridional flows both in the Sun and in other MainSequence stars.
2. A Cartesian model2.1. Model setup
As in Paper I, we first study the problem in a Cartesian geometry. Since our primary aimis to understand the behavior of the meridional flows generated (e.g. scaling of the solutions)in terms of the governing parameters, this approach is sufficient and vastly simplifies therequired algebra. In Section 3, we turn to numerical simulations to study the problem in aspherical geometry.In this Cartesian model section, distances are normalized to the solar radius R ⊙ , andvelocities to R ⊙ Ω ⊙ where Ω ⊙ is the mean solar angular velocity (the exact value is notparticularly relevant here). The coordinate system is ( x, y, z ), where x should be thought ofas the azimuthal coordinate φ , with x ∈ [0 , π ]; y represents minus the co-latitude and spansthe interval y ∈ [0 , π ] (the poles are at y = 0 and y = π while the equator is at y = π/ z -direction is the radial direction with z ∈ [0 , z = 0 is the interior, and z = 1 is the surface.In this framework, the system rotates with mean angular velocity Ω = (0 , , h torepresent the radiative–convective interface. Thus z ∈ [0 , h ] represents the “radiative zone” 4 –while z ∈ [ h,
1] represents the “convection zone”. From here on, h = 0 .
7. Figure 1 illustratesthe geometry of the Cartesian system. π/2 y= z r (cid:0)(cid:0)(cid:0)(cid:1)(cid:1)(cid:1) z=0z=hz=1y=0 Ω x −θ y φ "RADIATIVE" ZONE"CONVECTIVE" ZONE Fig. 1.— Cartesian model geometry and intended correspondence with the spherical case.The shaded area marks the convective region, where forcing is applied. The y = 0 and y = π/ π in the y − direction. For this simple Cartesian approach, we work with the Boussinesq approximation (thisassumption is dropped in the spherical model but doesn’t affect the predicted scalings).The background state is assumed to be stratified, steady, and in hydrostatic equilibrium.The background density and temperature profiles are denoted by ¯ T ( z ) and ¯ ρ ( z ) respectively.Density and temperature perturbations to this background state ( ρ and T ) are then assumedto be linearly related: ρ = − α T ( z ) T , where α T ( z ) is the coefficient of thermal expansion.For simplicity, we will assume that the background temperature gradient ¯ T z is constantthroughout the domain z ∈ [0 , α T → α T = 0 (see below). The alternative option of using a constant α T and a varying ¯ T z yields qualitatively equivalent scalings for the meridional flows (a statementwhich is verified in Section 3) although the algebra is trickier.The set of equations governing the system in this approximation are the momentum,mass and thermal energy conservation equations respectively: ∂ u ∂t + u · ∇ u + 2 e z × u = −∇ p + Ro ( z ) T e z + E ν ∇ u , ∇ · u = 0 , 5 – ∂T∂t + u · ∇ T + u · e z = E ν Pr ∇ T , (1)where u = ( u, v, z ) is the velocity field. In these equations, temperature perturbationshave been normalized to the background temperature difference across the box R ⊙ ¯ T z . Thegoverning non-dimensional parameters are:Ro( z ) = N ( z ) / Ω ⊙ the Rossby number ,E ν = νR ⊙ Ω ⊙ the Ekman number ,Pr = νκ T the Prandtl number , (2)where the dimensional quantity N ( z ) = α T ( z ) ¯ T z g is the square of the Brunt-V¨ais¨al¨a fre-quency ( g is the magnitude of gravity and is assumed to be constant). The microscopicdiffusion coefficients ν (the viscosity) and κ T (the thermal conductivity) are both assumedto be constant. In the Sun, near the radiative–convective interface, E ν ≃ × − and Pr ≃ × − (Gough, 2007).As mentioned earlier, the transition between the model radiative zone and convectionzone is measured by the behavior of α T ( z ), which goes from 0 for z > h to a finite value for z < h . This can be modeled in a non-dimensional way through the Rossby number, whichwe assume is of the form: Ro( z ) = Ro rz (cid:20) (cid:18) h − z ∆ (cid:19)(cid:21) , (3)where Ro rz is constant. The lengthscale ∆ may be thought of as the thickness of the “over-shoot” region near the base of the convection zone, but in practise is mostly used to ensurecontinuity and smoothness of the background state through the tanh function.In what follows, we further restrict our study to an “axially” symmetric ( ∂/∂x = 0),steady-state ( ∂/∂t = 0) problem. Within the radiative zone, the nonlinear terms u · ∇ u and u · ∇ T are assumed to be negligible. Within the convection zone on the other hand,anisotropic turbulent stresses are thought to drive the observed differential rotation . Wemodel this effect in the simplest possible way, replacing the divergence of the stresses by alinear relaxation towards the observed convection zone profile: u · ∇ u → F turb = u − u cz τ , (4) Note that for slowly rotating stars, such as the Sun, the direct generation of meridional flows byanisotropic stresses is a much weaker effect in the bulk of the convection zone. We neglect it here. u cz = u cz ( y, z ) e x and the function u cz ( y, z ) models the observed azimuthal velocityprofile in the solar convection zone. This is analogous to the prescription used by Spiegel& Bretherton (1968) in their study of the effect of a convection zone on solar spin-down,although in their model the convection zone was not differentially rotating. The dimension-less relaxation timescale τ can be thought of, for example, as being of the same order ofmagnitude as the convective turnover time divided by the rotation period. It is modeled as: τ − ( z ) = Λ2 (cid:20) (cid:18) z − h ∆ (cid:19)(cid:21) . (5)Note that in the real solar convection zone, τ varies by orders of magnitudes between thesurface ( τ ∼ − ) and the bottom of the convection zone ( τ ∼ u cz ( y, z ): u cz ( y, z ) = U ( z )2 (cid:20) (cid:18) z − h ∆ (cid:19)(cid:21) e iky = ˆ u cz ( z ) e iky , U ( z ) = U ( h ) + S ( z − h ) , (6)where k = 2 to match the equatorial symmetry of the observed solar rotation profile. Thetanh function once again is merely added to guarantee continuity of the forcing across theovershoot layer. The function U ( z ) describes the imposed “vertical shear”, and is for sim-plicity taken to be a linear function of z . If U ( h ) = 0, the forcing effectively vanishes atthe base of the convection zone. If U ( h ) = 0 on the other hand, a strong azimuthal shearis forced at the interface. The observed solar rotation profile appears to be consistent with U ( h ) and S both being non-zero (and of the order of 0.1, although since we are studying alinear problem, the amplitude of the forcing is somewhat irrelevant). Note that if S = 0 theforcing velocity u cz has zero vorticity.Finally, the observed asphericity in the temperature profile is negligible in the solarconvection zone; this is attributed to the fact that the turbulent convection very efficientlymixes heat both vertically and horizontally. We model this effect as: u · ∇ T → − D ( z ) ∇ T , (7)where the turbulent heat diffusion coefficient is modeled as D ( z ) = D (cid:20) (cid:18) z − h ∆ (cid:19)(cid:21) , (8)and thus vanishes beneath the overshoot layer. We will assume that the diffusion timescale1 /D (in non-dimensional units) is much smaller than any other typical timescale in thesystem ( D ≫ q ( y, z ) = ˆ q ( z ) e iky for each of the unknown quantities yields − v = E ν (cid:18) d ˆ u d z − k ˆ u (cid:19) − ˆ u − ˆ u cz τ ,2ˆ u = − ik ˆ p + E ν (cid:18) d ˆ v d z − k ˆ v (cid:19) − ˆ vτ ,0 = − dˆ p d z + Ro ( z ) ˆ T + E ν (cid:18) d ˆ w d z − k ˆ w (cid:19) − ˆ wτ , ik ˆ v + d ˆ w d z = 0 ,ˆ w = (cid:18) E ν Pr + D ( z ) (cid:19) d ˆ T d z − k ˆ T ! , (9)Note that as required, the imposed forcing term drags the fluid in the azimuthal direction:for τ →
0, ˆ u → ˆ u cz in the convection zone. The meridional flows ˆ v and ˆ w on the otherhand are generated by the y − component of the Coriolis force and by mass conservationrespectively (the essence of gyroscopic pumping, see McIntyre 2007).We now proceed to solve these equations to gain a better understanding of the meridionalflows and their degree of penetration into the radiative zone below. We use a dual approach,solving these equations first analytically under various limits, and then exactly using a simpleNewton-Raphson-Kantorovich (NRK) two-point boundary value algorithm. The analyticalapproximations yield predictions for the relevant scalings of the solutions in terms of thegoverning parameters (an in particular, the Ekman number and the Rossby number) whichare then confirmed by the exact numerical solutions. Although this limit is not a priori relevant to the physics of the solar interior, we beginby studying the case of an unstratified region, setting Ro rz = 0 (in this case, the thermalenergy equation can be discarded). This simpler problem, as we shall demonstrate, containsthe essence of the problem.In order to find analytical approximations to the solutions, we solve the governingequations separately in the convective zone and in the radiative zone. At this point, itmay be worth pointing out that in the unstratified case, the nomenclatures “convective”and “radiative” merely refer to regions which respectively are and are not subject to theadditional forcing. 8 –We assume that the transition region is very thin . In this case, τ − = Λ for z > h while τ − = 0 for z < h . Similarly, ˆ u cz ( z ) = U ( h ) + S ( z − h ) in the convection zonewhile ˆ u cz ( z ) = 0 in the radiative zone. Once obtained, the solutions are patched at theradiative–convective interface. In the convection zone, the equations reduce to − v = − Λ(ˆ u − ˆ u cz ) ,2ˆ u = − ik ˆ p − Λˆ v ,0 = − dˆ p d z − Λ ˆ w , ik ˆ v + d ˆ w d z = 0 , (10)where we have neglected the viscous dissipation terms in favor of the forcing terms since E ν ≪ Λ for all reasonable solar parameters. Combining them yieldsd ˆ w d z = k Λ ˆ w + 2 ik Λ4 + Λ dˆ u cz d z . (11)This second-order ordinary differential equation for ˆ w ( z ) suggests the introduction of a newlengthscale δ = √ k Λ , (12)so that the general solution to (10) isˆ w ( z ) = Ae z/δ + Be − z/δ − iSk Λ ,ˆ v ( z ) = − ikδ (cid:2) Ae z/δ − Be − z/δ (cid:3) ,ˆ u ( z ) = ˆ u cz ( z ) − ik Λ δ (cid:2) Ae z/δ − Be − z/δ (cid:3) ,ˆ p ( z ) = − ik ˆ u cz ( z ) − δ Λ (cid:2) Ae z/δ − Be − z/δ (cid:3) . (13) More precisely, ∆ ≪ E / ν , see Section 2.3.5. The original order of the system is much reduced in the convection zone since we ignored the effect ofviscous terms there. A and B are integration constants which must be determined by apply-ing boundary conditions (at z = 1) and matching conditions (at z = h ). Note from theˆ u -equation that the actual rotation profile approaches the imposed (observed) profile ˆ u cz provided A and B tend to 0, or when Λ ≫ δ → /k ). In the radiative region, the equations reduce to − v = E ν (cid:18) d ˆ u d z − k ˆ u (cid:19) ,2ˆ u = − ik ˆ p ,dˆ p d z = 0 , ik ˆ v + d ˆ w d z = 0 , (14)if we neglect viscous stresses in both y and z components of the momentum equation. Notethat viscous stresses in the x − equation cannot be dropped since they are the only forcebalancing the Coriolis force. These equations are easily solved:ˆ p ( z ) = p rz ,ˆ u ( z ) = − ik p rz ,ˆ v ( z ) = − iE ν k p rz ,ˆ w ( z ) = w rz − E ν k p rz z , (15)where p rz and w rz are two additional integration constants. Here, we recover the standardTaylor-Proudman constraint where in the absence diffusion or any other stresses, the velocitymust be constant along the rotation axis (here, e z ); in the limit E ν →
0, ˆ u ( z ) and ˆ w ( z )become independent of z , while ˆ v ( z ) → p rz and w rz form, together with A and B , a set of 4 unknownconstants which are determined by application of boundary and matching conditions. Sincewe have neglected viscous effects in the convection zone, we cannot require any boundary ormatching condition on the horizontal fluid motions. On the other hand, we are allowed toimpose impermeability ˆ w = 0 at the surface ( z = 1) and at the bottom ( z = 0). Moreover, 10 –we request the continuity of the radial (vertical) velocity and of the pressure at the interface( z = h ). Applying these conditions yields the set of equations w rz = 0 , Ae h/δ + Be − h/δ − iSk Λ = w rz − E ν k p rz h , − ik U ( h ) − δ Λ (cid:2) Ae h/δ − Be − h/δ (cid:3) = p rz , Ae /δ + Be − /δ − iSk Λ = 0 , (16)which have the following solution for A and B : A = E ν hk i U ( h ) + iSk Λ h − e (1 − h ) /δ h E ν k hδ Λ ii e h/δ (cid:2) − E ν k hδ Λ (cid:3) − e (2 − h ) /δ (cid:2) E ν k hδ Λ (cid:3) , B = 2 iSk Λ e /δ − Ae /δ . (17)These can be substituted back into (13) to obtain the meridional flow velocities in theconvection zone. While the exact form of A and B are not particularly informative, wenote that in the limit S = 0 (i.e. the forcing velocity has no azimuthal vorticity), both A and B scale as E ν . This implies that the amplitude of meridional flows everywhere in thesolar interior scales like E ν (even in the convection zone). The physical interpretation ofthis somewhat surprising limit is discussed in Section 2.3.4, but turns out to be of academicinterest only (Section 2.3.5).When S = 0 then A and B are of order S/k
Λ in the convection zone regardless of theEkman number, and respectively tend to A = 2 iSk Λ 1 − e (1 − h ) /δ e h/δ − e (2 − h ) /δ + O (E ν ) , B = 2 iSk Λ 1 − e ( h − /δ e − h/δ − e ( h − /δ + O (E ν ) , (18)as E ν →
0. This implies that ˆ w is of order S/k
Λ in the convection zone. Since significantflows are locally generated, one may reasonably expect a fraction of the forced mass flux topenetrate into the lower region, especially in this unstratified case.Using (17) in (16), solving for p rz , then plugging p rz into (15), we find that the generalexpression for ˆ w ( z ) in the radiative zone z ∈ [0 , h ] isˆ w ( z ) = − iE ν k (cid:18) U ( h ) + cosh((1 − h ) /δ ) − − h ) /δ ) δS (cid:19) z . (19) 11 –This implies that only a tiny fraction of the large mass flux circulating in the convectionzone actually enters the radiative zone. Instead, the system adjusts itself in such a way asto ensure that most of the meridional flows return above the base of the convection zone.We now compare this a priori counter-intuitive analytical result with exact numericalsolutions of the governing equations. The numerical solutions were obtained by solving(9) for Ro rz = 0 (unstratified case), and are uniformly calculated in the whole domain (i.e.there is nothing special about the interfacial point z = h ). The boundary conditions used areimpermeable boundary conditions at the top and bottom of the domain for ˆ w , and stress-freeboundary conditions for ˆ u and ˆ v .In Figure 2 we compare numerical and analytical solutions for ˆ w ( z ), in a case where theforcing function parameters are ∆ = 10 − , Λ = 10, U ( h ) = 0 and S = 1, for four valuesof the Ekman number. The analytical solution is described by equations (13) and (17) (forthe convection zone) and (15) (for the radiative zone). As E ν → w ( z ) ∝ z E ν in theradiative zone. The convection zone solution is also well-approximated in this case by theanalytical formula.A full 2D visualization of the flow for E ν = 10 − but otherwise the same governingparameters is shown in Figure 3. This figure illustrates more clearly the fact that themeridional flows are negligible below the interface, and mostly return within the convectionzone. Note that given our choice of the forcing function u cz ( y, z ) ∝ cos(2 y ), the inducedCoriolis force does not vanish at y = 0 or y = π (the “poles”). This explains why themeridional flows apparently cross the polar axis in this simple model. This is merely ageometric effect: in a true spherical geometry the forcing azimuthal velocity u cz ( r, θ ) would benull at the poles, and meridional flows cannot cross the polar axis. More realistic calculationsin spherical geometry are discussed in Section 3.Finally, it is interesting to note that the analytical solution for the azimuthal velocityˆ u ( z ) exhibits a “discontinuity” across the base of the convection zone, which tends toˆ u ( h + ) − ˆ u ( h − ) = (cid:18) ikδ Λ + ikδ Λ2 (cid:19) (cid:0) Ae h/δ − Be − h/δ (cid:1) (20)as E ν →
0. The numerical solutions of course are continuous, but the continuity is onlyassured by the viscosity in the system (in the y − direction) and the fact that the overshootlayer depth is finite. This is shown in Figure 4, together with a comparison of the numericalsolutions with the analytical solution, again confirming the analytical approximation derived. but a posteriori obvious, see Section 2.3.4
12 –Fig. 2.— Numerical (solid) and analytical (dashed) solutions for | ˆ w ( z ) | , in the case of astress-free bottom boundary. From the uppermost to lowermost curves, E ν = 10 − , 10 − ,10 − and 10 − respectively, confirming the analytical scaling that ˆ w ( z ) ∝ E ν z in the radiativezone while ˆ w ( z ) becomes independent of E ν in the convection zone. These solutions wereobtained with forcing defined by the parameters ∆ = 10 − , Λ = 10, U ( z ) = S ( z − h ) and S = 1.This highlights another and equally a priori counter-intuitive property of the system:the value of u rz in the radiative zone is markedly different from the imposed ˆ u cz ( h ) = U ( h )at the interface: u rz = U ( h ) + ikδ Λ2 (cid:0) Ae h/δ − Be − h/δ (cid:1) . (21)Hence, even if the imposed differential rotation is exactly 0 at the radiative–convectiveinterface (as it is the case in the simulation presented in Figure 4 since U ( z ) = S ( z − h )),a large-scale latitudinal shear measured by u rz may be present in the radiative zone, asillustrated in Figure 3. This shows that the propagation of the azimuthal shear into theradiative zone is non-local (i.e. does not rely on the presence of shear at the interface), andis instead communicated by the long-range pressure gradient. 13 –Fig. 3.— 2D visualization of the flow for E ν = 10 − , in the case of a stress-free bottomboundary. Shown as solid and dotted line respectively are linearly spaced streamlines ofcounter-clockwise and clockwise meridional flows. As predicted, the flows appear to returnentirely within the convection zone and carry a negligible mass flux into the radiative zone.Meanwhile the azimuthal velocity (ˆ u ) as displayed in the filled contours is constant alongthe rotation axis ( z − axis) below the interface ( z = h = 0 . − ,Λ = 10, U ( z ) = S ( z − h ) and S = 1, as in Figure 2. 14 –Fig. 4.— Numerical (solid) and analytical (dashed) solutions for ˆ u ( z ), in the case of astress-free bottom boundary. From the lowermost to uppermost curves, E ν = 10 − , 10 − ,10 − and 10 − respectively, confirming that ˆ u ( z ) tends to a constant in the radiative zone,while sustaining a finite discontinuity at the radiative-convective interface ( z = h = 0 . − , Λ = 10, U ( z ) = S ( z − h ) and S = 1, as in Figure 2. The stress-free bottom boundary conditions studied in the previous Section are at firstglance the closest to what one may expect in the real Sun, where the “bottom” boundarymerely represents the origin of the spherical coordinate system. However, let us now explorefor completeness (and for further reasons that will be clarified in the next Section) the caseof no-slip bottom boundary conditions.When the lower boundary is a no-slip boundary, the nature of the solution in the wholedomain changes. This change is induced by the presence of an Ekman boundary layer,which forms near z = 0. Just above the boundary layer, in the bulk of the radiative zone,the solution described in 2.3.2 remains valid. However, matching the bulk solution with theboundary conditions can no longer be done directly; one must first solve for the boundarylayer dynamics to match the bulk solution with the boundary conditions across the boundarylayer. This is a standard procedure (summarized in Appendix A for completeness), and leadsto the well-known “Ekman jump” relationship between the jump in ˆ w ( z ) and the jump in 15 –ˆ u ( z ) across the boundary layer:ˆ u bulk − ˆ u (0) = 2 ik E − / ν ( ˆ w bulk − ˆ w (0)) . (22)By impermeability, ˆ w (0) = 0. Moreover, by assuming that the total angular momentumof the lower boundary is the same as that of the convection zone, we require that ˆ u (0) = 0.Meanwhile, ˆ u bulk = u rz and ˆ w bulk = w rz in the notation of equation (15). So finally, forno-slip boundary conditions, we simply replace the impermeability condition ( w rz = 0) in(16) by u rz = 2 ik E − / ν w rz , (23)and solve for the unknown constants A , B , w rz and p rz as before.The exact expressions for the resulting integration constants A and B are now slightlydifferent from those given in (17), but are without particular interest. However, it can beshown that they have the same limit as in the stress-free case as E ν → S = 0). Thisimplies that the meridional flows driven within the convection zone , in the limit E ν → S = 0, are independent of the boundary condition selected at the bottom of theradiative zone. However, we now have the following expression for ˆ w ( z ) in the bulk of theradiative zone:ˆ w ( z ) = − (cid:18) k E / ν + E ν k z (cid:19) p rz ,= − ik E / ν (cid:18) U ( h ) + cosh((1 − h ) /δ ) − − h ) /δ ) δS (cid:19) + O (E ν ) z , (24)which has one fundamental consequence: the amplitude of the flows allowed to penetrate intothe radiative zone is now of order E / ν instead of being O (E ν ). This particular statement isactually true even if S = 0, although in that case both convection zone and radiative zoneflows scale with E / ν .Figure 5 shows a comparison between the approximate analytical formula and the nu-merical solution for the same simulations as in Figure 2, but now using no-slip bottomboundary conditions. For ease of comparison, the results from the stress-free numericalsimulations (for exactly the same parameters) have also been drawn, highlighting the muchlarger amplitude of the meridional flows down-welling into the radiative zone in the no-slipcase, and their scaling with E / ν . Figure 6 shows an equivalent 2D rendition of the solution,and illustrates the presence of large-scale mixing in the bulk of the radiative zone when thebottom-boundary is no-slip. 16 –Fig. 5.— Numerical (solid) and analytical (dashed) solutions for | ˆ w ( z ) | , in the case of ano-slip bottom boundary. From the uppermost to lowermost curves (as seen in the radiativezone), E ν = 10 − , 10 − , 10 − and 10 − respectively, confirming the analytical scaling thatˆ w ( z ) ∝ E / ν in the radiative zone while ˆ w ( z ) becomes independent of E ν in the convectionzone. These solutions were obtained with forcing defined by the parameters ∆ = 10 − ,Λ = 10, U ( z ) = S ( z − h ) and S = 1, as in Figure 2. For comparison, the previoussimulations with stress-free bottom boundary, for the same parameters, are shown as dottedlines. The various sets of solutions derived above can be physically understood in the followingway. Let us first discuss the solution in the convection zone. In the limit where u cz ( y, z ) isindependent of z (equivalently, S = 0), the azimuthal ( x − ) component of the vorticity ofthe forcing is zero. In that case there is no injection of x − vorticity into the system asidefrom that induced in the viscous boundary layers, and the amplitude of the meridional flowsgenerated in the convection zone scales with E ν . This limit is somewhat academic in thecase of the Sun, however given the observed rotation profile (see also Section 2.3.5). When S = 0, the amplitude of the induced meridional flows in the convection zone scales linearlywith S and is independent of viscosity.In the radiative zone, the Taylor-Proudman constraint enforces invariance of the flowvelocities along the rotation axis, except in regions where other forces balance the Coriolis 17 –Fig. 6.— The same as for Figure 3 but for no-slip boundary conditions. The Ekman layernear the lower boundary is clearly visible. For ease of comparison, the same streamlines areshown in the two plots. The two figures illustrates how the nature of the lower boundarycondition influences the mass flux through the radiative zone. 18 –force. In the non-magnetic, unstratified situation discussed in the two previous sections,the only agent capable of breaking the Taylor-Proudman constraint are viscous stresses,which are only significant in two thin boundary layers: one right below the convection zoneand the other one near the bottom boundary. These two layers are the only regions whereflows down-welling into the radiative zone are allowed to return to the convection zone . Thequestion then remains of what fraction of the mass flux entering the radiative zone returnswithin the upper Ekman layer, and what fraction returns within the lower Ekman layer. Thelatter, of course, permits large-scale mixing within the radiative zone.In the first case studied, the bottom boundary was chosen to be stress-free. This nat-urally suppresses the lower viscous boundary layer so that the only place where flows areallowed to return is at the radiative-convective interface. As a result, only a tiny fraction ofthe mass flux penetrates below z = h , and the turnover time of the remaining flows withinthe radiative zone is limited to a viscous timescale of the order of 1 / E ν Ω ⊙ .Following this reasoning, we expect and indeed find quite a different behavior when thebottom boundary is no-slip. In that case viscous stresses within the lower boundary layerbreak the Taylor-Proudman constraint and allow a non-zero mass flux (of order E / ν ) toreturn near z = 0. This flow then mixes the entire radiative zone as well, with an overallturnover time of order of 1 / E / ν Ω ⊙ (in dimensional units).To summarize, in this unstratified steady-state situation, the amount of mixing inducedwithin the radiative zone by convective zone flows depends (of course) on the amplitude ofthe convection zone forcing, but also on the existence of a mechanism to break the Taylor-Proudman constraint somewhere within the radiative zone. That mechanism is needed inorder to allow down-welling flows to return to the convection zone. But more crucially, thisphenomenon implies that the dynamics of the lower boundary layer entirely control the massflux through the system.Here, we studied the case of viscous stresses only. One can rightfully argue that thereare no expected “solid” boundaries in a stellar interior and that the overall behavior of thesystem should be closer to the one discussed in the stress-free case than the no-slip case. However, we chose here to study viscous stresses simply because they are the easiest availableexample.
In real stars viscous stresses are likely to be negligible compared with a varietyof other possible stresses: turbulent stresses at the interface with another convection zone,magnetic forces, etc. Nevertheless, these stresses will play a similar role in allowing flows tomix the radiative zone if they become comparable in amplitude with the Coriolis force, andhelp break the Taylor-Proudman constraint. This issue is discussed in more detail in Section4. 19 –
Before moving on to the more realistic stratified case, note that this unstratified systemholds one final subtlety. In all simulations presented earlier, the overshoot layer depth wasselected to be very small – and in particular, smaller than the Ekman layer thickness. Inthat case, the transition in the forcing at the base of the convection zone is indeed close tobeing a discontinuity, and the analytical solutions presented in Sections 2.3.1 and 2.3.2 area good fit to the true numerical solution.In the Sun, the overshoot layer depth ∆ is arguably always thicker than an Ekmanlengthscale E / ν . When this happens, the solution “knows” about the exact shape of theforcing function within the transition region, and therefore depends on it. This limit turnsout to be rather difficult to study analytically, and since in the case of the Sun we do notknow the actual profile of τ − ( z ), there is little point in exercise anyway.We can explore the behavior of the system numerically, however, for the profile τ − ( z )discussed in equation (5), in the limit where ∆ > E / ν . The example for which this effectmatters the most is the somewhat academic limit where S = 0 in the convection zone, but U ( h ) = 0. In this case, the asymptotic analysis predicts that the meridional flow amplitudesare O (E / ν ) in both the convection zone and in the radiative zone for the no-slip case. Wesee in Figure 7 that this is indeed the case in simulations where ∆ ≪ E / ν . However, whenthe overshoot thickness is progressively increased and becomes larger than the Ekman layerthickness, the amplitude of the meridional circulation in the convection zone is no longer O (E / ν ) but much larger. Meanwhile, the scaling of the radiative zone solutions with E ν remain qualitatively correct. The difference with the analytical solution in the convectionzone can simply be attributed to the fact that when the system knows about the shear within the overshoot layer the limit S = 0 is no longer relevant. While the previous section provides interesting insight into the problem, notably on therole of the Taylor-Proudman constraint, we now move to the more realistic situation wherestratification plays a role in the flow dynamics. In this section, we generalize our Cartesianstudy to take into account the stratification of the lower region (Ro rz = 0). For this purpose,we go back to studying the full system of equations (9). As before, we first find approximateanalytical solutions to derive the overall scaling of the solutions with governing parameters,and then compare them to the full numerical solutions of (9). The analytical solutions areobtained by solving the system in the convective region and radiative region separately, and 20 –Fig. 7.— Comparison of numerical simulations (solid lines) with analytical prediction(dashed line) for a no-slip bottom boundary, for E ν = 10 − , and forcing functions definedwith Λ = 10, S = 0 and U ( h ) = 1. The three numerical solutions are obtained for variousvalues of the overshoot layer depth: from lowermost to uppermost curves (as seen in the con-vective zone), ∆ = 10E / ν , E / ν , and 0 . / ν . The analytical solution assumes an infinitelythin overshoot layer and is therefore independent of ∆. Note that the analytical solutionin the convection zone is only a good approximation to the true solution if ∆ ≪ E / ν . Theoverall scalings in the radiative zone, however, are preserved.matching them at z = h . The equations in the convection zone are now given by − v = − Λ(ˆ u − ˆ u cz ) ,2ˆ u = − ik ˆ p − Λˆ v ,0 = − dˆ p d z − Λ ˆ w ,ˆ w = D d ˆ T d z − k ˆ T ! , 21 – ik ˆ v + d ˆ w d z = 0 , (25)where we have assumed that D ≫ E ν / Pr. Eliminating variables one by one yields the sameequation for ˆ w ( z ) in the convection zone as before (11), as well as a Poisson equation for ˆ T once ˆ w is known. The solutions are then as (13), together withˆ T ( z ) = T e kz + T e − kz + δ (cid:2) Ae z/δ + Be − z/δ (cid:3) D (1 − δ k ) + 2 iSk Λ D , (26)where the integration constants T and T remain to be determined. For the sake of analyticalsimplicity, we will assume that D ≫ The radiative zone equations are now − v = E ν (cid:18) d ˆ u d z − k ˆ u (cid:19) ,2ˆ u = − ik ˆ p ,0 = − dˆ p d z + Ro ( z ) ˆ T ,ˆ w = E ν Pr d ˆ T d z − k ˆ T ! , ik ˆ v + d ˆ w d z = 0 , (27)and can be combined to yieldd ˆ u d z − k (cid:18) rz2 (cid:19) d ˆ u d z + k PrRo rz2 u = 0 , (28)and similarly for ˆ T . The characteristic polynomial is( λ − k ) (cid:18) λ − PrRo rz2 k (cid:19) = 0 , (29)with solutions ± λ = ± k , ± λ = ±√ Pr Ro rz k . (30) 22 –These solutions are the same as those presented in Paper I, and will be referred to as the“global-scale” mode and the thermo-viscous mode respectively. Note that here, λ corre-sponds to k in Paper I.In this steady-state study, the quantity λ summarizes the effect of stratification. It isimportant to note that it contains information about the rotation rate of the star as well asthe Prandtl number, in addition to the buoyancy frequency. If λ ≪
1, then the thermo-viscous mode essentially spans the whole domain: the system appears to be “unstratified”,and is again dominated by the Taylor-Proudman constraint. On the other hand, if λ ≫ /λ as a result of the strong stratification of the system. The Taylor-Proudman constraint is irrelevant in this limit, since the magnitude of the buoyancy force ismuch larger than that of the Coriolis force.The calculation above was made in the limit where the viscous terms in the latitudi-nal and radial components of the momentum equation are discarded. Paper I shows thattwo additional Ekman modes are also present if they are instead kept. By analogy withthe unstratified case, we expect that these Ekman modes do not influence the solution forstress-free boundary conditions, but that additional care must be taken for no-slip boundaryconditions.Note that the equation for ˆ w instead simplifies tod ˆ w d z = k PrRo rz2 w , (31)and similarly for ˆ v (i.e. both equations are only second order in z , and only contain thethermo-viscous mode). The radiative zone ( z ∈ [0 , h ]) solutions are nowˆ u ( z ) = u e kz + u e − kz + u e λ z + u e − λ z ,ˆ v ( z ) = E ν k − λ ) (cid:2) u e λ z + u e − λ z (cid:3) ,ˆ w ( z ) = − ik E ν k − λ ) λ (cid:2) u e λ z − u e − λ z (cid:3) ,ˆ p ( z ) = − ik (cid:2) u e kz + u e − kz + u e λ z + u e − λ z (cid:3) ,ˆ T ( z ) = − ik Ro (cid:2) ku e kz − ku e − kz + λ u e λ z − λ u e − λ z (cid:3) . (32)where the 4 constants { u i } i =1 , are integration constants, to be determined. 23 – We now proceed to match the solutions in the two regions, assuming stress-free boundaryconditions near the lower boundary. Since there are in total 8 unknown constants (including A , B , T and T from the convection zone solution and { u i } i =1 , from the radiative zonesolution), we need a total of 8 matching and boundary conditions.At the lower boundary ( z = 0) we take ˆ w = dˆ u/ d z = 0; this condition in turn impliesthat ˆ T = 0. At the surface ( z = 1), we take as before ˆ w = 0, and select in addition ˆ T = 0.We then need 4 matching conditions across the interface: these are given by the continuityof ˆ w , ˆ p , ˆ T and d ˆ T / d z . Note that it is important to resist the temptation of requiring thecontinuity of ˆ v , since viscous stresses have been neglected in the analytical treatment of the y − component of the momentum equation in both radiative and convective zones. Moreover,we know that in the unstratified limit, ˆ u actually becomes discontinuous at the interface inthe limit E ν →
0. Since we expect the stratified solution to tend to the unstratified oneuniformly as Ro rz →
0, we cannot require the continuity of ˆ u at the interface .The equations and resulting solutions for the integration constants are fairly compli-cated. The most important ones are reported in the Appendix B for completeness, and areused to justify mathematically the following statements: • In the limit of Ro rz →
0, we find as expected that the solutions uniformly tend tothe unstratified solution summarized in equations (13), (15), and (17). Indeed, inthat case λ → • In the strongly stratified case (defined as λ ≫ k ), as described earlier, ˆ w in the radia-tive zone decays exponentially with depth on a lengthscale 1 /λ , with an amplitudewhich scales as E ν /λ . The flows are therefore very strongly suppressed, and return tothe convection zone within a small thermo-viscous layer. Note that E ν /λ = Ra − / where Ra is the usually defined Rayleigh number.The two limits are illustrated in Figure 8, which shows the numerical solution to (9) fortwo values of the Rossby number Ro rz , but otherwise identical parameters. In the stronglystratified limit ( λ = 10, using Pr = 0.01 and Ro rz = 10 ) we see that the solution decaysexponentially below the interface, with an amplitude which scales as E ν /λ as predictedanalytically. In the weakly stratified case ( λ = 0 .
1, using Pr = 0.01 and Ro rz = 1) thesolution tends to the unstratified limit and scales as E ν z . a fact which is again only obvious in hindsight
24 –Fig. 8.— Numerical solutions of (9) with the following parameters: ∆ = 0 .
01, Λ = 10, S = 0, U = 1, Pr = 0 .
01 and D = 10. Stress-free bottom boundary conditions are used.The solid lines correspond to the “strongly” stratified case with Ro rz = 10, with E ν = 10 − and 10 − for the top and bottom curves respectively. The dashed lines correspond to the“weakly” stratified case, with Ro rz = 1, with E ν = 10 − and 10 − for the top and bottomcurves respectively. Note that for k = 2, λ is simply equal to Pr / Ro rz . For comparison,the unstratified case (Ro rz = 0) is shown as dotted lines. At these parameters and with theseboundary conditions, Ro rz = 1 already belongs to the weakly stratified limit. By analogy with the previous section, we expect to recover the unstratified limit when λ →
0, so that ˆ w ( z ) ∝ E / ν in this no-slip case. In the strongly stratified limit on theother hand, the amplitude of the flows decays exponentially with depth below the interfaceas a result of the thermo-viscous mode and is negligible by the time they reach the lowerboundary. In that case, we do not expect the applied lower boundary conditions to affect thesolution, so that the scalings found in the strongly stratified limit with stress-free boundaryconditions should still apply: ˆ w ( z ) ∝ E ν /λ .These statements are verified in Figure 9. There, we show the results of a series ofnumerical experiments for no-slip boundary conditions where we extracted from the simu-lations the power α in the expression ˆ w ∝ E αν , and plotted it as a function of stratification( λ ). To do this, we integrated the solutions to equations (9) for the following parameters: 25 –∆ = 0 .
01, Λ = 10, S = 1, U ( h ) = 1, Pr = 0 .
01 and D = 10 and calculated ˆ w ( z = 0 .
5) for4 values of E ν : 10 − , 10 − , 10 − and 10 − . We estimated α by calculating the quantity α = log ˆ w ( z = 0 . , E ν = 10 − )ˆ w ( z = 0 . , E ν = 10 − ) (33)for the (E ν = 10 − , E ν = 10 − ) pair (diamond symbols) and similarly for the pairs (E ν =10 − , E ν = 10 − ) (triangular symbols) and (E ν = 10 − , E ν = 10 − ) (star symbols). In theweakly stratified limit ( λ → α → / λ ≫ α →
1, thus confirming our analysis. The transition between the two regimesappears to occur for slightly lower-than expected values of λ , namely 0.1 instead of 1.Fig. 9.— This figure shows the power α in the expression ˆ w ∝ E αν , as a function of λ (seemain text for detail). In the weakly stratified limit, α → / α → .
01, Λ = 10, S = 1, U ( h ) = 1, Pr = 0 .
01 and D = 10A final summary of our findings for the stratified case together with its implications formixing between the solar convection zone and the radiative interior, is deferred to Section4. There, we also discuss the consequences in terms of mixing in other stars. But first, wecomplete the study by releasing some of the simplifying assumptions made, and moving tomore realistic numerical solutions to confirm our simple Cartesian analysis. 26 –
3. A “solar” model
In this section we improve on the Cartesian analysis by moving to a spherical radiative–convective model. The calculations are thus performed in an axisymmetric spherical shell,with the outer radius r out selected to be near the solar surface, and the inner radius r in somewhere within the radiative interior. This enables us to gain a better understanding ofthe effects of the geometry of the system on the spatial structure of the flows generated. Inaddition, we use more realistic input physics in particular in terms of the background stratifi-cation, and no longer use the Boussinesq approximation for the equation of state. We expectthat the overall scalings derived in Section 2 still adequately describe the flow amplitudes inthis new calculation. However, the use of a more realistic background stratification adds anadditional complication to the problem: the background temperature/density gradients areno longer constant, so that the measure of stratification λ varies with radius (see Figure 10,for an estimate of λ in the Sun). This aim of this section is therefore to study the impactFig. 10.— Variation of λ /k =Pr . N/ Ω ⊙ in the Sun, as determined from Model S ofChristensen-Dalsgaard et al. (1996). The Prandtl number Pr is calculated using ModelS together with the formulae provided by Gough (2007) for the microscopic values of theviscosity ˆ ν and the thermal conductivity κ T (see also Garaud & Garaud, 2008). The insetzooms into the region near the base of the convection zone, which is the only region of theradiative zone where λ ≤ r → The spherical model used is analogous to the radiative-zone-only model presented inPaper I and described in detail (including the magnetic case) by Garaud & Garaud (2008).The salient points are repeated here for completeness, together with the added modificationsmade to include the “convective” region.We consider a spherical coordinate system ( r, θ, φ ) where the polar axis is aligned withrotation axis of the Sun. The background state is assumed to be spherically symmetric andin hydrostatic equilibrium. The background thermodynamical quantities such as density,pressure, temperature and entropy are denoted with bars (as ¯ ρ ( r ), ¯ p ( r ), ¯ T ( r ) and ¯ s ( r )respectively), and extracted from the standard solar model of Christensen-Dalsgaard et al. (1996). Perturbations to this background induced by the velocity field u = ( u r , u θ , u φ ) aredenoted with tildes. In the frame of reference rotating with angular velocity Ω ⊙ , in a steadystate, the linearized perturbation equations become ∇ · ( ¯ ρ u ) = 0 ,2 ¯ ρ Ω ⊙ e z × u = −∇ ˜ p + ˜ ρ g + f ∇ · Π − ¯ ρ Ω ⊙ u − u cz τ ( r ) ,¯ ρ ¯ c p ¯ T ¯ N g u r = ∇ · h ( f ¯ k T + R ⊙ Ω ⊙ D ( r )) ∇ ˜ T i ,˜ p ¯ p = ˜ ρ ¯ ρ + ˜ T ¯ T , (34)where ¯ c p is the specific heat at constant pressure, ¯ k T ( r ) = ¯ ρ ¯ c p ¯ κ T is the thermal conductivityin the solar interior, Π is the viscous stress tensor (which depends on the background viscosity¯ ν ) and g = − g ( r ) e r is gravity. Note that this set of equations is given in dimensional formhere, although the numerical algorithm used further casts them into a non-dimensionalform. Also note that both diffusion terms (viscous diffusion and heat diffusion) have beenmultiplied by the same factor f . This enables us to vary the effective Ekman numberE ν ( r ) = f E ⊙ ν = f ¯ ν ( r ) /R ⊙ Ω ⊙ while maintaining a solar Prandtl number at every radialposition. As a result, the quantity λ used in the simulations and represented in Figure 10is the true solar value (except where specifically mentioned).As in the Cartesian case, we model the dynamical effect of turbulent convection in theconvection zone through a relaxation to the observed profile in the momentum equation, anda turbulent diffusion in the thermal energy equation. The expressions for the non-dimensionalquantities τ ( r ) and D ( r ) are the same as in equations (5) and (8) with z replaced by r/R ⊙ ,and h = 0 .
713 instead of h = 0 . et al. u cz is selected to be u cz ( r, θ ) = r sin θ Ω cz ( θ ) e φ , (35)where Ω cz ( θ ) = Ω eq (cid:0) − a cos θ − a cos θ (cid:1) , (36)with a = 0 .
17 , a = 0 .
08 ,Ω eq π = 463 nHz , (37)which is a simple approximation to the helioseismically determined profile (Schou et al. eq is the observed equatorial rotation rate. As in Paper I, we finallyselect Ω ⊙ to be Ω ⊙ = Ω eq (cid:18) − a − a (cid:19) , (38)to ensure that the system has the same specific angular momentum as that of the imposedprofile u cz ( r, θ ).The computational domain is a spherical shell with the outer boundary located at r out = 0 . R ⊙ . It is chosen to be well-below the solar surface to avoid complications relatedto the very rapidly changing background in the region r > . R ⊙ . The position of the lowerboundary will be varied.The upper and lower boundaries are assumed to be impermeable. The upper boundaryis always stress-free, while the lower boundary is assumed to be either no-slip or stress-freedepending on the calculation. In the no-slip case, the rotation rate of the excluded core isan eigenvalue of the problem, calculated in such a way as to guarantee that the total torqueapplied to the core is zero. Finally, the boundary conditions on temperature are selected insuch a way as to guarantee that ∇ ˜ T = 0 outside of the computational domain, as in Garaud& Garaud (2008). We verified that the selection of the temperature boundary conditionsonly has a qualitative influence on the results, and doesn’t affect the scalings derived.The numerical method of solution is based on the expansion of the governing equationsonto the spherical coordinate system, followed by their projection onto Chebishev polynomi-als T n (cos θ ), and finally, solution of the resulting ODE system in r using a Newton-Raphson-Kantorovich algorithm. The typical solutions shown have 3000 meshpoints and 60-80 Fouriermodes. For more detail, see Garaud (2001) and Garaud & Garaud (2008). 29 – We first consider the artificial limit of weak stratification. In the following numericalexperiment, we use the available solar model background state, but divide the buoyancyfrequency ¯ N by 10 everywhere in the computational domain (all other background quantitiesremain unchanged). As a result, the new value of λ in the domain is artificially reducedfrom the one presented in Figure 10 by 10 , and is everywhere much smaller than one. Theposition of the lower boundary is arbitrarily chosen to be at r in = 0 . R ⊙ .Two sets of solutions are computed for no-slip lower boundary and for stress-free lowerboundary. Figure 11 is equivalent to Figure 5: it displays the radial velocity u r as a functionof radius near the poles (latitude of 80 ◦ ) for various values of f – in other words, E ν – andclearly illustrates the scalings of u r ∝ E / ν in the radiative zone for the no-slip case, and u r ∝ E ν for the stress-free case.Fig. 11.— Vertical velocity at 80 ◦ latitude in units of R ⊙ Ω ⊙ for an artificially weaklystratified simulation (where ¯ N was uniformly divided by 10 everywhere). The solid linesshow three simulations for f = 10 , f = 10 and f = 10 (from top to bottom) for theno-slip case. These correspond to E ν = 2 × − , E ν = 2 × − and E ν = 2 × − at the baseof the convection zone respectively, hence showing how u r ∝ E / ν . The dotted lines showsimulations with stress-free boundary conditions for the same parameters, showing u r ∝ E ν .In this calculation the overshoot depth ∆ was selected to be 0.01 R ⊙ , and Λ = 10. The valueof D is irrelevant in this very weakly stratified simulation. 30 –Figure 12 illustrates the geometry of the flow in both no-slip and stress-free cases for f = 10 (which corresponds to an Ekman number near the radiative–convective interfaceof about 2 × − ). The geometrical pattern of the flows observed within the convectionzone show a single-cell, with poleward flows near the surface and equatorward flows nearthe bottom of the convection zone. Below the convection zone we note the presence of threedistinct regions: the polar region, a Stewardson layer region (at the tangent cylinder) andan equatorial region. Flows within the equatorial region are weak regardless of the lowerboundary conditions. In the stress-free case, even in the tangent cylinder the flows tendto return mostly within the convection zone. If the lower boundary is no-slip on the otherhand, flows within the tangent cylinder are stronger, although the effect is not as obvious asin the Cartesian case because of the anelastic mass conservation equation used here.Fig. 12.— Normalized angular velocity ( ˜Ω / Ω eq ) and streamlines solutions to equations (34)for an artificially weak stratification (see Figure 11), and for f = 10 (corresponding toE ν = 2 × − at the base of the convection zone). On the left, we show the solutionwith no-slip lower boundary conditions, and on the right the stress-free solution. Dottedlines represent clockwise flows, and solid lines counter-clockwise flows. In this calculation,the overshoot layer depth ∆ was selected to be 0.01 R ⊙ , and Λ = 10. The value of D isirrelevant in this very weakly stratified simulation. 31 – Let us now consider the case of a true solar stratification. Since λ increases rapidlywith depth beneath the convection zone (from 0 to about 10 in the case of the Sun), theradius at which λ ≃ r for short) plays a special role: we expect the dynamics of thesystem to depend on the position of the lower boundary r in compared with r ≃ . R ⊙ .This is indeed observed in the simulations, as shown in Figure 13. If r in > r , then λ < u r ∝ E ν if the lower boundary is stress-free, and u r ∝ E / ν if the lower boundary is no-slip).On the other hand, if r in < r then the flows are strongly quenched by the stratificationbefore they reach the lower boundary. As a result, the radial velocities scale with E ν /λ regardless of the applied boundary conditions.The implications of this final result, namely the importance of the location of the stressesinvolved in breaking the Taylor-Proudman constraint in relation to the radius at which λ ≃
1, are discussed in Section 4.3.
4. Implications of this work for solar and stellar mixing4.1. Context for stellar mixing
The presence of mixing in stellar radiative zones has long been inferred from remainingdiscrepancies between models-without-mixing and observations (see Pinsonneault 1997 for areview). The most commonly used additional mixing source is convective overshoot, wherebystrong convective plumes travel beyond the radiative–convective interface and cause intensebut very localised (both in time and space) mixing events (Brummell, Clune & Toomre,2002). The typical depth of the layer thus mixed, the “overshoot layer”, is assumed to be asmall fraction of a pressure scaleheight in most stellar models.A related phenomenon is wave-induced mixing (Schatzman, 1996). While most of theenergy of the impact of a convective plume hitting the stably stratified fluid below is con-verted into local buoyancy mixing, a fraction goes into the excitation of a spectrum of gravitywaves, which may then propagate much further into the radiative interior. Where and whenthe waves eventually cause mixing (either through mutual interactions, thermal dissipationor by transferring momentum to the large-scale flow) depends on a variety of factors. It hasrecently been argued that the interaction of the gravity waves with the local azimuthal veloc-ity field (the differential rotation) would dominate the mixing process (Charbonnel & Talon 32 –2005), although this statement is not valid unless the wave-spectrum is near-monochromatic.For the typically flatter wave spectra self-consistently generated by convection, wave-inducedmixing has much more turbulent characteristics (Rogers, MacGregor & Glatzmaier 2008),and is again fairly localized below the convection zone.Mixing induced by large-scale flows comes in two forms: turbulent mixing resultingfrom instabilities of the large-scale flows, and direct transport by the large-scale flows them-selves. The former case is the dominant mechanism in the early stages of stellar evolutionwhen the star is undergoing rapid internal angular-momentum “reshuffling” caused by ex-ternal angular-momentum extraction (disk-locked and/or jet phase, early magnetic breakingphase). In these situations, regions of strong radial angular-velocity shear develop, whichthen become unstable and cause local turbulent mixing of both chemical species and an-gular momentum. Studies of these processes were initiated by the work of Endal & Sofia(1978). Later, Chaboyer & Zahn (1992) refined the analysis to consider the effect of stratifi-cation on the turbulence, and showed how this can differentially affect chemical mixing andangular-momentum mixing. Finally, Zahn (1992) proposed the first formalism which com-bines mixing by large-scale flows and mixing by (two-dimensional) turbulence. In additionto the flows driven by angular-momentum redistribution, he also considered flows driven bythe local baroclinicity of the rotating star, and showed that their effect can be representedas a hyperdiffusion term in the angular velocity evolution equation. In a quasi-steady, uni-formly rotating limit, the flows described are akin to a local Eddington-Sweet circulation.His formalism is used today in stellar evolution models with rotation (Maeder & Meynet,2000).In all cases described above, mixing is either localised near the base of the convectionzone (overshoot, gravity-wave mixing), or significant only in very rapidly rotating stars (localEddington-Sweet circulations) or stars which are undergoing major angular-momentum re-distribution (during phases of gravitational contraction, spin-down, mass loss, etc..). In thispaper, we have identified another potential cause of mixing, where the original energy sourceis the differential rotation in the stellar convective region: gyroscopic pumping (induced bythe Coriolis force associated with the differential rotation, see McIntyre 2007) drives large-scale meridional flows which may – under the right circumstances – penetrate the radiativeregion and cause a global circulation.This source of mixing is intrinsically non-local to the radiative zone. The simplest wayof seeing this is to consider a thought-experiment where the radiative–convective interface isimpermeable: as shown in Paper I, the amplitude of the meridional flows generated locally(i.e. below the interface) is then much smaller than the one calculated here. The originof the flows is also clearly independent of the baroclinicity (since the same phenomenon is 33 –observed in the unstratified limit), although the flows themselves can be influenced by thestratification. This implies that they are not related to Eddington-Sweet flows. Finally,contrary to some of the other mixing sources listed above, the one described here does notrely on the system being out-of-equilibrium: it is an inherently quasi-steady phenomenon,implying that this process is an ideal candidate for “deep mixing” for stars on the MainSequence.The process together with the conditions under which strong mixing might occur, arenow summarized and discussed.
The differential rotation observed in stellar convective envelopes (e.g. Barnes et al. mimics the effect of convective Reynolds stresses indriving the system towards a differentially rotating state. Using this model, we then derivethe expected mixing caused by large-scale meridional flows .We first found that large-scale flows are indeed self-consistently driven by gyroscopicpumping in the convection zone, as expected (McIntyre, 2007). The amplitude of these flows within the convection zone scales roughly as V cz ∼ τ R ⋆ (∆Ω) , (39)where (∆Ω) is the observed equator-to-pole differential rotation, R ⋆ is the stellar radius, and τ is, as discussed in Section 2.2, related to the ratio of the convective turnover time dividedby the rotation period. Note that for the Sun, with (∆Ω) ∼ . ⊙ , the typical amplitudeof the corresponding meridional flows would be of the order of 200 τ m/s – which doesn’tseem too unreasonable given the observations of subsurface flows (Giles et al. τ in the solar convection zone (see Section 2.2).Next, we studied how much mixing these flows might induce in the underlying radiativezone. In this quasi-steady formalism , we found that the magnitude of convection-zone- It is worth noting here that while we expect the details of the flow structure and amplitude to bedifferent when a more realistic forcing mechanism is taken into account, the overall scalings derived shouldnot be affected.
34 –driven flows decays exponentially with depth below the radiative–convective interface on thelengthscale l , where l = R ⊙ / Pr / Ro rz , as determined in Section 2.4.2 (see also Gilman &Miesch, 2004 and Garaud & Brummell, 2008). This penetration corresponds (in the linearregime) to a so-called “thermo-viscous” mode. The limit l ≪ R ⊙ corresponds to a stronglystratified limit, where the flow velocities are rapidly quenched beneath the convection zone.The limit l ≫ R ⊙ corresponds to the weakly stratified case, where the thermo-viscous modespans the whole radiative interior and the stratification has little effect on the flow. It isimportant to note that “weakly stratified” regions in this context can either correspond toregions with weak temperature stratification (small ¯ N ), or in rapid rotation, or with smallPrandtl number – this distinction will be used later.The amplitude of the flows upon entering the radiative zone V rz , together with l ,uniquely define the global circulation timescale in the interior (roughly speaking, l /V rz ).In the weakly stratified/rapidly rotating limit, we find that the fraction of the meridionalmass flux pumped in the convection zone which is allowed to enter the radiative zone isstrongly constrained by Taylor-Proudman’s theorem. This theorem, which holds when thepressure gradient and the Coriolis force are the two dominant forces and are therefore inbalance, enforces the invariance of all components of the flow velocities along the rotationaxis. Hence, flows which enter the radiative zone cannot return to the convection zone unlessthe Taylor-Proudman constraint is broken. However, additional stresses (such as Reynoldsstresses, viscous stresses, magnetic stresses) are needed to break this constraint. As a result,two regimes may exist. If the (weakly stratified/rapidly rotating) radiative zone is in pure Taylor-Proudman balance, then the system adjusts itself, by adjusting the pressure field, insuch a way as to ensure that the convection zone flows remain entirely within the convectionzone. On the other hand, if there are other sources of stresses somewhere in the radiativezone to break the Taylor-Proudman balance, then significant large-scale mixing is possiblesince flows entering the radiative zone are allowed to return to the convection zone. Fur-thermore, the resulting meridional mass flux in the radiative zone depends rather sensitivelyon the nature of the mechanism which breaks the Taylor-Proudman constraint (see nextsection).The strongly stratified/slowly rotating limit exhibits a very different behavior. Becauseof the strong buoyancy force, the Taylor-Proudman balance becomes irrelevant, the flowsare exponentially suppressed, and the induced radiative zone mixing is independent of thelower boundary conditions. However, note that since ¯ N tends to 0 at a radiative–convectiveinterface, there will always be a “weakly stratified” region in the vicinity of any convective more precisely, the perturbation to the pressure gradient around hydrostatic equilibrium
35 –zone. In that region the dynamics described in the previous paragraph apply.
In the illustrative model studied here, the only stresses available to break the Taylor-Proudman constraint are viscous stresses, which are only significant within the thin Ekmanlayer located near an artificial impermeable inner boundary. We do not advocate that thisis a particularly relevant mechanism for the Sun! However, it is a useful example of thesensitive dependence of the global circulation mass flux on the mechanism responsible forbreaking the Taylor-Proudman constraint.In the limit of weak stratification, we found that if the inner boundary is a stress-freeboundary then the global turnover time within the radiative zone is the viscous timescale.This is because stress-free boundary conditions effectively suppress the Ekman layer. On theother hand if the boundary layer is no-slip, then the global mass flux through the radiativezone is equal to the mass flux allowed to return through the Ekman layer . In that case, andaccording to well-known Ekman layer dynamics, the overall turnover time within the bulkof the domain is the geometric mean of the viscous timescale and the rotation timescale(1 / E / ν Ω ⊙ ), which correspond to a few million years only.Going beyond simple Ekman dynamics, a much more plausible related scenario for thesolar interior was studied by Gough & McIntyre (1998). They considered the same mech-anism for the generation of large-scale flows within the convection zone, studied how theseflows down-well into the radiative zone and interact with an embedded large-scale primordialmagnetic field. They showed that the field can prevent the flows from penetrating too deeplyinto the radiative zone, while the flows confine the field within the interior. In their model,this nonlinear interaction occurs in a thin thermo-magnetic diffusion layer, located somewhatbelow the radiative–convective interface. One can therefore see an emerging analogy withthe dynamics discussed here: in the Gough & McIntyre model, the field does act as a some-what impermeable barrier, and provides an efficient and elegant mechanism for breakingthe Taylor-Proudman constraint within the radiative zone. The only significant difference isthat the artificial Ekman layer is replaced by a more convincing thermo-magnetic diffusionlayer: the mass flux allowed to down-well into the radiative zone, and mix its upper regions,is now controlled by a balance between the Coriolis force and magnetic stresses (instead ofthe viscous stresses). With this new balance, they find that the global turnover time for thecirculation in the region between the base of the convection zone and the thermo-magneticdiffusion layer is of the order of a few tens of millions of years (which is still short comparedwith the nuclear evolution timescale). 36 –This mixed region is the solar tachocline. By relating their model with observations,Gough & McIntyre were able to identify the position of the magnetic diffusion layer to bejust at the base of the observed tachocline (around 0 . ± . R ⊙ , see Charbonneau et al. − cm/s) down to the magnetic diffusion layer. However,it is rather interesting to note that the Gough & McIntyre model could not have workedhad today’s tachocline been observed to be much thicker. It is also interesting to note thatfor younger, more rapidly rotating solar-type stars, a much larger region of the radiativezone can be considered “weakly stratified”, possibly leading to much deeper mixed regions ifthese stars also host a large-scale primordial field. The implications of these findings for Liburning, together with a few other interesting ideas, will be discussed in a future publication. Acknowledgments
This work was funded by NSF-AST-0607495, and the spherical domain simulations wereperformed on the UCSC Pleiades cluster purchased with an NSF-MRI grant. The authorsthank N. Brummell and S. Stellmach for many illuminating discussions, and D. Gough forthe original idea.
Appendix A: Ekman jump condition
Equation (15) provides the solution “far” from the lower boundary, in the bulk of thefluid. Let us refer to the limit of bulk solutions as z → u bulk (0 + ) (and similarly forthe other quantities). We now derive the Ekman solution close to the boundary, for theunstratified case. Let’s study the problem using the stream-function ψ with(ˆ u, ˆ v, ˆ w ) = ˆ u, d ˆ ψ d z , − ik ˆ ψ ! . (1)Moreover, let us assume that within the boundary layer, d ˆ ψ/ d z ≫ k ˆ ψ . The governingequations are then approximated by − ψ d z = E ν d ˆ u d z ,2ˆ u = − ik ˆ p + E ν d ˆ ψ d z , 37 –dˆ p d z = − ikE ν d ˆ ψ d z , (2)which simplify to d ˆ ψ d z = − ψ d z , (3)with solutions ˆ ψ ( z ) = ψ + ψ e λ + ψ e − λ z + ψ e λ z + ψ e − λ z ,ˆ u ( z ) = u − E ν (cid:20) λ ψ e λ z − λ ψ e − λ z + 1 λ ψ e λ z − λ ψ e − λ z (cid:21) , (4)where λ = (1 + i )E − / ν , λ = (1 − i )E − / ν . (5)The growing exponentials are ignored to match the solution far from the boundary layer;it then becomes clear that u = ˆ u bulk (0 + ), while − ikψ = ˆ w bulk (0 + ). Requiring no-slip,impermeable conditions at z = 0 implies ψ + ψ + ψ = 0 , λ ψ + λ ψ = 0 , u − E ν (cid:20) − λ ψ − λ ψ (cid:21) = 0 , (6)which in turn implies ψ = − λ λ ψ , ψ = λ λ − λ ψ , u = 2 E ν λ + λ λ λ ψ = 2 E − / ν ψ . (7)This last equation then uniquely relates the limit of the bulk solution ˆ u (0 + ) and ˆ w (0 + ) as z → u bulk (0 + ) = 2 ik E − / ν ˆ w bulk (0 + ) , (8)yielding the standard Ekman jump condition . Appendix B: Stratified stress-free solution
The boundary conditions discussed in Section 2.3.3 imply the following set of equations.At z = 0, ˆ w = 0 and ˆ u z = 0 (alternatively, ˆ T = 0):0 = u − u , 38 –0 = ku − ku + λ u − λ u . (9)At z = 1: ˆ w = 0 and ˆ T = 0: 0 = Ae /δ + Be − /δ − iSk Λ , (10)0 = T e k + T e − k . (11)Finally, matching conditions on ˆ w , ˆ p , ˆ T and d ˆ T / d z at z = h : − ikE ν k − λ λ u sinh( λ h ) = Ae h/δ + Be − h/δ − iSk Λ ,2 u cosh( kh ) + 2 u cosh( λ h ) = U ( h ) + ik δ Λ (cid:2) Ae h/δ − Be − h/δ (cid:3) , T e kh + T e − kh = − ik Ro [ ku sinh( kh ) + λ u sinh( λ h )] , T e kh − T e − kh = − ik Ro (cid:2) k u cosh( kh ) + λ u cosh( λ h ) (cid:3) , (12)where u and u were already eliminated using equations (11). We now proceed to eliminate A , B , T and T , which leaves two equations for u and u :2 G [ u cosh( kh ) + u cosh( λ h )] − δ Λ k E ν k − λ λ u sinh( λ h ) ,= δS (cid:0) e ( h − /δ (1 − G ) − (cid:1) + GU o ( h ) ,( F − ku sinh( kh ) + ku cosh( kh ) = − λ k u cosh( λ h ) − ( F − λ u sinh( λ h ) ,(13)where the functions F ( h, k ) and G ( h, h ) are geometric factors defined as F ( h, k ) = 21 − e k ( h − , G ( h, k ) = e ( h − /δ − e − h/δ e ( h − /δ + e − h/δ . (14)These equations form a linear system for u and u with u = − Hu , u = δS (cid:0) e ( h − /δ (1 − G ) − (cid:1) + GU o ( h )2 G [ − H cosh( kh ) + cosh( λ h )] − δ Λ k E ν k − λ λ sinh( λ h ) , (15)and where the function H ( h, k, λ ) is given as H ( h, k, λ ) = λ k λ k cosh( λ h ) + ( F −
1) sinh( λ h )cosh( kh ) + ( F −
1) sinh( kh ) . (16) 39 –These rather opaque solutions can be clarified a little by looking at the various relevantlimits. For weakly stratified fluids λ →
0. Then H ( h, k, λ ) = O ( λ ) →
0, and so u = δS (cid:0) e ( h − /δ (1 − G ) − (cid:1) + GU o ( h )2 G − δ Λ k E ν h + O ( λ ) . (17)In the limit E ν → u = O ( λ ) , u = 12 (cid:20) U o ( h ) − δS − cosh((1 − h ) /δ )sinh((1 − h ) /δ ) (cid:21) . (18)Folding this back into the original solution in the radiative zone then yieldsˆ w ( z ) = − ik E ν u z = − ik E ν (cid:20) U o ( h ) − δS − cosh((1 − h ) /δ )sinh((1 − h ) /δ ) (cid:21) z , (19)which is identical to equation (19).In the opposite, strongly stratified limit, λ ≫ k . Then we have instead H ( h, k, λ ) ≃ λ k cosh( λ h )cosh( kh ) + ( F −
1) sinh( kh ) , (20)so that this time u = O ( λ − ) →
0, and in the limit E ν → u = 12 cosh( kh ) (cid:20) U o ( h ) − δS − cosh((1 − h ) /δ )sinh((1 − h ) /δ ) (cid:21) . (21)Folding this back into the equation for ˆ w ( z ) in the radiative zone now yieldsˆ w ( z ) = − ik E ν λ sinh( λ z )cosh( λ h ) (cid:18) − tanh( kh )tanh( k (1 − h )) (cid:19) (cid:20) U o ( h ) − δS − cosh((1 − h ) /δ )sinh((1 − h ) /δ ) (cid:21) ,(22)therefore justifying the scaling discussed in 2.4.3. REFERENCES
Barnes, J.R., Collier Cameron, A., Donati, J.-F., James, D. J., Marsden, S. C., & Petit, P.,2005, MNRAS, 357, L1Brummell, N. H., Clune, T. L., & Toomre, J., 2002, ApJ, 570, 825Chaboyer, B. & Zahn, J.-P., 1992, A&A, 253, 173 40 –Charbonneau, P., 2005, LRSP, 2, 2Charbonnel, C. & Talon, S., 2005, Science, 309, 2189Christensen-Dalsgaard, J., et al. ∼ pgaraud/Garaud, P., 2007. in The Solar Tachocline , pp. 147–181, eds. Hughes, D. W., Rosner, R. &Weiss, CUP.Garaud, P. & Brummell, N. H., 2008, ApJ, 674, 498Garaud, P. & Garaud, J.-D., 2008, MNRAS, 391, 1239Giles, P. M., Duvall, T. L., Jr., Scherrer, P. H. & Bogart, R. S, 1997, Nature, 390, 52Gilman, P. A. & Miesch, M. S., 2004, ApJ, 611, 568Gough, D. O. & McIntyre, M. E., 1998, Nature, 394, 755Gough, D. O., 2007. in
The Solar Tachocline , pp. 3–30, eds. Hughes, D. W., Rosner, R. &Weiss, CUP.Kippenhahn, R., 1963, ApJ, 137, 664Kitchatinov L. L. & R¨udiger, G., 1993, A&A, 276, 96Kitchatinov L. L. & R¨udiger, G., 2005, Astr. Nachr., 326, 379LaBonte, B. J. & Howard, R. F., 1982, Solar Phys., 80, 361Maeder, A. & Meynet, G., 2000, ARA&A, 38, 143McIntyre, M. E., 2007. in
The Solar Tachocline , pp. 183-212, eds. Hughes, D. W., Rosner,R. & Weiss, CUP.Pinsonneault, M., 1997, ARA&A, 35, 557Rempel, M., 2005, ApJ, 622, 132 41 –Rogers, T. M., MacGregor, K. B. & Glatzmaier, G. A., 2008, MNRAS, 387, 616R¨udiger, G., 1989, in
Differential rotation and stellar convection. Sun and the solar stars ,Chapter 4. Publisher: Berlin, Akademie VerlagSchatzman, E., 1996, J. Fluid Mech., 322, 355Schou, J., et al. , 1998, ApJ, 505, 390Spiegel, E. A. & Bretherton, F. P., 1968, ApJ, 153, L77Zahn, J.-P., 1992, A&A, 265, 115
This preprint was prepared with the AAS L A TEX macros v5.2.
42 –Fig. 13.— Vertical velocity (in units of R ⊙ Ω ⊙ ) in nine different simulations, at latitude 80 ◦ .The background stratification in each case is solar, but the position of the lower boundaryis moved through the radiative zone from 0 . R ⊙ to 0 . R ⊙ and 0 . R ⊙ . The solid-line plotsare for no-slip lower boundary conditions while the dotted lines are for stress-free lowerboundary conditions. Three simulations are shown in each case: (from lowest to highestcurve) for f = 10 , f = 10 and f = 10 corresponding to E ν = 2 × − to E ν = 2 × − .The logarithmic scale clearly shows that u r scales with E ν in the radiative zone in the stress-free cases for all values of r in while in the no-slip case, r in scales with E ν if r in < ..