On the Dynamics of Overshooting Convection in Spherical Shells: Effect of Density Stratification and Rotation

There exists a large body of astrophysical literature on the topic of convective overshooting. Early studies were performed using modal expansions (e.g., Latour et al. 1976; Toomre et al. 1976; Latour et al. 1981) or within the context of mixing-length theory (see, for instance Shaviv & Salpeter 1973; Maeder 1975). However, as discussed in Renzini (1987), there were many inconsistencies among the results of the mixing-length theory studies due to their sensitivity to the underlying model assumptions. As a result, convective overshooting research came to rely on alternative, more sophisticated approaches, ranging from early modal-expansion models to modern, fully nonlinear computational models. These studies have been mainly focused on (1) the variation of the lengthscales associated with the overshooting motions in response to stellar parameters/properties and (2) the degree to which overshooting convection modifies the stable background stratification (thermal mixing) and the parameters under which it manages to turn part of it into an adiabatic layer (penetrative convection).

Early work focused on penetrative convection was performed by Zahn et al. (1982), who employed a modal-expansion approach to study overshooting in a Boussinseq system and found substantial penetration within the stable region. The effects of background density stratification were later included in the study of Massaguer et al. (1984), who demonstrated that increasing stratification in the CZ may lead to larger penetration lengthscales due to the enhanced flow asymmetries arising from compressibility effects. Scaling arguments have also been applied, as in Zahn (1991), who identified two dynamical regimes of interest below the base of the convective region: a penetrative layer for high Péclet numbers and a thermal adjustment (overshoot) layer for low Péclet numbers (where the Péclet number characterizes the strength of thermal advection relative to thermal diffusion). The existence of these two regimes was subsequently confirmed through the application of 2D, fully nonlinear simulations carried out in Cartesian and cylindrical geometries (e.g., Hurlburt et al. 1986, 1994; Rogers & Glatzmaier 2005; Rogers et al. 2006).

One of the first systematic studies of overshooting/penetrative convection via fully nonlinear simulations was that of Hurlburt et al. (1994; also see Hurlburt et al. 1986), who performed 2D Cartesian compressible calculations. That study explored the dependence of overshooting on the so-called stiffness parameter, a measure of the degree of relative stability of the RZ compared with the CZ. It also identified two scaling regimes, consistent with the findings of Zahn (1991), and showed that the penetrative regime was associated with lower stiffness parameters and larger penetration depths while the overshoot regime was related to higher stiffness parameters and possessed smaller overshoot lengthscales.

These results have largely been confirmed within the context of fully nonlinear 3D simulations. Studies carried out in both Cartesian and spherical geometries have demonstrated that the overshoot lengthscale decreases with increasing stiffness parameter (Brummell et al. 2002; Korre et al. 2019). However, as noted in Korre et al. (2019), the transition width between the CZ and the RZ also appears to play an important role in determining the amount of overshooting below the base of the convective region. Interestingly, these 3D studies did not identify pure penetration, namely the extension of the CZ further down into the stable RZ, but instead, they only observed partial thermal mixing in the RZ, associated with the stronger downflows.

This lack of pure penetration may be attributed to the lower filling factor of the plumes for 3D simulations relative to 2D simulations when run in otherwise similar parameter regimes (for more on this, see, e.g., Brummell et al. 2002; Rempel 2004). It may also result from the tendency of the overshoot lengthscale and degree of penetration to be established by the strongest downflows (see, e.g., Pratt et al. 2017; Korre et al. 2019). Both the strength of these downflows and the Péclet numbers achieved in simulations are explicitly dependent on the degree of thermal forcing as expressed through a Rayleigh number. The dependence of the overshoot lengthscale on the Rayleigh number has been explored in several studies (see Brummell et al. 2002; Rogers & Glatzmaier 2005; Korre et al. 2019), but different scalings were observed, possibly due to the different set-up configurations, geometries, and boundary conditions. Large density stratification, which enhances the asymmetry between upflows and downflows, may also lead to stronger overshooting and eventually an adiabatic penetrative layer within the RZ. This effect was discussed in Massaguer et al. (1984), but has yet to be investigated systematically via fully nonlinear, 3D global simulations.

Relatively few studies have explicitly examined the response of overshooting to the presence of rotation. While the details vary, these studies suggest that the overshooting depth δ decreases with decreasing Rossby number Ro (where Ro expresses the relative strength of inertial to Coriolis forces). Cartesian studies carried out in a rotating f-plane by Brummell et al. (2002) found that, for a limited range of Ro, δ ∝ Ro^0.15. Similar results were reported by Pal et al. (2007), who found that δ ∝ Ro^0.2 at the poles and δ ∝ Ro^0.4 at midlatitudes. Augustson & Mathis (2019) recently presented a model of rotating overshooting convection that employs the semi-analytical model of Zahn (1991) and estimated the scaling of the overshoot depth with the Rossby number. They also found that δ decreases with decreasing Rossby number and scales as Ro^0.3.

Overshooting in the Sun is thought to play a role in many aspects of the global flows and magnetic fields. For instance, temperature variations arising from rotating convective overshoot are thought to be responsible for inducing a thermal wind in the convection zone that leads to the observed tilt of differential rotation contours (Kitchatinov & Ruediger 1995; Rempel 2005; Miesch et al. 2006; Howe 2009). Moreover, overshooting may pump magnetic field into the solar tachocline region where it can be stretched into a large-scale toroidal field. As a result, the tachocline has been the presumed seat of the solar dynamo for some time, figuring particularly prominently in the flux-transport and interface class of dynamo models (e.g., Parker 1993; Dikpati & Charbonneau 1999; Rempel 2006).

For these reasons, a region of overshoot has been included in several 3D models of stellar convection. This includes studies within the solar context (e.g., Ghizaru et al. 2010; Brun et al. 2011; Racine et al. 2011; Guerrero et al. 2013, 2016; Augustson et al. 2015), as well as others related to more massive stars (e.g., Browning et al. 2004; Brun et al. 2005; Featherstone et al. 2009; Augustson et al. 2012, 2016). Throughout those studies, however, the region of overshoot was but one ingredient in numerical experiments targeting other aspects of the convection and its resulting dynamo. For instance, Brun et al. (2017) studied solar-like overshooting convection but mostly focused on the mean flows arising within different rotational regimes and the dependence of overshooting on stellar mass (also see Brun et al. 2011; Augustson et al. 2012). A systematic study of the effects of rotation on the convective overshooting dynamics, however, has yet to be performed in 3D global numerical calculations.

In this work, we use fully nonlinear 3D models of stratified convection in a spherical shell geometry to examine the dependence of convective overshooting on strong density stratification. We also develop scaling laws that describe the overshoot lengthscale with Rossby number and investigate the properties of angular momentum transport and meridional flow in the coupled CZ-RZ system. The paper is organized as follows. In Section 2, we describe our two-layered model formulation and provide the set of anelastic Navier–Stokes equations along with the initial conditions, and the boundary conditions used in our simulations. In Section 3, we present our results associated with the effect of the density stratification on the nonrotating overshooting dynamics. In Section 4, we study the dynamics related to the effect of rotation on the turbulent convective overshooting motions, as well as the dependence of the overshoot lengthscale on the Rossby number and on latitude. In Section 5, we focus on the rotational profiles arising at different Rossby numbers, the associated mean flows, and the angular momentum transport and balances within the convection zone and the overshoot region. Finally, in Section 6, we summarize our results, provide comparisons with previous work, and discuss the implications of these results in the solar context.

We aim to explore the dynamics associated with a two-zone system that consists of a convectively unstable region that lies above a stably stratified RZ. For all cases presented in this work, we use a fixed aspect ratio r_i /r_o = 0.45 where r_i is the inner radius and r_o is the outer radius of the spherical shell. Also, we assume that the convective region has a fixed depth L = r_o − r_c = 0.2408r_o , where r_c = 0.7592r_o is the radius at the bottom of the CZ. Our chosen depth L corresponds to the solar convection zone from ∼ 0.7187R_⊙ to ∼ 0.9467R_⊙, hence only accounting for its inner part and not the outer layers where radiative transfer processes take place. We are interested in the dynamics within deep stellar interiors where the flows are subsonic. As such, we employ the anelastic approximation that filters out the sound waves and assumes small perturbations of the thermodynamic variables compared with their mean (Gough 1969; Gilman & Glatzmaier 1981). Therefore, we solve the 3D Navier–Stokes equations under the anelastic and Lantz–Braginsky–Roberts approximations (Lantz 1992; Braginsky & Roberts 1995), the latter of which is exact in the convection zone where the reference state is adiabatic.

Then, the dimensional Navier–Stokes equations become

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {\boldsymbol{u}}}{\partial t}+{\boldsymbol{u}}\cdot {\rm{\nabla }}{\boldsymbol{u}}+2{{\rm{\Omega }}}_{o}\hat{z}\times {\boldsymbol{u}}\\ \ \ =\displaystyle \frac{g(r)}{{c}_{p}}S\hat{{\boldsymbol{r}}}-{\rm{\nabla }}(P/\bar{\rho })+\displaystyle \frac{1}{\bar{\rho }}{\rm{\nabla }}\cdot {\boldsymbol{D}},\end{array}\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot (\bar{\rho }{\boldsymbol{u}})=0,\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\begin{array}{l}\bar{\rho }\bar{T}\left(\displaystyle \frac{\partial S}{\partial t}+{\boldsymbol{u}}\cdot {\rm{\nabla }}S+{u}_{r}\displaystyle \frac{d\bar{S}}{{dr}}\right)\\ \ \ ={\rm{\nabla }}\cdot (\bar{\rho }\bar{T}\kappa {\rm{\nabla }}S)+Q+{\boldsymbol{\Phi }},\end{array}\end{eqnarray} \tag{ 3 }$$

where u = (u_r , u_θ , u_ϕ ) is the velocity field, P is the pressure, $\bar{\rho }$ is the reference density, $\bar{T}$ is the reference temperature, $d\bar{S}/{dr}$ is the background entropy gradient, S characterizes the entropy perturbations about the reference state, Ω_o is the frame rotation rate, g(r) is the gravity (where g ∝ 1/r²), and c_p is the specific heat at constant pressure. The viscosity ν and the thermal diffusivity κ are kept constant for simplicity. The viscous stress tensor D is given by

$$\begin{eqnarray}&&{\boldsymbol{D}}=2\bar{\rho }\nu ({e}_{{ij}}-\displaystyle \frac{1}{3}{\rm{\nabla }}\cdot {\boldsymbol{u}}),\end{eqnarray} \tag{ 4 }$$

and the viscous heating is denoted by

$$\begin{eqnarray}&&{\boldsymbol{\Phi }}=2\bar{\rho }\nu [{e}_{{ij}}{e}_{{ij}}-\displaystyle \frac{1}{3}{\left({\rm{\nabla }}\cdot {\boldsymbol{u}}\right)}^{2}],\end{eqnarray} \tag{ 5 }$$

where e_ij is the strain rate tensor. The internal heating term Q is proportional to the pressure such that ${L}_{\odot }=4\pi {\int }_{{r}_{i}}^{{r}_{o}}Q(r){r}^{2}{dr}$, where L_⊙ is the solar luminosity. As noted in Featherstone & Hindman (2016a), this profile is similar in structure to the one established in Model S (see, e.g., Christensen-Dalsgaard et al. 1996). We taper the profile using a tanh function below the base of the CZ such that

$$\begin{eqnarray}&&Q(r)\propto \displaystyle \frac{1}{2}\left(\tanh \left(\displaystyle \frac{r-{r}_{c}}{0.01{r}_{c}}\right)+1\right)(\bar{P}(r)-\bar{P}({r}_{o})),\end{eqnarray} \tag{ 6 }$$

where $\bar{P}$ is the reference pressure. It satisfies Q(r) = −∇ · F_rad, where F_rad is the radiative flux in the system given by

$$\begin{eqnarray}&&{F}_{\mathrm{rad}}(r)=\displaystyle \frac{1}{{r}^{2}}{\int }_{r}^{{r}_{o}}Q(r^{\prime} ){r}^{{\prime} 2}{dr}^{\prime} .\end{eqnarray} \tag{ 7 }$$

To close the set of Equations (1)–(3), we use an equation of state given by

$$\begin{eqnarray}&&\displaystyle \frac{\rho }{\bar{\rho }}=\displaystyle \frac{P}{\bar{P}}-\displaystyle \frac{T}{\bar{T}}=\displaystyle \frac{P}{\gamma \bar{P}}-\displaystyle \frac{S}{{c}_{p}},\end{eqnarray} \tag{ 8 }$$

assuming the ideal gas law

$$\begin{eqnarray}&&\bar{P}={\mathfrak{R}}\bar{\rho }\bar{T},\end{eqnarray} \tag{ 9 }$$

where ρ (T) are the density (temperature) perturbations about the reference state, ${\mathfrak{R}}$ is the gas constant, and γ = c_p /c_v is the heat capacity ratio (where c_v is the specific heat at constant volume).

To create our two-layered system, we can then select an appropriate $d\bar{S}/{dr}$ profile that satisfies an adiabatic polytropic solution with γ = 5/3 in the convective region r_c ≤ r ≤ r_o and that smoothly transitions to a stably stratified region for r_i ≤ r < r_c . We note, however, that there are more sophisticated approaches where the two-zone system can be created by assuming that the radiative heat flux depends on a varying (in density and temperature) thermal conductivity profile based on either computed and tabulated opacities or on approximations such as Kramers law (see, e.g., Käpylä 2019; Käpylä et al. 2019; Viviani & Käpylä 2021).

We also specify the number of density scale-heights in the convection zone through ${N}_{\rho }=\mathrm{ln}(\bar{\rho }({r}_{c}))/\bar{\rho }({r}_{o}))$; (see Appendix). Then, we integrate $d\bar{S}/{dr}$ from r_c inward, yielding a solution for $\bar{S}$. The integration constant is computed by matching the solution to the value of $\bar{S}$ at the base of the CZ.

Using the fact that the entropy for a monatomic ideal gas is $\bar{S}=\mathrm{ln}({\bar{P}}^{1/\gamma }/\bar{\rho })$, differentiating the latter with respect to the radius and applying the hydrostatic balance equation $\partial \bar{P}/\partial r=-\bar{\rho }g$ we obtain

$$\begin{eqnarray}&&\displaystyle \frac{\bar{\rho }}{{c}_{p}}\displaystyle \frac{d\bar{S}}{{dr}}+\displaystyle \frac{g}{\gamma }\exp (-\gamma \bar{S}/{c}_{p}){\bar{\rho }}^{(2-\gamma )}+\displaystyle \frac{d\bar{\rho }}{{dr}}=0,\end{eqnarray} \tag{ 10 }$$

which we can solve numerically for the reference density profile $\bar{\rho }(r)$ given any background entropy gradient profile.

We nondimensionalize Equations (1)−(3) using the depth of the convection zone [l] = L as the lengthscale, the viscous timescale [t] = L²/ν, and the velocity scale [u] = ν/L. Throughout the paper, we define the volume average of a quantity q(r, θ, ϕ) in the CZ as

$$\begin{eqnarray}&&{q}_{{cz}}=\displaystyle \frac{{\int }_{{r}_{c}}^{{r}_{o}}{\int }_{0}^{2\pi }{\int }_{0}^{\pi }{{qr}}^{2}\sin \theta d\theta d\phi {dr}}{{\int }_{{r}_{c}}^{{r}_{o}}{\int }_{0}^{2\pi }{\int }_{0}^{\pi }{r}^{2}\sin \theta d\theta d\phi {dr}}.\end{eqnarray} \tag{ 11 }$$

For the density and temperature scales, we choose $[\rho ]={\bar{\rho }}_{{cz}}$, and $[T]={\bar{T}}_{{cz}}$, respectively.

For the entropy perturbations S, we adopt a scaling related to the thermal energy flux F such that $[S]={{LF}}_{{cz}}/({\bar{\rho }}_{{cz}}{\bar{T}}_{{cz}}\kappa )$, where $F(r)={\int }_{{r}_{i}}^{r}Q(r^{\prime} ){r}^{{\prime} 2}{dr}^{\prime} /{r}^{2}$. We specify our chosen nondimensional $d\bar{S}/{dr}$ profile in the RZ such that

$$\begin{eqnarray}&&\displaystyle \frac{d\bar{S}}{{dr}}=\displaystyle \frac{A}{2}\left(1-\tanh \left(\displaystyle \frac{r-{r}_{b}}{d}\right)\right),\end{eqnarray} \tag{ 12 }$$

where A = 10, r_b = 0.6453r_o , and d = 0.0456r_o .

The nondimensional gravity (reference density $\bar{\rho }$, reference temperature $\bar{T}$) is equivalent to the dimensional gravity (reference density, reference temperature) divided by g_cz (${\bar{\rho }}_{{cz}}$, ${\bar{T}}_{{cz}}$). All of the quantities are now nondimensional, and the nondimensional anelastic Navier–Stokes equations become

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {\boldsymbol{u}}}{\partial t}+{\boldsymbol{u}}\cdot {\rm{\nabla }}{\boldsymbol{u}}+\displaystyle \frac{2}{\mathrm{Ek}}\hat{z}\times {\boldsymbol{u}}\\ \ \ =\displaystyle \frac{\mathrm{Ra}}{\Pr }g(r)S\hat{{\boldsymbol{r}}}-{\rm{\nabla }}(P/\bar{\rho })+\displaystyle \frac{1}{\bar{\rho }}{\rm{\nabla }}\cdot {\boldsymbol{D}},\end{array}\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot (\bar{\rho }{\boldsymbol{u}})=0,\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\begin{array}{l}\bar{\rho }\bar{T}\left(\displaystyle \frac{\partial S}{\partial t}+{\boldsymbol{u}}\cdot {\rm{\nabla }}S+{u}_{r}\displaystyle \frac{d\bar{S}}{{dr}}\right)\\ \ \ =\displaystyle \frac{1}{\Pr }{\rm{\nabla }}\cdot (\bar{\rho }\bar{T}{\rm{\nabla }}S)+\displaystyle \frac{1}{\Pr }{Q}_{{nd}}+\displaystyle \frac{\mathrm{Di}\ \Pr }{\mathrm{Ra}}{\boldsymbol{\Phi }},\end{array}\end{eqnarray} \tag{ 15 }$$

where Q_nd = LQ/F_cz is the nondimensional internal heating function. In Figure 1, we show the profile of Q_nd (r) against r/r_o for the run with N_ρ = 3.

We now have four nondimensional numbers: the flux Rayleigh number Ra, the Prandtl number Pr, the Ekman number Ek, and the dissipation number Di, which are defined respectively as

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{Ra} & = & \displaystyle \frac{{g}_{{cz}}{F}_{{cz}}{L}^{4}}{{c}_{p}{\bar{\rho }}_{{cz}}{\bar{T}}_{{cz}}{\kappa }^{2}\nu },\\ \Pr & = & \displaystyle \frac{\nu }{\kappa },\quad \mathrm{Ek}=\displaystyle \frac{\nu }{{{\rm{\Omega }}}_{o}{L}^{2}},\quad \mathrm{and}\quad \mathrm{Di}=\displaystyle \frac{{g}_{{cz}}L}{{c}_{p}{\bar{T}}_{{cz}}}.\end{array}\end{eqnarray} \tag{ 16 }$$

Note that, unlike Ra, Pr, and Ek, Di is not a free parameter. It depends on the chosen N_ρ and accounts for the fact that we have two energy scales—a thermal scale and a kinetic scale. Thus, Di appears in the thermal equation to account for the viscous heating term, which is kinetic in nature.

We have run a series of 3D numerical simulations solving Equations (13)–(15) using the Rayleigh code (Matsui et al. 2016; Featherstone & Hindman 2016a; Featherstone 2018) with a chosen nondimensional background entropy gradient $d\bar{S}/{dr}$ as described above and shown in Figure 2.

Zoom In Zoom Out Reset image size

We use impenetrable and stress-free boundary conditions for the velocity. For the entropy perturbations, we assume $\partial S/\partial r{| }_{{r}_{i}}=0$ at the inner boundary to account for the fixed flux of energy generated due to the nuclear burning at the stellar core, while we fix the entropy at the outer boundary such that $S{| }_{{r}_{o}}=0$. We note, however, that recent studies seem to suggest that a fixed flux outer boundary condition might also be a sensible choice in the stellar context (Matilsky et al. 2020). In all simulations, we assume a fixed Pr = 1 in order to avoid the onset of nonphysical modes known to arise in rapidly rotating anelastic convective systems for Pr < 1 (e.g., Calkins et al. 2015), while we vary N_ρ , Ra, and Ek (see Table 1). Each simulation is evolved from a zero initial velocity and small-amplitude perturbations in the entropy field until a statistically stationary and thermally relaxed state is reached.

Zoom In Zoom Out Reset image size

In Section 4, we focused on the fluctuating component of the kinetic energy in order to understand the dependence of the convective overshooting motions on Ro. We now concentrate on the differential rotational profiles and the associated dynamics within the different regimes arising at different Ro cases. We also investigate the mean flows resulting from velocity correlations due to the presence of rotation through the Coriolis force. In Figure 13, we present slices of the time- and azimuthally-averaged meridional circulation streamlines with the underlying contours of the mass flux (left panels) and differential rotation profiles ($\langle {u}_{\phi }\rangle /(r\sin \theta )$; right panels) for the four typical Ro cases. We see that for Ro ≈ 0.007, there are two large cells in the meridional circulation within the CZ and the overshoot region and four smaller ones close to the outer CZ boundary. All of these cells are localized at midlatitudes and the equator where the flow motions are the strongest (also see Figure 9(a)), similar to the differential rotation. We notice that the equator rotates faster than the poles, but convection in this case is inhibited by rotation, so it is in a rotationally constrained regime. At Ro ≈ 0.043, the profile exhibits a multicellular flow near the equator that propagates into the RZ particularly within the equatorial region. The differential rotation is again much stronger at the equator and decreases at midlatitudes and the poles in the CZ, while the stable region rotates almost uniformly. This rotational profile is qualitatively similar to the one observed in the solar CZ, although the latter is conical and not cylindrical as the one in this case. We will refer to this case as the "solar-like" case. At Ro ≈ 0.23, the meridional circulation has two distinct large cells within the CZ, while the differential rotation in the CZ is now characterized by faster rotating poles and a slower equator, and is more or less uniform in the bulk of the RZ. Hence, we will refer to this case as "anti-solar" (e.g., Gastine et al. 2014). Finally, in the largest Ro ≈ 1.7 case, the Coriolis force is very weak, and convection is not significantly affected by rotation. The meridional circulation consists of multiple cells and the convection zone now rotates uniformly at a slower rate than the RZ. Overall, in all of these cases (except for the rotationally constrained case of Ro ≈ 0.007), there is significant propagation of the meridional circulation into the stable layer, even beyond the overshoot region.

Zoom In Zoom Out Reset image size

To understand this behavior more quantitatively, in Figure 14, we plot the time- and spherically-averaged nondimensional kinetic energy related to the differential rotation ${\tilde{E}}_{\mathrm{DR}}$ and the meridional circulation ${\tilde{E}}_{\mathrm{MC}}$ as a function of radius where

$$\begin{eqnarray}&&{\tilde{E}}_{\mathrm{DR}}(r)=\displaystyle \frac{1}{2}\bar{\rho }\widetilde{{u}_{m,\phi }^{2}},\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{\tilde{E}}_{\mathrm{MC}}(r)=\displaystyle \frac{1}{2}\bar{\rho }(\widetilde{{u}_{m,r}^{2}}+\widetilde{{u}_{m,\theta }^{2}}).\end{eqnarray} \tag{ 29 }$$

Zoom In Zoom Out Reset image size

In Figure 14(a), we see that the kinetic energy associated with the differential rotation ${\tilde{E}}_{\mathrm{DR}}(r)$ is largest for the solar-like case and the anti-solar case both in the CZ and the RZ. It becomes very small at Ro ≈ 1.7, where the rotational influence is almost negligible, and in the lowest Ro ≈ 0.007 case where the motions are very weak overall and are unable to redistribute angular momentum as efficiently. Moreover, in all four cases, ${\tilde{E}}_{\mathrm{DR}}$ decreases with depth in the CZ, but there is still considerable kinetic energy below the base of the CZ (with the exception of Ro ≈ 0.007).

As shown in Figure 14(b), the kinetic energy associated with the meridional flows ${\tilde{E}}_{\mathrm{MC}}(r)$ decreases with decreasing Ro in the CZ and is very small for the lowest Ro ≈ 0.007 case (both in the CZ and in the RZ) where rotation dominates the dynamics. For the higher Ro cases, there is substantial kinetic energy in the meridional flows within the overshoot region and part of the stable zone. This indicates that these larger-scale mean flows, although mostly generated in the CZ, do not merely stop at the base of the convective region, but they travel deeper into the RZ.

We now focus on the solar-like case and the anti-solar case with Ro ≈ 0.043 and Ro ≈ 0.23, respectively. We look at the angular momentum transport and the main processes responsible for redistributing angular momentum within the two-zone spherical shell. Following Elliott et al. (2000) and Brun et al. (2004), by taking the ϕ component of the momentum equation and assuming a statistically stationary state, averaged in time and longitude, we obtain

$$\begin{eqnarray}&&\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial ({r}^{2}\langle {F}_{r}(r,\theta )\rangle )}{\partial r}+\displaystyle \frac{1}{r\sin \theta }\displaystyle \frac{\partial (\sin \theta \langle {F}_{\theta }(r,\theta )\rangle )}{\partial \theta }=0.\end{eqnarray} \tag{ 30 }$$

The radial flux F_r (r, θ) and the latitudinal flux F_θ (r, θ) can be expressed nondimensionally as

$$\begin{eqnarray}&&\begin{array}{l}{F}_{r}(r,\theta )=\bar{\rho }r\sin \theta \left({u}_{f,r}{u}_{f,\phi }+{u}_{m,r}{u}_{m,\phi }\right.\\ \ \ \left.+\,\displaystyle \frac{{u}_{m,r}r\sin \theta }{\mathrm{Ek}}+r\displaystyle \frac{\partial }{\partial r}\left(-\displaystyle \frac{{u}_{m,\phi }}{r}\right)\right),\end{array}\end{eqnarray} \tag{ 31 }$$

and

$$\begin{eqnarray}&&\begin{array}{l}{F}_{\theta }(r,\theta )=\bar{\rho }r\sin \theta \left({u}_{f,\theta }{u}_{f,\phi }+{u}_{m,\theta }{u}_{m,\phi }\right.\\ \ \ \left.+\,\displaystyle \frac{{u}_{m,\theta }r\sin \theta }{\mathrm{Ek}}+\displaystyle \frac{\sin \theta }{r}\displaystyle \frac{\partial }{\partial \theta }\left(-\displaystyle \frac{{u}_{m,\phi }}{\sin \theta }\right)\right),\end{array}\end{eqnarray} \tag{ 32 }$$

respectively. We can then consider the time and azimuthal average of the individual terms in Equation (31) such that

$$\begin{eqnarray}&&{F}_{r,R}=\langle \bar{\rho }r\sin \theta ({u}_{f,r}{u}_{f,\phi })\rangle ,\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray}&&{F}_{r,M}=\langle \bar{\rho }r\sin \theta ({u}_{m,r}{u}_{m,\phi })\rangle ,\end{eqnarray} \tag{ 34 }$$

$$\begin{eqnarray}&&{F}_{r,C}=\left\langle \bar{\rho }r\sin \theta \left(\displaystyle \frac{{u}_{m,r}r\sin \theta }{\mathrm{Ek}}\right)\right\rangle ,\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{F}_{r,V}=\left\langle \bar{\rho }r\sin \theta \left(-r\displaystyle \frac{\partial }{\partial r}\left(\displaystyle \frac{{u}_{m,\phi }}{r}\right)\right)\right\rangle .\end{eqnarray} \tag{ 36 }$$

In a similar manner, the individual terms in Equation (32) are

$$\begin{eqnarray}&&{F}_{\theta ,R}=\langle \bar{\rho }r\sin \theta ({u}_{f,\theta }{u}_{f,\phi })\rangle ,\end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray}&&{F}_{\theta ,M}=\langle \bar{\rho }r\sin \theta ({u}_{m,\theta }{u}_{m,\phi })\rangle ,\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&{F}_{\theta ,C}=\left\langle \bar{\rho }r\sin \theta \left(\displaystyle \frac{{u}_{m,\theta }r\sin \theta }{\mathrm{Ek}}\right)\right\rangle ,\end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray}&&{F}_{\theta ,V}=\left\langle \bar{\rho }r\sin \theta \left(-\displaystyle \frac{\sin \theta }{r}\displaystyle \frac{\partial }{\partial \theta }\left(\displaystyle \frac{{u}_{m,\phi }}{\sin \theta }\right)\right)\right\rangle .\end{eqnarray} \tag{ 40 }$$

Note that Equations (33) and (37) are associated with the Reynolds stresses, Equations (34) and (38) and Equations (35) and (39) correspond to the mean flows associated with the meridional circulation, and finally Equations (36) and (40) are related to the viscous stresses.

In Figure 15(a), (b), (c), and (d) we show the profiles of the latitudinal fluxes at Ro ≈ 0.043 (solar-like case). We observe that F_θ,C is large near the equator and at midlatitudes within the CZ. There is substantial penetration in the stable region where the distinct cells have opposite signs commensurate with the meridional circulation profile in Figure 13(c). Angular momentum in the convection zone and in the overshoot region seems to be dictated by the flux related to the action of the Coriolis force on the mean meridional flows, while all of the other fluxes are comparatively very small. In Figure 15(e), (f), (g), and (h), we notice that the radial flux F_r,C is again large within the CZ, while it also extends further down into the RZ transporting angular momentum mainly inward. The profile of F_r,R shows that transport by Reynolds stresses is weak compared to the transport by F_r,C ; however, it does contribute to the overall angular momentum transport with inward transport at midlatitudes and outward transport near the equator. The flux associated with the mean flows F_r,M is negligible, similar to F_θ,M , while F_r,V is also small overall although it does seem to become stronger close to the equator within the CZ where it transports angular momentum inward.

Zoom In Zoom Out Reset image size

For the anti-solar case of Ro ≈ 0.23, in Figure 16(c), we see that, similar to the solar-like case, F_θ,C is dominant within the CZ and the overshoot region, while F_θ,R , F_θ,M , and F_θ,V are all much weaker (Figure 16(a), (b), and (d)). In Figure 16(e), (f), (g), and (h), we notice that F_r,C is again the largest in magnitude compared with the other fluxes carrying angular momentum outward near the equator and inward at midlatitudes, while F_r,R , although not as large, is nonetheless present within the CZ carrying angular momentum inward (except at the poles where it is too weak), similar to the much weaker F_r,M . Finally, F_r,V is mostly larger close to the base of the CZ within the equatorial region transporting angular momentum in the same direction as F_r,C .

Zoom In Zoom Out Reset image size

We are interested in understanding what balances are achieved and how the mean flows are sustained in a steady state across the two-zone spherical shell. Two dominant terms that contribute to the meridional flows are the rotational shear associated with the differential rotation and the baroclinicity of the mean stratification. This balance, known as thermal-wind balance, occurs when the system is in both hydrostatic and geostrophic equilibrium, namely when there is a balance between the Coriolis force, the pressure, and the gravity (see, e.g., Ruediger 1989; Kitchatinov & Ruediger 1995; Rempel 2005; Miesch et al. 2006; Miesch & Hindman 2011; Acevedo-Arreguin et al. 2013). Within the anelastic approximation, the dimensional thermal-wind equation is then given by

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Omega }}}^{2}}{\partial z}\approx \displaystyle \frac{g}{\lambda {{rc}}_{p}}\displaystyle \frac{\partial \langle S\rangle }{\partial \theta },\end{eqnarray} \tag{ 41 }$$

where $z=r\cos \theta$, $\lambda =r\sin \theta$, and Ω = Ω_o + 〈u_ϕ 〉/λ is the total angular velocity. Nondimensionally, Equation (41) can be expressed as

$$\begin{eqnarray}&&\mathop{\underbrace{\displaystyle \frac{2}{\mathrm{Ek}}\left(r\cos \theta \displaystyle \frac{\partial \langle {u}_{\phi }\rangle }{\partial r}-\sin \theta \displaystyle \frac{\partial \langle {u}_{\phi }\rangle }{\partial \theta }\right)}}\limits_{\mathrm{LHS}}\approx \mathop{\underbrace{\displaystyle \frac{\mathrm{Ra}g(r)}{\Pr }\displaystyle \frac{\partial \langle S\rangle }{\partial \theta }}}\limits_{\mathrm{RHS}},\end{eqnarray} \tag{ 42 }$$

(where again 〈 · 〉 denotes a time- and azimuthally-averaged quantity). In Figure 17, we examine the thermal-wind balance for the solar-like case with Ro ≈ 0.043 and the anti-solar case with Ro ≈ 0.23. In Figure 17(a), we plot the left-hand-side (LHS) term and the right-hand-side (RHS) term of Equation (42) at different radii from the bottom of the CZ r_c down to ∼ r_c − δ_G for the solar-like case. The LHS term is associated with the Coriolis force and thus refers to the inertia of the differential rotation, while the RHS term accounts for the baroclinicity of the stratification through the latitudinal gradient of the entropy perturbations. The thermal-wind balance is satisfied within the overshoot region for this solar-like case where LHS ≈ RHS in that region. In Figure 17(b), we also show contours of the LHS and RHS terms with respect to θ and r. We observe that overall, the thermal-wind balance seems to be satisfied across the whole shell, and this is more pronounced within the equatorial region.

Similarly, in Figure 17(c), we plot the LHS and RHS terms of Equation (42) at all radii from the bottom of the CZ down to ∼ r_c − δ_G for the anti-solar case with Ro ≈ 0.23, and we see that the LHS term is somewhat larger than the RHS term within the overshoot region, especially at midlatitudes. Looking at Figure 17(d), we find that the LHS term is substantially different from the RHS term especially within the convection zone, thus indicating that the thermal-wind balance is not satisfied for the anti-solar case in the CZ, although it is somewhat sustained in the bulk of the stable region. We note that Equation (42) essentially shows that it is possible to maintain noncylindrical rotational profiles (i.e., when LHS ≠ 0) by accounting for latitudinal gradients of 〈S〉, but we observe no such behavior in this study.

Zoom In Zoom Out Reset image size

A well-known mechanism for the generation of meridional circulation via balances in the angular momentum equation, which had been originally extensively studied in the context of the Earth's atmospheric circulation and has been more recently revisited in studies of astrophysical fluid dynamics (McIntyre 2007; Garaud & Acevedo-Arreguin 2009; Garaud & Bodenheimer 2010; Wood et al. 2011; Miesch & Hindman 2011), is the so-called gyroscopic pumping. More specifically, gyroscopic pumping is associated with the advection of angular momentum by meridional flow due to the divergence of Reynolds stresses and viscous stresses. We now briefly explore the gyroscopic pumping balances in the solar-like case and the anti-solar case to see whether this mechanism takes place in both of these representative runs as well as which fluxes are dominant in maintaining it in the convective zone and the overshoot region.

Taking the zonal component ϕ of the momentum equation and multiplying it by $\lambda =r\sin \theta$, averaging over both time and longitude and finally assuming a steady state, we obtain the conservation of angular momentum equation of our system given by

$$\begin{eqnarray}&&\bar{\rho }\langle {{\boldsymbol{u}}}_{{mc}}\rangle \cdot {\rm{\nabla }}{ \mathcal L }=-({\rm{\nabla }}\cdot {{\boldsymbol{F}}}_{\mathrm{RS}}+{\rm{\nabla }}\cdot {{\boldsymbol{F}}}_{\mathrm{VS}})\end{eqnarray} \tag{ 43 }$$

where u _mc = (u_r , u_θ ) and where ${ \mathcal L }={\lambda }^{2}{\rm{\Omega }}$, which can be nondimensionally expressed as

$$\begin{eqnarray}&&{ \mathcal L }=r\sin \theta \left(\displaystyle \frac{r\sin \theta }{\mathrm{Ek}}+\langle {u}_{\phi }\rangle \right).\end{eqnarray} \tag{ 44 }$$

The nondimensional transport of angular momentum due to Reynolds stresses and due to viscous stresses are given respectively by

$$\begin{eqnarray}&&{{\boldsymbol{F}}}_{\mathrm{RS}}=\bar{\rho }r\sin \theta (\langle {u}_{f,r}{u}_{f,\phi }\rangle ,\langle {u}_{f,\theta }{u}_{f,\phi }\rangle ),\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}&&\begin{array}{l}{{\boldsymbol{F}}}_{\mathrm{VS}}=\bar{\rho }\left(\langle {u}_{m,\phi }\rangle \sin \theta -r\sin \theta \displaystyle \frac{\partial \langle {u}_{m,\phi }\rangle }{\partial r},\right.\\ \ \ \left.\,\langle {u}_{m,\phi }\rangle \cos \theta -\sin \theta \displaystyle \frac{\partial \langle {u}_{m,\phi }\rangle }{\partial \theta }\right).\end{array}\end{eqnarray} \tag{ 46 }$$

Following Featherstone & Miesch (2015), in Figure 18(a), we plot the left-hand-side term of Equation (43), i.e., $\bar{\rho }\langle {{\boldsymbol{u}}}_{{mc}}\rangle \cdot {\rm{\nabla }}{ \mathcal L }$, along with the right-hand-side term −(∇ · F _RS + ∇ · F _VS) at all radii from r_c down to ∼ r_c − δ_G for the solar-like case with Ro ≈ 0.043. We find that the balance is more or less achieved within the overshoot region with $\bar{\rho }\langle {{\boldsymbol{u}}}_{{mc}}\rangle \cdot {\rm{\nabla }}{ \mathcal L }$ being slightly larger than −(∇ · F _RS + ∇ · F _VS). In fact, by looking across the whole shell in Figure 18(b), we find that, overall, Equation (43) is satisfied within the CZ and the overshoot region where the differential rotation and the associated meridional flows are generated and advected. We note that we should not expect Equation (43) to be satisfied in the bulk of the RZ where the motions are not convectively driven and so are largely axisymmetric. Furthermore, we notice that −∇ · F _RS seems to be dominant in the bulk of the CZ while −∇ · F _VS is larger at lower latitudes and the equatorial region, especially close to and somewhat below the bottom of the CZ.

Zoom In Zoom Out Reset image size

In Figure 18(c), we show the left-hand-side and the right-hand-side terms of Equation (43) at different radii within the overshoot region for the anti-solar case with Ro ≈ 0.23. We observe that the gyroscopic pumping balance seems to be satisfied within that region although there are some discrepancies between the two terms mostly close to the equator with −(∇ · F _RS + ∇ · F _VS) being somewhat larger than $\bar{\rho }\langle {{\boldsymbol{u}}}_{{mc}}\rangle \cdot {\rm{\nabla }}{ \mathcal L }$. In Figure 18(d), we plot the individual terms across the whole shell and find that the gyroscopic balance is achieved. Similar to the solar-like case, the divergence of the flux associated with the Reynolds stresses is dominant throughout the shell while the viscous stresses play a less important role overall and are larger mainly at midlatitudes and the equatorial region of the CZ.

Stars operate in a regime where viscosity is very small (with typical Prandtl numbers Pr ∼ O(10⁻⁶)), so we would expect that the advection of the meridional flows would be due to only the divergence of the fluxes associated with the Reynolds stresses while viscous stresses should be negligible, unlike what we see in the runs presented here. Indeed, viscosity seems to have a fairly strong presence near the equatorial region for both the solar-like case and the anti-solar case, which would mean that we may need to consider even higher Ra runs (which are computationally much harder to achieve) so that the viscous stresses would not need to compensate for the weaker Reynolds stresses observed at this Ro and Ra.

We have presented the results of a series of numerical simulations in a two-zone spherical shell in order to examine the effect of rotation and density stratification on the overshooting dynamics associated with the interaction of a convectively unstable region overlying a stably stratified zone. Motivated by solar-type stars, we have assumed a fixed thermal gradient boundary condition at the inner boundary, as well as the existence of an internal heat source whose function mimics the radiative heating presented in more complicated solar models (Model S), while we have used a fixed aspect ratio that corresponds to the inner part of the solar convective region and a large portion of the RZ. We have conducted a systematic study by varying the density stratification in the CZ, the Rayleigh number, and the Ekman number, while we have assumed a fixed Prandtl number (Pr = 1). To create our two-layered configuration, we have considered a varying background entropy gradient $d\bar{S}/{dr}$, which leads to an adiabatic CZ that smoothly transitions into a subadiabatic RZ with a profile fixed for all of our runs—both nonrotating and rotating.

There have been many numerical studies of overshooting convection considering a convective region overlying a stably stratified zone, which, as we discussed in Section 1, have been mainly focused on the dependence of the overshooting on the relative stability and the transition width between the two layers as well as on the degree of the convective driving in the CZ. In all of these studies, similar to ours, they found that there is a substantial amount of overshooting δ below the base of the convective region that depends on the input parameters and the different set-ups and geometries, or whether they are 2D or 3D (see, e.g., Hurlburt et al. 1986, 1994; Singh et al. 1998; Brummell et al. 2002; Rogers & Glatzmaier 2005; Hotta 2017; Pratt et al. 2017; Käpylä 2019; Korre et al. 2019). For instance, Hurlburt et al. (1986) reported a value of δ = 0.67H_p , Singh et al. (1998) found that δ ranges from 0.4H_p to 0.88H_p , while Brummell et al. (2002) reported values of δ between 0.4H_p and 2H_p, where H_p is the pressure scale-height. In our runs the pressure scale-height at the bottom of the CZ defined as ${H}_{p}=-{((1/\bar{P}(r))d\bar{P}(r)/{dr})}^{-1}{| }_{r={r}_{c}}$ depends on N_ρ and ranges between ∼0.32 and ∼46 such that the overshoot lengthscale in terms of pressure scale-heights is 0.0074H_p ≤ δ = δ_G r_o ≤0.64H_p . So our findings are similar to those reported in previous studies, especially when we consider higher N_ρ cases.

We also found that the computed overshoot lengthscale δ_G for the nonrotating cases, which have fixed Pr and Ra, depends on the ratio of the density stratifications in the two zones given by the parameter R_ρ such that ${\delta }_{G}\propto {R}_{\rho }^{0.36}$. However, to the authors' knowledge, there is no previous systematic study on the dependence of overshooting on the density stratification as an input parameter. A notable exception is the work of Massaguer et al. (1984), who employed the anelastic approximation in a Cartesian geometry and reported the penetration lengthscale for two different density stratification cases and also compared these results to a Boussinesq case. In their work, they assumed stacked polytropes with a varying conductivity profile such that a convection zone is sandwiched between two stable layers and studied upward and downward penetration. However, they did not perform fully nonlinear 3D calculations, but instead they used anelastic modal equations whereby they expanded the horizontal structure of the flow using one- and two-mode hexagonal planforms in order to better resolve the vertical structure. As a result, their findings have a strong dependence on the chosen horizontal modes and cannot capture the fully nonlinear dynamics associated with overshooting convection. Nonetheless, they found that the increased density stratification leads to larger penetration and that the anelastic case predicts larger penetration depths than the Boussinesq case. In fact, they reported strong pure penetration whereby the convection zone extended further down into the stable region and part of the latter became adiabatic.

In our work, we do not observe pure penetration but rather overshooting, although there is a trend in the $d{\tilde{S}}_{T}/{dr}$ profiles (see Equation (21)) indicating that increasing the stratification along with the Rayleigh number even further might actually result in pure penetration. This can be seen in Figure 19, where we plot the total adjusted $d{\tilde{S}}_{T}/{dr}$ profile along with $d\bar{S}/{dr}$ for the case with N_ρ = 3 at three different Ra. We also notice that the slightly subadiabatic layer around the base of the CZ (discussed in Section 3) becomes even smaller at larger Ra, where convection is more vigorous.

Zoom In Zoom Out Reset image size

Massaguer et al. (1984) also showed that the local Pe decreased with decreasing stratification, which is not the case in our fully nonlinear 3D calculations where increased density stratification leads to smaller Reynolds (and Péclet numbers) due to somewhat weaker convection at the same Ra for larger N_ρ (see, e.g., Jones et al. 2009).

We postulate that, in a set-up where N_ρ in the CZ increases but somehow the density ratio in the RZ remains the same for all of the cases, R_ρ will be a monotonic function of N_ρ and as such it can lead to larger overshoot depths for larger values of N_ρ , but this remains to be verified by simulations that account for such a two-layered configuration.

There have been quite a few previous studies that have accounted for rotation in overshooting convective experiments. Using 3D compressible simulations, Brummell et al. (2002) showed that the overshoot depth decreases with decreasing Ro as a result of the braking of the vertical flows due to the horizontal interactions of the like-sign vortices as well as the fact that the downflows penetrate at an angle in the presence of strong rotation. They reported a scaling of δ ∝ Ro^0.15, while they also showed that it depends on the latitude with the poles and the equator having larger δ. However, their scaling was a result of only three different cases, and they also used the f-plane approximation in a Cartesian box that intrinsically cannot account for global flows. Similarly, Pal et al. (2007) found that δ decreases with Ro as δ ∝ Ro^0.2 at the poles and δ ∝ Ro^0.4 at midlatitudes. Augustson & Mathis (2019) studied overshooting convection under the effect of rotation by building on the linearized model of Zahn (1991), and they reported that the overshoot depth decreases as Ro^0.3.

In this work, we found that δ_G ∝ Ro^0.23, which is not that different from the predictions of these previous models, even though we operate in a spherical geometry, using a different set-up based on the anelastic approximation. Thus, the decrease in the amount of overshooting with increasing rotation rate seems to be a salient feature of the rotating overshooting dynamics. Additionally, we showed that the overshoot depth depends on latitude with the solar-like case (Ro ≈ 0.043) leading to larger overshoot depths closer to the poles and at midlatitudes where the flow is more weakly influenced by rotation and smaller depths near the equator where the convective motions penetrate at an angle due to their alignment with the axis of rotation. For the anti-solar case (Ro ≈ 0.23), again the overshoot depth is larger at higher latitudes, while it is just a bit larger close to the equator where the motions are not as much aligned with the axis of rotation as at midlatitudes. Possibly this is also a result of the topology of the plumes there, which appear to be more like the fly-wheeling motions described in Brummell et al. (2002). Overall when increasing Ro the latitudinal dependence becomes weaker due to the smaller effect of the rotation on convection.

We note that other studies of anelastic overshooting convection in spherical shells (see Browning et al. 2004; Brun et al. 2011; Augustson et al. 2012; Brun et al. 2017) showed that the overshoot lengthscale decreases with increasing Rossby number. However, these simulations were not focused on the dependence of the overshooting dynamics on Ro; the measure of the overshoot depth was based on the radial enthalpy flux rather than the full kinetic energy, and they were operating at lower Pe. Thus, it is not clear how to draw meaningful comparisons with the results of our study.

We showed that for sufficiently small Ro, we can obtain a solar-like differential rotation with a faster equator and slower poles in the CZ transitioning to a uniform rotation in the RZ. However, helioseismology has revealed that the Sun also possesses a conical profile, while our solutions lead to cylindrical ones.

Previous studies (e.g., Rempel 2005; Miesch et al. 2006) have shown that enforcing a strong latitudinal entropy gradient at the base of the CZ can result in a conical differential rotation similar to what we observe in the Sun. Looking at the entropy profile in Figure 3(b) of Miesch et al. (2006; also see Figure 8(d) from Matilsky et al. 2020 who accounted for a fixed flux outer boundary condition instead of a latitudinal entropy gradient at the inner boundary), we notice that, in the solar-like conical profiles achieved in their calculations, the entropy perturbations increase monotonically with latitude away from the equator.

In Figure 20, we show the time- and azimuthally-averaged entropy perturbations (where we have subtracted the spherically symmetric mean to enhance the latitudinal variations), volume averaged within the overshoot region 〈S〉_ov versus the latitude θ for the solar-like case of Ro ≈ 0.043. We observe that 〈S〉_ov is much weaker in magnitude and does not vary substantially in latitude compared with the entropy profiles corresponding to the conical differential rotation contours in Miesch et al. (2006). Also, despite the fact that the equatorial region is still cooler than the midlatitudes and the poles, the profile is nonmonotonic. These previous studies only considered a convection zone and thus they had to enforce a latitudinal entropy gradient at the inner boundary (or assume a fixed flux outer boundary condition) to achieve the tilt of the differential rotation contours. However, in our two-zone spherical shell model set-up, a strongly varying latitudinal entropy gradient might be obtained self-consistently with a suitable background $d\bar{S}/{dr}$ profile. Indeed, this could result in stronger entropy gradient perturbations ∂〈S〉/∂θ in the thermally equilibrated state, which could in turn lead to noncylindrical differential rotation profiles. This topic will be addressed in future work.

Zoom In Zoom Out Reset image size

We may attempt to provide some estimates on the solar value of the overshooting lengthscale using the results of our study. We note, however, that such extrapolations tend to be speculative to some degree, as values of many parameters used here are very different from the solar ones. For instance, Ra_⊙ ∼ O(10²⁰) and Pr_⊙ ∼ O(10⁻⁶); (e.g., Ossendrijver 2003), while our models operate in a much lower Ra and much higher Pr regime.

These caveats aside, we can use the scaling relationship between δ_G and the ratio of the densities R_ρ to estimate the solar overshoot depth. We consider the region of the convection zone spanning from the tachocline to the base of the granulation boundary layer (r_o ≈ 0.9983R_⊙), a region characterized by 11 e-foldings in density (e.g., Model S; Christensen-Dalsgaard et al.1996). Using this, and the values of the density at the base of the RZ at r_i ≈ 0.2R_⊙ and at the base of the CZ at r_c ≈ 0.71R_⊙, we find that R_ρ ≈ 0.003. Based on our scaling of ${\delta }_{G}=0.078{R}_{\rho }^{0.36}$ from Section 3, we estimate an overshoot lengthscale for the Sun of δ_G ≈ 0.0097, or, equivalently, δ = δ_G r_o ≈ 0.1H_p , where H_p ≈ 5.7 × 10⁹ cm is the pressure scale-height at the base of the solar CZ extracted from Model S. This estimated value for the solar overshoot lengthscale is similar to the estimates from other studies. For instance, Brummell et al. (2002) used 3D compressible simulations to estimate an overshoot depth in the range 0.02H_p –0.11H_p , while stellar evolution models typically assume overshoot depth values on the order of 0.1H_p . Helioseismic studies suggest that the upper limit for the solar overshoot depth is about 0.07H_p (see, e.g., Monteiro et al. 1994).

Carrying out a similar exercise based on Rossby number is more difficult, as there is still no clear consensus on the value of deep solar convection's Rossby number. Estimates of convective overturning time based on the B − V color index suggest that Ro_⊙ ∼ 0.04 (Corsaro et al. 2021), whereas Featherstone & Hindman (2016b) suggested that a value ≤0.01 was likely based on the observed spatial scale of supergranulation. In principle, this parameter could be inferred for the Sun through helioseismic measurements of deep convection, but attempts to do so have led to inconsistent results (e.g., Hanasoge et al. 2012; Greer et al. 2015; Nagashima et al. 2020). Given these present uncertainties, we refrain from using the scaling law of δ_G with Ro to obtain estimates of the solar overshoot lengthscale. Alternative observational approaches, however, such as those based on inertial mode observations (Gizon et al. 2021), may shed further light on this issue in the future.