© S. Daley-Yates et al. 2021

1 Introduction

Natural convection is a fundamental process governing behaviour at vastly differing length scales, from water boiling in a kettle to the interiors of stars (Busse 1978; Hurlburt et al. 1984; Bodenschatz et al. 2000). Convection leads to the dissipation of energy, the mixing of fluids, and the transfer of heat, all vital processes in astrophysics but especially so in the context of stellar interiors (Denissenkov 2010; Brown et al. 2013; Garaud et al. 2015) and (exo-)planetary and brown dwarf atmospheres (e.g. Saumon & Marley 2008; Morley et al. 2012; Tremblin et al. 2015, 2017).

Assessing the impact of thermal and compositional source terms (i.e. diabatic processes) has been a challenge, one that was recently tackled in Tremblin et al. (2019) as part of a general theory for thermo-compositional diabatic convection. In the above study, this theory was demonstrated to be applicable to thermohaline convection in Earth’s oceans (Stern 1960; Turner 1967; Gregg 1988), moist convection in Earth’s atmosphere (Arakawa & Jung 2011; Yano 2014), and, potentially, CO and CH₄ radiative convection in the atmospheres of giant exoplanets and brown dwarfs (Tremblin et al. 2015, 2016, 2017, 2019). The purpose of our study is to further establish the theoretical standing of thermo-compositional diabatic convection, proposed by Tremblin et al. (2019), by investigating the impact of dimensionality, resolution, and long timescales on convective numerical simulations forced by thermo-compositional diabatic processes. As an idealised test case, we chose a potential CO+O ↔ CO₂ transition in the secondary atmosphere of terrestrial exoplanets.

The composition and dynamics of rocky exoplanet atmosphere are currently receiving considerable attention (Ito et al. 2015; Madhusudhan 2019; Modirrousta-Galian et al. 2021). This is partially due to the imminent arrival of several next-generation telescopes, both space and ground based, which have the characterisation of exoplanet atmospheres in their remit. This is because of questions regarding the habitability and the imprint of the formation and evolutionary history of the planet (Stevenson et al. 2016; Tinetti et al. 2018).

The habitability of a planet’s atmosphere depends directly on the separate stages of its evolution; in the cases of Earth and Mars, for example, the early atmospheres were characterised by gaseous outflows from the forming planetary crust Elkins-Tanton (2008). In this type of environment there exists strong chemical gradients and high temperatures, due either to irradiation or to early evolutionary stages. This means that chemical reactions are abundant, for example the CO+O ↔ CO₂ transition. Such a transition is likely to have occurred early during the formation of the Earth’s atmosphere as CO₂ and CO are abundant via outgassing Forget & Leconte (2014). Understanding this mechanism is therefore not only interesting in the astrophysical context but also from the paleo-climatology prospective.

For hot secondary atmospheres of young or highly irradiated rocky exoplanets, we can assume that the vertical temperature gradient is negative (at least on the night side of tidally locked irradiated planets). This negative vertical temperature gradient, together with temperature-dependent chemical reactions, generally ensures that higher mean molecular weight species form above lower mean molecular weight species. In short, thermally driven chemistry often leads to heavier molecules above lighter molecules and, therefore, a positive vertical mean-molecular-weight gradient. In the context of the present study, this results in a CO to CO₂ transition from the deep to the upper atmosphere. In this physical situation, the atmosphere is unstable to Ledoux convection (otherwise known as natural convection). The denser CO₂ -dominated gas will fall and displace the lighter, less dense CO-dominated gas, which will in turn rise and replace the CO₂, resulting in an atmosphere inversion. This trivial situation would then remain, with heavier material below lighter material. However, as the mean molecular weight inversion proceeds, the chemical reaction mentioned above acts to replenish both the CO₂ at the top and the CO at the bottom of the atmosphere, restoring the unstable gradient. The two processes, inversion and chemical reaction, compete until an equilibrium is achieved and a continuous convective cycle is established.

Using this CO+O ↔ CO₂ transition as a test case, we have explored a range of atmosphere temperatures at low resolution in both 2D and 3D and one specific temperature at high resolution. We investigate both the susceptibility to convection at different CO-CO₂ mixtures and at different forcing timescales, as well as the system’s response to small-scale flow features, which are accessible at high resolution. In these idealised setups, we show that a temperature-gradient reduction can emerge when the compositional forcing timescale is sufficiently fast, potentially leading to the equivalent of the reddening seen in brown dwarf atmospheres. For more realistic situations, further detailed 1D studies are needed to verify that the compositional and thermal gradients, forcing timescales, and a realistic equation of state (EOS) will indeed lead to such a reduced gradient. Despite these additional checks, the possibility cannot be excluded a priori and has observational consequences that could be tested with, for example, the James Webb Space Telescope (JWST).

In the following sections we first describe the numerical setup used to study the relevant convective regimes. We then describe and discuss the simulation results and its limitations, before reaching our conclusions.

2 Numerical setup

In the following sections we follow the analytic modelling proposed by Tremblin et al. (2019), giving justification for the numerical regime and initial conditions that the simulations are conducted in.

2.1 Governing equations

The equations of hydrodynamics can be represented as the following hyperbolic partial differential equation: $${\partial }_{t}U+\nabla F=S,$$(1)

where U, F, and S are the state vector of conserved variables, flux vector, and vector of source terms, respectively, and are defined as: $$U=\left(\begin{array}{c}\rho \\ \rho u\\ \rho X\\ \rho \epsilon \end{array}\right),F=\left(\begin{array}{c}\rho \\ \rho u\otimes u+P\\ \rho Xu\\ u\left(\rho \epsilon +P\right)\end{array}\right),S=\left(\begin{array}{c}0\\ \rho \text{ pmbtextsl}g\\ \rho R(X,P,T)\\ \rho c_{p}H(X,P,T)\end{array}\right),$$(2)

where ρ, u, X, P, ε, c_p, and T are the density, velocity, mass-mixing ratio, pressure, total energy, and specific heat capacity at constant pressure and temperature, respectively. This system is closed using the ideal gas law, $$P=\frac{\rho k_{\text{B}}T}{\mu (X)m_{\text{H}}},$$(3)

where m_H is the mass of Hydrogen. The total energy is ε = e + u²∕2 + ϕ, where e is the internal energy and ϕ the gravitational potential, with g = −∇ϕ the gravitational acceleration vector. The mean molecular weight is related to the mass-mixing ratio via $$\frac{1}{\mu (X)}=\frac{X}{{\mu }_1}+\frac{1-X}{{\mu }_2},$$(4)

where μ₁ and μ₂ are the mean molecular weights of the two fluid species. We assume in this numerical study a constant adiabatic index γ in order to finely control the convective regime of the system and ensure that the simulations will be Schwarzschild stable and Ledoux unstable (i.e. will have a stabilising temperature gradient and de-stabilising mean-molecular-weight gradient). The purpose of this choice is to assess the impact of the mean-molecular-weight gradient alone in the simulations, we discuss in Sect. 3 the limitations of this choice for realistic test cases. With a constant adiabatic index, the EOS can be rewritten as P = ρe(γ − 1), with e = c_vT, and c_v the specific heat capacity at constant volume: c_v = (k_B∕μ(X)m_H(γ − 1)).

The ρR(X, P, T) is a source term for the equation governing the advection of the mass-mixing ratio. Physically it represents the conversion between chemical species. ρc_pH(X, P, T) is the source term for the energy equation and represents energy input into the system via numerous sources, for example, radiation, thermal diffusion or release of latent heat via chemical reactions. If both R(X, P, T) and H(X, P, T) are non-zero, then the system in general is non-conservative.

2.2 Diabatic and adiabatic convection

For a fluid obeying an ideal gas law and under the influence of an external plane parallel gravitational field and experiencing no other momentum, energy or mass source terms, the criterion for convection is: $${\nabla }_T-{\nabla }_{\text{ad}}>0,$$(5)

where ∇_T = − h_p∂T₀∕∂z and ∇_ad = (γ − 1)∕γ, with 1∕h_p = − ∂P₀∕∂z, P₀ is the initial vertical pressure profile, and T₀ is the initial temperature profile. This is known as Schwarzschild convection. In this regime, ∇_ad = constant (assuming an ideal gas law), a system’s susceptibility to convection is determined solely by ∇_T, the vertical temperature gradient. Therefore, we are free to set the system to be Schwarzschild stable simply by setting a temperature gradient such that ∇_T < ∇_ad at t = t₀. Such a system is then stable to convection for all t > t₀.

However, if there are multiple molecular species present in the fluid, the mean molecular weight need not be homogeneous. A gradient can be defined that represents a transition from one chemical species to another; ∇_μ. Modifying Eq. (5) to account for ∇_μ leads to the following expression: $${\nabla }_T-{\nabla }_{\text{ad}}-{\nabla }_{\mu }>0.$$(6)

This means that even if the Schwarzschild criterion for stability is met, the system can still be unstable. Equation (6) is known as the Ledoux convection criterion. This criterion is referred as ‘adiabatic’ convection in Tremblin et al. (2019) by opposition to the ‘diabatic’ criterion that involves the chemical and thermal source terms. (See the appendix in Tremblin et al. (2019) for a full derivation of the instability criteria.) As proposed in Tremblin et al. (2019), the conserved quantities in the presenceof source terms are not potential temperature and composition but linear (logarithmic) combinations of both, similar to moist potential temperature in Earth atmosphere. However, the situation is simpler in the context of Ledoux convection: the conserved quantity in the linear regime is simply proportional to the entropy. At saturation, the Ledoux criterion is equivalent to homogeneous entropy, hence the saturation is a nearly isentropic state. However, the chemical source term can still impact Ledoux convection in the non-linear regime by replenishing the mean-molecular-weight gradient. In that context, we can expect a reduction of the temperature gradient for a Ledoux unstable system susceptible to a chemical reaction with a nearly constant entropy profile (i.e. adiabatic profile) as shown with the 2D simulations presented in Tremblin et al. (2019). We emphasise that the temperature-gradient reduction in the Ledoux regime (with the adiabatic criterion) is then different from the temperature-gradient reduction with the diabatic criterion (moist or other source terms) since it does not happen in a sub-adiabatic regime here.

2.3 Physical system

As discussed in the introduction, the atmospheric chemical transition CO+O ↔ CO₂ is a representative process for exploring the role that diabatic convection plays in shaping thermodynamic structures in rocky exoplanet atmospheres. This chemical transition is consistent with the atmosphere of the early Earth and young or highly irradiated terrestrial exoplanets in general. As such, the magnitudes of the physical quantities used in this study are consistent with this context and have been estimated with the 1D atmospheric code ATMO (Tremblin et al. 2015). However, we emphasise that ATMO is not yet ready to perform precise models of secondary atmospheres due to the lack of the relevant collisionally induced absorptions. We therefore simply estimated the pressure and temperature ranges and timescales of the CO+O ↔ CO₂ in a high metallicity model ([M/H] ~ 3). As a consequence, the simulations presented in this paper should be seen as an idealised numerical experiment and we discuss further these limitations in Sect. 3.

2.3.1 Initial conditions

The vertical temperature profile, T(z), is specified ab initio as a linear transition between pre-specified deep (T₀) and upper (T₁) atmosphere temperatures; thus, $$T_{\text{eq}}(z)=T_0+z(T_1-T_0),$$(7)

resulting ina constant negative temperature gradient in the z direction (T₀ > T₁). A temperature gradient of ΔT = 10⁻⁴ K cm⁻¹ was kept the same for every simulated parameter set, only the absolute values were changed. In the following, the subscript_mid refers to the middle of the temperature range. So for the deep-atmosphere temperature, T_0,mid = 2500 K and upper-atmosphere temperature, T_1,mid = 2300 K. In combination with the initial temperature profile, T₀(z), we set the mass-mixing ratio profile as a parametrised transition from zero to unity while T_0,mid < T_eq(z) < T_1,mid, thus: $$X_{\text{eq}}\left(T_{\text{eq}}\right)=\{\begin{array}{ll}1\hfill & \text{if}{T}_{\text{eq}}(z)<T_{1,\text{mid}}\hfill \\ 0\hfill & \text{if}{T}_{\text{eq}}(z)>T_{0,\text{mid}}\hfill \\ \frac{T_{\text{eq}}(z)-T_{0,\text{mid}}}{T_{1,\text{mid}}-T_{0,\text{mid}}}.\hfill & \text{else}\hfill \end{array}$$(8)

If T_eq(z) is below (above) this range it is held at zero (unity). The result of this parametrisation is to impose a transition region in the mean molecular weight, which is calculated directly from the mass-mixing ratio, which in turn comes form the temperature profile. As we explore the deep-atmosphere temperature range, T₀ = 2300 K → 2700 K, the position of this transition region transits across the simulation domain from bottom to top. While this procedure may seam convoluted, the outcome is that at the extremities of the temperature range (2300 K or 2700 K) the mass-mixing ratio is either zero or unity everywhere within the simulation domain. This results in ∇μ(X_eq) = 0 and hence convective stability (if ∇_T is such that Eq. (6) < 0). Figure 1 illustrates this parametrisation and lack of convection in the regions where X_eq = constant.

All quantities, including X_eq, are homogeneous in the x-y plane. Hence, in the following, we have included the k (z direction) index but omitted the horizontal plane indices j (y direction) and i (x direction). Vertical hydrostatic equilibrium is achieved using a combination of the EOS (Eq. (3)) and a finite difference approximation to the hydrostatic equation: $$\frac{P_k-P_{k-1}}{\Delta z}=-\frac{{\rho }_k+{\rho }_{k-1}}{2}{g}_z.$$(9)

Rearrangingthis expression and substituting Eq. (3), leads to the following pressure profile, $$p_k=p_{k-1}\frac{1-\eta \frac{\mu {(X)}_{k-1}}{T_{k-1}}}{1+\eta \frac{\mu {(X)}_k}{T_k}},$$(10)

where η is a constant: $$\eta =\frac{\Delta z}{2}\frac{m_{\text{H}}{g}_z}{k_{\text{B}}}.$$(11)

The ρ_k is then derived from Eq. (3).

Even if ∇μ(X_eq)≠0, hydrostatic equilibrium means that an initial perturbation is required to seed convective motion. This is prescribed by: $$\begin{array}{ll}\Delta u\hfill & =A\left(\begin{array}{c}\sin (2\pi {\theta }_x)\cos (2\pi {\theta }_y)\sin (\pi {\theta }_z)\\ \cos (2\pi {\theta }_x)\sin (2\pi {\theta }_y)\sin (\pi {\theta }_z)\\ \cos (2\pi {\theta }_x)\cos (2\pi {\theta }_y)\cos (\pi {\theta }_z)\end{array}\right),\hfill \end{array}$$(12)

with ${\theta }_x=\frac{x-x_{\text{mid}}}{L_x}$, ${\theta }_y=\frac{y-y_{\text{mid}}}{L_y}$ and ${\theta }_x=\frac{z-z_{\text{mid}}}{L_z}$ being the distances from the box midpoint in units of box width. This sets up a perturbation of two convective roles in the x and y directions, both spanning the complete height of the box. This is illustrated in the column of Fig. 1.

The form of Eq. (10) requires the specification of the deep-atmosphere temperature and pressure (P₀), the initial perturbation also requiresan amplitude, A. These parameters are summarised in Table 1.

Fig. 1 Diagram illustrating the parametrisation of the mass-mixing ratio and the subsequent impact on the convectively unstable region of the simulation domain. This parametrisation directly leads to a mean-molecular-weight gradient that drives the convection. At the limits of the temperature range explored, there is convective stability in the domain, as ∇μ(Xeq) = 0.

2.3.2 Thermal and chemical source terms

Heating in the form of radiation transfer from the base of the atmosphere is mimicked by applying a Newtonian source term to the energy equation. This effectively forces the system towards an equilibrium state, T_eq (in this case the initial profile) on a timescale set by τ_rad via the equation: $$\Delta T=\frac{T^n+T_{\text{eq}}\frac{\Delta t}{{\tau }_{\text{rad}}}}{1+\frac{\Delta t}{{\tau }_{\text{rad}}}}-T^n.$$(13)

This change in temperature is then converted into a change in energy via the EOS and incorporated into the source term in Eq. (2), H(X, P, T).

In an identical manner, the mass-mixing ratio, X, is forced to equilibrium, X_eq, according to: $$\Delta X=\frac{X^n+X_{\text{eq}}(T^n)\frac{\Delta t}{{\tau }_{\text{chem}}}}{1+\frac{\Delta t}{{\tau }_{\text{chem}}}}-X^n,$$(14)

on the timescale τ_chem. This expression is equivalent to R(X, P, T). Tⁿ and Xⁿ are the states at the current time step, before the forcing is applied. The forced mass-mixing ratio is then used to compute the updated mean molecular weight. This process acts to maintain the spacial distribution of CO and CO₂ and is in competition with the advection of X. if no forcing terms were present, then mixing and diffusion would result in a homogeneous value of X across the simulation domain. The role of Eq. (14) is then to represent the capture, advection and subsequent release of energy through composition and the EOS in the system.

Both the heating of the system and the action of chemical reactions are constrained to occur on specific timescales, τ_rad and τ_chem, respectively. These are parameterised according to the results of 1D calculations conducted using the ATMO code (Tremblin et al. 2015) and have the values τ_rad = 100 s and τ_chem = 10 s (these values are also in Table 1). The estimations of these timescales follow the same procedure used in Tremblin et al. (2015) to estimate the fingering instability criterion. However, we should point out here that the chemical network used in ATMO is targeted for hydrogen-dominated planets (Venot et al. 2012). These timescales are therefore rough estimates, and more detailed 1D modelling is needed to assess them properly for secondary atmospheres.

The radiative timescale, τ_rad, can be evaluated in a diffusion approximation as done in Tremblin et al. (2015): it typically varies from orders of magnitude in the atmosphere from 10⁷ s at 100 bars down to 0.1 s at 0.01 bar. Close to 1 bar τ_rad is of the order of 100 s. Using the ATMO code coupled to the out-of-equilibrium network of Venot et al. (2012), we explored the chemical timescale, τ_chem, by adjusting K_zz to get a quenching point at 1 bar. The timescale deduced K_zz is of the order of 10 s but we point out that it can also vary by orders of magnitude in the atmosphere. since we use a chemical network andnot a parametrisation we do not have easily access to the reaction timescales. As an indication, one can use the timescale for CO2 in the case of hydrogen-dominated planets that was computed in Zahnle & Marley (2014); it varies from 10⁻⁸ s at 100 bars to 10⁷ s at 0.01 bars for the pressure-temperature profile modelled with ATMO. These numbers cannot be used directly for our setup since we do not simulate a hydrogen-dominated planet but it gives an idea of how strongly the chemical timescale can vary. Based on the fact that both the chemical and radiative timescales vary strongly in the atmosphere and will be dependent on many details of the planet under consideration, we explore different cases in the paper, with a chemical timescale faster than the radiative timescale (base model) and vice versa.

To determine the impact of chemical timescale, τ_chem, on the convergence of the simulation results; several simulations of the central deep-atmosphere temperature T₀ = 2522 K were conducted. Four simulations were run with τ_chem = 1, 10, 100, 1000 s, the results of which are discussed in Sect. 3.3. However, the simulation suite that covered the full T₀ range was limited to τ_chem = 10.

2.3.3 Boundary conditions

The physical extent of the simulation domain was bounded by several types of conditions. Periodic boundary condition were imposed in the two horizontal directions, x and y. In the case of the vertical z direction, a custom condition was used for the lower and upper boundaries. At the upper and lower atmosphere’s boundaries, both the temperature and the passive scalar, X, was allowed to float, meaning that both quantities in the boundary were extrapolated from the computationally active region via a backwards or forwards (lower and upper boundary, respectively) finite difference approximation, assuming a constant gradient. We use Neumann boundary conditions for the transverse velocities, and reflexive boundary conditions for the normal velocity at the lower and upper boundary. Then the conservative variables are computed from the temperature and mean molecular weight, μ, (from the mass-mixing ratio) and Eq. (10) and the ideal gas law. This results in hydrostatic equilibrium enforced both below and above the simulation domain. This setup results in a simulation domain where there is stable non-convective atmosphere above and below the simulated region. These boundary conditions ensure that the normal velocity at the interface of the domain are exactly zero, hence the simulations conserve the total mass at machine precision.

2.4 Simulations

The numerical simulations were conducted using the ARK¹ (all-regime Kokkos) code (Padioleau et al. 2019), which implements the all-regime method of Chalons et al. (2016). This approach allows for both low and high Mach number flows to be captured in the same simulation. A secondary feature of the code is the ability to capture hydrostatic balance and ensures it to machine precision. This property is why the method is described as being well-balanced. The codes ability to reproduce convection in both the Rayleigh-Benard and Schwarzschild regimes has been demonstrated by Padioleau et al. (2019), making the code particularly well suited to exploring convection in the Ledoux regime.

ARK leverages the Kokkos framework for shared memory parallelism Edwards et al. (2014) in the C++ programming language. Inessence, Kokkos allows for code to be written in a generic manner but be executed on multiple hardware architectures, while maintaining high numerical performance. This programming paradigm is known as performance portability, a property critical for the present study as it allows for the use of accelerators such as Intel KNLs and GPGPUs. Convergence of the simulations, at the resolution employed in a reasonable time frame, could only be achieved via the use of such accelerators.

2.4.1 Computational setup

The simulations are broken down into two distinct categories. The first are relatively low resolution (relative to the maximum resolution used) and are designed specifically to not be numerically intensive, allowing for a parameter study in temperature to be conducted. The second category concentrates exclusively on a single simulation, allowing for higher resolution, and thus serves as a convergence study to verify the validity of the low resolution simulations. In addition, this high resolution simulation acts to test the scalability of ARK on large numbers of GPGPUs.

This second category of simulation was carried out on the Jean Zay supercomputer of IDRIS under a GENCI allocation for ‘grand challenge’. This machine allowed us to run ARK on up to 1000 Nvidia V100 GPGPUs. A weak scaling study to test ARK on this architecture is presented in Fig. 2 and illustrates the efficiency of ARK at the scales afforded by Jean Zay. This scaling study was conducted with 496³ numerical cells per GPGPU with each run lasting 1000 time steps.

Fig. 2 Weak scaling of the ARK code. Good efficiency is maintained for all numbers of GPGPUs, dropping to slightly lower than 80% for 600 or greater GPGPUs. This scaling test was conducted with a numerical grid of 4963 cells and used the same convection simulation as described in the text, but with a maximum iteration count of 1000.

2.4.2 Numerical algorithm and grid

The well-balanced all-regime method employed by ARK is a conservative, first-order accurate split scheme. This splitting allows the transport and acoustic terms to be solved separately with the gravitational acceleration added to the acoustic part, allowing for machine-accurate hydrostatic balance; full details of the method and its implementation in ARK are available in Chalons et al. (2016) and Padioleau et al. (2019).

The study was executed in two distinct phases. First, a parameter study in which a deep-atmosphere temperature range of 2300−2700 K was covered, with ten separate temperature steps. This increment number was chosen to be a balance between computation cost and parameter coverage necessary to illustrate the diabatic convection’s dependence on the deep-atmosphere temperature. The numerical grid used for this parameter study was 256³ and was kept constant for all values of deep-atmosphere temperature.

In the second phase, convergence was tested. This was achieved in separate resolution jumps where the resolution was doubled each jump from an initial low resolution of 320³. This was done to relax the simulation into higher resolutions without introducing blocking effects, where one cell is replaced by more than eight new cells at a time. The resolution jumps are as follows: 310³, 620³, 1240³, 2480³, and 4960³.

Due to the computationally intensive nature of the study, it is worth noting several steps that were taken to overcome the unique challenges encountered when performing simulations at this scale. At full resolution, the simulation grid was 5.2 TB in size (4960³ grid cells, each containing 6 state variables at double precision accuracy). As such, the number of snapshots possible was restricted and, therefore, time series analysis required us to implement routines for the on-the-fly analysis of global variables. In addition, for visual inspection, a slicing routine was added that outputs the x-y, x-z, and y-z 2D planes only, significantly reducing storage overhead and enabling high frame rate videos of the simulation to be produced.

Next we present the simulation results beginning with a stability analysis followed by the time series analysis with a discussion on the observable implication of the results, finally the resolution study is presented.

Fig. 3 Pseudo-2D profiles constructed from 1D profiles of the criterion for Ledoux (left) and Schwarzschild (right) convection as a function of both deep-atmosphere temperature and altitude. The profiles are averaged in both the x and y directions. Every deep-atmosphere temperature value represents a different simulation. The dark blue contour indicates unity and hence the threshold for the convective instability. Where the criterion is greater than unity, the system is susceptible to convection. We can see from the position of this contour that the system is unstable to Ledoux convection for all T0, both in the initial conditions and in steady state. This is not the case for Schwarzschild convection, with the initial conditions stable for all T0, both initially and once the system is in equilibrium.

3 Results and discussion

First we determine the system’s susceptibility to convection (both Schwarzschild and Ledoux) in the initial and steady state conditions. We then analyse the time series results and how they impact possible observables relevant to future space missions. We conclude by examining the system’s convective properties at high resolution.

3.1 Convective stability

Applying the stability criterion of Eqs. (6) and (5) to the initial conditions and to the final steady state system, produces the results shown in Fig. 3. At all temperatures and altitudes the initial conditions are stable to Schwarzschild convection. This is by design, as the system needs to be Schwarzschild stable to exclude this type of convection from the simulations such that our analysis can be done purely in the Ledoux regime. Lower temperatures and higher altitudes are displaying greater stability.

In contrast, the system is susceptible to Ledoux convection for a portion of the simulation domain for every temperature. The maximum altitude of the convectively susceptible region increases as the deep-atmosphere temperature increase. This dependence is due to the mass-mixing ratio being a direct function of the deep-atmosphere temperature and temperature range of the study. A pattern directly mirrored from the Schwarzschild criterion case above. However, in the case of the Ledoux criteria, the system exceeds zero at every deep-atmosphere temperate value. This observation makes clear the impact of the addition of the mean-molecular-weight gradient on the system’s convective stability.

This global behaviour is expected due to the manner in which the mean molecular weight is prescribed, as a linearly increasing profile that transitions from low to high altitude as the imposed deep-atmosphere temperature is progressively increased with each subsequent simulation. As the mean-molecular-weight gradient is the defining difference between the two convective criteria (Eqs. (6) and (5)) and as the system is Schwarzschild stable both initial and for all t > 0, we can conclude that any convection established is solely due to the imposed mean-molecular-weight gradient and maintained by the chemical and thermal source terms.

The same criterion is calculated for the system at t = t_max and is shown in the bottom row of Fig. 3. We can see that the susceptibility to both Schwarzschild and Ledoux convection, specified in the initial conditions, is maintained through to the end of the simulation. This illustrates the lack of stabilising or destabilising influences on the system beyond those we impose through the initial conditions and the source Eqs. (13) and (14).

As an aside, when running prototype simulations in 2D, the system was found to readily transition to shear; under the same conditions where the 3D simulations do not. Figure 4 illustrates this 2D shear profile with two distinct higher velocity bands running perpendicular to the direction of gravity and opposite to each other. A narrow boundary layer exists between the two flow streams where the velocity falls by an order of magnitude. The horizontal and vertical velocity components are also plotted below the velocity map in Fig. 4 and highlight the horizontal velocity’s growth and the degree to which shear has dominated the simulation at late times. These 2D shear modes have already been observed in the context of double-diffusive convection (Garaud & Brummell 2015) and seems to be a result of the dimensionality of the setup, an effect that we also identify in this paper.

From a physical perspective, this system should be susceptible to mixing via the Kelvin-Helmholtz instability. However, from a numerical point of view, the simulation is overall first-order accurate. The Riemann solver employed was similar to the Harten-Lax-van Leer-Contact (HLLC) type, known to be necessary for capturing this type of instability; however, as the spatial reconstruction and time integration is first order, the degree of dissipation in the model may simply be too great for the instability to establish.This is, however, purely speculative and we leave further discussion of the role of shear to a future publication in which the general case of shear in 3D simulations will be studied (see also Garaud & Brummell 2015; Garaud et al. 2019).

With the stability of the system and the convective regime established we will analyse the time-dependent nature of the simulations.

3.2 Time series

A time series of the parameter study is shown in Fig. 5. To avoid confusion, here it is worth reiterating the difference between what is meant by average bottom temperature, ${\overline{T}}_{\text{bottom}}$, a time-dependent quantity measured from the simulations; and the forced deep-atmosphere temperature, T₀, which is a parameter. Each separate coloured line in the legend of Fig. 5 refers to a different forced deep temperature, T₀.

The velocity magnitude, top left, was used as a steady state measure. Once this quantity becomes approximately constant with time, we assume steady state has been achieved. As the initial conditions are seeded by a convective perturbation, if the system itselfis unstable to convection (see Sect. 3.1), then this initial perturbation will grow rapidly until saturated ata maximum value. Once reached, the velocity magnitude will slightly oscillate about this maximum. However, if the system is stable to convection, then the initial perturbation will decay away steadily until swamped by numerical noise.

Both of the above behaviours are seen in Fig. 5, with both extremes in the deep-atmosphere temperature (2300 and 2700 K) displaying a decaying velocity magnitude. The velocity magnitude grew in every other simulation until saturation, which occurred at ~1000 s for all convectively unstable simulations.

The right plot of the top row of Fig. 5 shows the total energy evolution. Departure from constant energy may at first seem unexpected; however, it is physically justified as the simulation includes a source term and boundary conditions are not reciprocal on all faces. Indeed the z-direction boundary conditions are customised to produce the simulation conditions we seek. If all the boundary conditions were reciprocal and source terms set to zero, then the total energy would remain constant as the numerical scheme is conservative.

The top atmosphere temperature, ${\overline{T}}_{\text{top}}$, exhibits an increase across the T₀ range (except in the extremes 2300 and 2700 K). This increase, together with the decrease in ${\overline{T}}_{\text{bottom}}$, means a decrease in the vertical temperature gradient. This result indicates that the presence of convection directly leads to a reduction inthe upper atmosphere temperature with respect to that specified initially and that which would be maintained through the forcing profile.

As the temperatures T₀ = 2300 K and T₀ = 2700 K do not experience this temperature-gradient reduction – instead maintaining the hydrostatic equilibrium imposed initially but nonetheless allowed to evolve dynamically – we can conclude that convection in this model atmosphere is responsible for modifying the temperature structure. In the case of bulk planetary atmospheres, this will result in an observable upper-atmosphere temperature corresponding to a deep-atmosphere temperature (if hydrostatic equilibrium is assumed), which is larger than that actually present. This notion will be explored in greater depth in Sect. 3.5.

The division between the convectively stable (2300 and 2700 K) and unstable (2344 K → 2655 K) T₀ simulations is consistent in all the bulk quantities tracked throughout the evolution. The mid temperatures show the greatest degree of energy increase into the system. This is evident through the largest temperature changes at the top and bottom of the simulation domain.

Fig. 4 Illustration of the development of shear in 2D simulaitons. Top: velocity map from a 2D simulation illustrating the transition to shear. The quivers indicate the velocity direction and magnitude. There are two distinct regions of oppositely directed flow separated by a transition of approximately equal width. Bottom: time series of $\overline{\left|v_x\right|}$ and $\overline{\left|v_z\right|}$ showing the decay of the vertical component and the growth and saturation of the horizontal component. After a sufficient length of time (twice that of the 3D simulations), steady state has been reached and any convective motion is suppressed in favour of shear.

Fig. 5 Time evolution of bulk fluid quantities over the course of the parameter study simulations. For each of the simulated deep-atmosphere temperatures (coloured from yellow to deep blue according to decreasing temperature), the average velocity magnitude, total energy, average deep-atmosphere temperature, and average temperature of the upper atmosphere are plotted. The behaviour of the average velocity magnitude is the primary metric for determining whether steady state has been reached. This occurs at 7000 s.

3.3 Thermal forcing versus deep-atmosphere temperature

In the absence of convection, the temperature profile of the model would be completely imposed by the Newtonian forcing profile, for example in the radiative part of the atmosphere. The deep-atmosphere temperature would be simply given by T₀ (z_bottom). If we plot this forcing deep-atmosphere temperature against the final deep-atmosphere temperature, we can analyse the impact that convection has on the temperature profile. This is displayed in the left hand plot of Fig. 6 for initial (t = 0 s), intermediate (t = 500 s) and steady state (t = 10 000 s) time states.

If we assume that the region we are probing is stable but close to the radiative-convective boundary in the planet, we can assume that the energy flux is well approximated by ${\varphi }_{\text{base}}={\sigma }_{\text{SB}}{T}_0{(z_{\text{bottom}})}^4$. As a consequence, this energy flux is also a measure of the equivalent forcing convective flux in the deep atmosphere in a model following radiative-convective equilibrium. We are not explicitly forcing an energy flux at the bottom of the domain; however, the deep-atmosphere forcing temperature can, to some extend, be seen as a proxy for a forcing energy flux in our setup.

The dotted (initial condition adiabatic profile) and solid curves illustrate the difference in behaviour for the adiabatic and diabatic convecting systems, respectively. The system effectively undergoes a bifurcation, as the convective motion builds and the instability saturates. The deep-atmosphere temperature is forced away from the Newtonian (sub-adiabatic) profile towards the diabatic profile in the mid temperatures, where the mean-molecular-weight gradient is greatest. This effect is most apparent early in the simulation while the instability is still growing and before saturation. The dashed curve (t = 500 s) in Fig. 6 shows this initial phase and has a much more pronounced bifurcation than the steady state solution.

The behaviour described above means that the net effect of the convective instability (and therefore the presence of a mean-molecular-weight gradient), is that the deep-atmosphere temperature is cooler than the setup without mean-molecular-weight gradient. It is clear then that the upper atmosphere possesses a greater temperature than that allowed by thermal forcing and the deep-atmosphere temperature.

The bifurcation in the temperature plot in Fig. 6 is similar to the bifurcation proposed in Tremblin et al. (2019) in the case of the O-CH₄ transition inthe atmosphere of brown dwarfs and extrasolar giant exoplanets. It therefore suggests that the analogy to the Nukiyama curve for the boiling crisis experiment (Nukiyama 1934) is relevant and can be used to deduce bifurcation in the behaviour of convective systems in general.

The ratio of the forcing timescales, τ_chem∕τ_rad, determines the convergence of the deep-atmosphere temperature perturbation away from the adiabatic profile. To test this convergence, multiple iterations of the T₀ = 2522 K simulation were conducted with τ_chem = 1, 10, 100 and 1000 s. The right plot of Fig. 6 shows the deep-atmosphere temperature as a function of τ_chem; in the interest of brevity τ_rad = 100 s for all the simulations and only τ_chem was varied. As can be seen, the deep-atmosphere temperature does not undergo any further reduction once τ_chem∕τ_rad < 1, and converges at ~ 70 K less than the initial condition temperature of 2522 K. It also shows that the reduced temperature gradient is not present if the chemical timescale is too slow.

Fig. 6 Deep-atmosphere forcing (equilibrium) temperature and the ratio of timescales explored, illustrating convergence. Left: relationship between the deep-atmosphere forcing (equilibrium) temperature and the resulting deep-atmosphere temperature. The dotted curve represents the initial adiabatic conditions. The solid curve illustrates the deviation from this due to the convective motion. The final profile shows a departure from the adiabatic branch towards the diabatic branch. The dashed curve highlights a transient early stage when the perturbation towards the diabatic branch is most pronounced. Right: ratio of the forcing timescales τchem ∕τrad for T0 = 2522 K with τchem = 1, 10, 100, and 1000 s; τrad = 100 s for all values. Convergence is reached once τchem∕τrad < 1.

3.4 Convergence

To test the convergent properties of our computational setup, we selected a deep-atmosphere temperature in the middle of the temperature range, T₀ = 2522 K. We performed a set of simulations, varying the resolution between 310³ → 4960³. The lowest resolution is different from the resolution of the time series study simply to provide a value that can be recursively doubled-up to a sufficiently high resolution while simultaneously allowing for parallel decomposition of the computational domain. A snapshot of the final resolution (4960³) is displayed in Fig. 7, illustrating the turbulent fluid present in the simulation beyond that of the bulk convective motion.

To asses the impact of the resolution jumps and to ensure convergence at each intermediate stage, the average velocity magnitude and deep-atmosphere temperature were tracked. The results are displayed in Fig. 8 with each colour indicative of a different resolution. It can clearly be seen that the resolution jumps have no impact on these bulk quantities. With the velocity magnitude maintaining the convectively unstable profile. Indeed every resolution leaves the maximum velocity magnitude unchanged. This narrative persists in the case of the deep-atmosphere temperature, where the reduction due to the action of convection remains consistent across every resolution jump. It is also worth noting that the amplitudes of the fluctuations are around ~20 K, that is, approximately 25% of the overall temperature reduction. These fluctuations in the case of CO-CH₄ radiative convection (Tremblin et al. 2020) could potentially explain rotational spectral modulations in the atmosphere of brown dwarfs (Buenzli et al. 2012; Apai et al. 2013; Biller 2017; Artigau 2018).

Finally, to illustrate this separation of length and energy scales in relation to the convective instability and turbulent higher order behaviour, we have plotted in Fig. 9 the kinetic energy power spectrum, E(||k||). These data are taken from the x-y plane at the midpoint in the z direction, for each resolution in the convergence study. Overlaid is also the Kolmogorov spectrum, E(k) ∝ k^−5∕3, for comparison.

At large length scales (small ||k||), each resolution results in approximately the same contribution to E(||k||); however, there is an anomalous dip for all resolutions were 10⁻⁵ < ||k|| (x ≲ 6.3 km), At these scales, the bounding box of the simulation starts to be probed. As the reciprocal boundaries used in the simulation are not reflected in the calculation of E(||k||), when the smallest vales of ||k|| are probed there is only contributions from the box corners and hence an incomplete picture is obtained.

For 10⁻⁵ < ||k|| < 10⁻⁴ (0.63 km ≲x ≲ 6.3 km), we initially see the same E(||k||) for every resolution. This result is expected as the largest cell separation is ~ 0.065 km (resolution of 310³), meaning that wave numbers ||k|| ≈ 10⁻⁵ are spanned by ~10 numerical cells.

Only when 10⁻⁴ < ||k|| (x ≲ 0.63 km) do we see the impact of higher resolution and the additional turbulent modes that this allows for. The viscous dissipation at the lower limit is numerical in nature, this means that as kinetic energy cascades from large to small eddies, the cell width sets the lower limit on the length scale of the turbulent eddies.

Temperature stratification of the atmosphere is the focus of this study and as this is impacted by the large-scale convective motion, it is tempting to dismiss the role of small-scale structures. However, turbulence introduces instabilities that have the potential to disrupt large-scale behaviour, for example the transition to shear modes (see e.g. Garaud & Brummell 2015). This does not occur, and the unlocking of this kinetic energy cascade with higher resolution does not result in disruption of the convective behaviour. This indicates that small-scale dynamics do not inhibit bulk convective motion and inhabit separate regions of the kinetic energy power spectrum.

Fig. 7 Rendering of the density of the high resolution simulation’s entire computational volume. Units are in cgs. Both large- and small-scale features are evident, including Rayleigh-Taylor-like and Kelvin-Helmholtz-like instabilities. Despite the high degree of mixing via these turbulent processes, the large-scale vertical density gradient is still apparent and the convective motion is persistent.

Fig. 8 Time evolution of bulk fluid quantities for the resolution study. One can see that each successive resolution jump leaves the bulk velocity magnitude and deep-atmosphere temperature unchanged. This demonstrates that convergence is achieved at the low resolutions used in the parameter study.

3.5 Observational implications

In order to place the above results into the context of observational interpretation we compare the expectations from thermal forcing alone to the simulated curve, whose data comes directly from the simulation domain and assumes nothing about the thermal profile. we can see there is considerable difference over the 20 km of simulated atmosphere with the most pronounced difference being ~ 160 K for T₀ ≈ 2400 K in Fig. 6. This is illustrated more concisely in the left-hand plot of Fig. 10 where we can see the temperaturegradient is modified by the presence of a mean-molecular-weight gradient. This modification moves the system towards a constant entropy profile. This can be seen in the right hand-side of Fig. 10, were we plot the profile of the ratio of potential temperature to mean molecular weight, $$\frac{\theta (z)}{\mu (X)}=\frac{T(z)}{\mu (X)}{\left(\frac{P_0}{P(z)}\right)}^{\frac{\gamma -1}{\gamma }}.$$(15)

The point isthat the reduced temperature gradient occurs even in the situation where the entropy is homogenised (a point argued in Tremblin et al. 2019).

For the parameters used in this study, the scale height is H_p ~ 70 km. Therefore, over approximately a third of the scale height, the temperature deviation from the adiabatic profile is up to ~ 6% in the most pronounce case. It should be noted though that this effect is only present in the portion of the atmosphere where there is sufficient mean-molecular-weight gradient to promote convective motion. This ~ 6% temperature deviation will only impact this portion of the atmosphere depending on the chemical timescale used in the simulation, the overall temperature reduction across other sections of the atmosphere, at different altitudes, will depend on this chemical timescale and could be higher or lower.

The implication for observations of planetary atmospheres, is that the chemical composition and therefore the potential for molecular weight gradients need to be considered before assumptions can be made about the deep-atmosphere and overall temperature profiles. Convective instabilities such as Ledoux convection, driven by chemical and thermal source terms, as investigated here, result directly from these composition gradients.

Turning the above analysis upside down, the temperature of the upper atmosphere is often placed in the context of that which would be produced by the thermal transfer of energy from the deep to the upper atmosphere. Therefore, the presence of source terms and mean-molecular-weight gradients allows the upper atmosphere to be considerably hotter than that what would be afforded by a deep-atmosphere temperature according to the thermal forcing alone.

This potential for a reduced temperature gradient directly mirrors the so called reddening observed in the spectra of brown dwarfs (Tsuji et al. 1996; Chabrier et al. 2000; Allard et al. 2001; Marley et al. 2010, traditionally explained with silicate clouds), where the deep-atmosphere infrared flux is found to be smaller than expected based on adiabatically extrapolating the observed upper-atmosphere temperature down to the deep atmosphere. It has been proposed by Tremblin et al. (2016, 2019) that these theories overestimate the deep-atmosphere temperature by assuming an adiabatic convective profile. Diabatic convection provides therefore a compelling mechanism for accounting for this reddening in brown dwarfs. This type of behaviour could be at the origin of the need for a reduced temperature gradient in the atmosphere of cold brown dwarfs of spectral types T and Y (induced there by the transition between N₂ and NH₃, Leggett et al. 2017, 2019), with growing evidence that the increased mid-infrared flux cannot be explained by clouds for these objects. However, we emphasise again that the temperature-gradient reduction in the Ledoux regime is of a different nature since it happens even if the simulation is nearly adiabatic at saturation. The common point is that we can obtain reduced temperature gradients even when the temperature gradient is stabilising and this effect can have observational implications.

We have seen from the above simulations that the chemical transition of CO + O ↔ CO₂, and their separation into atmospheric layers with the associated boundary between them, is capable of providing the mean-molecular-weight gradient required for the establishment of Ledoux convection. Both CO₂ and CO are readily available in the atmospheres of terrestrial exoplanets, especially in the early stages of atmospheric formation. We therefore speculate that the mechanism responsible for the reddening observed in the spectra of brown dwarfs can also be present in the spectra of terrestrial exoplanetary atmospheres in the form of reduced temperature gradients in the Ledoux regime. This type ofmechanism can also be activated in the presence of other chemical transitions, or other chemical and thermal source terms in general, following the general theory developed in Tremblin et al. (2019) and should not be excluded a priori. The setup used in this paper is, however, idealised: we list in the next subsection the limitations of this work for more realistic applications.

The spectra and temperature of exoplanetary atmospheres will feature in the future observational objectives of instruments such as JWST and the European Extremely Large Telescope. It is therefore vital that the results are interpreted with robust theoretical models. We therefore make the case that the role of chemically and thermally driven convection needs to be incorporated into any analysis of observations concerned with planetary atmospheres.

Fig. 9 Kinetic energy power spectrum for each resolution jump used in the convergence study. The Kolmogorov spectrum, E(k) ∝ k−5∕3, is shown as the dashed black line for illustrative purposes. The E(||k||) is calculated in 2D at the midplane in the z direction. All resolutions converge at small ||k|| (large length scales), indicating that all resolutions capture the large-scale convective behaviour. At large ||k|| (small length scales), we see the cascading of kinetic energy characteristic of turbulence. At this point the different simulations begin to diverge as the additional resolution resolves the smaller scales to differing degrees in the simulations.

3.6 Limitations

We recap here the approximations and the three limitations of the idealised setup used in this paper. First, the temperature and compositional gradients need to be evaluated from detailed 1D models with all the physics (chemistry, opacities) relevant for rocky-exoplanet secondary atmospheres. Second, the chemical timescale needs to be evaluated with chemical networks validated for the conditions of rocky-exoplanet secondary atmospheres. Third, the EOS is oversimplified here, keeping a constant adiabatic index in order to finely control the convective behaviour: a more realistic EOS needs to be used.

For the last point, we highlight that this is not only a question of varying the adiabatic index with composition: since we are in a reactive fluid, the heat capacity has a reactive contribution well-defined when the system is close to chemical equilibrium. This leads to the following equation for c_v: $$\begin{array}{lllll}{c}_v\hfill & =\hfill & {\left(\frac{\partial e}{\partial T}\right)}_V\text{cr}\hfill & =\hfill & X_{\text{eq}}{c}_{v_1}+(1-X_{\text{eq}})c_{v_2}+\frac{\text{d}{X}_{\text{eq}}(T)}{\text{d}T}(e_1-e_2).\hfill \end{array}$$(16)

The resulting adiabatic index can then be deduced from the relation γ − 1 = k_b∕(c_vμ). This form should be relatively valid in our setup since we use a fast chemical timescale; however, we point out that these thermodynamic quantities are likely challenging to define properly when the system is out of chemical equilibrium (see e.g. Lebon & Jou 2008).

The setup used in this paper is therefore relatively idealised. It, however, shows that we cannot exclude a priori the occurrence of a reduced temperature gradient in these atmospheres.

Fig. 10 Temperature and potential temperature profiles illustrating the deviation from the adiabatic behaviour. Left: vertical temperature-gradient profiles comparing the adiabatic (solid line) to the diabatic (dashed line), corresponding to the dotted and solid curves, respectively, in the left hand plot. The colouring indicates a separate simulation, each with a different initial deep-atmosphere temperature. Each of these profiles correspond to a data point on the curves in the left hand plot. The central temperatures show the largest deviation from the adiabatic profile. Right: vertical profiles of potential temperature for both the initial and final conditions of the T0 = 2522 K simulation. The impact of the mean-molecular-weight gradient is to homogenise entropy.

4 Conclusions

We have numerically confirmed the ability of compositional and thermal source terms to induce a reduction in an atmospheric temperaturegradient in the Ledoux regime that could happen in the context of rocky exoplanet atmospheres. This has been done via both time series analysis and a convergence study that spanned an order of magnitude in resolution. We have also explored the impact of 2D versus 3D in these setups, showing that the shear modes emerging in 2D do not impact the 3D simulations.

We have studied an idealised test case based on a CO + O ↔ CO₂ transition, which has been used as a prototype chemical transition that is readily present in the atmospheres of hot rocky exoplanets. This makes it relevant not only for the study of present-day young or irradiated terrestrial exoplanet atmospheres but also paleo-climatology, as the thermal profiles and parameters used in this study are comparable to the early Earth’s climate. We have used this transition to demonstrate that a mean-molecular-weight gradient reduces the atmosphere’s temperature gradient, a result that does not converge away with resolution. We find that the above transition leads to a bifurcation of the atmosphere’s temperatureprofile, with the largest deviation coinciding with a maximum in the mean-molecular-weight gradient, a result that directly impacts how the deep-atmosphere temperature should be deduced from that of the upper atmosphere.

Although more detailed analyses with 1D models are needed, our results indicate the possibility of an analogy of the reddening (currently observed in the spectra of brown dwarfs) in the spectra of terrestrial exoplanet atmospheres. This is a new prediction for rocky exoplanet atmospheres and strongly motivates the use of diabatic models when interpreting atmospheric spectra. These types of observations will feature in future campaigns of instruments designed for the characterisation of the atmospheresof exoplanets, such as JWST and the European Extremely Large Telescope, making robust theoretical models vital for accurately interpreting observations.

Acknowledgments

We thank the referee for their helpful comments. S.D.Y. and P.T. acknowledge support by the European Research Council under grant agreement ATMO 757858. This work was granted access to the HPC resources of the supercomputer Joliot-Curie@TGCC under the allocation 2019-100925 made by GENCI and of the supercomputer Jean-Zay@IDRIS under the “grand challenge” 2019-100951.

References

All Tables

All Figures

	Fig. 1 Diagram illustrating the parametrisation of the mass-mixing ratio and the subsequent impact on the convectively unstable region of the simulation domain. This parametrisation directly leads to a mean-molecular-weight gradient that drives the convection. At the limits of the temperature range explored, there is convective stability in the domain, as ∇μ(Xeq) = 0.
In the text

	Fig. 2 Weak scaling of the ARK code. Good efficiency is maintained for all numbers of GPGPUs, dropping to slightly lower than 80% for 600 or greater GPGPUs. This scaling test was conducted with a numerical grid of 4963 cells and used the same convection simulation as described in the text, but with a maximum iteration count of 1000.
In the text

Fig. 3 Pseudo-2D profiles constructed from 1D profiles of the criterion for Ledoux (left) and Schwarzschild (right) convection as a function of both deep-atmosphere temperature and altitude. The profiles are averaged in both the x and y directions. Every deep-atmosphere temperature value represents a different simulation. The dark blue contour indicates unity and hence the threshold for the convective instability. Where the criterion is greater than unity, the system is susceptible to convection. We can see from the position of this contour that the system is unstable to Ledoux convection for all T0, both in the initial conditions and in steady state. This is not the case for Schwarzschild convection, with the initial conditions stable for all T0, both initially and once the system is in equilibrium.

In the text

Fig. 4 Illustration of the development of shear in 2D simulaitons. Top: velocity map from a 2D simulation illustrating the transition to shear. The quivers indicate the velocity direction and magnitude. There are two distinct regions of oppositely directed flow separated by a transition of approximately equal width. Bottom: time series of $\overline{\left|v_x\right|}$ and $\overline{\left|v_z\right|}$ showing the decay of the vertical component and the growth and saturation of the horizontal component. After a sufficient length of time (twice that of the 3D simulations), steady state has been reached and any convective motion is suppressed in favour of shear.

In the text

Fig. 5 Time evolution of bulk fluid quantities over the course of the parameter study simulations. For each of the simulated deep-atmosphere temperatures (coloured from yellow to deep blue according to decreasing temperature), the average velocity magnitude, total energy, average deep-atmosphere temperature, and average temperature of the upper atmosphere are plotted. The behaviour of the average velocity magnitude is the primary metric for determining whether steady state has been reached. This occurs at 7000 s.

In the text

Fig. 6 Deep-atmosphere forcing (equilibrium) temperature and the ratio of timescales explored, illustrating convergence. Left: relationship between the deep-atmosphere forcing (equilibrium) temperature and the resulting deep-atmosphere temperature. The dotted curve represents the initial adiabatic conditions. The solid curve illustrates the deviation from this due to the convective motion. The final profile shows a departure from the adiabatic branch towards the diabatic branch. The dashed curve highlights a transient early stage when the perturbation towards the diabatic branch is most pronounced. Right: ratio of the forcing timescales τchem ∕τrad for T0 = 2522 K with τchem = 1, 10, 100, and 1000 s; τrad = 100 s for all values. Convergence is reached once τchem∕τrad < 1.

In the text

Fig. 7 Rendering of the density of the high resolution simulation’s entire computational volume. Units are in cgs. Both large- and small-scale features are evident, including Rayleigh-Taylor-like and Kelvin-Helmholtz-like instabilities. Despite the high degree of mixing via these turbulent processes, the large-scale vertical density gradient is still apparent and the convective motion is persistent.

In the text

	Fig. 8 Time evolution of bulk fluid quantities for the resolution study. One can see that each successive resolution jump leaves the bulk velocity magnitude and deep-atmosphere temperature unchanged. This demonstrates that convergence is achieved at the low resolutions used in the parameter study.
In the text

Fig. 9 Kinetic energy power spectrum for each resolution jump used in the convergence study. The Kolmogorov spectrum, E(k) ∝ k−5∕3, is shown as the dashed black line for illustrative purposes. The E(||k||) is calculated in 2D at the midplane in the z direction. All resolutions converge at small ||k|| (large length scales), indicating that all resolutions capture the large-scale convective behaviour. At large ||k|| (small length scales), we see the cascading of kinetic energy characteristic of turbulence. At this point the different simulations begin to diverge as the additional resolution resolves the smaller scales to differing degrees in the simulations.

In the text

Fig. 10 Temperature and potential temperature profiles illustrating the deviation from the adiabatic behaviour. Left: vertical temperature-gradient profiles comparing the adiabatic (solid line) to the diabatic (dashed line), corresponding to the dotted and solid curves, respectively, in the left hand plot. The colouring indicates a separate simulation, each with a different initial deep-atmosphere temperature. Each of these profiles correspond to a data point on the curves in the left hand plot. The central temperatures show the largest deviation from the adiabatic profile. Right: vertical profiles of potential temperature for both the initial and final conditions of the T0 = 2522 K simulation. The impact of the mean-molecular-weight gradient is to homogenise entropy.

In the text