A model of rotating convection in stellar and planetary interiors: I - convective penetration

TL;DR

This paper develops a model for stellar and planetary convection that predicts how convection properties and penetration depth depend on rotation rate, diffusivities, and stellar structure, providing insights into convective boundary dynamics.

Contribution

It introduces a monomodal convection model linking velocity, superadiabaticity, and length scale to rotation and diffusivities, and relates convective penetration depth to local flow parameters.

Findings

Convection velocity and superadiabaticity scale with rotation rate and diffusivities.

Penetration depth depends on the Rossby number, diffusivities, and pressure scale height.

Upward and downward penetration exhibit similar but distinct scaling behaviors.

Abstract

A monomodal model for stellar and planetary convection is derived for the magnitude of the rms velocity, degree of superadiabaticity, and characteristic length scale as a function of rotation rate as well as with thermal and viscous diffusivities. The convection model is used as a boundary condition for a linearization of the equations of motion in the transition region between convectively unstable and stably-stratified regions, yielding the depth to which convection penetrates into the stable region and establishing a relationship between that depth and the local convective Rossby number, diffusivity, and pressure scale height of those flows. Upward and downward penetrative convection have a similar scaling with rotation rate and diffusivities, but they depend differently upon the pressure scale height due to the differing energetic processes occurring in convective cores of early-type stars versus convective envelopes of late-type stars.

Tables

Table 1. Table 1: Frequently used symbols in the convection model.

$a=k_x/k_z$	Maximizing horizontal wavevector
$k_z=\pi /{\ell }_0$	Maximizing vertical wavevector
$K_0=\kappa k_z^2/N_{\ast ,0}$	Normalized thermal diffusivity
$O_0=\sqrt{6}/\left(5\pi {\mathrm{Ro}}_c\right)$	Normalized Coriolis coefficient
${\mathrm{Ro}}_c=\sqrt{6}{N}_{\ast ,0}/\left(10\pi {\Omega }_0\right)$	Convective Rossby number
$\hat{s}=s/N_{\ast }$	Normalized growth rate
$v_0=\sqrt{6}{N}_{\ast ,0}{\ell }_0/\left(5\pi \right)$	Velocity of the nonrotating case
$V_0=\nu k_z^2/N_{\ast ,0}$	Normalized viscosity
$q=N_{\ast ,0}/N_{\ast }$	Ratio of buoyancy timescales
$z^3=k^2/k_z^2$	Normalized wavevector

Equations

$$\left({\partial }_t-\nu {\nabla }^2\right)\mathrm{v}$$

$$\left({\partial }_t-\kappa {\nabla }^2\right)T^{\prime }$$

$$\left({\partial }_t-\kappa {\nabla }^2\right)T^{\prime }=\beta {\mathrm{v}}_z,$$

$$\left({\partial }_t-\nu {\nabla }^2\right)\zeta =2\mathbf{\Omega }\cdot \nabla {\mathrm{v}}_z,$$

$$\left({\partial }_t-\nu {\nabla }^2\right){\nabla }^2{\mathrm{v}}_z=g{\alpha }_T{\nabla }_{\perp }^2{T}^{\prime }-2\mathbf{\Omega }\cdot \nabla \zeta ,$$

$$\left({\partial }_t-\kappa {\nabla }^2\right){\left({\partial }_t-\nu {\nabla }^2\right)}^2{\nabla }^2{\mathrm{v}}_z+g{\alpha }_T\beta {\nabla }_{\perp }^2\left({\partial }_t-\nu {\nabla }^2\right){\mathrm{v}}_z$$

$$+4\mathbf{\Omega }\cdot \nabla \left[\mathbf{\Omega }\cdot \nabla \left({\partial }_t-\kappa {\nabla }^2\right){\mathrm{v}}_z\right]=0.$$

$$[\left(s+\kappa k^2\right){\left(s+\nu k^2\right)}^2{k}^2+g{\alpha }_T\beta k_{\perp }^2\left(s+\nu k^2\right)$$

$$+4{\left(\mathbf{\Omega }\cdot \mathbf{k}\right)}^2\left(s+\kappa k^2\right)]{\mathrm{v}}_z=0.$$

$$\hat{s}$$

$$K$$

$$\left(\hat{s}+K{z}^3\right)\left(z^3{\left(\hat{s}+V{z}^3\right)}^2+O^2{\left(\cos \theta +a_y\sin \theta \right)}^2\right)$$

$$-\left(z^3-1\right)\left(\hat{s}+V{z}^3\right)=0,$$

$$O^2$$

$$s\mathrm{v}=-i\mathbf{k}{P}^{\prime }/⟨\rho ⟩+g{\alpha }_T{T}^{\prime }\hat{\mathbf{z}}-2\mathbf{\Omega }\times \mathrm{v}-\nu k^2\mathrm{v},$$

$$i\mathbf{k}\cdot \mathrm{v}=0.$$

$$\left(s+\nu k^2\right)⟨{\mathrm{v}}^2⟩=g{\alpha }_T⟨T^{\prime }{\mathrm{v}}_z⟩.$$

$$F=\frac{\rho c_P\Re [\left(s+\nu k^2\right)⟨{\mathrm{v}}^2⟩]}{g{\alpha }_T}.$$

$$F$$

$$=\frac{⟨\rho ⟩c_P{N}_{\ast }^3}{g{\alpha }_T{k}_z^2}\left[\frac{\Re [\hat{s}]}{z^3}+V\right]{\left|\hat{s}\right|}^2,$$

$${\hat{s}}^2=\frac{a^2-O^2{\left(\cos \theta +\sin \theta a_y\right)}^2}{1+a^2},$$

$$\left[\hat{s}+V{z}^3\right]\left[z^3\left(\hat{s}+K{z}^3\right)\left(\hat{s}+V{z}^3\right)-a^2\right]=0.$$

$$\hat{s}$$

$$-\frac{1}{2}\left(V+K\right)z^3.$$

$$\hat{s}$$

$${\mathrm{v}}_0=\frac{s_0}{k_0}=\frac{\sqrt{6}}{5}\frac{N_{\ast ,0}}{k_z}=\frac{\sqrt{6}}{5\pi }{N}_{\ast ,0}{\ell }_0.$$

$$\mathrm{R}{\mathrm{o}}_{\mathrm{c}}=\frac{{\mathrm{v}}_0}{2{\Omega }_0{\ell }_0}=\frac{\sqrt{6}{N}_{\ast ,0}}{10\pi {\Omega }_0},$$

$$O$$

$$q$$

$$O$$

$$K$$

$$V$$

$$\left(\hat{s}+K_{0}q{z}^3\right)\left(z^3{\left(\hat{s}+V_{0}q{z}^3\right)}^2+O_0^2{q}^2{\cos }^2\theta \right)$$

$$-\left(z^3-1\right)\left(\hat{s}+V_{0}q{z}^3\right)=0,$$

$$F=\frac{F_0}{q^3}\left[\frac{{\hat{s}}^3}{z^3}+V_{0}q{\hat{s}}^2\right],$$

$$\frac{d\hat{s}}{dz}=\frac{3{\hat{s}}^2}{z\left(3\hat{s}+2V_{0}q{z}^3\right)}.$$

$$\frac{d\hat{s}}{dz}=-3z^2[V_{0}q\left(1-2z^3\right)+K_0{q}^3\left(O_0^2{\cos }^2\theta +4V_0^2{z}^9\right)$$

Full text

Convection, Instabilities, Turbulence – Stars: evolution, rotation

1 Introduction

In the context of stellar and planetary physics, convection driven through buoyancy and doubly-diffusive instabilities plays a critical role in mixing and transport processes in convectively unstable regions of stars (e.g., Miesch & Toomre, 2009; Kupka & Muthsam, 2017; Garaud, 2018). It primarily serves to transport the energy released deep within the star or planet through regions where radiative energy transport is inefficient. Moreover, convection directly and nonlocally transports heat and chemicals through advection and it acts diffusively through entrainment and dissipative processes. In the presence of rotation, convection can also transport angular momentum through Reynolds stresses and meridional flows to establish and maintain a differential rotation (e.g., Glatzmaier & Gilman, 1982; Kichatinov, 1986; Kichatinov & Rudiger, 1993; Brummell et al., 1996; Busse, 2002; Brun & Toomre, 2002; Miesch & Toomre, 2009; Augustson et al., 2012; Brun et al., 2017). However, parameterizing the impact of rotation on the convection and its related transport properties over evolutionary time scales remains a relatively open problem.

Convective flows cause mixing not only in regions of superadiabatic temperature gradients but in neighboring subadiabatic regions as well, as motions from the convective region contain sufficient inertia to extend into those regions before being buoyantly braked or turbulently eroded (e.g., Massager, 1990; Hurlburt et al., 1994; Miesch, 2005; Lecoanet et al., 2015; Viallet et al., 2015). This convective penetration and turbulence can thus alter the chemical composition and thermodynamic properties in those regions, softening the transition between convectively stable and unstable regions, with the consequence being that the differential rotation, opacity, and thermodynamic gradients are modified (e.g., Spiegel, 1963; Zahn, 1991; Canuto, 1992; Freytag et al., 1996; Brummell et al., 2002; Browning et al., 2004; Augustson et al., 2012; Brun et al., 2017; Pratt et al., 2017a). Such processes have an asteroseismic signature as has been observed in many kinds of stars (e.g., Aerts et al., 2003; Degroote et al., 2010; Neiner et al., 2012; Briquet et al., 2012; Zhang et al., 2013; Moravveji et al., 2016; Pedersen et al., 2018). Indeed, penetrating convection can modify the depth of the convection zone and any extant tachocline leading to a frequency shift in asteroseismically detected modes (e.g., Dziembowski & Pamyatnykh, 1991; Monteiro et al., 2000; Christensen-Dalsgaard et al., 2011; Montalbán et al., 2013) and advect magnetic field and angular momentum into the stable region leading to changes in dynamo action and the differential rotation (e.g., Miesch, 2005; Browning et al., 2006, 2007; Jones et al., 2010; Featherstone et al., 2009; Masada et al., 2013; Augustson, 2013; Augustson et al., 2016). In stars with a convective core, convective overshoot and penetration can lead to a greater amount of time spent on stable burning phases as fresh fuel is mixed into the burning region (e.g., Mowlavi & Forestini, 1994; Meakin & Arnett, 2007; Arnett et al., 2009; Maeder, 2009; Viallet et al., 2013; Jin et al., 2015). Penetrative convection induces chemical mixing near the transition between convectively stable and unstable regions of stars, leading to changes in the spectral characteristics of the atmosphere (e.g., Schatzman, 1977; Vauclair et al., 1978; Freytag et al., 2010; Baraffe et al., 2017). Hence, from the standpoint of stellar evolution, an open problem is to understand how the depth of penetration and the character of the convection in this region change with rotation, magnetism, and diffusion.

In this work, a generalized version of a heuristic model for convection for rotating systems is developed following Stevenson (1979). This model of convection is then employed to estimate the convective depth of penetration above and below a convection zone. The linearized Boussinesq equations of motion that yield the growth rate used in the model are given in §2.1. The heat-flux maximization principle is discussed in §2.2. The model of convection is derived for a diffusion-free setting in §2.4, when including thermal diffusion in §2.5, and when taking into account both thermal and viscous diffusion in §2.6. This convection model is then used in conjunction with the linearized model of convective penetration derived in Zahn (1991) to estimate the depth of convective penetration as shown in §3. A summary of the results and perspectives are presented in §4.

2 Heat-Flux Maximized Convection Model

As a first step, the heuristic model will be considered to be local such that the length scales of the flow are much smaller than either density or pressure scale heights. This is equivalent to ignoring the global dynamics and assuming that the convection can be approximated as local at each radius and colatitude in a star or planet. As such, one may consider the dynamics to be in the Boussinesq limit (Section 2.1). The linearized dynamics is used to construct the heat-flux maximization principle.

2.1 Linearized Boussinesq Dynamics

Boussinesq devised one of the simplest convective systems, consisting of an infinite layer of a nearly incompressible fluid with a small thermal expansion coefficient $\alpha_{T}=-\partial\ln{\rho}/\partial T|_{P}$ that is confined between two infinite impenetrable plates differing in temperature by $\Delta T=T(z_{2})-T(z_{c})$, where it is assumed for this model that $T(z_{2})<T(z_{c})$, and separated by a distance $\ell_{0}=z_{2}-z_{c}$, as in Figure 1. The thermodynamic variables are further expanded about their averages as $q=\langle q\rangle+\bar{q}(z)+q^{\prime}$, with $\langle q\rangle$ being the volumetric mean, $\bar{q}$ being the horizontal average, and $q^{\prime}$ being the dynamical perturbation. Following Spiegel & Veronis (1960) and Chandrasekhar (1961), the dynamics of a rotating Boussinesq system in Cartesian coordinates may be ascertained from a parametric expansion of the Navier-Stokes equation under the assumption that quantities are scaled relative to and expanded in the parameter $\alpha_{T}\Delta T$. This is equivalent to requiring that the maximum density fluctuation be small in the sense that $\alpha_{T}\Delta T\ll 1$. Under these assumptions, the equations of motion become

[TABLE]

where $t$ is the time coordinate and $\partial_{t}$ is the partial derivative with respect to it, $\bm{\nabla}$ is the spatial coordinate gradient. To first order in $\alpha_{T}\Delta T$, the continuity equation requires that the velocity field $\bm{\mathrm{v}}$ be solenoidal ($\bm{\nabla}\bm{\cdot}{\bm{\mathrm{v}}}=0$). The vertical direction $\bm{\hat{z}}$ is anti-aligned to the local direction of the constant effective gravity $g$, e.g. $\bm{\hat{z}}=-\mathbf{g}/g$, so that fluid with $T^{\prime}>0$ rises when $\Delta T<0$. As the Cartesian f-plane domain assumed here is meant to correspond to a small, local region in the global spherical geometry, the $x$ direction corresponds to the azimuthal direction and the $y$ direction to the latitudinal direction (See Figure 1). Here $\langle\rho\rangle$ is the mean density, $P^{\prime}$ is the pressure perturbation, $T^{\prime}$ is the temperature perturbation, $\mathrm{v}_{z}$ is the vertical component of the velocity field, $\bm{\Omega}$ is the local rotational vector that is taken to be in the y-z plane, and $\kappa$ and $\nu$ are the thermal and viscous diffusivities. The vertically-aligned thermal gradient to first order in the expansion parameter $\alpha_{T}\Delta T$ is defined as $\beta=d\bar{T}/dz+g/c_{P}$, where $c_{P}$ is the specific heat capacity at constant pressure. The adiabatic lapse rate $g/c_{p}$ arises from the first order expansion of the $PdV$ work in the thermal energy equation.

To construct a dispersion relationship for later use in the heat-flux maximization turbulence model, one may take the $z$ component of both the curl and the double curl of Equation 1. Dropping the nonlinear terms, these operations yield the following set of linearized Boussinesq equations

[TABLE]

with $\zeta$ being the vertical component of the vorticity and $\nabla_{\perp}$ being the gradient transverse to the vertical direction as in Equations 79, 84, and 85 of Chandrasekhar (1961). As seen in many papers regarding Boussinesq dynamics (e.g., Chandrasekhar, 1961), this set of equations can be reduced to a single third-order in time and eighth-order in space equation for the vertical velocity as

[TABLE]

For impenetrable and stress-free boundary conditions the solutions of Equation 6 are periodic in the horizontal, sinusoidal in the vertical, and exponential in time, e.g. $\mathrm{v}_{z}=\mathrm{v}\sin{\left[k_{z}\left(z-z_{c}\right)\right]}\exp{(i\bm{k}_{\perp}\bm{\cdot}\bm{r}+st)}$, where $\bm{k}_{\perp}$ is the horizontal wavevector, $s$ is the growth rate, $\bm{r}$ is the local coordinate vector, and $\mathrm{v}$ is a to be determined constant velocity amplitude. To satisfy the impenetrable, stress-free, and fixed temperature boundary conditions, it is required that the vertical wavenumber be $k_{z}=n\pi/\ell_{0}$. This can also be considered to be equivalent to modeling only the portion of the spectrum of convective eddies or motions whose Mach numbers are small and those that are much smaller than the integral scale of the system. The introduction of this solution into Equation 6 yields the following dispersion relationship that relates $s$ to the wavevector $\bm{k}$ as

[TABLE]

This equation can be simplified and made non-dimensional by dividing through by the appropriate powers of $N_{*}$ and $k_{z}$, where $N_{*}^{2}=|g\alpha_{T}\beta|$ is the absolute value of the square of the Brunt-Väisälä frequency as in Barker et al. (2014), which is otherwise negative in a convective region. Utilizing the following auxiliary quantities

[TABLE]

the dispersion relationship becomes

[TABLE]

where $4(\bm{\Omega}\cdot\bm{k})^{2}/N_{*}^{2}=k_{z}^{2}O^{2}\left(\cos{\theta}+a_{y}\sin{\theta}\right)^{2}$ with

[TABLE]

where $\Omega_{0}$ is the bulk rotation rate of the system.

In assessing the behavior of the Boussinesq heat flux, both with and without diffusion, there are several limits of the diffusive parameters that can be considered, each with its own physical justification. In highly turbulent regimes, it is often assumed that both diffusivities can be neglected, which is the case considered in S79, and expanded upon below in §2.4. Within stellar interiors, the molecular value of the thermal Prandtl number ($\mathrm{Pr}=\nu/\kappa$) is typically very small, being of the order of $10^{-5}$ or smaller, so the limit $\nu\rightarrow 0$ may be considered as in §2.5. The case where thermal diffusion can be ignored ($\kappa\rightarrow 0$) can be relevant to very high Prandtl number systems such as those found in geophysical contexts, but it is not directly considered here. However, since the general case with both diffusivities can be treated here (§2.6), turbulent diffusivities may also be incorporated into the model, where typically it is assumed that $\mathrm{Pr}$ is approximately unity. Moreover, to make contact with numerical simulations and laboratory experiments, it is useful to retain diffusive effects.

2.2 Maximum Heat Transport

The secular impacts of rotation and magnetic fields on stellar and planetary evolution are of keen interest within the astrophysical community (e.g., Maeder, 2009; Mathis, 2013). Accordingly, extensions to MLT and its ilk have been introduced (e.g., Stevenson, 1979; Canuto & Hartke, 1986; Canuto & Mazzitelli, 1991; Zhou, 1995; Xiong et al., 1997). As expounded upon in Stevenson (1979, hereafter S79), a surprisingly effective approach to including rotation in MLT is to hypothesize a convection model where the convective length-scale, degree of superadiabaticity, and velocity are governed by the linear mode that maximizes the convective heat flux. This heuristic method was made more rigorous when incorporated into the spectrally nonlocal large-scale turbulence model of Canuto & Hartke (1986), where it is used to set the scale of energy injection of a Heisenberg-Kolmogorov spectrum of turbulent energy cascade (Kolmogorov, 1941; Heisenberg, 1948). Subsequently, it sets the local values of the turbulent diffusivities. Though, these theories have not yet widely been employed in stellar evolution computations, with one exception being Ireland & Browning (2018). One can also incorporate approximate 2D dynamics through models of the Reynolds and Maxwell-stresses, which has been significantly developed over the last several decades (e.g., Durney & Spruit, 1979; Hathaway, 1984; Ruediger, 1989; Garaud et al., 2010).

The S79 model of rotating convection has its origins in the principle of maximum heat transport proposed by Malkus (1954). In that principle, an upper limit for a boundary condition dependent turbulent heat flux is established that depends upon the smallest Rayleigh unstable convective eddy. The size of this eddy is determined with a variational technique that is similar to that developed in Chandrasekhar (1961) for the determination of the Rayleigh number, which is the ratio of the buoyancy force to the viscous force multiplied by the ratio of the thermal to viscous diffusion timescales. This technique then permits the independent computation of the rms values of the fluctuating temperature and velocity amplitudes. Yet the heat-flux maximization principle is not fully rigorous as there is no apriori theoretical justification as to why the turbulence would arrange itself so as to maximize the heat flux (e.g., Howard, 1963). It is instead built upon three assumptions: first the mean-temperature gradient must be everywhere negative, second there is a finite range of wavenumbers that are effective at transporting heat, and third that the highest vertical wavenumber contributing to the heat transport is marginally stable with respect to a given mean temperature profile. As Howard (1963) has shown, however, the first assumption can be relaxed if the power integrals of the Boussinesq equations are considered in place of the Boussinesq equations. It was also shown that there will be a hierarchy of such constraint integrals that close the theory for each order of expansion. So, it may be more practical to appeal to a more mathematically complete means to compute the most probable equilibrium states, which is the variational technique established in Prigogine & Glansdorff (1965). Indeed, as Spiegel (1962); Townsend (1962) and Howard (1963) all point out, these optimal solutions are formed from solutions to the linear equations and hence they are not true solutions to the full nonlinear Boussinesq equations. However, this flaw can be mollified if, as in S79, it is conjectured that the convective modes are amplitude limited by instabilities in a nonlinearly saturated state.

The Rayleigh-Bénard experiments carried out in Townsend (1959) aim to directly test the heat-flux maximization principle of Malkus (1954), where the temperature and its gradient are measured to compute the heat flux and the moments of the temperature distribution. These values are then compared against the mean temperature gradient predicted in Malkus (1954). Townsend (1959) finds that indeed the measured mean temperature gradient is well-described by those expected from the heat-flux maximization principle, being inversely proportional to the vertical coordinate outside of the boundary layer. This is in contrast to the self-similarity solution of Priestley (1954), which predicts a scaling proportional to the inverse cube root of the vertical coordinate. Howard (1963) reexamines those experimental data with an emphasis on testing wether or not both the mean-fields and the fluctuating fields predicted in Malkus (1954) correspond to the measured flow structures, which are both found to match the theoretical values to within a factor of order unity. Additionally, Howard (1963) constructs a variational technique to assess the Rayleigh number ($Ra=-g\alpha_{T}\Delta T\ell_{0}^{3}/\nu\kappa$) and the Nusselt number (which measures the ratio of the total heat flux to the rate of thermal diffusion) of the heat-flux maximizing flows. It is also shown that, when the continuity equation is imposed in addition to the power integrals, the experimentally measured Nusselt number is asymptotically well-approximated by the variationally obtained value although its magnitude is overestimated by about a factor of two. This approach can also be taken to compare Malkus’s theory to more recent laboratory and numerical experiments of Rayleigh-Bénard convection (e.g., Funfschilling et al., 2005; Ahlers et al., 2006; Anders & Brown, 2017). Furthermore, similar experiments of Rayleigh-Bénard convection in rotating cylinders have measured the Nusselt number scaling with respect to the Rayleigh number and with the Ekman number [$\mathrm{Ek}=\nu/(2\Omega\ell_{0}^{2})$] in the rotating cases (Zhong et al., 2009; King et al., 2009; King & Aurnou, 2012; King et al., 2012, 2013; Cheng et al., 2015; Aurnou et al., 2015; Weiss et al., 2016). Thus, a Malkus-like turbulence theory that includes rotation, as in S79 and as shown here, can in principle also be compared to those experiments, which will be a focus of later work.

Laboratory experiments (Coates & Ivey, 1997; Weiss et al., 2016) and numerical simulations have lent some credence to this convection model (Käpylä et al., 2005; Barker et al., 2014; Sondak et al., 2015). In particular, those simulations indicate that the low convective Rossby number scaling regime established in S79 appears to hold up well for three decades in convective Rossby number and for about one decade in Nusselt number. Moreover, a recent implementation of Stevenson’s asymptotic scaling regimes in convective Rossby number have been used to provide a local model of a rotationally-dependent mixing-length theory (Ireland & Browning, 2018). There, it is shown that the entropy gradient of fully convective stars can be significantly modified by rapid rotation in the bulk of their interiors, excluding the surface region where the convective Rossby number is large. Therefore, the adiabat of the star is modified, leading to structural and evolutionary changes. What remains to be shown is how such a model of convection can impact the mixing for more modest rotation rates where the convective Rossby numbers are closer to unity, the depth of convective penetration, as well as the amplitude of wave-driven and shear-induced transport mechanisms.

2.3 General Framework

Following S79, the nonlinear saturation conjecture used in conjunction with the maximum heat transport principle requires that the primary contribution to each scale of the temperature and velocity fluctuation is from the convective term at the corresponding scale. This implies that its growth rate is $s=\left|\bm{\mathrm{v}}\right|\left|\bm{k}\right|$, which also provides a way to estimate the heat flux $F=\rho c_{P}\Re{[\mathrm{v}_{z}T^{\prime}]}$ with linear dynamics. Moreover, in the theory of Malkus and Howard, the mode that carries the largest fraction of the heat flux is the fundamental vertical mode, e.g. $n=1$ so that $k_{z}=\pi/\ell_{0}$, but whose horizontal structure is undetermined. It is this mode that Stevenson’s turbulence model is built around. Nominally, however, the maximization should be carried out scale-wise for each of the modes that contribute to the heat flux. This task will be left for later work that examines non-Boussinesq turbulence using Malkus’s theory as a point of comparison. Thus, Equation 1, with the assumptions used above to formulate the dispersion relationship, yields

[TABLE]

Taking the dot product of $\bm{\mathrm{v}}$ and Equation 11, using Equation 12 to eliminate the pressure, and taking a volume average it can be seen that

[TABLE]

Therefore,

[TABLE]

Employing the hypothesis that the velocity in the quasi-stationary turbulent state scales as $\langle\mathrm{v}^{2}\rangle=ss^{*}/k^{2}$, one has the heat flux

[TABLE]

where the definitions given in Equations 8 have been used to simplify the latter expression.

It is useful to eliminate some of the possible solutions to the maximization of this heat flux. First, consider the nondiffusive case, where

[TABLE]

from which it is clear that any value of $a_{y}$ will reduce the growth rate. This also places a lower limit on the value of $a$ for which the solutions are real with a growing branch, with $a_{x}>O\cos{\theta}$, which is shown to be satisfied by the heat-flux maximization in Stevenson (1979) and in the next subsection. Likewise, in the nonrotating but diffusive case, one has that

[TABLE]

This dispersion relationship has a growing real solution for sufficiently small $V$ and $K$ as

[TABLE]

When linearized with respect to $V$ and $K$, this yields

[TABLE]

which is maximized when $a_{y}=0$. So, in all the instances of parameter regimes considered below, the maximization will be carried out over $z^{3}=1+a_{x}^{2}$ with $a_{y}=0$. Note that this is equivalent to considering only the two-dimensional rolls aligned with the rotation axis, which in the rotating case is due largely to the Taylor-Proudman constraint. These modes also have a larger growth rate than three-dimensional modes for non-rotating Rayleigh-Bénard convection, but nevertheless the 2D modes can approximate many of the behaviors of the 3D modes (e.g., van der Poel et al., 2013).

It is also useful to recast some of the parameters so that the system can be more readily compared to the S79 results. The characteristic velocity $\mathrm{v}_{0}$ is derived from the growth rate and maximizing wavevector in the nonrotating and nondiffusive case, which are $s_{0}^{2}=3/5|g_{0}\alpha_{T}\beta_{0}|=3/5N^{2}_{*,0}$, with $\beta_{0}$ being the thermal gradient and $g_{0}$ being the effective gravity in the nonrotating case, and $k_{0}^{2}=5/2k_{z}^{2}$ as in (Stevenson, 1979) Equation 36, leading to

[TABLE]

This permits the definition of the convective Rossby number $\mathrm{Ro_{c}}$, which for this Boussinesq system is

[TABLE]

which implies that

[TABLE]

Next consider the variation of the superadiabaticity, which for this system is given by $\epsilon=H_{P}\beta/T$, meaning that $N_{*}^{2}=|g\alpha_{T}T\epsilon/H_{P}|$, where $H_{P}$ is the pressure scale height. The comparison made throughout this paper is between cases with additional included physical effects and the base case that is nondiffusive and nonrotating. The potential temperature gradient in this latter case is ascertained from the Malkus-Howard turbulence model, which yields a value of $N_{*,0}$. So far, all quantities have been normalized with respect to $N_{*}$. Instead, it is useful to compare them relative to $N_{*,0}$ and to introduce the ratio of superadiabaticities as an additional unknown as

[TABLE]

Therefore, all parametric quantities have the following equivalencies

[TABLE]

So, the dispersion relationship (Equation 9) and the heat flux (Equation 16) may be written as

[TABLE]

where $F_{0}=\langle\rho\rangle c_{P}N_{*,0}^{3}/\left(g\alpha_{T}k_{z}^{2}\right)$. When maximized for a positive real value of $\hat{s}$, Equation 27 implies that

[TABLE]

Depending upon the method of finding $\hat{s}$, it can be quite complicated to compute directly as a function of $z$. Instead, it is useful to compute the implicit derivative of the dispersion relationship (Equation 26) to obtain $d\hat{s}/dz$, which yields

[TABLE]

Finally, to assess the scaling of the superadiabaticity, the velocity, and the horizontal wavevector, a further assumption must be made in which the maximum heat flux is invariant to any parameters, namely that $\max{\left[F\right]}=\max{\left[F\right]}_{0}$ so the heat flux is equal to the maximum value $\max{\left[F\right]}_{0}$ obtained in the Malkus-Howard turbulence model for the nonrotating case. This is based on the assumption that the total heat flux should not change with rotation, which is not true in general as the centripetal acceleration of the star can lead to a variation in the central temperature and density of the star and thus to a variation of the total luminosity of the star at a fixed mass (e.g., Rieutord et al., 2016). Furthermore, the flux will be latitude dependent due to the rotational variation of the local gravity (e.g., Maeder, 1999; Wang et al., 2016). Indeed, an analysis of 3D nonlinear global-scale convection simulations of rapidly rotating oblate stars has been carried out in Wang et al. (2016). These simulations indicate that the Von Zeipel theorem (von Zeipel, 1924) holds well even for the emergent flux at the outer boundary, with the exception of a region near the equator where there is a flux enhancement. Thus, while acknowledging that the flux derived from the Von Zeipel theorem is not rigorously representative of the heat flux in a convection zone, it provides a very good approximation. Following Maeder (1999) and the von Zeipel theorem, one may show that the heat flux $F_{V}$ on an isobar with $P=\langle P\rangle$ varies as

[TABLE]

where $L\!\left(\langle P\rangle\right)$ is the total luminosity on the isobar and $M_{*}\!\left(\langle P\rangle\right)$ is the total mass integrated to the isobar. Approximating the effective gravity and the shift in the mass as in Maeder (1999) (Equation 2.3), this becomes

[TABLE]

where $r$ and $\theta$ correspond to the origin of the local coordinate system as in Figure 1. Taking the ratio of the rotating and non-rotating fluxes, one has that

[TABLE]

which in terms of the convective Rossby number and the buoyancy frequency is

[TABLE]

Noting the definition of the buoyancy frequency amplitude and the local density, one can see that

[TABLE]

Therefore, since the linearized change in flux is proportional to the Boussinesq expansion parameter $|\alpha_{T}\Delta T|$ and the inverse square of the convective Rossby number, this requires $\mathrm{Ro_{c}}>5^{-1}\pi^{-1}|\alpha_{T}\Delta T|^{1/2}$ for the approximations of the model to hold. Thus, for now, these effects will be ignored as they are sensitive to the global-scale dynamics, microphysics of the nuclear energy generation rate, and the equation of state. Therefore, building this convection model consists of three steps: defining a dispersion relationship that links $\hat{s}$ to $q$ and $z$, maximizing the heat flux with respect to $z$, and assuming an invariant maximum heat flux that then closes this three variable system.

In the case of planetary and stellar interiors, the viscous damping timescale is generally longer than the convective overturning timescale (e.g., $V_{0}\ll N_{*,0}$). Thus, the maximized heat flux invariance is much simpler to treat. In particular, the heat flux invariance condition under this assumption is then

[TABLE]

implying that

[TABLE]

where $\tilde{s}=2^{1/3}3^{1/2}5^{-5/6}$ and $\max{\left[F\right]}_{0}=6/25\sqrt{3/5}F_{0}$ follows from the definition of the flux and the maximizing wavevector used to define $\mathrm{v}_{0}$ above in Equation 21.

The assumption of this convection model is that the magnitude of the velocity is defined as the ratio of the maximizing growth rate and wavevector. With the above approximation, the velocity amplitude can be defined relative to the nondiffusive and nonrotating case scales without a loss of generality as

[TABLE]

So only the maximizing wavevector needs to be found in order to ascertain the relative velocity amplitude. For reference, the symbols that will be frequently used from this section are listed in Table 1.

2.4 Diffusion-free Models

As considered in S79, the simplest case to consider is one in which all diffusive mechanisms are neglected. Given Equation 28, the maximization condition requires

[TABLE]

where the derivative of the growth rate can be determined from the nondiffusive limit of Equation 29 as

[TABLE]

when equated this yields the following equation for the growth rate

[TABLE]

However, the dispersion relationship (Equation 26) for the nondiffusive case requires

[TABLE]

Using the heat flux invariance condition (Equation 37) so that $\hat{s}=\tilde{s}qz$, the two above equations can be manipulated to isolate a single equation for the maximizing value of $z$ and a dependent equation for $q$ as

[TABLE]

Together with Equations 25, these two equations yield the following scalings with convective Rossby number and maximizing wavevector as

[TABLE]

With the quintic form of Equation 46, its solutions are not representable as radicals or other functions. So, it is left in its current form to be solved numerically for a given colatitude and convective Rossby number. However, at low convective Rossby number, these equations can be asymptotically approximated: for $\mathrm{Ro}_{c}\ll 1$ it is clear that $z\propto\mathrm{Ro}_{\mathrm{c}}^{-2/5}$ since $z^{2}\ll z^{5}$. It then follows that $\epsilon/\epsilon_{0}\propto\mathrm{Ro}_{\mathrm{c}}^{-4/5}$ and $\mathrm{v}/\mathrm{v}_{0}\propto\mathrm{Ro}_{\mathrm{c}}^{1/5}$, all of which match the scalings given in S79, the direct numerical simulations of Käpylä et al. (2005) and Barker et al. (2014), and can be seen in Figure 2(a). In contrast, at very large convective Rossby number, these quantities converge to unity as expected.

At a first glance, one can see from Equations 45 and 46 that at the equator there is no rotational scaling of any of these quantities. However, there is a rapid transition from this behavior just a few degrees away from the equator. At a fixed convective Rossby number, the superadiabaticity has a latitudinal dependence of $\cos^{2}{\theta}$, meaning that its value at a larger colatitude is reduced with respect to its maximum value at the pole. Therefore, in order to maintain a constant heat flux, the velocity must increase toward the equator. In contrast, the horizontal wavevector is also a monotonically decreasing function of colatitude, with a maximum at the pole and minimum at the equator. Hence, the heat carrying motions are of a larger scale at larger colatitudes than at the pole as expected in theoretical considerations (e.g., Busse, 2002; Dormy et al., 2004) and as seen in numerical simulations of spherical rotating convection in both fully convective domains and in spherical shells (e.g., Glatzmaier, 1985; Miesch, 2005; Brun et al., 2005; Browning, 2008; Augustson et al., 2016). These behaviors are illustrated in Figure 2(b).

2.5 Adding Thermal Diffusion

The procedure developed above can be applied again in the case with thermal diffusion as the heat flux maximization condition remains unchanged from Equation 39. Taking the limit that $V\rightarrow 0$ in Equation 29, one finds that

[TABLE]

which implies the following heat flux maximization condition

[TABLE]

The dispersion relationship (Equation 26) with $V_{0}\rightarrow 0$ is

[TABLE]

Again utilizing the maximized heat flux constraint $\hat{s}=\tilde{s}qz$ and the definitions in Equations 25, the two previous equations can be shown to be equivalent to

[TABLE]

So the connection to the nondiffusive case is clear and can be recovered directly when taking the limit as $K_{0}\rightarrow 0$. Employing the relationship between $O$ and the convective Rossby number given in Equation 23, the superadiabaticity and wavevector can be defined as

[TABLE]

Immediately, one can see that the superadiabaticity increases as the thermal diffusion is increased because the thermal gradient must increase to drive convection. The nondiffusive component is unchanged and so it also grows with lower convective Rossby numbers. The form of Equation 53 indicates that the wavevector must increase with decreasing convective Rossby number, whereas its behavior with varying the thermal diffusion rate ($K_{0}$) is not obvious and is illustrated in Figure 4(a) for the fully diffusive case. There it is seen that thermal diffusion induces a shift in the value of the wavevector ratio that is a factor of $(7/10)^{1/3}$, or approximately 13%, lower for convective Rossby numbers below the value of $K_{0}$. For $K_{0}>1$, the asymptotic value of $z$ at large convective Rossby number is reduced by the same factor, from $(5/2)^{1/3}$ to $(7/4)^{1/3}$, indicating that the horizontal scale of the flows becomes larger. The velocity scales inversely with the horizontal wavevector; so, as before, the velocity must decrease with increasing rotation rate. These behaviors are shown for a typical value of the thermal diffusion in Figure 3(a). The velocity decreases with lower convective Rossby number given that the wavevector increases, with the same nontrivial behavior with the thermal diffusivity.

2.6 Including Viscosity

When including viscosity within this heuristic framework, the maximum heat flux should formally be derived from the algebraic system formed by the dispersion relationship for $\hat{s}$ and from

[TABLE]

In particular, the above equation, Equation 9, and Equation 28 give

[TABLE]

While this set of equations can be solved numerically, they can be simplified in the case of planetary and stellar interiors to a relatively simple scaling behavior with viscosity. In particular, the viscous component of the heat flux may be neglected if it is again assumed that $V_{0}\ll N_{*,0}$), so that as above $\hat{s}\approx\tilde{s}qz+\mathcal{O}(V_{0})$. Subsequently, following the method shown for the case of thermal diffusion, one may find the implicit wavevector derivative of the growth rate $\hat{s}$ (Equation 29), which implies that the constraining dispersion relationship and the equation resulting from the heat flux maximization are

[TABLE]

So one may solve for $q$ from the former equation, whose solution upon substitution into the latter equation yields an equation solely for the wavevector $z$:

[TABLE]

The horizontal wavevector is then recovered from the roots of the fourteenth-order Equation 60, whereas the superadiabaticity is defined as

[TABLE]

Clearly, the inclusion of viscosity increases the degree of superadiabaticity beyond that due to the thermal diffusion. Moreover, as before, it also grows with lower convective Rossby numbers. Given the form of Equation 60, the wavevector must increase with decreasing convective Rossby number, whereas again its behavior with varying $K_{0}$ is not obvious, but is relatively benign as is shown in Figure 4(a). In contrast, viscosity has a much larger impact on the amplitude of the change in the wavevector as exhibited in Figure 4(b). The velocity decreases with lower convective Rossby number given that the wavevector increases, with the same nontrivial behavior seen earlier in the case with thermal diffusion.

However, for values of $\left\{K_{0},V_{0}\right\}<N_{*,0}$, one may find a heuristic fit to the horizontal wavevector for the diffusive case given by

[TABLE]

which can be used more directly to estimate the scaling of the relative superadiabaticity and velocity amplitudes. The form of the fit follows from the observation of the qualitative features of the scaling of the horizontal wavevector given successively in Figures 2, 3, and 4, where it is well-described by hyperbolic tangents to within a few percent even in the transitional region of $\mathrm{Ro}_{\mathrm{c}}\approx 1/10$. The accuracy of the fit is illustrated in Figure 4(b).

2.7 Diffusion Approximation for Turbulent Transport

In stellar models, a diffusive treatment of transport processes can be adopted for thermodynamics, chemicals, angular momentum, and magnetic field (e.g., Heger et al., 2000; Maeder, 2009; Mathis, 2013). Within a convection zone, the turbulent mixing of those quantities can be approximated through a parameterized vertical diffusion $D_{\mathrm{v}}=1/3\,\alpha H_{P}\mathrm{v}_{\mathrm{MLT}}$. Since this depends upon the pressure scale height $H_{P}$, it is worth noting that it impacts the superadiabaticity as well. The definition of the Brunt-Väisälä frequency is $N^{2}=g\alpha_{T}T\epsilon/H_{P}$. Thus, allowing $H_{P}$ to vary, one has that $N^{2}/N_{0}^{2}=xH_{P,0}/H_{P}=x/h$. To maintain the constancy of the heat flux, this then requires $x=\epsilon/\epsilon_{0}$ to depend linearly on $h$, with $x\rightarrow xh$.

Within the context of this convection model and for a constant mixing-length parameter $\alpha$, the diffusion coefficient ratio scales as

[TABLE]

where $D_{0}$ is the local diffusion coefficient in the absence of rotation. The scaling of $D_{\mathrm{v}}$ is shown in Figure 6 and described in detail in the following section, as it suffices to note that the diffusion will increase with scale-height ratio $h$ and decrease as the rotation rate is increased. This convection model does however exclude any horizontal diffusion associated with the increased horizontal shear typically present in rotating convection.

3 Convective Penetration

3.1 Models Excluding Rotation

Standard MLT is not able to address convective penetration directly, rather it must be modified to account for this nonlocal effect. Such nonlocalities arise naturally in the generalizations of MLT as they directly model nonlinear interactions. One of the first and simplest such nonlocal generalizations of MLT can be found in Spiegel (1963). There the convective heat flux is treated as an unknown distribution in position and velocity space and as such it follows the dynamics described by a collisional Boltzmann equation. The collisional processes are a stochastic source arising from nonlinear processes and a damping term. The latter is the connection to MLT as it limits the path length of a convective eddy to the mixing length $\ell_{\mathrm{MLT}}$, but it is not required to be small in some sense as it needs to be in standard mixing length theory. Another approach is to include the kinetic energy flux, which leads to a convection theory that depends upon the entropy profile in the star and thus upon global integral constraints upon the energy flux in the star (e.g., Roxburgh, 1978). Other models have also been developed that variously extend MLT to include nonlinear processes through finite amplitude analyses (e.g., Veronis, 1963; Sparrow et al., 1964). Subsequently, modal expansions of the Boussinesq and anelastic equations provided solutions that established that nonlinear penetration was substantially deeper than what linear theory predicted (e.g., Musman, 1968; Moore & Weiss, 1973). In anelastic convection, the density stratification induces asymmetries between upflows and downflows that further increases the penetration depth (e.g., Zahn et al., 1982; Massaguer et al., 1984). Yet as is often the case, the penetration depth determined in these models depends sensitively on their assumptions, with a wide disparity in the calculated penetration depths. The problems and inconsistencies with many of these models were reviewed in Renzini (1987).

Three pioneering studies analytically estimate the convective penetration depth by prescribing the convective flows as a boundary condition and examining the impact on the region of penetration. One of the first models to attempt to estimate the importance of overshooting above a convective core is from Roxburgh (1965), in which ergodicity and isotropy are assumed, leaving a monodimensional and time-independent equation of motion. A linearization of the thermodynamics with respect to the vertical displacement then permits the equation of motion to be integrated to yield an estimate of the depth of penetration of a fluid element. Building upon the work of Roxburgh (1978), Schmitt et al. (1984) uses analytical models of turbulent plume ensembles and examines their statistical properties, finding that the penetration depth scales as $\mathrm{v}_{p}^{3/2}f^{1/2}$, where $\mathrm{v}_{p}$ is the plume velocity and $f$ is the fractional area filled by the plumes. Likewise, following Roxburgh (1965) and Roxburgh (1978) and by directly parameterizing the filling factor $f$, Zahn (1991, hereafter Z91) found a similar scaling result through a spatial linearization of the equations of motion and the energy transport. Both of these studies predict penetration depths as fairly large fractions of a pressure scale height. Z91 also confirmed the modal results of Massaguer et al. (1984) where the flow asymmetry due to the density stratification causes downflows to penetrate farther into an underlying stable layer than upflows penetrate into an overlying layer. As a first estimate of the convective penetration depth, one may consider the Z91 linearized model for convective penetration. This model neglects the effects of rotation in both the physical description of the model as well as in its parameterization of the convective dynamics. The former modification will be saved for later work. However, the extended version of the S79 model derived above provides a means of estimating the latter alteration. In particular, one can harness the linearized framework of the Z91 model to obtain an order of magnitude estimate of the depth of penetration while allowing for the effects of rotation to be included through a modified value of the superadiabaticity and mixing length velocity within the convective region. It is also a formalism that may be easily implemented in stellar evolution codes.

3.2 Convective Penetration with Rotation

3.2.1 Theory

The Z91 model has four phenomenologically distinct zones: a convective zone that ends when the radiative energy flux becomes equal to the convective energy flux, a subadiabatic penetrative zone where flows render the region nearly adiabatic, a transitional thermal boundary layer with overshooting convection where the Péclet number is reduced to below unity, and a radiative zone where waves and conduction carry the total energy flux. There are a few basic assumptions underlying the model. The first assumption is that the velocity and the temperature share the same planform, namely that they are perfectly correlated. The second assumption is that only the downflows for downward penetrating flows and upflows for upward penetrating flows are effective at carrying enthalpy, which is parameterized through their filling factor $f$, which results from a horizontal average over their planform $\mathrm{v}_{z}(x,y,z)=\mathrm{v}(z)h(x,y)$ such that $f=\overline{h^{2}(x,y)}$. From simulations, the value of $f$ appears to be about $1/3$, with slight variations of order of 10% that are sensitive to the level of turbulence, rotational, and magnetic effects (Brummell et al., 2002; Stein et al., 2009a, b). However, rotation does impact the structure of the convection as a function of colatitude, whose structure in turn depends upon the degree of supercriticality of the system. For rotating systems, the convective structures should possess a modicum of alignment with the rotational axis, which implies that the mechanism of penetration will vary with colatitude. For instance, consider that columnar flow structures will mostly conform to the classical paradigm of radially-aligned penetrating flows near the pole, but they will increasingly become horizontally-shearing, yet still localized, flows at lower colatitudes. This will be combined with any shear-induced mixing by large-scale flows such as differential rotation (e.g., Zahn, 1992; Maeder, 2003; Mathis et al., 2004; Mathis, 2013). Yet, for the purposes of this model, the filling factor $f$ shall be considered to be a fixed parameter of the theory. Additionally, it shall be assumed that the local magnitude of the gravity $g$, the local conductivity $\kappa$, the thermal expansion coefficient $\alpha_{T}$, the nuclear energy generation rate $Q$, and the specific heat capacity at constant pressure $c_{P}$ are unaffected by rotation, namely that the centripetal acceleration may be ignored to first order.

A further assumption is that the convective enthalpy flux can be linearized in the region of penetration such that

[TABLE]

where $M(r)$ is the mass, $L(r)$ is the luminosity, and $\chi(r)$ is the radiative conductivity such that $\chi=\rho c_{P}\kappa$ as seen in Z91 (Equation 4.4). Finally, another assumption is that the vertical inertial force balances the buoyancy force on average such that the pressure contribution vanishes and so that the density perturbations can be linearized with $\rho^{\prime}=\alpha_{T}T^{\prime}$. To follow Z91 as closely as possible, the system is considered only at the pole so that the direct effects of the local Coriolis acceleration $2\Omega_{0}\sin{\theta}v_{x}$ may be neglected, which is equivalent to having a buoyant braking timescale that is shorter than the rotational timescale. Instead, the Coriolis effect implicitly influences the penetration depth by modifying the upper boundary value of the velocity, which is taken from the convection zone. Thus, when horizontally-averaged and a pressure equilibrium is assumed as in Z91 (Equation 3.8), the force balance becomes

[TABLE]

where $c=\overline{h^{3}}/\overline{h^{2}}$. This equation can then be linked to the convective energy flux through Equation 64 and integrated across the penetrative region to yield an estimate for the penetration depth $L_{P}$ relative to the pressure scale height $H_{P}$.

In Z91, there are two cases for the scaling of the convective penetration depth with velocity and spherically-symmetric thermodynamic quantities: one for penetrative convection into a stable region below a convection zone, such as takes place near the base of the solar convection zone, and one for penetrative convection into a stable region above the convection zone, such as takes place near the convective cores of intermediate mass and high mass stars. The first of these two regimes neglects the variation in the total luminosity and mass and so only retains the term related to changes in the radiative conductivity. It scales as

[TABLE]

where $\nabla_{\mathrm{ad}}$ is the adiabatic temperature gradient and $\chi_{P}=\partial\ln{\kappa}/\partial\ln{P}|_{S}$ is the adiabatic logarithmic derivative of the radiative conductivity with respect to pressure. The opposite case of convective penetration into a stable region above a convection zone scales as

[TABLE]

where $r_{0}$ is the radius at the edge of the core. Note that $\nabla_{\mathrm{ad}}=\nabla+\epsilon$, with $\nabla$ being the temperature gradient and $\epsilon$ being the superadiabatic gradient as above. However, the basic assumption of the model is that $\epsilon$ does not grow large enough to modify the background temperature gradient in a steady state, and so its variation with rotation rate does not strongly influence the depth of penetration. Thus, the ratio of the penetration depth with rotation and diffusion to the nonrotating inviscid value for convective penetration into a stable layer either above or below a convection therefore scales as

[TABLE]

As seen in the previous section, the velocity amplitude of the mode that maximizes the heat flux decreases with lower diffusivities and lower convective Rossby numbers. Therefore, the penetration depth necessarily must decrease when the convective Rossby number is decreased. This behavior follows intuitively given that the reduced vertical momentum of the flows implies that the temperature perturbations are also reduced (via Equation 13). Thus, due to the decreased buoyant thermal equilibration time and the reduced inertia of the flow the penetration depth must decrease. In contrast, the velocity and the horizontal scale of the flow increase with greater diffusivities in order to offset the reduced temperature perturbations in the case of a larger thermal conductivity. In the case of a larger viscosity, the horizontal scale of the velocity field is increased, whereas, for a fixed thermal conductivity, the thermal perturbations are of a smaller scale. Thus, to maintain the heat flux, the amplitude of the velocity must increase in order to compensate for the reduced correlations between the two fields. The scaling behaviors of the penetration depth are illustrated as a function of diffusivities and convective Rossby number in Figure 5.

3.2.2 Comparison with Penetrative Convection Experiments

Penetrative convection has been studied extensively in fully compressible and anelastic numerical experiments in 2D within the context of the solar convection zone (Hurlburt et al., 1986, 1994; Rogers & Glatzmaier, 2005; Rogers et al., 2006), for the convective cores of intermediate mass stars (Rogers et al., 2013), and in the context of modeling laboratory experiments for temperature stratified water (Lecoanet et al., 2015; Couston et al., 2017; Toppaladoddi & Wettlaufer, 2017). These studies all tend to find that penetration into stable layers below the convection zone is more extensive than into overlying stable layers and that the depth of penetration depends sensitively on the stiffness of the interface between them and upon the strength of the convective driving. Moreover, in rotating convection models of equatorial planes of intermediate mass stars, it is directly seen that rotation reduces the overshooting depth as the Coriolis force leads to an azimuthal deflection of convective plumes and to less available radial kinetic energy (Rogers et al., 2013). While 2D simulations are a useful first step to understanding convective penetration, they are limited in assessing the convective properties of stellar interiors given the 3D nature of the convection there, which leads to much different convective structures and spectral energy transfer properties.

In constrast to the behavior of penetrative convection seen in local simulations and from scaling arguments such as those described above, the depth of penetration in fully nonlinear global-scale 3D simulations tends to increase with decreasing convective Rossby number and in spherical geometry (e.g., Miesch et al., 2000; Miesch, 2005; Browning et al., 2004; Brun et al., 2011; Augustson et al., 2012; Augustson, 2013; Alvan et al., 2014; Augustson et al., 2016; Brun et al., 2017). However, these simulations are typically in a low Péclet number regime $Pe\approx\mathcal{O}(10)$ and thus are still influenced by thermal diffusion. The stiffness of the interface often has been a priori softened so that the simulation is more computationally feasible, meaning that the degree of convective overshoot relative to penetration is larger than in a stellar context and that the excitation mechanisms are also more dominated by Reynolds stresses rather than buoyantly driven. So, it may be that the stiffness of the transition region is too weak in current global-scale modelling approaches, permitting large-scale entraining and overshooting flows rather than restricting them to plume-like dynamics. Moreover, the sweeping motions of the relatively laminar globally connected flows induce mixing on a scale of the correlation length of those flows (Viallet et al., 2015), which enhances the depth of penetration and is something that will be assessed in later work.

Three-dimensional simulations have been conducted that can simultaneously resolve gravity waves as well as the mechanisms that excite them, namely shear and convective penetration in both nonrotating Cartesian settings (e.g., Hurlburt et al., 1994; Lecoanet et al., 2015; Couston et al., 2018) and f-planes (e.g., Julien et al., 1996; Brummell et al., 2002; Pal et al., 2007, 2008). As in 2D simulations, these 3D simulations also find the following: penetration into stable layers below the convection zone is more extensive than into overlying stable layers, the depth of penetration depends sensitively on the stiffness of the interface between them, and the penetration depth is latitudinally dependent.

In the 3D f-plane simulations of rotating convection described in Julien et al. (1996) and Brummell et al. (2002), it is found that the penetration depth into a stable layer below a convective region scales as $L_{P}\propto\mathrm{Ro}_{\mathrm{c}}^{0.15}$, that is it decreases with increasing rotational influence, due primarily to a reduction in the flow amplitude. Likewise, in a similar suite of f-plane simulations examined in Pal et al. (2007), it is found that there is a decrease in the penetration depth with increasing rotation rate that scales as $L_{P}\propto\mathrm{Ro}_{\mathrm{c}}^{0.2}$ at the pole and to $L_{P}\propto\mathrm{Ro}_{\mathrm{c}}^{0.4}$ at mid-latitude. For penetration into a stable layer above a convection zone, on the other hand, f-plane simulations from Pal et al. (2008) indicate that there is a much weaker rotational scaling, being statistically consistent with no scaling. However, the range of parameters examined is quite restricted due to the computational requirements to resolve structures as well as maintain highly supercritical flows at much lower convective Rossby numbers. Nevertheless, the depth of convective penetration as assessed in those numerical simulations appears to be roughly consistent with the heuristic model derived above, where $L_{P}/L_{P,0}\propto\mathrm{Ro}_{\mathrm{c}}^{3/10}$, which follows from $\mathrm{v}/\mathrm{v}_{0}\propto\mathrm{Ro}_{\mathrm{c}}^{1/5}$ in the nondiffusive and low convective Rossby number limit of the convection model.

Couston et al. (2017) have examined the influence of the stiffness of the convective-radiative transition in numerical analogs of a laboratory experiments involving the buoyancy transition of water that occurs near $4^{\circ}$ C, which is also a decent analog for convective stellar cores. In these 2D simulations, there is a smooth transition between plume dominated and entrainment dominated dynamical regimes that is independent of the Peclet number. Instead, it is sensitive to the stiffness of the convective-radiative transition region and the Rayleigh number. According to those local domain numerical simulations, and to the theory developed in Hurlburt et al. (1994), the interface stiffness provides a further distinction drawn between the regimes of penetrative convection. Particularly, the penetration depth depends upon the relative stability $S=(2n_{s}-3)/(3-2n_{c})$ between the stably-stratified and convective regions, where $n_{s}$ is the polytropic index in the stable region and $n_{c}$ is the index in the convective region. Indeed, Hurlburt et al. (1994) employs a linearized convection model similar to the one derived in Z91, which has a continuous thermal conductivity, but instead the conductivity is taken to be piecewise constant. Using this model, it can be shown that the total penetration depth depends upon the depth of the adiabatic layer and the depth of the thermal adjustment layer. However, the depth of the two regions scale differently with the relative stability parameter $S$. In particular, the depth of the adiabatic layer scales as $L_{a}\propto S^{-1}$ and the thermal adjustment layer scales as $L_{t}\propto S^{-1/4}$. Therefore, as $S$ increases the depth of the adiabatic layer decreases significantly, leaving the thermal adjustment layer as the primary means of defining the depth of penetration. Note that as stated in Zahn (1991) (Equation 3.18), the penetration model developed there and expanded on above scales proportionally to the inverse of the interface stiffness. For the simulations of penetration into a stably-stratified layer below a convection zone, Brummell et al. (2002) and Pal et al. (2007) find that the depth of penetration scales primarily as $L_{P}\propto S^{-1/4}$, for values of $S$ between about unity and around 30. Yet for simulations of convective penetration into a stably-stratified region above a convection zone, the opposite behavior is found, namely that $L_{P}\propto S^{-1}$ for $S$ between unity and ten. Given the similarity in the parameters used in those simulations, one might conclude that the penetration into a stably-stratified layer below a convection zone has physics more akin to a thermal adjustment layer, whereas penetration into a layer above is more akin to an adiabatic penetrative layer.

3.3 A Diffusive Approach

This model of penetrative convection can be introduced into stellar models utilizing a diffusive approach. Such a parameterization of mixing processes has been extensively examined (e.g., Freytag et al., 1996; Herwig et al., 1999; Denissenkov et al., 2013; Viallet et al., 2015; Lecoanet et al., 2016). A complementary model has been established through an extreme-value statistical analysis of 3D penetrative convection simulations (Pratt et al., 2017b). A frequent assumption made of turbulent flows is to assume that the underlying distribution of a flow quantity is normal. However, it is shown through the simulations and analysis of Pratt et al. (2017b) and Pratt et al. (2017a) that the more rare events occurring in the tails of the actual distribution function, e.g. the extremal penetrative flows with a higher velocity and a greater entropy deficit than the average flow, have a larger impact than is expected in Gaussian turbulence models on the turbulent mixing coefficients. This analysis was carried out with a novel method to quantify this non-Gaussianity of the dynamics. Their method permits a statistical examination of the turbulent particle dispersion and subsequently the construction of a model for a turbulent diffusion based upon the Gumbel distribution (Pratt et al., 2017a). This model for diffusion and overshoot has been applied to the problem of lithium depletion occurring in rotating low-mass stars (Baraffe et al., 2017), where it was found to mimic the observed trends well when the effects of rotation on the diffusion treatment were introduced. There the rotational effects are included as an empirically defined threshold rotation rate above which the diffusion is quenched ostensibly due to the action of rotation on the convective flows. Also, the penetrative convective depth was a parameter and was not dynamically estimated. However, the initial estimate of Baraffe et al. (2017) may be improved. Using the above extension of the Z91 model, one can estimate both the penetration depth and the level of diffusion. Taking the parameters of the Gumbel distribution as in Pratt et al. (2017a), yields the following description of the radial dependence of the diffusion coefficient

[TABLE]

where $r_{c}$ is the base of the convection zone and where $\mu$ and $\lambda$ are the empirically determined parameters from Baraffe et al. (2017). An illustration of the scaling behavior due both to the variation of $D_{\mathrm{v}}$ with the mixing length velocity $\mathrm{v}$ and with the depth of penetration $L_{P}$ are shown in Figure 6, where the parameters are taken to make contact with Sun-like stars where the transition region begins around $r\approx 0.7R_{\sun}$. Note, however, the convective velocity will vary in the convection zone. So Figure 6 is only meant to illustrate the rough behaviors of this diffusive model near the convectively stable-to-unstable transition point. From Equation 70 and as seen in Figure 6, the radial structure of the diffusion coefficient follows from the scaling of the velocity, namely the diffusion will globally decrease with decreasing convective Rossby number, and will further vary if the local angular velocity is taken into account. The depth of penetration is perhaps most notable, in that its strong rotational dependence can lead to severe restrictions on the region in which the diffusion acts. This is evidenced through the transition between a long diffusive tail at high convective Rossby number to a step-like function at low convective Rossby number. Thus, as discussed earlier, the amount of and the depth over which mixing induced above or below a convective region may be severely reduced in stars whose convection possesses a low convective Rossby number.

In Lecoanet et al. (2016), the mixing induced below a convection zone is characterized as a means of determining whether or not carbon flames can be disrupted by overshooting convection in evolved massive stars of between about 7-11 solar masses. To estimate the amount of mixing, the results of an idealized model of these flames are assessed with 3D hydrodynamic simulations of the Boussinesq equations. In particular, a passive scalar field is evolved that permits the measurement and fitting of a height-dependent turbulent chemical diffusivity that mimics the effect of the overshooting convection. This diffusion coefficient is modeled as a composition of two error functions, which has a super-exponential tail in the mixing region somewhat similar to the model constructed above.

What is currently lacking, however, are extensive parameter studies using 3D simulations of convective overshoot that also include rotation. Such simulations would at least provide numerical evidence regarding the hypothesized dependence of the mixing depth and its shape. For instance, the simulations carried out in Korre et al. (2018) have also addressed the shape of the overshooting region in a Boussinesq system and its associated potential mixing properties in 3D simulations. As in Lecoanet et al. (2016), it is found that the horizontally-averaged kinetic energy is super-exponentially reduced with depth below the convection zone. Specifically, in those simulations, the average kinetic energy is best modelled with a Gaussian. Moreover, the radial correlation length scale downflows is larger than the length scale associated with the averaged kinetic energy profile. This length scale may be interpreted as depth to which the strongest downflows travel before stopping. This behavior is at least roughly consistent with the model proposed here, and with the results of Pratt et al. (2017b) and Pratt et al. (2017a), although it is in contrast to the models proposed in Freytag et al. (1996) and Herwig (2000). For the time being, the examinations in Julien et al. (1996) and Brummell et al. (2002) will have to suffice as evidence that the overshooting is reduced with rotation rate, while the shape of the mixing region likely depend upon the model setup and equations solved in similar simulations.

4 Summary and Discussion

A model of rotating convection originating with Stevenson (1979) has been extended to include thermal and viscous diffusion for any convective Rossby number. Moreover, a systematic means of developing such models for an arbitrary dispersion relationship have also been shown. An explicit expression is given for the scaling of the horizontal wavenumber in terms of the convective Rossby number and diffusion coefficients under the constraint that the values of the diffusive time scale are less than the convective time scale (Equations 46, 53, and 60). The scalings of the velocity and superadiabaticity in terms of that wavenumber are also given. Asymptotically at low convective Rossby number and without diffusion, these match the expressions given in Stevenson (1979) (see Figure 2), as well as the numerical results found in the 3D simulations of Käpylä et al. (2005) and Barker et al. (2014).

In §3, the model of rotating convection is employed to assess the convective Rossby number scaling of the depth of convective penetration, utilizing the linearized model of Zahn (1991). Due to the reduced velocity and increased superadiabaticity at lower convective Rossby number, the penetration depth decreases proportionally to $\mathrm{Ro}_{\mathrm{c}}^{3/10}$ when diffusive processes are neglected. This estimate of the penetration depth is then employed to construct an estimate for the diffusive mixing coefficient in the region of penetration based upon the numerical experiments of Pratt et al. (2017b) and Pratt et al. (2017a). In the 3D f-plane simulations of rotating convection of Julien et al. (1996) and Brummell et al. (2002) it is found that $L_{P}\propto\mathrm{Ro}_{\mathrm{c}}^{0.15}$. Similarly, in the simulations of Pal et al. (2007) and Pal et al. (2008) it is found that $L_{P}$ varies in latitude, with $L_{P}\propto\mathrm{Ro}_{\mathrm{c}}^{0.2-0.4}$ from the pole to mid-latitudes. Hence, the depth of convective penetration as assessed in those numerical simulations is roughly consistent with those that follow from the coupling of the convection model with the results of Zahn (1991).

This convection model will soon be exploited to begin to study the impact of reduced convective penetration in stellar and planetary evolution models when rotation is included. Yet its impact on the structure of the convection can already be anticipated. From Figure (3), since this convection model is local, a qualitative picture of those aforementioned impacts can be constructed as follows: in regions where the local convective Rossby number is less than unity, the mixing-length velocity will be reduced and the superadiabatic temperature gradient will be enhanced. Therefore, depending upon the value of the superadiabatic gradient in the nonrotating case, the full temperature gradient may increase sufficiently to modify the location of the ionization zones as well as regions of large opacity changes. However, most of those thermodynamic changes will be felt most keenly deeper within the convection zone where the convective Rossby number can be smaller than near the photosphere for stars, where the convective Rossby number is typically of order unity or larger. Near the boundary of convective regions the typical mixing-length velocity is small due to the increasing importance of the radiative transport of energy. Thus, in those regions, the local convective Rossby number can be quite small and yet the thermal diffusivity can become larger relative to the dynamical time scale. With that in mind, and again appealing to Figure (3), it is clear that the velocity will be even further reduced and the superadiabatic gradient further increased. This implies that the transition to the region of convective penetration and the radiative zone itself will be shallower (or deeper for a convective core) when compared to the nonrotating case. Hence, in addition to the reduced depth of convective penetration with decreasing convective Rossby number seen in §3, the convection zone depth will be reduced when considering convection in the presence of rotation with this model. Due to decreased diffusion and transport, sharper thermodynamic and chemical gradients may be present in convectively stable regions.

Being a first step, there are several improvements that can be made to this convection model. One clear omission in this work is the impact of the magnetic field. However, magnetism has been considered in both Stevenson (1979) and Canuto & Hartke (1986). An alternative model for the impact of the magnetic field on the superadiabaticity, derived from the variational principle of Bernstein et al. (1958), can be found in Gough & Tayler (1966). Thus, convection and magnetic field will be the focus of a paper in this series. Second, the conjecture of Stevenson (1979) gives the scaling of only the dominant heat-carrying mode. To build a theory of the turbulent spectrum and to assess the changing cascade of energy and helicity, one could extend the theory in a manner similar to that shown in Malkus (1954), with a mode by mode analysis of the heat flux to construct its accompanying spectrum. One attempt at generalizing this turbulence model has already been undertaken utilizing the Heisenberg-Kolmogorov turbulence model, as shown in Canuto & Hartke (1986). It would also be prudent to further numerically assess the validity of these scaling laws with 3D simulations with a larger range of Nusselt number and to examine the impact of diffusion. The theory can be further improved by considering more sophisticated models of the structure of the flows, such as applying the results obtained for rotating plumes considered in Pedley (1968) or Grooms et al. (2010). Additionally, one relatively simple improvement to the treatment of convective penetration will be to include the effects of rotation in the linearized equations of motion in an f-plane, expanding further upon the Zahn (1991) model.

These waves, which are internal gravity waves modified by rotation through the Coriolis acceleration (e.g., Dintrans & Rieutord, 2000), propagate in stably stratified stellar radiation zones. They provide a means for transporting angular momentum and mixing chemical species in these regions (e.g., Pantillon et al., 2007; Mathis et al., 2008; Mathis, 2009; Lee et al., 2014). They are one of the mechanisms that may explain the strong extraction of angular momentum in the Sun (e.g., Talon & Charbonnel, 2005), in subgiant and red giant stars (e.g., Fuller et al., 2014; Belkacem et al., 2015a, b; Pinçon et al., 2017), and in intermediate-mass and massive stars (e.g., Rogers, 2015) as revealed by the weak differential rotation observed thanks to helio- and asteroseismology (e.g., García et al., 2007; Beck et al., 2012; Mosser et al., 2012; Deheuvels et al., 2012; Kurtz et al., 2014; Saio et al., 2015; Murphy et al., 2016; Aerts et al., 2017; Gehan et al., 2018). However, their amplitude and frequency spectrum strongly depend on the properties of the turbulent convective flows that excite them stochastically at the radiation/convection interface and in the bulk of convective regions (e.g., Schatzman, 1993; Zahn et al., 1997; Rogers et al., 2006; Belkacem et al., 2009b; Lecoanet & Quataert, 2013; Rogers et al., 2013; Alvan et al., 2014). Therefore, one should consider the impact of the modification of their convective excitation source by rotation. Although some first attempts have been made to do this (Belkacem et al., 2009a; Mathis et al., 2014; Rogers, 2015), a systematic study should be done. This will be adressed in the second article of this series.

Another possible use of the convection model developed here is for tidal dissipation. Turbulent friction is one of the crucial physical mechanisms driving the dissipation of tidal flows in stellar and planetary convective regions (e.g., Zahn, 1966, 1989; Ogilvie & Lesur, 2012). This friction acts both on the equilibrium tide and on tidal inertial waves in these layers. Recently, in Mathis et al. (2016), the asymptotic version of Stevenson’s theoretical scaling laws was used to construct a corresponding local model of tidal waves in order to understand the consequences of rotating convective turbulence for linear tidal dissipation. With this convective model, the turbulent friction acting on the tides was significantly reduced in rapidly rotating objects. Thus, to better capture the behavior of the tidal friction in the convective regions in stars and planets at modest convective Rossby number, it will be useful to apply the model of convection developed here.

Acknowledgments

The authors thank the anonymous referee for their constructive comments, which have improved the description of the convection model and its implications. K. C. Augustson and S. Mathis acknowledge support from the ERC SPIRE 647383 grant and PLATO CNES grant at CEA/DAp-AIM. The authors also thank Q. André, A. Astoul, M. Browning, A. S. Brun, V. Prat, R. Raynaud, and J. Toomre for fruitful conversations.