Validity of Sound-Proof Approaches in Rapidly-Rotating Compressible Convection: Marginal Stability vs. Turbulence

TL;DR

This study evaluates the validity of sound-proof approximations, specifically the anelastic and pseudo-incompressible models, in rapidly-rotating compressible convection, highlighting their regimes of accuracy and limitations through linear and nonlinear simulations.

Contribution

It clarifies the conditions under which the anelastic and pseudo-incompressible approximations remain valid in rapidly-rotating compressible convection regimes.

Findings

Pseudo-incompressible approximation remains valid at marginal stability.

Anelastic approximation converges to compressible results at high Rayleigh numbers.

Both models show good agreement in highly supercritical turbulent regimes.

Abstract

The validity of the anelastic approximation has recently been questioned in the regime of rapidly-rotating compressible convection in low Prandtl number fluids (Calkins et al. 2015). Given the broad usage and the high computational efficiency of sound-proof approaches in this astrophysically relevant regime, this paper clarifies the conditions for a safe application. The potential of the alternative pseudo-incompressible approximation is investigated, which in contrast to the anelastic approximation is shown to never break down for predicting the point of marginal stability. Its accuracy, however, decreases close to the parameters corresponding to the failure of the anelastic approach, which is shown to occur when the sound-crossing time of the domain exceeds a rotation time scale, i.e. for rotational Mach numbers greater than one. Concerning the supercritical case, which is naturally characterised by smaller rotational Mach numbers, we find that the anelastic approximation does not show unphysical behaviour. Growth rates computed with the linearised anelastic equations converge toward the corresponding fully compressible values as the Rayleigh number increases. Likewise, our fully nonlinear turbulent simulations, produced with our fully compressible and anelastic models and carried out in a highly supercritical, rotating, compressible, low Prandtl number regime show good agreement. However, this nonlinear test example is for only a moderately low convective Rossby number of 0.14.

Equations

$$\mu =\rho \nu ,k_t=c_p\rho \kappa ,$$

$${\text{upartial}}_t{\rho }_1=$$

$${\rho }_0{\text{upartial}}_t\mathbf{v}=$$

$${\rho }_0{T}_0{\text{upartial}}_t{s}_1=$$

$$p_0=(c_p-c_v){\rho }_0{T}_0,\frac{p_1}{p_0}=\frac{{\rho }_1}{{\rho }_0}+\frac{T_1}{T_0},$$

$$\frac{{\text{upartial}}_z{s}_0}{c_p}=\frac{{\text{upartial}}_z{T}_0}{T_0}-\frac{c_p-c_v}{c_p}\frac{{\text{upartial}}_z{p}_0}{p_0},\frac{s_1}{c_p}=\frac{T_1}{T_0}-\frac{c_p-c_v}{c_p}\frac{p_1}{p_0}.$$

$${\text{upartial}}_z{p}_0=-{\rho }_{0}g,{\text{upartial}}_z\left(k_t{\text{upartial}}_z{T}_0\right)=-S,$$

$${\text{upartial}}_z{T}_0=-\left(\frac{g}{c_p}+\frac{\Delta T}{d}\right)$$

$$\displaystyle D\,=\,{gd}\big{/}({c_{p}T_{r}})\,,$$

$$\displaystyle\epsilon\,=\,{\Delta T}\big{/}{T_{r}}\,.$$

$$n=\frac{\gamma }{\gamma -1}\frac{D}{D+ϵ}-1,$$

$$\displaystyle\gamma\,=\,{c_{p}}\big{/}{c_{v}}\,=\,{5}/{3}$$

$$n_{ad}=\frac{1}{\gamma -1}=\mathrm{1.5.}$$

$$\frac{ϵ}{D}=\frac{n_{ad}-n}{n+1}.$$

$$\displaystyle N_{\rho}\,=\,\ln\bigl{(}{\rho_{\mathrm{bot}}}\big{/}{\rho_{\mathrm{top}}}\bigr{)}\,=\,\,-\,n\,\ln\left[1-(D+\epsilon)\right]$$

$$\displaystyle Pr\,=\,{\nu_{r}}\big{/}{\kappa_{r}}\,,$$

$$\displaystyle Ra\,=\,{gd^{3}\Delta T}\big{/}\bigl{(}{\kappa_{r}\nu_{r}T_{r}}\bigr{)}\,,$$

$$\displaystyle Ek\,=\,{\nu_{r}}\big{/}\bigl{(}{2{\varOmega}d^{2}}\bigr{)}\,,$$

$${\hat{T}}_0(\hat{z})=$$

$${\hat{p}}_0(\hat{z})=$$

$$\frac{ϵD}{\gamma -1}\frac{1}{{\hat{T}}_0}{\widehat{\text{upartial}}}_t{\hat{p}}_1-ϵ{\hat{\rho }}_0{\widehat{\text{upartial}}}_t{\hat{s}}_1=$$

$${\hat{\rho }}_0{\widehat{\text{upartial}}}_t\hat{\mathbf{v}}=$$

$${\hat{\rho }}_0{\hat{T}}_0{\widehat{\text{upartial}}}_t{\hat{s}}_1=$$

$${\hat{\rho }}_1=\frac{D}{\gamma -1}\frac{1}{{\hat{T}}_0}{\hat{p}}_1-{\hat{\rho }}_0{\hat{s}}_1,{\hat{s}}_1=\frac{{\hat{T}}_1}{{\hat{T}}_0}-D\frac{{\hat{p}}_1}{{\hat{p}}_0},$$

$${\rho }_{\mathrm{cont}}=$$

$$p_{\mathrm{cont}}=$$

$$M_{\mathrm{phase}}^{\mathrm{top}}$$

$$M_{\mathrm{rot}}$$

$${\hat{v}}_{\mathrm{rms}}=\sqrt{\frac{1}{\hat{V}}\int \mathrm{d}\hat{V}{\hat{\mathbf{v}}}^2},$$

$$Nu\,=\,-\,\upartial_{z}\langle\hat{T}_{1}\rangle\bigr{|}_{\hat{z}=0}\,,$$

$${\hat{\omega }}_z={\widehat{\text{upartial}}}_x{\hat{v}}_y-{\widehat{\text{upartial}}}_y{\hat{v}}_x$$

$$ϵ{\widehat{\text{upartial}}}_t{\hat{\rho }}_1+\hat{\nabla }\cdot \left[({\hat{\rho }}_0+ϵ{\hat{\rho }}_1)\hat{\mathbf{v}}\right]=$$

$$\displaystyle(\hat{\rho}_{0}+\epsilon\hat{\rho}_{1})\bigl{[}\hat{\upartial}_{t}\hat{\bm v}+(\hat{\bm v}{\bm\cdot}\hat{\bm\nabla})\hat{\bm v}\bigr{]}\,=\,$$

$$-\frac{1}{Ek}\sqrt{\frac{Pr}{Ra}}({\hat{\rho }}_0+ϵ{\hat{\rho }}_1)\hat{\mathbf{z}}\times \hat{\mathbf{v}},$$

Full text

1 Introduction

The interiors of stellar and many planetary bodies can be characterised as compressible rapidly-rotating fluids that have small viscous compared to thermal diffusivities, i.e., Prandtl numbers much less than unity. These astrophysical fluid dynamical systems can be modeled with the fully compressible equations that follow from first principles of physics. Their generality comprises a broad range of temporal and spacial scales corresponding to physical phenomena reaching from sound waves over buoyancy induced flows up to atmospheric jets. Accounting for each of these processes, however, implies high computational costs for numerical simulations. In order to more adequately understand the structure and evolution of planetary and astrophysical objects, it is often favourable to channel all available resources into only the relevant physical phenomena, which has motivated the development of reduced forms of the fully compressible equations. For the purpose of modeling planetary and stellar convection, it is believed that dynamically unimportant sound waves can safely be neglected, which, due to the reduced spectrum of time-scales that needs to be resolved, effectively decreases the necessary computational effort.

Several so-called sound-proof approaches retaining the major compressibility effects are in use for modeling these kind of systems. The most widely employed ones are the anelastic approximation (e.g. Batchelor, 1953; Ogura and Phillips, 1962; Gough, 1969; Gilman and Glatzmaier, 1981; Braginsky and Roberts, 1995; Lantz and Fan, 1999; Jones et al., 2009; Glatzmaier, 2014; Verhoeven et al., 2015), the pseudo-incompressible approximation (Durran, 1989, 2008; Klein and Pauluis, 2012; Vasil et al., 2013; Wood and Bushby, 2016), and the closely related low Mach number approach (Majda and Sethian, 1985; Bell et al., 2004; Almgren et al., 2006a, b, 2014). Two necessary conditions for all these approaches to be valid are that the Mach number, which is defined as the ratio of the fluid velocity to the local sound speed, is much smaller than unity and that pressure perturbations are small relative to the depth-dependent background pressure. For the pseudo-incompressible and low Mach number approaches this is also sufficient, but for the anelastic approximation all thermodynamic perturbations from a depth-dependent background state must also be small.

Only few studies are available that verify these theoretical predictions quantitatively for thermal convection, with most of them confirming a good agreement between anelastic and fully compressible computations and paying less attention to the pseudo-incompressible approach. Whereas Berkoff et al. (2010) addressed linear magneto-convection, Lecoanet et al. (2014) investigated differences between thermal conductivity and large-eddy entropy diffusion (Glatzmaier, 1984) and additionally considered the pseudo-incompressible case. This was followed by Calkins et al. (2015a) verifying the accuracy of the anelastic approximation at the point of marginal stability under the influence of rotation in Prandtl number unity fluids. However, in a follow-up study Calkins et al. (2015b) showed that the linear anelastic equations fail to produce physically meaningful results for marginally stable, rapidly-rotating, low Prandtl number systems. Against their theoretical expectations, they report that this breakdown of the anelastic approximation is caused by the temporal derivative of the density perturbation becoming an important component in the continuity equation. Furthermore, they suspect related problems in the astrophysically relevant nonlinear turbulence regime. So far, Verhoeven et al. (2015) carried out the only systematic one-to-one comparison for fully compressible and anelastic turbulent convection, which again proved the functionality of the anelastic equations. They, however, did not consider the problematic regime of rapid-rotation and low Prandtl number.

Although there has been no such study that proves the breakdown of the anelastic approximation for the nonlinear case, the findings of Calkins et al. (2015b) are alarming, as they implicitly question the validity of works that investigate astrophysical objects by means of the anelastic equations. There are two possible ways to deal with this problem: Firstly, alternative sound-proof approaches must be studied concerning their applicability at the point of marginal stability. Secondly, the possible breakdown of the anelastic approximation in the astrophysically more relevant turbulence regime must be either confirmed or disproved. Here we extend the work of Calkins et al. (2015b) on rapidly-rotating, low Prandtl number compressible convection and present a study on the accuracy of the linear pseudo-incompressible equations at the point of marginal stability. This approach is promising as it retains parts of the problematic temporal derivative of the perturbational density term and thus may provide a more accurate alternative to the anelastic approximation. Furthermore, we present the first one-to-one comparison between fully compressible and anelastic numerical simulations of rapidly-rotating convection in the fully nonlinear turbulent regime at low Prandtl number.

The following questions are addressed:

• •

How do the anelastic and the pseudo-incompressible approaches vary with respect to their accuracy in predicting the point of marginal stability for rapidly-rotating compressible convection in low Prandtl number fluids?

• •

Can we identify the physical process constraining their accuracy, which causes the problems found by Calkins et al. (2015b)?

• •

Is it safe to apply the anelastic approximation to the astrophysically relevant fully nonlinear turbulent regime of rapidly-rotating compressible convection at low Prandtl number?

The remainder of this paper is organised as follows. In section 2, we start with defining our idealised model and discussing the differences among the anelastic, pseudo-incompressible and fully compressible approaches. Then in section 3 we describe our numerical results concerning the marginal point of stability, supercritical linear convection and the fully nonlinear turbulence regime. Finally, general conclusions are drawn in section 4.

2 Model

In this section the fully compressible, pseudo-incompressible and anelastic linear equations are discussed along with the numerical approach to solve them.

2.1 Governing equations

Linear compressible convection within a Newtonian ideal gas is investigated in a plane layer geometry rotating with angular velocity $\bm{\varOmega}={\varOmega}\hat{\bm z}$. The rotation axis is aligned with the unit vector $\hat{\bm z}$ and antiparallel to the constant gravity $\bm g=-g\hat{\bm z}$. The fluid is characterised by constant specific heat capacities at fixed volume and pressure, $c_{v}$ and $c_{p}$. The dynamic viscosity and thermal conductivity,

[TABLE]

are taken as fixed functions of $z$ and relate to the kinematic quantities $\nu$ and $\kappa$, which are the viscous diffusivity and the thermal diffusivity, respectively.

The linear approximation to the governing equations for fully compressible convection describing the temporal evolution of density $\rho$, fluid velocity $\bm v$ and entropy $s$ within a reference frame rotating at angular velocity $\bm{\varOmega}$ are

[TABLE]

with $t$, $p$ and $T$ specifying time, pressure and temperature, and the indices [math] and $1$ denoting the z-dependent background state and its time and 3-D space dependent perturbation, respectively. ${\varPi}_{ij}=\mu\left(\upartial_{j}v_{i}+\upartial_{i}v_{j}-2/3(\bm\nabla{\bm\cdot}\bm v)\delta_{ij}\right)$ is the viscous stress tensor for a Newtonian fluid and $S$ denotes a heat source or sink. Note that we have neglected centrifugal forces in (2b) that express the conservation of mass, momentum and energy. They can be closed with an equation of state, i.e. the ideal gas law

[TABLE]

and a thermodynamic expression relating entropy, temperature and pressure,

[TABLE]

The background state is chosen to satisfy hydrostatic and thermal equilibrium

[TABLE]

which results in these background terms dropping out of equations (2b) and (2c), respectively. Furthermore, the background temperature gradient

[TABLE]

is assumed constant in $z$ and defined by the sum of the adiabatic temperature gradient $-g/c_{p}$ and the ratio of superadiabatic temperature drop, $\Delta T$, prescribed by the boundary conditions and the domain depth $d$. The background state results from solving equations (2a), (3f,ga), (4f,ga) and (6). Therefore, for a constant thermal conductivity, $S$ vanishes; but, for a constant thermal diffusivity, $S$ is negative, i.e., a heat sink.

2.2 Non-dimensionalisation and parameters

The background state can be non-dimensionalised by using the bottom temperature $T_{r}$, bottom density $\rho_{r}$, bottom pressure $p_{r}=(c_{p}-c_{v})\rho_{r}T_{r}$ and $c_{p}$ for entropy. The superadiabatic temperature difference $\Delta T$ prescribed by the boundary conditions of the system is chosen to be the scale for the temperature perturbations. As temperature, density and entropy perturbations are usually assumed to be closely correlated in low Mach number flows (see e.g. Clayton, 1968), the density and entropy perturbations are scaled correspondingly with $\Delta\rho=\rho_{r}\Delta T/T_{r}$ and $\Delta s=c_{p}\Delta T/T_{r}$. The domain depth $d$ is used as reference length scale. The velocity is non-dimensionalised with a convective free-fall velocity $v_{r}=\sqrt{\Delta\rho gd/\rho_{r}}$ and correspondingly time is scaled with the free-fall time $t_{r}=d/v_{r}=\sqrt{\rho_{r}d/(\Delta\rho g)}$. The pressure perturbation scale $\Delta p=\rho_{r}v_{r}^{2}=\Delta\rho gd$ is inferred from the fact that pressure extracts kinetic energy from the vertical flows to drive horizontal motions. The scales for viscous and thermal diffusivities $\nu_{r}$ and $\kappa_{r}$ are defined as their respective values at the bottom boundary of the domain.

When applying these scales, the following non-dimensional parameters emerge. The dissipation number,

[TABLE]

defines the absolute value of the non-dimensional global adiabatic temperature gradient. With $-(D+\epsilon)$ being the total dimensionless temperature gradient, see (2.3a), the superadiabaticity of the system is given by

[TABLE]

The polytropic index

[TABLE]

with the ratio of specific heats

[TABLE]

being chosen to represent a monatomic ideal gas, is an alternative (dependent) parameter for the superadiabaticity $\epsilon$, with a superadiabatic polytropic index satisfying $n<1.5$. The adiabatic polytropic index for an ideal monatomic ideal gas (with $\epsilon=0$) consistently results in

[TABLE]

The ratio of superadiabatic to adiabatic temperature gradient $\epsilon/D$ often given by 1-D solar models111Christensen-Dalsgaard et al. (1996) for example specify values for $\epsilon/D$ of ${\mathrm{O}}(10^{-4})$ in the bulk of the solar convection zone. follow from the derivatives of equations (2.3a,b).

[TABLE]

Another useful and widely used parameter is the number of density scale heights

[TABLE]

specifying the system’s density contrast $\rho_{\mathrm{bot}}/\rho_{\mathrm{top}}$ between bottom and top boundary that can be used instead of the Dissipation number $D$. The Prandtl number,

[TABLE]

is the ratio of momentum to thermal diffusivity. The Rayleigh number,

[TABLE]

controls the vigour of convection with large $g$, $d$, and $\epsilon=\Delta T/T_{r}$ enhancing and large diffusivities $\nu$ and $\kappa$ reducing the convective vigour. More formally $Ra$ is the ratio of the product of the diffusive timescales $d^{2}/\kappa_{r}~{}d^{2}/\nu_{r}$ to the square of the free-fall timescale $t_{r}^{2}$. The Ekman number,

[TABLE]

is the ratio of the rotation time scale to the viscous diffusion time scale.

Note that all non-dimensional parameters given in this subsection are defined at the bottom boundary and may vary strongly over the whole domain.

2.3 Non-dimensional equations

The background state results in

[TABLE]

and the governing non-dimensional equations read

[TABLE]

[TABLE]

with the hat $\hat{}$ denoting non-dimensional quantities, e.g. $T_{0}=T_{r}\hat{T}_{0}$.

2.4 Differences in the anelastic, pseudo-incompressible and fully compressible equations

The density perturbation term in the continuity equation (2.3a) has been expressed in terms of pressure and entropy by using (2.3d) in order to illustrate the approximations carried out in the anelastic and pseudo-incompressible equations in the following.

The anelastic approximation neglects the time derivative of density $\epsilon\hat{\upartial}_{t}\hat{\rho}_{1}=\epsilon D/(\gamma-1)/\hat{T}_{0}\hat{\upartial}_{t}\hat{p}_{1}-\epsilon\hat{\rho}_{0}\hat{\upartial}_{t}\hat{s}_{1}$ in the continuity equation (2.3a), which is reasonable as long as all perturbations remain small with $\epsilon\ll 1$. The pseudo-incompressible approximation is less restrictive and only neglects the time derivative of the perturbational pressure term. This approach is justified for low Mach number flows but allows for temperature, density and entropy perturbations to be large222Note firstly that the prefactor in the perturbational pressure term equals the squared Mach number (see Appendix B) and secondly that $-\rho_{0}/c_{p}\upartial_{t}s_{1}$, which is the dimensional form of the $-\epsilon\hat{\rho}_{0}\hat{\upartial}_{t}\hat{s}_{1}$ term in the continuity equation (2.3a), equals the temporal derivative of the pseudo-incompressible density $\upartial_{t}\rho^{*}$ as first defined in Durran (1989).. In order to check for the validity of both sound-proof approaches the relative magnitudes of the terms being neglected in each approximation, i.e. the perturbational density and pressure terms in the continuity equation

[TABLE]

will be quantified for various parameters in this paper. We have chosen to focus on the values at the top boundary as the sound speed is the lowest at this location.

For simplicity we will assume $\epsilon\rightarrow 0\Longleftrightarrow n\rightarrow 1.5$ for each anelastic simulation carried out making the $\epsilon$ parameter obsolete. This yields the same equations as when assuming a small but finite $\epsilon$ and neglecting all terms involving $\epsilon$ in (2.3a-d) and (2.3a) resulting in an adiabatic background state, which is perturbed by a constant superadiabatic temperature drop prescribed by the boundary conditions. Similar approaches are common practice when applying the anelastic approximation for nonlinear turbulent convection simulations (see, e.g. Gastine et al., 2014; Heimpel et al., 2015) and do not increase the error of the anelastic approximation (Lantz and Fan, 1999). Concerning the pseudo-incompressible approach, however, the parameter $\epsilon$ is required, as it still appears in the continuity equation (2.3a) in the entropy term.

2.5 Numerical approach for the linear stability problem

In order to solve the linear equations (2.3a-e) numerically, each variable is represented by the typical normal mode ansatz, e.g., \hat{T}_{1}=\hat{T}(\hat{z})\exp\bigl{[}\hat{r}\hat{t}+{\mathrm{i}}\bigl{(}\hat{\omega}\hat{t}+\hat{k}_{x}\hat{x}+\hat{k}_{y}\hat{y}\bigr{)}\bigr{]}. Here, $\hat{r}$ denotes the growth rate, $\hat{\omega}$ is the oscillation frequency and $\hat{k}_{x}$ and $\hat{k}_{y}$ are the horizontal wavenumbers with $\hat{k}=\sqrt{\hat{k}_{x}^{2}+\hat{k}_{y}^{2}}$ (not to be confused with the thermal conductivity $\hat{k}_{t}$). The equations, as they are solved numerically, are given in Appendix E, with the critical frequency $\hat{\omega}_{c}$ and the critical Rayleigh number $Ra_{c}$ being eigenvalues that depend on the critical wavenumber $\hat{k}_{c}$, which is determined by using a nested intervals scheme for $\hat{r}=0$. Growth rates $\hat{r}$ can be determined by prescribing the Rayleigh number $Ra$ and the wavenumber $\hat{k}$. A positive $\hat{r}$ with a non-zero $\hat{\omega}$ means that the onset of convection is a temporal oscillation of all variables with their mean amplitudes increasing exponentially in time, whereas the instability is non-oscillatory for $\hat{\omega}=0$. As the solutions corresponding to critical modes $\hat{k}_{c}$ are independent of the individual choice of $\hat{k}_{x}$ and $\hat{k}_{y}$ and manifest themselves in two-dimensional convection rolls, we restrict our linear analysis to the $k_{y}=0$ modes for simplicity without losing generality (see Calkins et al., 2015b).

The numerical framework is based on a second order accurate Newton-Raphson-Kantorovich (NRK) method for solving eigenvalue problems. While the computational domain is periodic in horizontal direction, stress free, impermeable and fixed temperature (i.e. $\hat{T}_{1}=0$) top and bottom boundaries are applied.

The numerical approaches for solving the nonlinear fully compressible and anelastic equations given in Appendix A are described in Verhoeven and Stellmach (2014); Verhoeven et al. (2015).

3 Results

In this section results from a suite of anelastic, pseudo-incompressible and fully compressible computations are presented in order to test their accuracy and limitations. All linear computations shown are carried out at fixed Ekman number $Ek=10^{-6}$ and reside in the rapidly-rotating geostrophic limit at low Rossby number, i.e. the Coriolis forces are mostly balanced by the pressure gradient and therefore the results will be essentially the same for computations with smaller Ekman numbers. This was shown by Calkins et al. (2015b) for the marginal stability cases and will be checked for the supercritical simulations presented in the following. We start with studying the point of marginally stable convection, which is followed by investigating the onset of convection in a supercritical setup and eventually finish with the fully nonlinear turbulence regime.

3.1 Marginal stability

Although linear marginal modes at the onset of convection are not relevant to turbulent convection in stars and planets, we start with the investigation of this topic. This will pave the way for understanding possible problems in sound-proof approaches for supercritical cases. Accordingly, some of the parameter values considered in this subsection do not correspond to an astrophysically achievable system.

Figure 1 shows plots of critical wavenumbers $\hat{k}_{c}=\sqrt{\hat{k}_{x}^{2}+\hat{k}_{y}^{2}}$, frequencies $\hat{\omega}_{c}$ and Rayleigh numbers $Ra_{c}$ against the number of density scale heights $0\leq N_{\rho}\leq 5$ for all approaches with various polytropic indices $1\leq n\leq 1.5$ at low Prandtl numbers $0.01\leq Pr\leq 0.1$. The results from Calkins et al. (2015b) were reproduced for the anelastic and the fully compressible equations for these constant dynamic viscosity $\mu$ and thermal conductivity $k_{t}$ cases.333Compare figures (1a,b,c,d) to Calkins et al. (2015b) figures (3a,c,b,d), respectively. Note that they plot slightly different values than we do. While they plot $\hat{k}_{c}Ek^{1/3}$, we plot $\hat{k}_{c}$, which effectively makes our values $100$ times larger for the $Ek=10^{-6}$ cases shown. Moreover they scale their critical frequency $\hat{\omega}_{c}^{\mathrm{Calkins}}$ with a timescale inferred from the free-fall velocity based on values at the top boundary, see their equation (2.6). We, however, use reference values at the bottom boundary. The relation between our critical frequency and theirs is given by \hat{\omega}_{c}=\exp\bigl{[}(N_{\rho})^{1/(2n)}\bigr{]}\hat{\omega}_{c}^{\mathrm{Calkins}}. See also their supplementary material for critical Rayleigh numbers. In contrast to the anelastic approximation that produces infinitely small $Ra_{c}$ and $\hat{k}_{c}$ at finite $N_{\rho}$, which Calkins et al. (2015b) considered ”unphysical”, the pseudo-incompressible approach never fails and always produces physically meaningful results with $Ra_{c}>0$ and $\hat{k}_{c}>0$ in our computations. The results obtained with the pseudo-incompressible equations, however, deviate from those of the fully compressible equations when the relative magnitude of the perturbational pressure term $p_{\mathrm{cont}}$ (see (19b)) gets larger than ${\mathrm{O}}(0.1)$. These transitions are marked with circles in figure 1 and generally correspond to lower $N_{\rho}$ the lower the Prandtl number. These results suggest that the anelastic and the pseudo-incompressible approximation both should be applied carefully in marginally stable, rapidly-rotating, low Prandtl number convection systems.

Note that in our simplified model the failure of the anelastic approximation generally involves a drastic decrease in critical wavenumbers $k_{c}$ corresponding to very large aspect ratios ${\varGamma}$, which are given by the horizontal wavenumber $k_{c}$ and the depth of the domain $d$ with ${\varGamma}=2\pi/(k_{c}d)=2\pi/\hat{k}_{c}$. Instead of solving for ${\varGamma}$, as we do here, one may choose to model marginally stable convection in, for example, a particular spherical shell geometry by prescribing the radii of both the inner and outer boundaries, which sets the value of ${\varGamma}$ and so places a lower limit on the possible values of the horizontal wavenumber $\hat{k}_{c}$. How this affects the validity of the anelastic approximation for rapidly-rotating, low Prandtl number, marginally stable convection is not clear. Two previous studies of this linear stability problem for a 3-D spherical shell Glatzmaier and Gilman (1981a, b); Drew et al. (1995) agree for cases with $Pr=1$ and $10$ but disagree for $Pr=0.1$. It is worth noting that many anelastic studies of convection in stars and giant planets have modeled the transport of heat by turbulent eddies as a diffusion process based on a turbulent thermal diffusivity and the gradient of entropy, instead of the gradient of temperature. The argument (Glatzmaier, 1984) is that the heat transport by eddies, that are too small to numerically resolve in global models but contribute significantly more heat transport than that by radiation or conduction, should be proportional to the local entropy gradient since turbulence tends to maintain a constant entropy. Jones et al. (2009) has shown that using entropy diffusive heat flux in an anelastic linear stability model of rotating marginally stable convection cannot have negative critical Rayleigh numbers as Drew et al. (1995) found for temperature diffusive heat flux at $Pr=0.1$.

The constant $\mu$ and $k_{t}$ setup and the corresponding parameter range investigated so far, which has been adopted from Calkins et al. (2015b), involves some problems that might explain the observed breakdown of the anelastic approximation and the inaccuracies in the pseudo-incompressible approach: First, the diffusivities $\nu$ and $\kappa$ vary strongly with height for large $N_{\rho}$ (see (2.1a,b)) partly resulting in diffusion velocities close to the sound speed. Second, the phase velocity $\hat{v}_{\mathrm{phase}}=\hat{\omega}/\hat{k}$, which is the velocity of the pattern of the dominant perturbation corresponding to the oscillatory instability through the domain, exceeds the sound speed for specific parameters. Third, for other parameters the typical rotation time $t_{\mathrm{rot}}=1/(2{\varOmega})$ falls below the sound-crossing time of the the domain $t_{\mathrm{sound}}=d/v_{\mathrm{sound}}$. All of these issues potentially introduce timescales shorter than the free-fall time, which, although $\epsilon$ is small, hinders the justification of the approximations carried out in sound-proof approaches.

In order to exclude the possibility that the accuracy problems of the anelastic and the pseudo-incompressible approach are an artefact due to the very high diffusivities near the top boundary, figure 2 shows the results from a second series of computations using the same parameters as before but for constant diffusivities $\nu$ and $\kappa$ throughout the domain preventing larger diffusion velocities near the top boundary. Although the plots show individual differences to figure 1, the overall result is similar. Alike the anelastic equations that fail to produce physically meaningful outcomes at similar $N_{\rho}$ as in figure 1, the pseudo-incompressible approach deviates from the fully compressible method at nearby $N_{\rho}$ compared to the constant $\mu$ and $k_{t}$ case.

Interestingly, our fully compressible computation with constant diffusivities, $Pr=0.01$ and $n=1$ yields strongly decreasing critical wavenumbers with $\hat{k}_{c}\ll 1$ and Rayleigh numbers of ${\mathrm{O}}(10^{3})$ as $N_{\rho}$ approaches a value of approximately $3.2$. We speculate that this rather surprising result is caused by the influence of the heat sink that is necessary in order to maintain the background state, see (4f,gb). The investigation of this topic is beyond the scope of this paper and left for future studies.

According to all our results the anelastic and the pseudo-incompressible approximation both seem to lose accuracy when the perturbational density term and the perturbational pressure term, respectively, have magnitudes that cannot be neglected in the continuity equation anymore, see (19a-19b). In sound-proof approaches it is typically assumed that the magnitude of the neglected perturbational terms relative to the magnitude of the other terms in the continuity equation scale with the square of the Mach number $M^{2}$, with $M$ being defined as the ratio of a typical fluid velocity $v_{r}$ and the speed of sound $v_{\mathrm{sound}}=\sqrt{c_{p}(c_{p}-c_{v})T/c_{v}}$, see Appendix B. For linear stability calculations the mean amplitudes of the fluid velocity and thermodynamic perturbations are arbitrary and increase exponentially with time while the background state does not change. Therefore, the regular Mach number $M$ is not a useful diagnostic for linear calculations. However, the phase velocity $\hat{v}_{\mathrm{phase}}=\hat{\omega}/\hat{k}$ is independent of time in linear computations. As the sound speed decreases for smaller temperatures, the phase Mach number is the largest at the top boundary where the fluid is cold. Such a phase top Mach number is derived in the non-dimensionalisation used here in Appendix C and reads

[TABLE]

Additionally, in order for sound-proof approaches to be valid the rotation time scale must be larger than the sound-crossing time of the domain (Braginsky and Roberts, 1995). A corresponding rotational Mach number can be defined as in appendix D, resulting in

[TABLE]

Please note that for the linear marginal stability case we are considering here, the Rayleigh number in (21) must be consistently replaced by $Ra_{c}$, which is calculated and not a prescribed parameter as it would be for the supercritical onset of convection or a nonlinear simulation. As both, the phase and the rotational Mach numbers, depend on the superadiabaticity $\epsilon$, no finite values can be determined for our anelastic computations with $\epsilon=0$ directly. Furthermore, in the anelastic limit, $\epsilon\rightarrow 0$, the critical Rayleigh number can become infinitely small, $Ra_{c}\rightarrow 0$, prohibiting the calculation of $M_{\mathrm{rot}}$. This can be compensated by using fully compressible computations close to the anelastic limit with $\epsilon\ll 1\Longleftrightarrow 1.5-n\ll 1$. For $Ek=10^{-6}$, $Pr=0.01$ and $Pr=0.1$, and constant $\mu$ and $k_{t}$ and constant diffusivities $\nu$ and $\kappa$, figure 3 shows plots of $N_{\rho}$, for which the afore mentioned Mach numbers exceed a specific limit, against $1.5-n$ being proportional to $\epsilon$ for $\epsilon\ll 1$. The fully compressible simulations are closest to the anelastic limit for small $1.5-n$, for which both Mach number cases show convergence against a constant value of $N_{\rho}$. Concerning the rotational Mach number case this plateau is caused by the constancy of $\epsilon/Ra_{c}$ in this parameter range; see figure 1. Figure 3 shows that the breakdown of the anelastic approximation (marked by the black line) coincides with the $N_{\rho}$ corresponding to the rotational Mach number $M_{\mathrm{rot}}$ exceeding $1$ (red). For comparison the values of $N_{\rho}$ correlating with the phase Mach number $M_{\mathrm{phase}}=0.1$ are displayed in blue. The limit of $0.1$ was chosen, as the limit of $1$ was not exceeded for the computations with small $1.5-n$.

In summary, the accuracy of the anelastic approximation for marginally stable convection is controlled by the rotational Mach number $M_{\mathrm{rot}}$ in our simulations. Our results suggest that Calkins et al. (2015b) found the anelastic equations to fail, as the typical rotation time was smaller than the sound-crossing time of the domain in their respective simulations. We speculate that the pseudo-incompressible approximation has a more gentle changeover to being imprecise than the anelastic approximation to failure because the temporal derivative term in the continuity equations compensates for the missing pressure term and allows for numerically stable but nevertheless inaccurate solutions.

3.2 Supercritical onset of convection

So far, we focused on the marginal point of stability only. Natural systems, however, are often characterised by strongly supercritical Rayleigh numbers, $Ra\gg Ra_{c}$, which in combination with small superadiabaticity leads to small rotational Mach numbers, see (21). As the supercriticality increases, non-oscillatory (i.e., $\hat{\omega}=0$) convective instabilities emerge. The critical wavenumbers $\hat{k}_{c}^{\mathrm{non-osc}}$ and Rayleigh numbers $Ra_{c}^{\mathrm{non-osc}}$ corresponding to this kind of instability do not depend on the Prandtl number and are given in figure 4 for the anelastic ($n=1.5$) and the fully compressible case ($n=1.49$) with $Ek=10^{-6}$ and constant $\mu$ and $k_{t}$. In contrast to the oscillatory instability, both cases closely match for all $N_{\rho}$ considered.

The relative magnitudes of the time derivative term in the continuity equation at the top boundary $\rho_{\mathrm{cont}}$, see (19a), is plotted against the Rayleigh number for oscillatory and non-oscillatory fully compressible linear convection in figure 5 for $n=1.49$, $Pr=0.1$ and different Ekman numbers and numbers of density scale heights. The outcome of figure 5 is twofold: First, the results shown do not vary for $Ek\leq{\mathrm{O}}(10^{-6})$ for $N_{\rho}=2$ in the left panel, and for $Ek\leq{\mathrm{O}}(10^{-7})$ for $N_{\rho}=5$ in the right panel as the geostrophic limit is reached for all corresponding Rayleigh numbers displayed. Second, $\rho_{\mathrm{cont}}$ decreases for increasing $Ra$, with $\rho_{\mathrm{cont}}\ll 1$ for $Ra>Ra_{c}^{\mathrm{non-osc}}$. This strongly suggests the functionality of the anelastic approximation for large Rayleigh numbers for rapidly-rotating convection, although the Prandtl number is low, the density contrast is large and the geostrophic limit is reached.

In order to further ascertain this, the growth rates $\hat{r}$ corresponding to both non-oscillatory and oscillatory instabilities for anelastic and fully compressible convection are displayed in figure 6 for $Ek=10^{-6}$ and $Pr=0.1$. For given Rayleigh numbers $Ra$, numbers of density scale heights $N_{\rho}$ and polytropic indices $n$, the growth rates $\hat{r}$ are plotted against the wavenumber $\hat{k}$ . The anelastic approximation neither fails nor gives rise to spurious growth for any parameters investigated in the supercritical regime. In fact fully compressible growth rates generally converge to the anelastic limit case for $n\rightarrow 1.5$. While $\hat{r}$ strongly varies between different polytropic indices for large $N_{\rho}$ and low $Ra$, the match increases for lower $N_{\rho}$ and larger $Ra$, with the latter seeming to be the most relevant parameter. For moderately high $Ra=2Ra_{c}^{\mathrm{non-osc}}$ the growth rates of non-oscillatory and oscillatory instabilities are indistinguishable for $n\geq 1.49$ and just differ to the ones for $n=1.2$ within a few percent.

In summary, the problems of the anelastic approximation found by Calkins et al. (2015b) only happen close to marginal stability and then only for rotational Mach numbers of order one or larger. They do not occur for Rayleigh numbers more than twice the critical Rayleigh number corresponding to the fully compressible framework. Given the very high Rayleigh numbers required for nonlinear turbulent convection, this strongly suggests the insignificance of this phenomenon for the turbulence regime relevant to planets and stars.

3.3 Nonlinear turbulent convection

In order to test the functionality of the anelastic approximation in rapidly-rotating turbulent convection in low Prandtl number fluids, a one-to-one comparison with constant $\mu$ and $k_{t}$ is carried out in the parameter range where Calkins et al. (2015b) suspect the breakdown of the anelastic equations.

The employed parameters for the fully nonlinear anelastic simulation are $Pr=0.1$, $Ek=10^{-6}$, $N_{\rho}=2$ and the Rayleigh number is defined by two times the critical value for non-oscillatory anelastic convection given in figure 4, which yields $Ra\approx 1.98\times 10^{9}$. The choice of these parameters ensures a major role of Coriolis forces with a convective Rossby number (Gilman, 1977) of $Ro_{c}=Ek\sqrt{Ra/Pr}\approx 0.14$ and correspond to the geostrophic limit according to the linear results shown in figure 5. For the corresponding fully compressible computation we have to additionally provide the superadiabaticity, which is chosen to be $\epsilon=0.1$. This choice of the superadiabaticity is a trade-off between a low value, which is predicted for planetary and stellar convection zones, to computational feasibility. We further keep the adiabatic number of density scale heights $N_{\rho}^{ad}=2$, which yields $N_{\rho}\approx 2.17$ and $n\approx 1.20$. The aspect ratio ${\varGamma}=l_{x}/d=l_{y}/d$, i.e. the horizontal width divided by the domain height, is inspired by the critical wavenumber associated with the non-oscillatory instability. The choice of ${\varGamma}\approx 0.39$ ensures that several wavelengths of the critical non-oscillatory ($\sim 6$) and oscillatory instabilities ($\sim 2$) fit into the domain in each horizontal direction, see figures 1 and 4.

The anelastic and the fully compressible simulations turn out to be very similar.

Figure 7 shows the time evolution of the root-mean-square velocity

[TABLE]

with $\hat{V}$ being the volume of the domain, for the fully compressible and the anelastic simulation in red and blue, respectively. Panel (a) displays the exponential growths of $v_{\mathrm{rms}}$, which match the theoretical predictions for the fully compressible (orange) and anelastic (green) non-oscillatory instabilities. Panel (b) shows the long-time evolution of $v_{\mathrm{rms}}$ being qualitatively the same for both cases. After the first five free-fall times characterised by exponential growth, one large-scale cyclonic vortex emerges in each simulation, which accumulates kinetic energy until about $\hat{t}=200$ when the computations reach statistical equilibrium. The time-averaged fully compressible value for $v_{\mathrm{rms}}$ matches the corresponding anelastic one within $2\%$ in this final state. Also the net heat transport in terms of the Nusselt number

[TABLE]

where the angle brackets imply temporal and horizontal averaging, turns out to be comparable in the final state, with $Nu=13.2$ for the anelastic and $Nu=12.2$ for the fully compressible case representing a difference of roughly $8\%$. This conformity in the output parameters agrees with direct comparisons of anelastic and fully compressible non-rotating convection (Verhoeven et al., 2015) and shows (albeit indirectly) that the time derivative of the density fluctuation in the continuity equation is not important in our nonlinear simulations.

In order to give a visual impression of the fluid dynamics, figure 8 shows snapshots of the anelastic simulation in statistical equilibrium displaying (a) the vertical vorticity

[TABLE]

and (b) the superadiabatic temperature $\hat{T}_{1}$. Panel (a) clearly shows the large-scale cyclonic vortex (red), whose angular momentum is balanced by many small-scale anticyclones (blue). In panel (b) cold material (blue) is located near the top and hot fluid (red) close to the bottom boundary.

Such symmetry breaking single large-scale cyclonic vortices typically emerge on the verge of the regime of nonlinear rapid rotation (see, e.g., Guervilly et al., 2014; Favier et al., 2014). Ekman numbers of ${\mathrm{O}}(10^{-7})$ and below are necessary in order to reach the nonlinear geostrophic turbulence regime, in which strong cyclonic and anticyclonic vortices are generated with the symmetry being restored (Stellmach et al., 2014; Rubio et al., 2014). Such parameter ranges are clearly desirable, albeit not accessible with our fully compressible code, which prevents us from carrying out a direct comparison with anelastic results for these extreme parameters. A broader parameter range needs to be tested in future studies.

4 Conclusions

The validity of the anelastic approximation in the astrophysically relevant regime of rapidly-rotating compressible convection has recently been questioned by Calkins et al. (2015b). The high computational efficiency and broad application of sound-proof approaches provided the motivation for this paper, that extends their work by further constraining and reviewing the applicability of different sound-proof models in this regime.

As a starting point we focused on the pseudo-incompressible approximation (Durran, 1989), which comprises more discreet simplifications in comparison to the anelastic approach. This is indeed reflected in our computations and in contrast to the anelastic approximation, the pseudo-incompressible approach does not fail in linear stability calculations at the point of marginal stability. Instead its results slowly deviate from the fully compressible solution as the perturbational pressure term in the continuity equation becomes more and more non-negligible, which is usually the case for an increasing number of density scale heights and decreasing Prandtl number shortly after the anelastic approximation breaks down. Our results confirm what Calkins et al. (2015b) conclude about anelastic models of marginally stable convection; i.e., that the anelastic approximation for compressible convection in the rapidly-rotating low Prandtl number regime can be inaccurate at marginal stability. We find that the anelastic approach breaks down for simulations, in which the sound-crossing time of the computational domain exceeds the rotation time scale. Correspondingly, a rotational Mach number is defined as the ratio of both values, which in analogy to the classical Mach number needs to be small in order for sound-proof approaches to be valid.

As the rotational Mach number is inversely proportional the square-root of the Rayleigh number, the situation is much different for higher supercriticality, where the anelastic approximation neither fails nor gives rise to spurious growth for any parameters investigated. Instead our computed fully compressible growth rates generally converge against the ones of the anelastic limit case as the superadiabaticity decreases. In this regime, we found the Rayleigh number to be the most constraining parameter in the sense that low values—close to marginal stability—show distinct differences between anelastic and fully compressible results, whereas moderately high values—of the order of the critical value of non-oscillatory convection—were sufficient to yield a match within a few percent for moderately low superadiabaticities. As Rayleigh numbers in planetary and stellar convection systems are often much higher and superadiabaticities lower than this, we expect the problems causing the breakdown of the anelastic approximation at marginal stability to be insignificant in such systems and predict the agreement between the approximated and the full equations to get even closer.

These findings could be confirmed by fully nonlinear turbulent convection simulations. The computations in the regime where Calkins et al. (2015b) suspected the breakdown of the anelastic approximations, showed close qualitative and quantitative agreement between the anelastic and the fully compressible case. However, further work needs to be done in order to exclude effects other than the one suggested by Calkins et al. (2015b), which may cause problems in the anelastic approximation in the geostrophic turbulence regime.

To sum up, the anelastic approximation breaks down for rapidly-rotating compressible convection in low Prandtl number fluids at the point of marginal stability when the rotational Mach number is greater than one. According to our results, however, these problems disappear in the astrophysically more relevant regime of turbulence typically characterised by small rotational Mach numbers. Although we did not specifically test the pseudo-incompressible approximation for supercritical convection, we do not see any reason why the same should not be true for this less invasive approach.

Acknowledgments

We would like to thank two anonymous reviewers for comments that greatly improved the manuscript. We gratefully acknowledge that Prof. Pascale Garaud from UCSC provided her NRK routines, which were the basis of our linear convection code. Support was provided for J.V. by the National Aeronautics and Space Administration under grants OPR NNX13AK94G and PGG NNX14AN70G. The simulations were run at UCSC on the Hyades supercomputer obtained with an MRI grant from the National Science Foundation.

Appendix A Nonlinear equations of compressible convection

Nonlinear equations for fully compressible and anelastic convection for the constant $\mu$ and $k_{t}$ case as given by Verhoeven et al. (2015) with additionally considering rotation are displayed below. The background profiles $\hat{\rho}_{0}$, $\hat{T}_{0}$ and $\hat{p}_{0}$ are assumed to be adiabatic.

A.1 Fully compressible equations

The non-dimensional fully compressible equations read

[TABLE]

with \hat{e}_{ij}=\frac{1}{2}\bigl{(}\hat{\upartial}_{j}\hat{v}_{i}+\hat{\upartial}_{i}\hat{v}_{j}\bigr{)} being the strain rate tensor.

A.2 Anelastic equations

The non-dimensional anelastic equations result from (25a-d) in the limit of $\epsilon\rightarrow 0$:

[TABLE]

Appendix B Mach number

The Mach number $M$ is defined as the ratio of a reference velocity, which we assume to be the free-fall velocity

[TABLE]

to the sound speed

[TABLE]

resulting in

[TABLE]

Exploiting $\hat{\rho}_{0}^{\mathrm{top}}=\exp{\left(-N_{\rho}\right)}$ yields the Mach number at the top boundary,

[TABLE]

where the sound speed is the lowest and thus $M$ has a maximum.

Appendix C Phase Mach number

The phase Mach number $M_{\mathrm{phase}}$ will be defined as the ratio of the phase velocity corresponding to the oscillatory instability

[TABLE]

to the sound speed

[TABLE]

resulting in

[TABLE]

Exploiting $\hat{\rho}_{0}^{\mathrm{top}}=\exp{\left(-N_{\rho}\right)}$ and equation (30) yields the phase Mach number at the top boundary,

[TABLE]

where the sound speed is the lowest and thus $M_{\mathrm{phase}}$ has a maximum.

Appendix D Rotational Mach number

The rotational Mach number $M_{\mathrm{rot}}$ will be defined as the ratio of the velocity corresponding to the domain depth $d$ and the angular frequency ${\varOmega}$, where

[TABLE]

to the sound speed at the bottom

[TABLE]

resulting in

[TABLE]

This is equivalent to defining the rotational Mach number as the ratio of the sound-crossing time of the domain to a rotation time scale,

[TABLE]

based on the sound speed $v^{\mathrm{bot}}_{\mathrm{sound}}$ at the bottom boundary. Note that the value of $t^{\mathrm{bot}}_{\mathrm{sound}}$ is similar to the one of the real sound-crossing time $t_{\mathrm{sound}}=\int_{0}^{d}{\mathrm{d}}z/v_{\mathrm{sound}}$.

Appendix E Numerically solved linear equations

This section summarises the equations as they are solved numerically in the linear convection code, which is based on a Newton-Raphson-Kantorovich (NRK) method. As it only deals with non-dimensional quantities, the hat $\hat{}$ is left out for clarity. Equations (2.3a-e) can be simplified by expressing $\rho_{1}$ and $s_{1}$ in terms of ${T}_{1}$ and ${p}_{1}$, which reduces the number of equations from five to three,

[TABLE]

The parameters $\alpha$ and $\beta$ have been introduced in order to account for the three different approaches under investigation. While in the anelastic approximation the perturbational density term is neglected in the continuity equation (i.e., $\alpha=\beta=0$), the pseudo-incompressible approximation just neglects the perturbational pressure term resulting in $\alpha=0$ and $\beta=1$. The fully compressible equations do not contain any simplifications, thus $\alpha=\beta=1$. We further define $P_{1}=\sqrt{Ra/Pr}\,p_{1}$ in order to simplify the numerical method, compare equation (39b) with (42f), (43f).

For the constant dynamic viscosity $\mu$ and thermal conductivity $k_{t}$ case we assume

[TABLE]

whereas for constant diffusivities $\nu$ and $\kappa$ the corresponding dynamic viscosity and thermal conductivity vary with depth (see (2.1a,b) and (2.3b)):

[TABLE]

Note that the ′ sign denotes first derivatives with respect to $z$, i.e., the numerically independent variable $\mu^{\prime}(z)$ is the analytical derivative of $\mu(z)$.

When using the typical normal mode ansatz, e.g., ${T}_{1}={T}({z})\exp\left[{r}{t}+{\mathrm{i}}\left({\omega}{t}+{k}_{x}{x}+{k}_{y}{y}\right)\right]$, with the wavenumber $k_{x}$ and $k_{y}$ not to be confused with the thermal conductivity $k_{t}$, equations (39a-c) can be solved separately for the real part (index re),

[TABLE]

and the imaginary part (index im)

[TABLE]

The equations for the eigenvalues

[TABLE]

guarantee that they are constants, independent of $z$.