1. Introduction
One of the long-lasting problems in stellar astrophysics is the nature of the mechanisms responsible for the angular momentum and chemicals transport within the stably stratified radiative envelope of rotating massive stars. Such transport can result from various physical processes such as internal gravity waves (e.g. Rogers et al. Reference Rogers, Lin, McElwaine and Lau2013; Lee, Neiner & Mathis Reference Lee, Neiner and Mathis2014), turbulence from shear instabilities due to differential rotation (e.g. Zahn Reference Zahn1974, Reference Zahn1992; Prat & Lignières Reference Prat and Lignières2013, Reference Prat and Lignières2014; Prat et al. Reference Prat, Guilet, Viallet and Müller2016; Garaud, Gagnier & Verhoeven Reference Garaud, Gagnier and Verhoeven2017; Gagnier & Garaud Reference Gagnier and Garaud2018; Kulenthirarajah & Garaud Reference Kulenthirarajah and Garaud2018) and centrifugal driving of meridional circulations. The latter is usually associated with the steady baroclinic state of the radiative envelope of rotating massive stars (e.g. Garaud Reference Garaud2002; Rieutord Reference Rieutord2006; Espinosa Lara & Rieutord Reference Espinosa Lara and Rieutord2013; Hypolite & Rieutord Reference Hypolite and Rieutord2014; Rieutord & Beth Reference Rieutord and Beth2014; Hypolite, Mathis & Rieutord Reference Hypolite, Mathis and Rieutord2018). In such a state, isobar and isopycnic lines (or equipotentials) are different and result in a baroclinic torque $(\boldsymbol {\nabla } p \times \boldsymbol {\nabla } \rho )/\rho ^2$ where p is the pressure and ρ is the density. Because of the strength of the pressure and density gradients in stellar interiors, a slight misalignment of the two vectors is sufficient to drive a sizeable baroclinic flow. In fact, the misalignment remains less than a degree, even for the most rapidly rotating stars (Espinosa Lara & Rieutord Reference Espinosa Lara and Rieutord2011). The effects of baroclinicity can, however, be supplemented by the meridional circulation induced by spin-up or spin-down flows resulting from stellar contraction or expansion (Hypolite & Rieutord Reference Hypolite and Rieutord2014). Furthermore, in massive stars, rapid rotation usually combines with strong radiation-driven winds that may lead to significant latitude-dependent outward mass flow and angular momentum flux (e.g. Maeder & Meynet Reference Maeder and Meynet2000; Pelupessy, Lamers & Vink Reference Pelupessy, Lamers and Vink2000; Georgy, Meynet & Maeder Reference Georgy, Meynet and Maeder2011; Gagnier et al. Reference Gagnier, Rieutord, Charbonnel, Putigny and Espinosa Lara2019a,Reference Gagnier, Rieutord, Charbonnel, Putigny and Espinosa Larab). When this is the case, typically for stars more massive than seven solar masses, the baroclinic flow is supplemented by a meridional circulation resulting from the surface stress induced by the mass loss (Zahn Reference Zahn1992). The transport of products of nuclear reactions in the stellar core can then enrich the surface layers in chemical elements, locally increasing opacity and thus enhancing the radiation-driven outward mass and angular momentum fluxes (Kudritzki, Pauldrach & Puls Reference Kudritzki, Pauldrach and Puls1987; Puls, Springmann & Lennon Reference Puls, Springmann and Lennon2000; Puls, Vink & Najarro Reference Puls, Vink and Najarro2008). The physical process of wind-driven spin-down flow is therefore of crucial importance for the understanding of the secular evolution of rotating massive stars.
So far, however, there is still no consensus on the appropriate way to model the mixing induced by such fluid flows in stellar evolution codes. Indeed, while the transport of chemicals is always accounted for as a diffusive process (as justified by Chaboyer & Zahn Reference Chaboyer and Zahn1992), the angular momentum transport is either treated as an advection–diffusion process following Zahn (Reference Zahn1992), Meynet & Maeder (Reference Meynet and Maeder1997) and Maeder & Zahn (Reference Maeder and Zahn1998) or as a purely diffusive process (Paxton et al. Reference Paxton, Bildsten, Dotter, Herwig, Lesaffre and Timmes2011). Moreover, when facing fast rotation, developments beyond the current one-dimensional (1-D) model approximations become necessary. In this context, the achievement of the first self-consistent 2-D models of rapidly rotating early-type stars, worked out by Espinosa Lara and Rieutord (e.g. Espinosa Lara & Rieutord Reference Espinosa Lara and Rieutord2013; Rieutord, Espinosa Lara & Putigny Reference Rieutord, Espinosa Lara and Putigny2016) and in which the differential rotation as well as the meridional circulation arising from baroclinicity are computed self-consistently, opens the door to the exploration of the evolution of fast stellar rotators. Because the sources of rotation-induced mixing are multiple, a meticulous study of each transport mechanism appears to be necessary for a better understanding of stellar evolution. In this paper we address this issue and investigate the properties of the primary and secondary flows driven by radiation-driven outward mass and angular momentum fluxes at the surface of rotating massive stars.
However, the complete modelling of astrophysical rotation-induced mixing is a complex problem. Hence, to understand its different facets, it is useful to study simplified set-ups which incorporate, step by step, the various physical phenomena that contribute to the whole realistic model. For instance, we shall be particularly interested in the scaling laws which control the viscous effects.
When a star loses mass, the associated wind extracts angular momentum. At some place inside the star, but close to the surface, this extraction generates a radial differential rotation that further extracts angular momentum from the deeper layers. The upper layers therefore impose a torque to the interior of the star, which slowly spins down. As is well known (see Greenspan Reference Greenspan1968), the spin-down of an incompressible fluid inside a rigid container with no-slip boundaries occurs on a time scale $P_{rot}/\sqrt {E}$, where
$P_{rot}$ is the rotation period of the fluid and
$E$ is the Ekman number (see below for its definition). Since in most situations
$E\ll 1$, the spin-down time scale is much shorter that the viscous diffusion time scale
$P_{rot}/E$. In stars, however, boundary conditions are not rigid. Rather, the wind imposes a global angular momentum loss which amounts to a torque applied to the inner layers of the star. Hence, at some depth, the fluid is spun-down by a (turbulent) viscous stress. As in the no-slip case, secondary meridional flows arise as shown by Friedlander (Reference Friedlander1976).
The set-up is further complicated by the convective core of the massive stars at hand. Compared to the envelope, the core may be viewed as a very viscous region, since thermal convection is highly turbulent there. To simplify, we shall represent it with the limiting case of a solid ball. A solid boundary or an important jump in viscosity both lead to the formation of Stewartson shear layers staying along the tangent cylinder of the core, parallel to the rotation axis (e.g. Stewartson Reference Stewartson1966; Rieutord Reference Rieutord2006). Hence, the secondary flows are certainly more complex that those arising in a full sphere. Moreover, we wish to know how long it takes for chemical elements generated in the core by nucleosynthesis to reach the stellar surface, where they can be observed. Hence, it is important to know the meridional circulation turnover time and its dependence with viscosity. In addition thermal stratification and density variations between the core and the surface also influence the spin-down flow.
Our simplified model takes into account all these features of the problem so as to examine their particular role. Let us summarise this model: we consider a spherical shell containing a viscous fluid where the rotation of the core is imposed and where an imposed large-scale stress is applied on the outer surface. We first consider a fluid of constant density and temperature, then discuss the case with a stable thermal stratification (as requested by radiative envelopes) and, finally, we allow for radial density variations using a polytropic envelope within the anelastic approximation (e.g. Jones, Kuzanyan & Mitchell Reference Jones, Kuzanyan and Mitchell2009). We further assume that the angular momentum extraction and the associated surface stress are weak enough that linearised equations can be used for the axisymmetric flow determinations.
The paper is organised as follows: in § 2 we introduce our model and describe the numerical method used for solving the equations of motion. We then consider the case of an incompressible flow and we discuss the time scales governing the transient phase and the asymptotic properties of the stationary primary (differential rotation) and secondary (meridional circulation) flows in the limit of small Ekman numbers in § 3. In § 4, we discuss the role of thermal stratification using the Boussinesq approximation, and we finally include density variations of the background using the anelastic approximation in § 5. Discussion and conclusions follow in § 6.
2. Formulation of the problem
2.1. Description
A fluid of constant kinematic viscosity $\nu$ is enclosed between two spheres (see below). The inner one is rigidly rotating with a constant angular velocity
$\varOmega _c$ while the outer shell supports a prescribed tangential surface stress
$\tau _{*}(\theta )$;
$R$ and
$\eta R$ are the radii of the outer and inner shells, respectively, and
$\theta$ is the colatitude. We sketch out this model in figure 1. We consider the driving stress to be sufficiently weak that the rotation period of the system is much shorter than the typical turnover times associated with the flow in the rotating frame of reference. We thus consider the nonlinear terms to be negligible, which we justify a posteriori in appendix A. The dimensionless equation of vorticity and mass conservation equation in the rotating frame with angular velocity
$\varOmega _c$, and in the context of an ‘anelastic’ flow, read
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn1.png?pub-status=live)
where we have used $R$ as the length scale,
$(2\varOmega _c)^{-1}$ as the time scale and
$\rho _c$, the density at the inner shell, as the density scale. Here,
$\boldsymbol{e}_z$ is the unit vector in the z-direction,
$\boldsymbol{u}$ is the non-dimensionalised velocity vector,
$E$ is the Ekman number defined as
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn2.png?pub-status=live)
and
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn3.png?pub-status=live)
is the dimensionless viscous force. This expression of the viscous force is meant to represent the case where the envelope is pervaded by some small-scale turbulence with constant diffusive properties as often used in stellar physics (Brandenburg et al. Reference Brandenburg, Jennings, Nordlund, Rieutord, Stein and Tuominen1996; Käpylä, Mantere & Brandenburg Reference Käpylä, Mantere and Brandenburg2012). Centrifugal effects are neglected altogether, and the flows are considered axisymmetric.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig1.png?pub-status=live)
Figure 1. Schematic view of the system: the inner shell of radius $\eta R$ rotates with an angular velocity
$\varOmega _c$ while the outer one of radius
$R$ rotates differentially at
$\varOmega _s(\theta )$ by means of a prescribed tangential surface stress
$\tau _{*}(\theta )$. The dashed lines correspond to the edges of the tangent cylinder
$\mathcal {C}$, circumscribing the inner sphere.
This system is then completed by boundary conditions. On the outer shell, we impose a specified stress, which represents the torque resulting from the stellar wind. The simplest dimensionless expression of such a stress, which we take equatorially symmetric, is
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn4.png?pub-status=live)
where $[\sigma ]$ is the dimensionless stress tensor and
$A$ is a positive constant that sets the amplitude of the prescribed (braking) stress. Although the wind imposes a small radial flow at the surface, we shall ignore it in this first approach, and impose a vanishing normal velocity on the outer shell. Finally, we impose no-slip boundary conditions at the inner boundary, as justified by the important viscosity jump at the stellar envelope/core interface. Boundary conditions thus read
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn5.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn6.png?pub-status=live)
2.2. Numerical method
We use a spectral decomposition and expand the fields in spherical harmonics for the angular part and Chebyshev polynomials for the radial part (see Rieutord & Valdettaro Reference Rieutord and Valdettaro1997). We write
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn7.png?pub-status=live)
where
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn8.png?pub-status=live)
with $Y^{m}_l$ being the usual normalised scalar spherical harmonic function. Because the
${\boldsymbol {q}}$-flow is divergenceless both in the incompressible and anelastic models,
$v_{m}^{l}$ can be expressed as a function of
$u_m^l$ only (Rieutord Reference Rieutord1987). Projecting the vorticity equation (2.1) on
$\boldsymbol {R}^{m}_l$ and
$\boldsymbol {T}^{m}_l$ for an axisymmetric flow (
$m=0$), we obtain the following system of equations for the radial parts
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn9.png?pub-status=live)
where $\varLambda =l(l+1)$. Here,
$A_{l-1}^{l}$,
$A_{l+1}^{l}$,
$B_{l-1}^{l}$ and
$B_{l+1}^{l}$ are the coupling coefficients defined as
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn10.png?pub-status=live)
The viscous force is projected on the spherical harmonics basis as well, namely
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn11.png?pub-status=live)
where
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn12.png?pub-status=live)
with
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn13.png?pub-status=live)
and where $b(r)=1/\rho (r)$ is the inverse density function.
Finally, the boundary conditions in the rotating frame of reference are projected on the spherical harmonics and read
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn14.png?pub-status=live)
where $\delta _{ij}$ is the Kronecker symbol.
3. The incompressible flow
In this section, and as a first step, we solve (2.1) for constant density throughout the domain, namely for $\rho (r)= 1$, with the boundary conditions (2.5), and with the initial condition
${\boldsymbol {u}}= \boldsymbol {0}$. Of course, before the steady state is reached, the flow is in a transient state during which the driving surface stress has not yet been fully communicated to the interior of the fluid. We now briefly analyse this transitional stage so as to estimate the transient time scale.
3.1. The transient phase
As to assess the relevance of the steady-state approximation for geophysical and astrophysical applications, it is important to determine the time scales the flow characteristics are governed by. To do so, we measure the time at which the evolution of the total angular momentum in the rotating frame may be considered to its end. We note that the change in angular momentum of the fluid is due to the difference between the torque exerted on the fluid by the outer boundary condition and the torque that the fluid exerts on the steady inner sphere. Therefore, a stationary state is reached when
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn15.png?pub-status=live)
where $\varGamma (1)$ and
$\varGamma (\eta )$ are the torques about the
$z$-axis exerted on the outer and inner boundary surfaces, respectively. The torque about the
$z$-axis exerted on a layer located at some radius
$r$ can be written
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn16.png?pub-status=live)
where $\boldsymbol {\sigma }$ is the local stress applied on the sphere of radius
$r$. We expand
$\boldsymbol {\sigma }$ on the spherical harmonics basis, hence the torque about the
$z$-axis reads
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn17.png?pub-status=live)
where
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn18.png?pub-status=live)
Using Legendre polynomials recurrence relations, (3.3) actually reads
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn19.png?pub-status=live)
Applying boundary conditions (2.13), (3.1) can finally be written
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn20.png?pub-status=live)
where the non-dimensional surface density $\rho _s=\rho (1)=1$ for the considered incompressible model. We monitor the evolution of the relative difference between the torques
$\Delta \varGamma /\varGamma (1)$ for various Ekman numbers and show the result in figure 2. It shows that on a time scale
$O(E^{-1})$ the stress is completely communicated to the interior flow, and a steady state is reached.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig2.png?pub-status=live)
Figure 2. Relative difference between the torque applied on the outer sphere and the torque exerted by the fluid on the stationary inner sphere, as a function of the reduced time $E t$. We find the steady state to be reached on a
$O(E^{-1})$ time scale, that is on a viscous time scale, in the asymptotic regime of small Ekman numbers.
Of course in the astrophysical or geophysical context, it may become necessary to relax the constant rotation of the inner sphere. The torque exerted by the fluid on this boundary can make the inner sphere spin-down. According to the angular momentum theorem
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn21.png?pub-status=live)
where $\varGamma _*$ is the dimensional torque,
$I_c$ is the moment of inertia of the core assumed to be in solid-body rotation and
$L^{core}_z$ is its total angular momentum. Starred quantities are dimensional.
The core spin-down then engenders an additional Euler force $\dot {\varOmega _c}{\boldsymbol {e}}_z \times {\boldsymbol {r}}$ since our frame is attached to the core. This force can, however, be neglected if the time scale associated with the core spin-down is larger than the (viscous) time scale during which the angular momentum is redistributed. In that case, the flow can reach a quasi-steady state. In the present work, for simplicity, we shall use this assumption (
$\dot {\varOmega _c}=0$). Astrophysically, we justify its adoption by the Roche approximation, often used for massive stars, where the whole mass of the star is assumed to be in the core.
Another important time scale is that of the rise of the Stewartson layer. This layer appears as the result of the equatorial singularity of the inner Ekman boundary layer. It is well known that within an equatorial band of latitude $O(E^{1/5})$, the thickness of the Ekman boundary layer is
$\delta _E = O(E^{2/5})$ (Roberts & Stewartson Reference Roberts and Stewartson1963). Hence, this singular equatorial viscous boundary layer is expected to be fully developed, and to initiate the development of the Stewartson layer, on the
$O(\delta _E^2 /E)=O(E^{-1/5})$ time scale, that is much shorter than the
$O(E^{-1})$ viscous time scale. Likewise, the central Stewartson layer of thickness
$O(E^{1/3})$ is fully developed on the
$O(E^{-1/3})$ time scale, again, much shorter than the viscous time scale. We therefore expect the Stewartson layer to start developing within a few revolutions of the inner shell and to be fully developed on a time scale that is much shorter than the viscous one on which we expect the shear flow inside it to evolve towards the steady state. We verify this in figure 3 showing the amplitude of the streamfunction
$\psi _h/\psi _{h,st}$, taken at cylindrical radius
$s=\eta$ and
$z=1/2$, and normalised by its value at the steady state, as a function of the reduced time
$Et$ for various Ekman numbers.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig3.png?pub-status=live)
Figure 3. Amplitude of the streamfunction taken at $s=\eta$ and
$z=1/2$ normalised by its value at the steady state
$\psi _h/\psi _{h,st}$, as a function of the reduced time
$Et$ for various Ekman numbers. We find the Stewartson to appear within a few inner sphere revolutions and to evolve in a viscous time scale.
The analysis of the transient phase preceding the settling of a steady state thus underlines that the entire meridional flow evolves on a viscous time scale of order $E^{-1}$. In the astrophysical context characterised by extremely small Ekman numbers (typically
$E<10^{-10}$), such a scaling implies that a steady state of the radiative envelope of massive stars is not reached during their lifetime.
3.2. The steady flow
Even if it is not reached during the lifetime of the star, the steady state is worth studying since it owns many simple features and its structure is very similar to that of the transient flow.
It is well known that the linear and stationary vorticity equation of an inviscid incompressible flow verifies Taylor–Proudman theorem (Proudman Reference Proudman1916; Taylor Reference Taylor1921), namely
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn22.png?pub-status=live)
For quasi-inviscid interior flows, that is for sufficiently small Ekman numbers, we thus expect a quasi-geostrophic solution for the velocity field. Figure 4 shows the normalised angular velocity as a function of the cylindrical radial coordinate $s$ for
$E=10^{-7}$,
$\eta =0.35$ and for various radial distances
$r$. Figure 5 shows a 2-D view of the meridional circulation as well as the differential rotation in the rotating frame for the same model. The resulting flow is characterised by a nested Stewartson shear layer along the tangent cylinder
$\mathcal {C}$ where the meridional circulation is essentially concentrated. This narrow shearing region separates two regions of quasi-geostrophic flow: the volume inside and outside
$\mathcal {C}$. Inside
$\mathcal {C}$, the no-slip conditions on the inner core impose an almost rigid rotation, while outside
$\mathcal {C}$ a columnar differential rotation appears as a consequence of the surface stress.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig4.png?pub-status=live)
Figure 4. Normalised angular velocity as a function of the cylindrical radial coordinate $s$ for
$E=10^{-7}$,
$\eta =0.35$,
$A=0.01$ and for various
$r$.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig5.png?pub-status=live)
Figure 5. Meridional view of the streamfunction $\psi$
$(a)$ and of the differential rotation in the rotating frame of reference
$\delta \varOmega = (\varOmega -\varOmega _c)/2\varOmega _c$
$(b)$ for
$E=10^{-7}$,
$\eta =0.35$ and
$A=0.01$. Three regions can be distinguished: a Stewartson layer on the outer boundary of a tangent cylinder
$\mathcal {C}$ of radius
$\eta$ separating two regions of weak amplitude circulation. The region inside the cylinder
$\mathcal {C}$ corotates with the inner sphere and the region outside
$\mathcal {C}$ rotates differentially with a columnar profile.
3.2.1. Flow outside the core tangent cylinder
$\mathcal {C}$
We now seek for an analytical solution for both the primary and secondary flows outside $\mathcal {C}$ (
$s > \eta$). In the limit of small Ekman numbers, the divergenceless velocity satisfies the Taylor–Proudman theorem, we thus seek a geostrophic solution for both primary and secondary flows outside the tangent cylinder
$\mathcal {C}$. Such a solution usually does not satisfy the viscous boundary conditions, however. We thus decompose the dynamical variables as
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn23.png?pub-status=live)
where overlined variable correspond to the interior geostrophic fields, and tilded variables are their boundary layer corrections. Introducing the $O(1)$ stretched boundary layer coordinate
$\zeta = (1-r)/\sqrt {E}$, and using the geostrophic equilibrium relation
${\boldsymbol {e}}_z \times \bar {{\boldsymbol {u}}} = - \boldsymbol {\nabla } \bar {p}$, the equation of motions, or Ekman boundary layer equation, reads
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn24.png?pub-status=live)
where $\boldsymbol {n} \equiv {\boldsymbol {e}}_r$ is the outwardly directed normal vector at the outer shell. Integrating the azimuthal component of (3.10) from
$0$ to
$\infty$ yields, at the lowest order,
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn25.png?pub-status=live)
Here, $\tilde {Q}$ is the outer Ekman layer
$\theta$-directed volume flux.
Using the boundary condition for the prescribed tangential stress on the outer shell (2.5), (3.11) can be rewritten, at the lowest order
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn26.png?pub-status=live)
where $F(s)= \bar {u}_{\phi }$ is the geostrophic solution for the azimuthal velocity in the interior. To ensure mass conservation, the divergence of
$\tilde {{\boldsymbol {u}}}$ in the boundary layer generates further interior motion by establishing a weak amplitude secondary normal flow
$\tilde {u}_r$ known as Ekman pumping. This radial velocity can be determined by integrating the mass conservation equation in the outer Ekman boundary layer
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn27.png?pub-status=live)
This equation shows that, near the outer shell, and in the Ekman layer, for the boundary condition on the tangential stress to be ensured, there must exist a $O(\sqrt {E})$ tangential flow that, in turn, generates a
$O(E)$ circulation in the interior. We now express the geostrophic azimuthal velocity
$F(s)$ as a function of the prescribed tangential surface stress
$\tau (\theta )$ driving differential rotation. To do so, we ensure that the no-penetration boundary condition
$u_r=0$ at
$r=1$ is enforced, namely
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn28.png?pub-status=live)
where $u_s$ and
$u_z$ are obtained from the azimuthal component of the momentum equation and from the mass conservation equation, respectively. Equation (3.14) finally yields a third-order ordinary differential equation for
$F(s)$
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn29.png?pub-status=live)
where $K=A \sqrt {3/4{\rm \pi} }$. The solutions of (3.15) can be expressed as an integral functional of the prescribed stress (Friedlander Reference Friedlander1976), that is
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn30.png?pub-status=live)
The solution for the geostrophic azimuthal velocity outside the tangent cylinder $\mathcal {C}$, avoiding singularities of the vorticity at
$s = 1$, and accounting for the no-slip boundary condition at the inner shell finally reads
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn31.png?pub-status=live)
with
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn32.png?pub-status=live)
and for which
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn33.png?pub-status=live)
We compare the angular velocity profiles from full numerical solutions to the analytical expression $\delta \varOmega = F(s)/s$, for various Ekman numbers and two inner core radii (
$\eta =0.1$ and
$0.35$) in figure 6. The analytical asymptotic solution reproduces rather well our numerical results when the Ekman number is below
$10^{-9}$. Our boundary layer analysis also predicts the secondary meridional flow amplitude to be
$O(E)$ outside the tangent cylinder
$\mathcal {C}$. In figure 7 we show the streamfunction, defined as
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn34.png?pub-status=live)
as a function of the cylindrical radial coordinate $s$ for various Ekman numbers. We note that
$\psi$ indeed scales as
$E$ outside the Stewartson layer, which appears as oscillations of
$\psi$ near
$s=\eta$. Furthermore, using (3.19a,b), the streamfunction associated with the geostrophic solution for the velocity field reads
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn35.png?pub-status=live)
implying that streamlines are $z=\textrm {Cst}$ lines in a meridional plane. Hence, the outer Ekman boundary layer expels fluid in the
$s$-direction, towards the Stewartson shear layer. The streamfunction
$\psi$ as well as the corresponding streamlines, from (3.21) and from full numerical solutions are represented in figure 8, where we have masked the Stewartson layer. We see that the analytical and numerical solutions nicely match. Interestingly, the flow (3.19a,b) reconnects with the Stewartson layer and partly sources the upwelling flow in this region (see below § 3.2.3).
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig6.png?pub-status=live)
Figure 6. Angular velocity in the rotating frame of reference as a function of the cylindrical radial coordinate $s$ for various Ekman numbers,
$r=0.7$,
$\eta =0.1$
$(a)$ and
$\eta =0.35$
$(b)$, and
$A=0.01$. The black dashed line corresponds to the analytical solution
$\delta \varOmega = F(s)/s$ (see (3.17)).
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig7.png?pub-status=live)
Figure 7. Value of $\psi E^{-1}$ along a meridian at
$r=0.7$ as a function of the cylindrical radial coordinate
$s$ for various Ekman numbers,
$\eta =0.35$ and
$A=0.01$. The predicted
$O(E)$ scaling of the secondary flow amplitude is verified outside the Stewartson layer, which is located at
$s=0.35$.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig8.png?pub-status=live)
Figure 8. Meridional view of the streamfunction and streamlines from full numerical simulations $(a)$, and from the boundary layer analysis
$(b)$, for
$E=10^{-8}$,
$\eta =0.1$ and
$A=0.01$. The tangent cylinder
$\mathcal {C}$ as well as the adjacent reconnecting layer are masked.
The foregoing solution shows that the stress-driven spin-down flow is rather different from the spherical Taylor–Couette flow outside the tangent cylinder. For such a flow, the amplitude of the meridional circulation in this region is vanishing exponentially away from the Stewartson layer (Proudman Reference Proudman1956; Dormy, Cardin & Jault Reference Dormy, Cardin and Jault1998). In our case there is a residual flow, exactly perpendicular to the rotation axis, directed to it, and which scales as the Ekman number.
3.2.2. Flow inside the tangent cylinder
$\mathcal {C}$
Let us now investigate the geostrophic flow inside the tangent cylinder $\mathcal {C}$. Writing
$u_{\phi } = \ell /(r \sin \theta )$ where
$\ell$ is the specific angular momentum, the azimuthal component of the momentum equation and vorticity equation respectively read (e.g. Goldstein Reference Goldstein1938, Reference Goldstein1965; Proudman Reference Proudman1956; Stewartson Reference Stewartson1966)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn36.png?pub-status=live)
where
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn37.png?pub-status=live)
The boundary conditions (2.5) in terms of streamfunction and specific angular momentum in turn read
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn38.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn39.png?pub-status=live)
In the Ekman boundary layers, (3.22) may be rewritten
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn40.png?pub-status=live)
We integrate the first equation over $r$ to obtain the fourth-order differential equation on
$\psi$
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn41.png?pub-status=live)
where $\bar {\psi }$ is a function of
$\theta$ only that can be determined by the boundary conditions (Proudman Reference Proudman1956). Let us first focus on the outer Ekman layer localised near
$r=1$. The solution of (3.26) satisfying the boundary conditions (3.24) as well as the condition
$\psi (\xi \to \infty ) = \bar {\psi }$, where
$\xi =(1-r)\sqrt {\cos \theta / 2E}$, reads
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn42.png?pub-status=live)
From the first equation of (3.25), we may further write
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn43.png?pub-status=live)
where $\alpha =\sqrt {\cos \theta /(2E)} \gg 1$ in the asymptotic limit of small Ekman numbers. The boundary condition (3.24b) then implies
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn44.png?pub-status=live)
where $\bar {\ell }$ is also a function of
$\theta$ only that can be determined by the boundary conditions. Hence, the azimuthal and meridional components of the quasi-geostrophic velocity outside the Ekman and Stewartson shear layers, are not independent and satisfy the outer Ekman jump condition at the lowest order
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn45.png?pub-status=live)
In a similar way, the study of the inner Ekman layer with boundary conditions (3.24a) yields
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn46.png?pub-status=live)
where $\xi '=(r-\eta )\alpha$, and the Ekman jump condition at the inner shell reads
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn47.png?pub-status=live)
Hence (3.32) implies $\bar {\psi }= O(\bar {\ell }\sqrt {E})$ and thus (3.30) implies
$\bar {\psi }= O(E)$ and
$\bar {\ell }= O(\sqrt {E})$. We may therefore rewrite (3.30) as
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn48.png?pub-status=live)
Finally, the geostrophic solution for the velocity inside the tangent cylinder $\mathcal {C}$ (
$s < \eta$) reads, at the lowest order,
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn49.png?pub-status=live)
and the meridional components of the geostrophic velocity at the lowest order
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn50.png?pub-status=live)
We compare these analytical $z$-directed velocity and angular velocity profiles with full numerical solutions in figure 9, for various
$E$. We find our analytical expression to be in good agreement with the numerical solutions. Remarkably, the meridional flow inside the tangent cylinder is parallel to the rotation axis and directed towards the inner core, which is quite different from the outer-
$\mathcal {C}$ meridional flow. Once inside the inner Ekman layer, the flow heads towards the equator where the Ekman layer thickens and eventually changes scale at the equatorial singularity (e.g. Roberts & Stewartson Reference Roberts and Stewartson1963; Stewartson Reference Stewartson1966; Hollerbach Reference Hollerbach1994; Dormy et al. Reference Dormy, Cardin and Jault1998; Marcotte, Dormy & Soward Reference Marcotte, Dormy and Soward2016). The fluid then returns to the outer Ekman layer following the Stewartson layer. Figure 10 shows the meridional view of two arbitrarily selected streamlines, one inside and one outside the tangent cylinder
$\mathcal {C}$. It illustrates the different shape of the secondary flow in the two regions, as well as the reconnecting shear layer redirecting the
$s$-direction flow outside
$\mathcal {C}$ to the Stewartson layer where it flows parallel to the rotation axis towards the outer boundary.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig9.png?pub-status=live)
Figure 9. The $z$-directed velocity
$u_z E^{-1}$
$(a)$ and angular velocity in the rotating frame of reference
$(b)$ inside the tangent cylinder
$\mathcal {C}$ as a function of the cylindrical radial coordinate
$s$, for
$E=10^{-7}$ and
$E=10^{-8}$,
$r=0.7$ and
$\eta =0.35$. The black dashed line corresponds to the analytical geostrophic solutions (3.34a,b) and (3.35a,b).
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig10.png?pub-status=live)
Figure 10. Meridional view of two arbitrarily selected streamlines, inside and outside the tangent cylinder $\mathcal {C}$. The arrows indicate the direction of propagation of the fluid along the streamlines.
3.2.3. Flow in the nested Stewartson shear layer
In figures 7 and 9 we see that the meridional flow strongly deviates from the analytical solution as one nears the Stewartson layer. From the numerical solution, the amplitude of the meridional circulation turns out to be much stronger in this layer than outside it. This is obvious if we compute the total meridional kinetic energy of the secondary flow
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn51.png?pub-status=live)
where $V_m^2= u_r^2 + u_{\theta }^2$. In figure 11(a) we show
$E_{k, tot}$ as a function of the Ekman number. We find
$E_{k, tot}$ to scale as
$E$ in the asymptotic regime
$E\to 0$, whereas the contribution of the flows outside the Stewartson layer remains
$O(E^2)$. Hence, the Stewartson layer deserves some investigation.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig11.png?pub-status=live)
Figure 11. $(a)$ Total meridional kinetic energy of the stationary flow as a function of the Ekman number, for
$A=0.01$ and
$\eta =0.1$. The black dashed line corresponds to the fit in the asymptotic regime (
$E \lesssim 10^{-9}$) and yields
$E_{k, tot} \propto E^{1.074}$.
$(b)$
$E_{k, tot} E^{-1.074}$ as a function of the reduced time
$Et$ for various Ekman numbers,
$A=0.01$ and
$\eta =0.35$.
We first recall that in the Taylor–Couette flow the Stewartson layer is a nested shear layer that can be split into three layers of width $E^{2/7}$,
$E^{1/3}$ and
$E^{1/4}$ (Stewartson Reference Stewartson1966). The
$E^{1/3}$-layer is the central layer, while the
$E^{2/7}$-layer is on the inner side and the
$E^{1/4}$ on the outer side. We now wish to find out if these three nested layers are present in our models. Let us start by noticing that, according to (3.13), the Ekman pumping may be written
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn52.png?pub-status=live)
Hence, introducing the adimensional Cartesian stretched coordinate of order $O(1)$ in a shear layer parallel to the rotation axis, outside the tangent cylinder, and of thickness
$O(E^{\gamma })$,
$\xi =E^{-\gamma }(s-\eta )$, we get
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn53.png?pub-status=live)
where $F=\bar {u}_{\phi }$. We conclude that, in such layers, the radial velocity resulting from the Ekman pumping/suction is of order
$O(E^{1-2\gamma })$. The mass conservation equation, in turn, reads
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn54.png?pub-status=live)
implying that in such a shear layer $u_z = O(E^{-\gamma }u_s)$, that is
$u_z \gg u_s$ in the asymptotic regime of small Ekman numbers.
Simple manipulations of the momentum and vorticity equations allow us to eliminate the velocity $\boldsymbol {u}$ and lead to the following sixth-order partial differential equation for the pressure (Greenspan Reference Greenspan1968)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn55.png?pub-status=live)
representing a balance between Coriolis and viscous forces. The equality of these two terms thus requires $\gamma =1/3$, and indicates a shear layer of thickness
$O(E^{1/3})$ as expected. Hence, in this central ageostrophic layer, we expect
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn56.png?pub-status=live)
Let us now recall that the quasi-geostrophic $E^{1/4}$ layer is associated with the equilibrium between the Ekman pumping flow, and the internal friction of the
$O(1)$ azimuthal flow (e.g. Stewartson Reference Stewartson1966; Barcilon Reference Barcilon1968). This equilibrium is, in fact, exactly that which is in place in the entire region outside the tangent cylinder
$\mathcal {C}$, and the
$E ^{1/4}$-layer is therefore not needed. Indeed, the azimuthal component of the momentum equation in the boundary layer of thickness
$O(E^{\gamma })$ reads
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn57.png?pub-status=live)
that is $u_s=O(E^{1-2\gamma })$ if we assume
$u_{\phi }=O(1)$. In addition, (3.39) further implies
$u_z=O(E^{1-3\gamma })$. However, the Ekman pumping still demands
$u_z=O(E^{1-2\gamma })$ (see (3.37)). Hence,
$\gamma =0$ and the equilibrium between the Ekman pumping flow and the internal friction of the
$O(1)$ azimuthal flow does indeed occur in a
$O(1)$ ‘layer’ that is the entire
$s>\eta$ region. Note that if we had considered no-slip boundary conditions at the outer shell, that is the case of a spherical Taylor–Couette flow, the
$z$-directed velocity resulting from the Ekman pumping would be of order
$O(\sqrt {E} \boldsymbol {\nabla } \times \bar {u}_{\phi } {\boldsymbol {e}}_{\phi } )= O(E^{1/2 - \gamma })$. The matching of the
$u_z$-flux and the Ekman pumping implies that
$1/2 - \gamma =1-3\gamma$, and thus
$\gamma =1/4$ as expected.
Let us finally consider the case of a quasi-geostrophic $z$-directed free shear layer of thickness
$O(E^{\gamma '})$ located inside the tangent cylinder
$\mathcal {C}$ and near
$s= \eta$. It follows that in such a region,
$\ell$ is a function of
$s$ only. Note that in the case of a spherical Taylor–Couette flow,
$\gamma '=2/7$ and
$u_{\phi }(s=\eta ) = O(E^{1/28})$ (e.g. Stewartson Reference Stewartson1966; Marcotte et al. Reference Marcotte, Dormy and Soward2016). In this region, the azimuthal component of the momentum equation (3.22) integrated over
$z$ reads
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn58.png?pub-status=live)
where $f(s)$ is determined by applying the inner Ekman layer jump condition (3.32) at
$z=0$ and
$s - \eta \ll 1$ (Stewartson Reference Stewartson1966), yielding
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn59.png?pub-status=live)
In addition, at $z=\sqrt {1-\eta ^2}$, that is at
$r=1$ and
$s=\eta$, (3.30) reads
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn60.png?pub-status=live)
and (3.43) finally becomes
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn61.png?pub-status=live)
where we have introduced the stretched $O(1)$ shear layer coordinate
$x=(1-s/\eta )/E^{\gamma '}$. At this point, it is interesting to determine the correct balancing in (3.46). We note that the only balance implying the existence of an asymptotically narrow shear layer is that of the first two terms, the third one being of higher order. This balance enforces
$\gamma ' = 2/7$, that is a standard inner Stewartson shear layer of thickness
$O(E^{2/7})$. Hence, the second-order ordinary differential equation for the specific angular momentum at the lowest order and in the
$E^{2/7}$ shear layer reads
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn62.png?pub-status=live)
where
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn63.png?pub-status=live)
This equation can be transformed into a Bessel equation, and the solution asymptotically matching the solution (3.34a,b), that is decaying to zero as $x \to + \infty$, reads (see also Stewartson Reference Stewartson1966; Moore & Saffman Reference Moore and Saffman1968; Marcotte et al. Reference Marcotte, Dormy and Soward2016)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn64.png?pub-status=live)
where $K_{\nu }(z)$ is the modified Bessel function of second kind and
$C$ is a constant to be determined. We note that
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn65.png?pub-status=live)
and thus, at the lowest order
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn66.png?pub-status=live)
Finally, $\textrm {d} \ell /\textrm {d}s$ is rendered continuous across
$s= \eta$ in the two quasi-geostrophic regions, that is
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn67.png?pub-status=live)
which yields the specific angular momentum at $s=\eta$
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn68.png?pub-status=live)
We verify the shear layer thickness as well as the analytical azimuthal velocity profile (3.49) comparing it with full numerical solutions in figure 12(a). The lowest-order solution for the azimuthal velocity in the entire $s \leq \eta$ domain is shown in figure 12(b). Finally, using (3.43), the corresponding streamfunction reads
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn69.png?pub-status=live)
and thus
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn70.png?pub-status=live)
in the $E^{2/7}$-layer, where
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn71.png?pub-status=live)
and
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn72.png?pub-status=live)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig12.png?pub-status=live)
Figure 12. $(a)$ Value of
$u_{\phi } E^{-2/7}$ as a function of the stretched cylindrical radial coordinate
$(s-\eta )E^{-2/7}$ for various Ekman numbers,
$\eta =0.35$,
$z=0.7$ and
$A=0.01$. The black dashed line corresponds to the analytical geostrophic solution (3.49).
$(b)$ Angular velocity in the rotating frame of reference as a function of the cylindrical radial coordinate
$s$. The dashed lines correspond to the analytical solutions (3.34a,b) and (3.49).
We verify the validity of the solution (3.54) in the $E^{2/7}$-layer away from the
$E^{1/3}$-layer in figure 13. We note that
$\psi$ and its derivatives remain discontinuous across
$s=\eta$, hence, contrary to the differential rotation (3.49), the solution (3.54) is not valid in the close vicinity of
$s=\eta$, and it is the ageostrophic
$E^{1/3}$-layer that smooths them out. Furthermore, since
$3/7>1/3$, the
$z$-component of the velocity is expected to be maximum in the
$E^{1/3}$ layer.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig13.png?pub-status=live)
Figure 13. Value of $\psi E^{-5/7}$ as a function of the cylindrical radial coordinate
$s$ for various Ekman numbers,
$\eta =0.35$,
$z=0.7$ and
$A=0.01$. The dashed lines correspond to the analytical geostrophic solutions (3.54).
We further verify the scalings for the thickness of the two Stewartson shear layers, as well as that of the amplitude of their associated $s$- and
$z$-directed velocity numerically, in figures 14 and 15. Although the amplitude of the
$z$-directed velocity is, as expected, maximum in the central Stewartson layer, we find a slight discrepancy with the expected power laws for the velocity near
$s=\eta$. This discrepancy is, in fact, not surprising as these two Stewartson shear layers are nested. The velocity in the region where the two coexist, that is the
$E^{1/3}$-layer, is thus expected to follow both (3.41a,b) and (3.55). In practice, we measure
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn73.png?pub-status=live)
which lies in between the velocity scalings in the two layers. For simplicity's sake, we may consider the amplitude scaling of meridional velocity in the $E^{1/3}$-layer to be that of (3.41a,b), that is
$u_s=O(E^{2/3})$ and
$u_z=O(E^{1/3})$.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig14.png?pub-status=live)
Figure 14. Value of $u_s E^{-2/3}$
$(a)$ and
$u_z E^{-1/3}$
$(b)$ as a function of the stretched cylindrical radial coordinate
$(s-\eta )E^{-1/3}$ for various Ekman numbers,
$\eta =0.35$,
$z=0.7$ and
$A=0.01$.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig15.png?pub-status=live)
Figure 15. Value of $u_s E^{-5/7}$
$(a)$ and
$u_z E^{-3/7}$
$(b)$ as a function of the stretched cylindrical radial coordinate
$(s-\eta )E^{-2/7}$ for various Ekman numbers,
$\eta =0.35$,
$z=0.7$ and
$A=0.01$. The black dashed lines correspond to the solution (3.55).
Ultimately, we find our results to be consistent with the scaling of the kinetic energy $E_{k, tot}\propto E$, which suggests that the main contribution comes from a
$E^{1/3}$-velocity field spread over a free shear layer of width
$E^{1/3}$. This
$E^{1/3}$ amplitude scaling of the velocity field in the Stewartson layer is actually verified throughout the transient, as can be seen in figure 11(b). This implies that despite a steady state seemingly not being reached in the lifetime of massive stars, the time relevant for the advective transport of angular momentum and of chemicals within their radiative envelope scales as
$E^{-2/3}$, since the whole meridional velocity grows on a viscous time scale (e.g. figure 3).
3.3. Conclusions
The picture which results from the foregoing discussion is rather neat. The spherical shell is split into three domains:
(i) The domain outside the tangent cylinder
$\mathcal {C}$,
$s>\eta$ where at first order, and outside any layer,
(3.59a–c)Thus for a braking torque,\begin{equation} u_s = -\frac{2KE}{3s}, \quad u_{\phi} = -\frac{K}{3} \left( s \ln s/\eta - \frac{1}{s} + \frac{s}{\eta^2} \right), \quad u_z=0. \end{equation}
$u_{\phi }<0$ and
$u_s<0$ so that there is a weak radial flow towards the Stewartson layer.
(ii) The domain inside the tangent cylinder,
$s<\eta$ where we find that
(3.60a–c)\begin{equation} u_s=0, \quad u_{\phi}=- \frac{ \sqrt{2E}Ks}{\sqrt{1-s^2}} \left(1-\frac{s^2}{\eta^2}\right)^{1/4}, \quad u_z = -\frac{E K}{\sqrt{1-s^2}} \left( 1 + \frac{1}{1-s^2}\right). \end{equation}
(iii) In between, the Stewartson layer is composed of two nested shear layers of thicknesses
$O(E^{2/7})$ and
$O(E^{1/3})$, and is dominated by a flow parallel to the rotation axis and directed towards the outer shell. In the
$E^{2/7}$-layer,
(3.61a–c)and in the\begin{equation} u_s\sim E^{5/7}, \quad u_{\phi} \sim E^{2/7}, \quad u_z \sim E^{3/7} \end{equation}
$E^{1/3}$-layer,
(3.62a–c)Hence the meridional circulation is completely dominated by the shear flow of the\begin{equation} u_s\sim E^{2/3}, \quad u_{\phi} \sim E^{1/3}, \quad u_z \sim E^{1/3}. \end{equation}
$E^{1/3}$ Stewartson layer.
4. The role of thermal stratification
4.1. Description
Stable stratification introduces a restoring buoyancy force that acts to inhibit vertical motions and may partially or entirely eliminate the control of the flow exercised by the Ekman layers. Indeed, the tendency for viscous layers to generate secondary circulation is counteracted by the stable stratification. Thus, we may wonder to what extent the flow inside the stably stratified radiative envelope of a massive star losing angular momentum can remain properly modelled with a geostrophic solution. To investigate this question we introduce a radially stable stratification and use the Boussinesq approximation. We still neglect the centrifugal acceleration. The linearised and dimensionless vorticity, mass conservation and heat equations read
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn78.png?pub-status=live)
where we have used $\Delta T (2\varOmega _c)^2/ \mathcal {N}^2$ as the temperature scale;
$\Delta T$ is the temperature difference between outer and inner shells, and
$\mathcal {N}$ is the Brunt–Väisälä frequency defined as
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn79.png?pub-status=live)
where $\alpha$ is the coefficient of thermal expansion and
$\boldsymbol {g}= -g {\boldsymbol {r}}$ is the gravitational field. We have assumed the equilibrium temperature gradient to be produced by a uniform distribution of heat source (Chandrasekhar Reference Chandrasekhar1961), which implies that
$\boldsymbol {\nabla } T_{eq, *} = (\Delta T / R^2) r {\boldsymbol {e}}_r$;
$\delta T$ is the temperature perturbation from such thermal equilibrium, and
$Pr= \nu /\kappa _T$ is the Prandtl number, which only enters the equations of motion combined with
$S= \mathcal {N}^2/(2 \varOmega _c)^2$, as actually noted by Barcilon & Pedlosky (Reference Barcilon and Pedlosky1967). Here,
$\kappa _T$ is the heat diffusivity that we have assumed to be constant. In what follows, we use the dimensionless parameter introduced by Garaud (Reference Garaud2002), namely
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn80.png?pub-status=live)
as the parameter characterising the effective strength of thermal stratification.
We complete these equations with the boundary conditions (2.5) for the velocity field, and we impose $\delta T=0$ at the boundaries. The initial conditions are
${\boldsymbol {u}}=\boldsymbol {0}$, and
$\delta T=0$. Projecting (4.1) on spherical harmonics yields the following system of equations for radial parts
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn81.png?pub-status=live)
where
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn82.png?pub-status=live)
Figure 16 shows the differential rotation in the rotating frame of reference and the meridional circulation for $E=10^{-7}$ and various
$\lambda$ once the steady state is settled. We see that the flow departs from the quasi-geostrophic equilibrium as
$\lambda$ increases. At high
$\lambda$ the angular velocity profile thus becomes shellular, namely it only depends on the radial distance.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig16.png?pub-status=live)
Figure 16. Meridional view of the streamfunction $\psi$ (a,c,e,g) and of the differential rotation in the rotating frame of reference
$\delta \varOmega = (\varOmega -\varOmega _c)/ (2\varOmega _c$) (b,d,f,h) for
$E=10^{-7}$,
$\eta =0.35$,
$A=0.01$ and various
$\lambda$.
4.2. Horizontal boundary layers
In the case of a thermally homogeneous fluid, we have seen that horizontal boundary layers, and specifically Ekman layers play a crucial role on the interior dynamics of the flow. Hence, one may wonder how such layers impact the interior flow in a thermally stratified fluid. We first consider the condition of existence of such horizontal boundary layers and determine their thickness scaling with the relevant adimensional parameters.
Inside boundary layers, horizontal length scales are much larger than the transverse one. Radial derivatives are therefore expected to prevail. We thus write the simplified equations, assuming $\partial _r \gg \partial _{\theta }$ in horizontal boundary layers
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn83.png?pub-status=live)
Combining these equations, only keeping the highest-order derivatives for each coefficient ($E^2$,
$\lambda$ and 1) yields the general horizontal boundary layer equation (see also Raze, Lignières & Mimoun Reference Raze, Lignières and Mimoun2017)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn84.png?pub-status=live)
Introducing the $O(1)$ stretched radial coordinate
$\zeta = (1-r)/\sqrt {E}$, (4.7) can be re-written, in the outer boundary layer and in the asymptotic regime of small Ekman numbers,
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn85.png?pub-status=live)
We see that if $\lambda \ll 1/E$, we get the actual Ekman layer where the Coriolis force is balanced by the viscous shear. As remarked by Barcilon & Pedlosky (Reference Barcilon and Pedlosky1967), we note that, whenever Ekman layers are present, their structure is independent of the stratification strength because the forces in balance in such layers essentially involve horizontal motions. Similarly, if
$\lambda \gg 1$, introducing the stretched radial coordinate
$\gamma = (1-r)\sqrt {\lambda }$, (4.7) can be re-written, in the outer boundary layer,
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn86.png?pub-status=live)
Hence, if $1 \ll \lambda ^2 \ll 1/E^2$, the viscous term drops out and buoyancy balances the Coriolis force. This is typical of a thermal boundary layer of width
$\delta _{\lambda }= O(1/\sqrt {\lambda })$. We note that in the parameter regime
$\lambda \ll 1$, only Ekman layers are present at the boundaries, while in the regime
$1 \ll \lambda \ll 1/E$ both thermal and Ekman boundary layers coexist. In this latter case thermal boundary layers are always much thicker than Ekman layers.
4.3. Massive stars interior flows
From the foregoing discussion, it turns out that the dynamics of the flows may be quite different whether $\lambda \gg 1$ or
$\lambda \ll 1$. Let us now place the case of rapidly rotating massive stars.
We note that, in the envelope of massive stars, the radiative viscosity and radiative heat diffusion largely dominate diffusion of collisional origin (Espinosa Lara & Rieutord Reference Espinosa Lara and Rieutord2013). These two quantities read
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn87.png?pub-status=live)
where $T$ is the temperature,
$\kappa _R$ is the Rosseland mean opacity,
$a = 4 \sigma /c$ is the radiation density constant with
$\sigma$ the Stefan–Boltzmann constant,
$c$ the speed of light and
$c_p$ is the specific heat capacity at constant pressure. In such a radiation-dominated system, the radiative Prandtl number reads
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn88.png?pub-status=live)
where $c_s$ is the adiabatic sound speed, showing that, naturally,
$Pr_{rad}\ll 1$.
However, the differential rotation of the radiative envelope is expected to drive some small-scale turbulence. Zahn (Reference Zahn1992) proposes that the turbulent vertical kinematic viscosity associated with marginal shear stability reads as
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn89.png?pub-status=live)
where $Ri_c \simeq 1/4$ is the critical Richardson number. In figure 17, we plot the radiative and turbulent Prandtl numbers radial profiles at the equator of a 15 solar mass star as given by a 2-D ESTER model (Rieutord et al. Reference Rieutord, Espinosa Lara and Putigny2016) for various rotation rates defined as the ratio
$\omega$ between the equatorial angular velocity and the equatorial Keplerian angular velocity. The associated radial profiles of
$\lambda$ are also shown. Although the turbulent Prandtl number is a few orders of magnitude larger than the radiative one in fast rotating stars, turbulent
$\lambda$ is just one order of magnitude larger.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig17.png?pub-status=live)
Figure 17. The $Pr$ and
$\lambda$ radial profiles at the equator measured from ESTER 2-D stellar models of a
$15\ M_{\odot }$ star, for various angular velocity ratios
$\omega$. The radiative Prandtl number
$Pr_{rad}$ and the corresponding
$\lambda$-parameter are represented in full lines, and their turbulent counterparts are represented in dashed lines.
We find that the turbulent $\lambda$-parameter is roughly independent of
$\omega$, and never exceeds
$10^{-3}$ in the considered models. The radiative
$\lambda$-parameter, on the other hand, increases for decreasing
$\omega$ to the point it may locally exceed
$10^{-4}$ when
$\omega \lesssim 0.1$. As massive stars are often considered to be fast rotators, we conclude that the thermal stratification regime of their radiative envelope corresponds to
$\lambda \ll 1$.
4.4. Asymptotically weak thermal stratification regime
We now study the case of weak temperature stratification, relevant to the radiative envelope of rotating massive stars. In other contexts the $\lambda \ll 1$ regime may also be reached if the Brunt–Väisälä frequency is just small. We have seen that in this regime, only Ekman layers characterised by an
$O(E)$ radial velocity at their edge, are present at the boundaries. We thus assume
$\lambda \ll 1$ for the weak stratification, but also suppose that
$E\ll \lambda$ as expected in stars.
4.4.1. The steady flow
We first focus on the steady flow. The radial component of the steady momentum equation reads
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn90.png?pub-status=live)
where differential rotation is driven by the $O(1)$ surface stress, that is
$u_{\phi }=O(1)$. In the asymptotic regime of small Ekman number, the temperature deviation from equilibrium is therefore, at most,
$O(1)$. If that is the case, namely if
$\delta T = O(1)$, the heat equation (4.1) implies that the interior radial velocity
$u_r$ is
$O(E/\lambda )$, namely larger than the
$O(E)$ Ekman pumping (3.13), in the considered stratification regime. Because of the boundary conditions imposed by Ekman layers, an
$O(E/\lambda )$ radial velocity must vanish at
$r=\eta$ and
$r=1$. Furthermore, the radial component of the vorticity equation reads
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn91.png?pub-status=live)
implying that $u_r$ is
$z$-independent up to
$O(E)$ corrections (see Barcilon & Pedlosky Reference Barcilon and Pedlosky1967). Hence, a consistent solution for the interior flow is
$u_r = O(E)$ and consequently
$\delta T = O(\lambda )$, which is completed by the classical Ekman boundary layer with
$\tilde {u}_r=O(E)$ and
$\tilde {u}_{\theta }=O(\sqrt {E})$.
Let us now write the variables of the interior flow, as a power expansion of the small parameter $\lambda$, in the asymptotic regime of small Ekman numbers. They read
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn92.png?pub-status=live)
Injecting (4.15) into the momentum equation then yields, at zeroth and first orders
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn93.png?pub-status=live)
and
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn94.png?pub-status=live)
Hence, at the lowest order, and identically to the homogeneous case, the Taylor–Proudman theorem is satisfied and the $O(1)$ geostrophic flow dynamics is consequently entirely controlled by the Ekman pumping. This is indeed verified in figure 16. On the other hand, the
$O(\lambda )$ flow is affected by temperature deviation from equilibrium produced by the
$O(E)$ radial velocity. We conclude that the thermal stratification regime in the radiative envelope of rotating massive stars has a negligible influence on the primary and secondary quasi-geostrophic and stationary flows.
4.4.2. The transient flow
As in the unstratified case, we now wish to determine the scaling of the transient time with the relevant non-dimensional parameters, that is, the Ekman number and the $\lambda$-parameter, in the
$E\ll \lambda \ll 1$ limit. We write the radial component of the vorticity equation
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn95.png?pub-status=live)
as well as the $z$-derivative of the heat equation
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn96.png?pub-status=live)
Adding (4.18) and (4.19), and using the $\phi$-component of the momentum equation, we get
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn97.png?pub-status=live)
Hence, if $\delta T$ remains
$O(\lambda )$ during the transient, we expect the quasi-geostrophic steady state to be reached on the
$O(E^{-1})$ viscous time scale in the limit
$\lambda \ll 1$.
We verify this conclusion with our numerical solution and measure the transient time scale $\tau _t$ as the time required for the relative difference between the torques at the boundaries to be less than
$0.01\,\%$. We show the scaled values
$E\tau _t$ as a function of
$\lambda$ for
$E=10^{-6}$ in figure 18. We find that indeed, as for the thermally homogeneous case, the transient time scale is
$O(E^{-1})$ in the weakly stratified limit.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig18.png?pub-status=live)
Figure 18. The scaled transient time scale $E\tau _t$ given by the numerical solution as a function of
$\lambda$ for
$E=10^{-6}$. The steady state is reached on a
$O(E^{-1})$ time scale, both in the strongly and weakly stratified limits.
The thermal stratification regime in the interior of slowly rotating stars may, on the other hand, lie in the $\lambda \gg 1$ limit. In this regime, it is well known that the Eddington–Sweet time scale associated with the angular momentum redistribution by meridional circulation is so long that it is unlikely that the system would ever relax to a steady state. Hence, the study of the time-dependent interior dynamics of slowly rotating stars, which may largely depend on initial conditions and on the damping rate of the baroclinic modes, is outside the scope of this work. However, for the sake of completeness and to appreciate some effects of a strong thermal stratification in the stress-driven barotropic spin-down flow, we give a short account of the
$\lambda \gg 1$ regime in appendix B.
5. The flow in a polytropic envelope
5.1. The background
The next step towards a realistic model of the radiative envelope of a massive star is to include the strong density variations of the fluid between the convective core and the stellar surface. Typically density varies by a factor of $10^9$ in the envelope of a star of 15 solar masses. To take this density distribution into account in a simple way we assume that the gas can be described by a polytropic equation of state as usually done in stellar physics (e.g. Maeder Reference Maeder2009). Hence the hydrostatic background state verifies
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn98.png?pub-status=live)
where $p_{*}$,
$\rho _*$ and
$r_*$ are the dimensional pressure, density and radial coordinate, respectively. As previously, we also use the Roche approximation and assume that the core gathers all the mass
$M$ of the star. In (5.1a,b)
$G$ is the gravitational constant,
$\kappa$ is a constant related to the thermal conditions at the core boundary and
$\varGamma =1+1/n$, where
$n$ is the so-called polytropic index. In the following we shall set
$n=3$, which is typical for radiative envelopes. We note that for such an index the fluid is stably stratified (e.g. Dintrans & Rieutord Reference Dintrans and Rieutord2001; Rieutord & Dintrans Reference Rieutord and Dintrans2002, for instance), but as concluded from the previous section the small value of the Prandtl number allows us to neglect buoyancy effects. The following results therefore assume no buoyancy force.
Equation (5.1a,b) can be easily solved and gives the density profile
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn99.png?pub-status=live)
where
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn100.png?pub-status=live)
Here $\rho _c = \rho _{*}(\eta R)$ is the density at the core boundary and
$\rho _s=\rho _{*}(R)/\rho _c$ is the adimensional surface density. Note that such a background is also used in numerical simulations of convection in stellar envelope (Raynaud et al. Reference Raynaud, Rieutord, Petitdemange, Gastine and Putigny2018) or planetary atmospheres (Gastine & Wicht Reference Gastine and Wicht2012). Hence, in this section, we solve (2.1) with the background polytropic density profile (5.2), the boundary conditions (2.5) and with the initial condition
${\boldsymbol {u}}= \boldsymbol {0}$.
5.2. The transient phase
Again, we first focus on the transient phase preceding the settling of a steady state. In particular, we aim at determining whether the scaling of the governing time scales are modified by density variations of the background. We measure the transient time as the time for which relative difference between the torque exerted by the fluid on the stationary inner sphere and the torque applied on the outer sphere, is less than $0.01 \%$. We monitor the evolution of the relative difference between the torques
$\Delta \varGamma / \varGamma (1)$ given by (3.6), for various surface density ratios. We show the resulting steady-state times in figure 19 for the polytropic indexes
$n=3$ and
$n=3/2$. The latter index corresponds to an isentropic monatomic gas, hence to a neutral thermal stratification. We find the steady-state time to depend on the density ratio between the core and the surface. For the polytropic index
$n=3$, relevant to the envelope of massive stars, we find
$E \tau _t \propto \rho _s^{0.3}$ for
$10^{-3} \leq \rho _s$, and
$E \tau _t \propto \rho _s^{0.09}$ for
$\rho _s < 10^{-4}$. Of course, these power laws quantitatively depend on the chosen density profile, as can be seen with the
$n=3/2$ case. Unlike thermal stratification, density stratification thus (mildly) influences the time scale required to reach a steady state in a stellar envelope.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig19.png?pub-status=live)
Figure 19. Scaled steady-state time scale $E\tau _t$ as a function of the surface density
$\rho _s$, for
$E=10^{-6}$ models with various values of
$\rho _s$,
$n=3/2$ and
$n=3$. The dashed and dotted lines correspond to the steady-state time scalings with
$\rho _s$.
5.3. The steady flow
We now focus on the steady flow. We may observe that if viscosity, nonlinearity and buoyancy are neglected then the steady flow obeys a Taylor–Proudman theorem applied to the momentum $\rho {\boldsymbol {u}}$ instead of the velocity, namely
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn101.png?pub-status=live)
Unfortunately, we cannot reiterate the boundary layer analysis of the incompressible case since the viscous force (2.3) now includes terms that depend on $r=\sqrt {s^2+z^2}$ and make the partial differential equation not separable. We may, however, observe that since the density variations do not add any new length scale, the viscous balance in the shear layers remains similar and we can still expect the presence of a Stewartson layer along the tangent cylinder.
We now revert to numerical solutions to make progress. We thus solve (2.1) with boundary condition (2.5). As already mentioned, we neglect the effects of buoyancy. We first focus on the differential rotation and the meridional circulation. It is convenient to introduce the streamfunction $\chi$ associated with the meridional momentum, namely
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn102.png?pub-status=live)
Figure 20 shows the differential rotation and this new streamfunction in the rotating frame of reference for $E=10^{-7}$.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig20.png?pub-status=live)
Figure 20. Meridional view of the streamfunction $\chi =\rho \psi$ (a) and of the differential rotation in the rotating frame of reference
$\delta \varOmega = (\varOmega -\varOmega _c)/2\varOmega _c$ (b) for
$E=10^{-7}$,
$\eta =0.35$,
$\rho _s=10^{-4}$ and
$A=0.01$.
We note that as for the incompressible flow, the amplitude of the secondary flow is dominated by the Stewartson layer located at the edge of the tangent cylinder $\mathcal {C}$. The differential rotation, on the other hand, and as expected, is no longer cylindrical and very
$z$-dependent, with a maximum value that is close to the outer shell and the tangent cylinder. Figure 21(a) shows the differential rotation as a function of the cylindrical radial coordinate for
$E=10^{-7}$ at different radii
$r$, illustrating the
$z$-dependence of
$u_{\varphi }$, while figure 21(b) shows that
$q_{\varphi }=\rho (r)u_{\varphi }$ verifies the Taylor–Proudman constraint.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig21.png?pub-status=live)
Figure 21. Normalised angular velocity (a) and normalised angular velocity multiplied by the density $\rho$ (b), and as a function of the cylindrical radial coordinate
$s$ for
$E=10^{-7}$,
$\eta =0.35$,
$\rho _s=10^{-4}$,
$A=0.01$ and various radii.
As far as meridional circulation is concerned, this is still an ${O}(E)$ flow outside the Stewartson layer as shown by figures 22(a) and 23(a). We note that outside the tangent cylinder, streamlines are no longer straight lines, even for the
${\boldsymbol {q}}$ momentum field. In the inner cylinder and near the rotation axis we still get streamlines parallel to the rotation axis for the momentum
${\boldsymbol {q}}$. This latter feature comes from the fact that the Ekman layer has the same structure as in the constant density case and drives an Ekman circulation which has a unique component along the
$z$-axis. From, (3.35a,b) we deduce that
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn103.png?pub-status=live)
in this region.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig22.png?pub-status=live)
Figure 22. $(a)$ Value of
$\chi E^{-1}$ as a function of the cylindrical radial coordinate
$s$ for various Ekman numbers,
$\eta =0.35$,
$r=0.7$,
$\rho _s=10^{-4}$ and
$A=0.01$. The predicted
$O(E)$ scaling of the secondary flow is verified away from the Stewartson nested layers.
$(b)$ Meridional view of the streamfunction
$\chi E^{-1}$ outside
$\mathcal {C}$, for the
$E=10^{-9}$ model.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig23.png?pub-status=live)
Figure 23. $(a)$ Value of
$u_z E^{-1}$ as a function of the cylindrical radial coordinate
$s$ for various Ekman numbers,
$\eta =0.35$,
$r=0.7$,
$\rho _s=10^{-4}$ and
$A=0.01$.
$(b)$ Meridional view of the streamfunction
$\chi$ inside
$\mathcal {C}$ and away from the
$O(E^{2/7})$ Stewartson layer, for the
$E=10^{-9}$ model.
Finally, the Stewartson layer seems to conserve its structure if we focus on the momentum, as shown by figures 24 and 25. Hence, the $E^{1/3}$ scale is still the dominating scale.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig24.png?pub-status=live)
Figure 24. Value of $q_s E^{-2/3}$ (a) and
$q_z E^{-1/3}$ (b) as a function of the stretched cylindrical radial coordinate
$(s-\eta )E^{-1/3}$ for various Ekman numbers,
$\eta =0.35$,
$z=0.7$,
$\rho _s=10^{-4}$ and
$A=0.01$.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig25.png?pub-status=live)
Figure 25. Value of $q_s E^{-5/7}$ (a) and
$q_z E^{-3/7}$ (b) as a function of the stretched cylindrical radial coordinate
$(s-\eta )E^{-2/7}$ for various Ekman numbers,
$\eta =0.35$,
$z=0.7$,
$\rho _s=10^{-4}$ and
$A=0.01$.
6. Summary and conclusions
Aiming at a better description of the dynamics of the radiative envelopes of massive stars, which lose both mass and angular momentum through radiative winds, we have considered the problem of the spin-down flow of a viscous fluid inside a spherical shell. The spin-down is assumed to be driven by a tangential stress prescribed on the outer shell. This problem is quite close to the classical spherical Couette flow (Zikanov Reference Zikanov1996; Rieutord et al. Reference Rieutord, Triana, Zimmerman and Lathrop2012), but includes new features like density and thermal stratification, or stress driving, that needed new investigations.
To start with, we therefore considered this problem in the case of a mild driving so as to be able to deal with linear equations. After examining the case of constant density, we investigated the role of both density and thermal stable stratification. Since stars are large-size bodies, the Ekman number is very small. We therefore restricted our analysis to small Ekman numbers.
Our numerical solutions showed that the meridional kinetic energy is concentrated in the Stewartson shear layer that is tangent to the inner shell, both for incompressible and anelastic stationary flows assuming no thermal stratification. Outside this layer, a boundary layer analysis of the Ekman layers allowed us to exhibit analytical solutions for the quasi-geostrophic and incompressible primary and secondary stationary flows. We found the latter flow to be essentially perpendicular to the rotation axis outside the tangent cylinder, flowing from the outer Ekman layer towards the Stewartson layer. Inside the tangent cylinder, we found this poloidal flow to be parallel to the rotation axis, being pumped into the Ekman layer attached to the inner core. The mass flux is then returned towards the outer Ekman layer through the Stewartson layer. Outside of the Stewartson layer, the amplitude of the meridional flow (Ekman circulation) scales like $E$, as a consequence of the surface stress driving. In our model, the Stewartson layer is composed of two nested free shear layers of thicknesses
$O(E^{2/7})$ and
$O(E^{1/3})$. The analysis of the quasi-geostrophic
$E^{2/7}$-layer located inside the tangent cylinder allowed us to derive asymptotic solutions for the primary and secondary flows, and a simple analysis of the ageostrophic
$E^{1/3}$-layer indicates that it dominates the entire meridional circulation with a maximum amplitude of the
$z$-directed velocity scaling as
$E^{1/3}$.
We then accounted for a stable thermal stratification by introducing a radial temperature gradient using the Boussinesq approximation. Two limits of the parameter $\lambda =Pr \mathcal {N}^2/(2\varOmega _c)^2$ show up. In the limit of strong thermal stratification,
$\lambda \gg 1$, the angular velocity profile becomes shellular (only radially dependent), while the circulation is concentrated in thermal boundary layers and the radial motion is strongly inhibited outside of them. As a consequence, the Stewartson layer is suppressed. A more in-depth consideration of this regime, although certainly relevant to slowly rotating stars, has been deliberately left aside. Indeed, such slow rotators may never relax to a steady state and thus might crucially depend on initial conditions. The study of the
$\lambda \gg 1$ asymptotic regime, including the time-dependent baroclinic flow, certainly calls for a separate investigation. However, for rapidly rotating stars the relevant limit is
$\lambda \ll 1$. In this case both the structure and amplitude of the incompressible and stationary flow are unaffected by the thermal stable stratification.
Radiative envelopes of massive stars exhibit, however, strong density variations between the convective core and the surface, making the incompressible and Boussinesq approximations too restrictive. But since thermal stratification has little impact for our stars, we can neglect buoyancy in the dynamics and concentrate on the effects of density stratification through the anelastic approximation. We chose to describe the gas of the stellar envelope by a polytropic equation of state with polytropic index $n=3$. We found that, assuming no further thermal stratification, the amplitude scaling of the stationary flow remains
$O(E)$ outside the Stewartson layer and
$O(E^{1/3})$ inside. However, outside the tangent cylinder, the streamlines are no longer straight lines and the flow becomes very dependent on the coordinate parallel to the rotation axis.
Because of the weakness of the meridional circulation, the stationary state of the flow is not quite relevant because it implies very long transients that might actually exceed the lifetime of the star. Massive stars are indeed short lived (a few million years). We therefore checked the time scale associated with the transient phase and found it to be governed by the viscous diffusion time, that is $O(E^{-1})$ indeed. The density variations of the hydrostatic background seem to shorten it slightly. This implies that such a steady flow is typically reached on a time that is longer than the lifetime of the star. Fortunately, though the flow structure outside the tangent cylinder evolves during the transient, the dominating
$E^{1/3}$ amplitude scaling of the meridional velocity stands. Hence, for all the models considered, the time relevant to the advective transport of chemicals within the radiative envelope of massive stars scales as
$E^{-2/3}$, which is much shorter than the advection of angular momentum acting on
$O(E^{-1})$ time scale.
The conclusion of the foregoing study is that the Stewartson layer is a key feature for the transport of chemical elements between the core and the surface of a massive star. The advective time scale is indeed $O(E^{-1/3})$ shorter than the spin-down time scale. In our case indeed, spin-down is driven by a stress and the steady state is reached on a viscous diffusion time, which is longer than the star's life. The much shorter time associated with the meridional current of the Stewartson layer allows chemical elements produced in the core to be transported to the surface of the star and be possibly observable. We have shown that neither the stable stratification of the envelope, nor its strong density variations inhibit the rise of the Stewartson layer. The stable stratification is bypassed thanks to the low Prandtl number, while density variations have little influence on the mass flux
$\rho {\boldsymbol {u}}$ basically because they occur on a large scale.
These conclusions are not the end of the story of course since other effects might complicate the scenario. The first effect one might think of is the local anisotropic turbulence of the envelope. Indeed, the differential rotation induces a local shear that is unstable. This shear instability is reduced by the stable stratification (Richardson criterion) but eased by the large heat diffusion. Zahn (Reference Zahn1992) and Maeder & Zahn (Reference Maeder and Zahn1998) have argued that such instabilities will lead to a strongly anisotropic turbulence, making horizontal transport much more efficient than the vertical one. It is therefore an open question whether the Stewartson layer can resist to an anisotropic turbulent diffusion and in which circumstances.
If we still wish more realism, the interface between the core and the envelope will need a more detailed description. In this region, stable chemical stratification builds up in the course of stellar evolution. This so-called $\mu$-barrier might isolate the core from the envelope. However, the thermal gradient is still unstable and overstable double-diffusive convection is suspected to develop there (Garaud Reference Garaud, Rieutord, Baraffe and Lebreton2020). The impact of rotation on the dynamics of this region is almost unknown and will also motivate future studies. Finally, magnetic fields, of fossil or core dynamo origin, may have a dramatic impact on a massive star's interior dynamics. Indeed, depending on the field geometry, it can, for instance, amplify (Hollerbach Reference Hollerbach1997) or suppress the Stewartson layer and yield a super-rotating jet (Kleeorin et al. Reference Kleeorin, Rogachevskii, Ruzmaikin, Soward and Starchenko1997; Dormy et al. Reference Dormy, Cardin and Jault1998). Their consideration is also a matter for future work.
Acknowledgements
We are very grateful to the referees for their insightful remarks which allowed us to improve the presentation and the depth of our work. The numerical computations have been performed using the HPC resources from CALMIP (grant 2016-P0107), which is gratefully acknowledged.
Declaration of interests
The authors report no conflict of interest.
Appendix A. The linear approximation
A.1. Validity condition
This appendix presents an a posteriori justification for the use of linear approximation. For simplicity, we focus on the incompressible flow. In spherical coordinates and in the inertial frame, the nonlinear dimensional momentum equation reads
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn104.png?pub-status=live)
where
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn105.png?pub-status=live)
is the dimensional viscous force, and $\mu =\rho _{*}(r) \nu$ is the dynamical viscosity. Let us decompose the velocity field as the sum of the bulk component
$\varOmega _c {\boldsymbol {e}}_z \times {\boldsymbol {r}}$ and a residual velocity field
${\boldsymbol {u}}_{*}^r$
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn106.png?pub-status=live)
For an axisymmetric flow, the nonlinear term can be rewritten
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn107.png?pub-status=live)
We further note that the centrifugal term derives from a potential that can be gathered with the pressure into $\varPi$. Finally, we write the nonlinear adimensional momentum equation in the rotating frame
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn108.png?pub-status=live)
From now on the adimensional velocity field in the rotation frame ${\boldsymbol {u}}^r$ will be written
${\boldsymbol {u}}$ for simplicity. We now make explicit the nonlinear term as
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn109.png?pub-status=live)
The Ekman boundary layer analysis presented in § 3.2.1 revealed $u_r$ and
$u_{\theta }$ to be of order
$E$ (except in the Stewartson layer). For an axisymmetric flow, and neglecting
$O(E^2)$ terms, (A 6) can be simplified as
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn110.png?pub-status=live)
The Coriolis acceleration, in turn, reads
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn111.png?pub-status=live)
Since we work with the vorticity equation, we need
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn112.png?pub-status=live)
and
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn113.png?pub-status=live)
where
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn114.png?pub-status=live)
and
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn115.png?pub-status=live)
We note that, in the steady state, outside boundary layers $\partial u_{\phi } / \partial z = O(E)$ according to the Taylor–Proudman theorem, we can thus write the amplitude scaling of both terms
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn116.png?pub-status=live)
Hence, if $u_{\phi } = O(1)$, nonlinear terms in the vorticity equation are of the same order as the curl of the Coriolis acceleration and therefore cannot be neglected. However, inside the tangent cylinder, we have found, in § 3.2.2, that
$u_{\phi }= O(\sqrt {E})$. In that region the nonlinear terms can be neglected. Outside the tangent cylinder, the amplitude of the azimuthal velocity scales as
$A$ (see (3.17)), and the nonlinear terms are therefore
$O(A E)$. Hence, for sufficiently small
$A$, we expect the linear approximation to be relevant to our problem. We show the ratio between the Euclidian norm of (A 9) and that of (A 10) calculated a posteriori, for
$A=0.01$ and
$A=1$, in figure 26. We find, as expected, this ratio to be of order
$A$ outside the tangent cylinder. The linear solution can therefore be considered as a good approximation to the nonlinear and incompressible vorticity equation provided
$A \ll 1$.Similarly, for the anelastic flow, we find
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn117.png?pub-status=live)
where this time $u_{\phi } = O(A \rho _s)$ outside the tangent cylinder. Since
$\rho _s\ll 1$, the conditions of linearity are more easily met.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig26.png?pub-status=live)
Figure 26. Meridional view of the ratio between the Euclidian norm of the vorticity equation nonlinear term and that of the curl of the Coriolis acceleration, for $A=0.01$ (a) and
$A=1$ (b), and for an incompressible flow.
A.2. Application to massive stars
Let us now determine a typical value for $A$ when considering the differential rotation to be driven by a radiative wind at the surface of a massive star. We consider the outwardly accelerated outer layers to exert a viscous stress on the underlying layers of the star, spinning them down. Hence, we assume the local vertical angular momentum flux from the outwardly accelerated flow at the stellar surface to amount for the torque applied to the inner layers of the stars. That is
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn118.png?pub-status=live)
where $\dot {\ell }_z \equiv \dot {m} \varOmega _s (R \sin \theta )^2$ is the local angular momentum flux,
$\dot {m}$ is the associated local mass flux that is assumed isotropic,
$R$ is the stellar radius and
$[\sigma _*]$ is the dimensional stress tensor. Combining the expression for the imposed azimuthal stress (2.5) with (A 15) yields the adimensional amplitude of the surface stress resulting from the outward angular momentum flux
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn119.png?pub-status=live)
where $\rho _{s,*} = \rho _s \rho _c$ is the dimensional surface density. We compute all quantities with the ESTER 2-D code (Rieutord et al. Reference Rieutord, Espinosa Lara and Putigny2016; Gagnier et al. Reference Gagnier, Rieutord, Charbonnel, Putigny and Espinosa Lara2019a) for a
$15\ M_{\odot }$ stellar model with a rotation period of one day and we assume a (turbulent) kinematic viscosity at the surface
$\nu = 10^{12}\ \textrm {cm}\,^{2}\,\textrm {s}^{-1}$ (Zahn Reference Zahn1992; Espinosa Lara & Rieutord Reference Espinosa Lara and Rieutord2013). We plot the amplitude of the surface stress weighed by surface density
$A\rho _s$ as a function of co-latitude in figure 27. We further note that rotating stars are not strictly spherically symmetric because of the centrifugal force, in particular when rotation is rapid. Hence, the stellar radius as well as the surface density have latitudinal dependencies. Besides the resulting slight variation on the stellar surface, we see that
$A \rho _s$ is of the order
$5 \times 10^{-8} \ll 1$ for this model. Hence, according to the results of appendix 1, the linear approximation can be used to model massive stars losing angular momentum, provided the kinematic viscosity at the surface is no less than
${\sim }10^7\ \textrm {cm}^{2}\,\textrm {s}^{-1}$, a typical value in radiation-dominated surface layers of massive stars (Espinosa Lara & Rieutord Reference Espinosa Lara and Rieutord2013).
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig27.png?pub-status=live)
Figure 27. Amplitude of the surface stress weighed by the adimensional surface density $A\rho _s$ resulting from the angular momentum outward flux, as a function of co-latitude, for a
$15\ M_{\odot }$ 2-D ESTER stellar model rotating with a period of one day.
Appendix B. Asymptotically strong thermal stratification regime
In this appendix, we study the case of strong temperature stratification, that is in the $\lambda \gg 1$ asymptotic regime. We have seen, in § 4.2, that in this regime, both Ekman and thermal horizontal boundary layers coexist, with respective thickness of order
$\sqrt {E}$ and
$1/\sqrt {\lambda }$.
B.1. The steady flow
In the asymptotic regime of small Ekman numbers, we have seen that because $u_{\phi } = O(1)$, the temperature deviation from equilibrium is at most
$O(1)$. Unlike the weakly stratified case, however, the associated
$O(E/\lambda )$ interior radial velocity is less than the
$O(E)$ stratification independent Ekman pumping, and is therefore a consistent solution for the interior flow. Let us now write the dynamical variables in the interior, as a power expansion of the small parameter
$1/\sqrt {\lambda }$ corresponding to the width of the thermal boundary layer, in the asymptotic regime of small Ekman numbers. They read
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn120.png?pub-status=live)
Injecting (B 1) in the momentum equation then yields the $O(1)$ and
$O(1/\sqrt {\lambda })$ interior equations, that is
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn121.png?pub-status=live)
and
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn122.png?pub-status=live)
Hence, for large values of $\lambda$, thermal stratification inhibits vertical motions in the interior, and the Ekman layer pumping/suction no longer controls the dynamics of the interior flow. The Taylor–Proudman is not verified and is replaced by the thermal wind equation
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn123.png?pub-status=live)
to $O(E {\boldsymbol {u}})$, and the
$O(1)$ azimuthal velocity can be obtained solving
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn124.png?pub-status=live)
This yields, using boundary conditions (2.13)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn125.png?pub-status=live)
Hence, at the lowest order, $\delta \varOmega =u_{\phi }/(r \sin \theta )$ is shellular in the asymptotic regime of large
$\lambda$. Figure 28 shows the angular velocity radial profiles at
$\theta = {\rm \pi}/2$, for
$E=10^{-7}$ and for various
$\lambda$. We see that the asymptotic limit where the
$O(1/\sqrt {\lambda })$ terms can be neglected corresponds to
$\lambda \gtrsim 10^4$.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig28.png?pub-status=live)
Figure 28. Angular velocity radial profiles in the rotating frame $\delta \varOmega$, at
$\theta ={\rm \pi} /2$,
$E=10^{-7}$ and for various
$\lambda$. The black dashed line corresponds to the asymptotic solution (B 6).
Another interesting feature of the $\lambda \gg 1$ regime is that the circulation is concentrated in the boundary layers (see also Rieutord Reference Rieutord2006; Raze et al. Reference Raze, Lignières and Mimoun2017). This can be seen in figure 16 for
$\lambda =10^2$. In the interior, the radial motion is so strongly inhibited by thermal stratification that it prevents the existence of the Stewartson layer, thus undermining rotation-induced advection.
In the thermal layers, radial gradients are increased by a factor $\sqrt {\lambda }$ which gives
$\hat {u}_r= O(E)$, and from the continuity equation
$\hat {u}_{\theta }= O(\sqrt {\lambda }E)$. In figure 29, we plot
$u_{\theta }/(E \sqrt {\lambda })$ as a function of the stretched radial coordinate
$(r-\eta )\sqrt {\lambda }$, for various combinations of Ekman numbers and
$\lambda$ parameters. We find that, indeed, the latitudinal velocity scales as
$\sqrt {\lambda E}$ in the outer part of the layer of thickness
$\delta _{\lambda }$, that is in the thermal boundary layer region, outside the Ekman layer (typically for
$(r-\eta )\sqrt {\lambda } \gtrsim 0.5)$. Additionally, we verify the
$O(\sqrt {E})$ amplitude scaling of the latitudinal velocity in the Ekman boundary layer. Indeed, figure 29 shows that provided fixed values of
$\delta _E / \delta _{\lambda }$, that is fixed values of
$\sqrt {\lambda E}$, the latitudinal velocity scales as
$\sqrt {\lambda }E$ in the Ekman boundary layer as well. However, taking
$\sqrt {\lambda E} = D$, where
$D$ is some constant, implies that
$\sqrt {\lambda }E = D \sqrt {E}$, hence
$\hat {u}_{\theta } \propto \tilde {u}_{\theta } = O(\sqrt {E})$.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_fig29.png?pub-status=live)
Figure 29. Value of $u_{\theta } ( \sqrt {\lambda }E)^{-1}$ as a function of the stretched radial coordinate
$(r-\eta )\sqrt {\lambda }$, for various combinations of Ekman numbers and
$\lambda$ parameters,
$\eta =0.35$,
$\theta ={\rm \pi} /8$ and
$A=0.01$.
B.2. The transient flow
We finally determine the scaling of the transient time with the Ekman number and the $\lambda$-parameter, in the
$\lambda \gg 1$ limit. Assuming that in this regime,
$\delta T$ remains
$O(1)$ during the transient, (4.18) simplifies to
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20201012175618520-0670:S0022112020007120:S0022112020007120_eqn126.png?pub-status=live)
indicating that the (non-geostrophic) steady state is reached on a $O(E^{-1})$ time scale as well. This is verified in figure 18. We note a weak
$\lambda$ dependence of this time scale in the intermediate stratification regime.