Issue 
A&A
Volume 653, September 2021



Article Number  A54  
Number of page(s)  13  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202040120  
Published online  07 September 2021 
Idealised 3D simulations of diabatically forced Ledoux convection
Application to the atmosphere of hot rocky exoplanets
Université ParisSaclay, UVSQ, CNRS, CEA, Maison de la Simulation,
91191
GifsurYvette,
France
email: sddyates@gmail.com
Received:
11
December
2020
Accepted:
15
June
2021
Aims. We investigate the impact of dimensionality, resolution, and long timescales on convective numerical simulations forced by thermocompositional diabatic processes. We focus our study on simulations that are stable to the Schwarzschild criterion but unstable to the Ledoux one (i.e. simulations with a stabilising temperature gradient and a destabilising meanmolecularweight gradient). We aim to establish the possibility of a reduced temperature gradient in such setups.
Methods. A suite of 3D simulations incorporating both time series and convergence studies were conducted using a highperformance numerical hydrodynamic code. We used, as a simplified and idealised test case, a sample region of the secondary atmosphere of a hot rocky exoplanet, of the order of the scale height of the system, within which the chemical transition CO + O ↔ CO_{2} could occur. Newtonian cooling was employed to force an equilibrium temperature, and a chemical source term was used to maintain a negative meanmolecularweight gradient in the vertical direction.
Results. Our results demonstrate that a meanmolecularweight gradient and a chemical source term can reduce the atmosphere temperature gradient, a result that does not converge away with resolution and is stable when exploring long timescales. Simulations in two dimensions are prone to the development of shear modes, as already seen in the literature for doublediffusive convection. The 3D convective steady state is not impacted by these shear modes, suggesting that this phenomenon is linked to the dimensionality of the problem. We also show that the presence of the reduced temperature gradient is a function of the forcing timescales, disappearing if the chemical forcing is too slow. We find that the above transition leads to a bifurcation of the atmosphere’s temperature profile when the chemical forcing is fast. Such a bifurcation is reminiscent of the bifurcation seen in the boiling crisis for steam or liquid convection.
Conclusions. With the reduced temperature gradient in these idealised setups, there exists the possibility of an analogy of the reddening (currently observed in the spectra of brown dwarfs) in the spectra of rocky exoplanet atmospheres. This possibility needs, however, to be checked with detailed 1D models in order to precisely characterise the equilibrium thermal and compositional gradients, the thermal and compositional forcing timescales, and the impact of a realistic equation of state to, in turn, assess if the regime identified here will develop in realistic situations. However, the possibility of this reddening cannot be excluded a priori. This prediction is new for terrestrial atmospheres and represents strong motivation for the use of diabatic models when analysing the atmospheric spectra of rocky exoplanets that will be observed with, for example, the James Webb Space Telescope.
Key words: planets and satellites: atmospheres / planets and satellites: terrestrial planets
© S. DaleyYates et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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 thermocompositional 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_{4} 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 thermocompositional diabatic convection, proposed by Tremblin et al. (2019), by investigating the impact of dimensionality, resolution, and long timescales on convective numerical simulations forced by thermocompositional diabatic processes. As an idealised test case, we chose a potential CO+O ↔ CO_{2} 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; ModirroustaGalian et al. 2021). This is partially due to the imminent arrival of several nextgeneration 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 ElkinsTanton (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_{2} transition. Such a transition is likely to have occurred early during the formation of the Earth’s atmosphere as CO_{2} 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 paleoclimatology 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 temperaturedependent 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 meanmolecularweight gradient. In the context of the present study, this results in a CO to CO_{2} 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_{2} dominated gas will fall and displace the lighter, less dense COdominated gas, which will in turn rise and replace the CO_{2}, 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_{2} 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_{2} 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 COCO_{2} mixtures and at different forcing timescales, as well as the system’s response to smallscale flow features, which are accessible at high resolution. In these idealised setups, we show that a temperaturegradient 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: (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: (2)
where ρ, u, X, P, ε, c_{p}, and T are the density, velocity, massmixing ratio, pressure, total energy, and specific heat capacity at constant pressure and temperature, respectively. This system is closed using the ideal gas law, (3)
where m_{H} is the mass of Hydrogen. The total energy is ε = e + u^{2}∕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 massmixing ratio via (4)
where μ_{1} and μ_{2} 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 destabilising meanmolecularweight gradient). The purpose of this choice is to assess the impact of the meanmolecularweight 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_{v}T, 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 massmixing ratio. Physically it represents the conversion between chemical species. ρc_{p}H(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 nonzero, then the system in general is nonconservative.
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: (5)
where ∇_{T} = − h_{p}∂T_{0}∕∂z and ∇_{ad} = (γ − 1)∕γ, with 1∕h_{p} = − ∂P_{0}∕∂z, P_{0} is the initial vertical pressure profile, and T_{0} 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_{0}. Such a system is then stable to convection for all t > t_{0}.
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: (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 nonlinear regime by replenishing the meanmolecularweight 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 temperaturegradient reduction in the Ledoux regime (with the adiabatic criterion) is then different from the temperaturegradient reduction with the diabatic criterion (moist or other source terms) since it does not happen in a subadiabatic regime here.
2.3 Physical system
As discussed in the introduction, the atmospheric chemical transition CO+O ↔ CO_{2} 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_{2} 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 prespecified deep (T_{0}) and upper (T_{1}) atmosphere temperatures; thus, (7)
resulting ina constant negative temperature gradient in the z direction (T_{0} > T_{1}). A temperature gradient of ΔT = 10^{−4} K cm^{−1} 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 deepatmosphere temperature, T_{0,mid} = 2500 K and upperatmosphere temperature, T_{1,mid} = 2300 K. In combination with the initial temperature profile, T_{0}(z), we set the massmixing ratio profile as a parametrised transition from zero to unity while T_{0,mid} < T_{eq}(z) < T_{1,mid}, thus: (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 massmixing ratio, which in turn comes form the temperature profile. As we explore the deepatmosphere temperature range, T_{0} = 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 massmixing 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 xy 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: (9)
Rearrangingthis expression and substituting Eq. (3), leads to the following pressure profile, (10)
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: (12)
with , and 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 deepatmosphere temperature and pressure (P_{0}), the initial perturbation also requiresan amplitude, A. These parameters are summarised in Table 1.
Fig. 1 Diagram illustrating the parametrisation of the massmixing ratio and the subsequent impact on the convectively unstable region of the simulation domain. This parametrisation directly leads to a meanmolecularweight gradient that drives the convection. At the limits of the temperature range explored, there is convective stability in the domain, as ∇μ(X_{eq}) = 0. 
Simulation parameters.
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: (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 massmixing ratio, X, is forced to equilibrium, X_{eq}, according to: (14)
on the timescale τ_{chem}. This expression is equivalent to R(X, P, T). T^{n} and X^{n} are the states at the current time step, before the forcing is applied. The forced massmixing ratio is then used to compute the updated mean molecular weight. This process acts to maintain the spacial distribution of CO and CO_{2} 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 hydrogendominated 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^{7} 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 outofequilibrium 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 hydrogendominated planets that was computed in Zahnle & Marley (2014); it varies from 10^{−8} s at 100 bars to 10^{7} s at 0.01 bars for the pressuretemperature profile modelled with ATMO. These numbers cannot be used directly for our setup since we do not simulate a hydrogendominated 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 deepatmosphere temperature T_{0} = 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_{0} 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 massmixing 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 nonconvective 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^{1} (allregime Kokkos) code (Padioleau et al. 2019), which implements the allregime 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 wellbalanced. The codes ability to reproduce convection in both the RayleighBenard 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^{3} 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 496^{3} 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 wellbalanced allregime method employed by ARK is a conservative, firstorder 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 machineaccurate 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 deepatmosphere 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 deepatmosphere temperature. The numerical grid used for this parameter study was 256^{3} and was kept constant for all values of deepatmosphere 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^{3}. 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^{3}, 620^{3}, 1240^{3}, 2480^{3}, and 4960^{3}.
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^{3} 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 onthefly analysis of global variables. In addition, for visual inspection, a slicing routine was added that outputs the xy, xz, and yz 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 Pseudo2D profiles constructed from 1D profiles of the criterion for Ledoux (left) and Schwarzschild (right) convection as a function of both deepatmosphere temperature and altitude. The profiles are averaged in both the x and y directions. Every deepatmosphere 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 T_{0}, both in the initial conditions and in steady state. This is not the case for Schwarzschild convection, with the initial conditions stable for all T_{0}, 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 deepatmosphere temperature increase. This dependence is due to the massmixing ratio being a direct function of the deepatmosphere 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 deepatmosphere temperate value. This observation makes clear the impact of the addition of the meanmolecularweight 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 deepatmosphere temperature is progressively increased with each subsequent simulation. As the meanmolecularweight 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 meanmolecularweight 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 doublediffusive 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 KelvinHelmholtz instability. However, from a numerical point of view, the simulation is overall firstorder accurate. The Riemann solver employed was similar to the HartenLaxvan LeerContact (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 timedependent 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, , a timedependent quantity measured from the simulations; and the forced deepatmosphere temperature, T_{0}, which is a parameter. Each separate coloured line in the legend of Fig. 5 refers to a different forced deep temperature, T_{0}.
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 deepatmosphere 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 zdirection 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, , exhibits an increase across the T_{0} range (except in the extremes 2300 and 2700 K). This increase, together with the decrease in , 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_{0} = 2300 K and T_{0} = 2700 K do not experience this temperaturegradient 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 upperatmosphere temperature corresponding to a deepatmosphere 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_{0} 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 and 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 deepatmosphere temperatures (coloured from yellow to deep blue according to decreasing temperature), the average velocity magnitude, total energy, average deepatmosphere 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 deepatmosphere 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 deepatmosphere temperature would be simply given by T_{0} (z_{bottom}). If we plot this forcing deepatmosphere temperature against the final deepatmosphere 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 radiativeconvective boundary in the planet, we can assume that the energy flux is well approximated by . As a consequence, this energy flux is also a measure of the equivalent forcing convective flux in the deep atmosphere in a model following radiativeconvective equilibrium. We are not explicitly forcing an energy flux at the bottom of the domain; however, the deepatmosphere 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 deepatmosphere temperature is forced away from the Newtonian (subadiabatic) profile towards the diabatic profile in the mid temperatures, where the meanmolecularweight 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 meanmolecularweight gradient), is that the deepatmosphere temperature is cooler than the setup without meanmolecularweight gradient. It is clear then that the upper atmosphere possesses a greater temperature than that allowed by thermal forcing and the deepatmosphere 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 OCH_{4} 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 deepatmosphere temperature perturbation away from the adiabatic profile. To test this convergence, multiple iterations of the T_{0} = 2522 K simulation were conducted with τ_{chem} = 1, 10, 100 and 1000 s. The right plot of Fig. 6 shows the deepatmosphere 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 deepatmosphere 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 Deepatmosphere forcing (equilibrium) temperature and the ratio of timescales explored, illustrating convergence. Left: relationship between the deepatmosphere forcing (equilibrium) temperature and the resulting deepatmosphere 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 T_{0} = 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 deepatmosphere temperature in the middle of the temperature range, T_{0} = 2522 K. We performed a set of simulations, varying the resolution between 310^{3} → 4960^{3}. The lowest resolution is different from the resolution of the time series study simply to provide a value that can be recursively doubledup to a sufficiently high resolution while simultaneously allowing for parallel decomposition of the computational domain. A snapshot of the final resolution (4960^{3}) 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 deepatmosphere 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 deepatmosphere 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 COCH_{4} 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 xy 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^{−5} < 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^{−5} < k < 10^{−4} (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^{3}), meaning that wave numbers k ≈ 10^{−5} are spanned by ~10 numerical cells.
Only when 10^{−4} < 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 largescale convective motion, it is tempting to dismiss the role of smallscale structures. However, turbulence introduces instabilities that have the potential to disrupt largescale 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 smallscale 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 smallscale features are evident, including RayleighTaylorlike and KelvinHelmholtzlike instabilities. Despite the high degree of mixing via these turbulent processes, the largescale 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 deepatmosphere 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_{0} ≈ 2400 K in Fig. 6. This is illustrated more concisely in the lefthand plot of Fig. 10 where we can see the temperaturegradient is modified by the presence of a meanmolecularweight gradient. This modification moves the system towards a constant entropy profile. This can be seen in the right handside of Fig. 10, were we plot the profile of the ratio of potential temperature to mean molecular weight, (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 meanmolecularweight 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 deepatmosphere 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 meanmolecularweight gradients allows the upper atmosphere to be considerably hotter than that what would be afforded by a deepatmosphere 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 deepatmosphere infrared flux is found to be smaller than expected based on adiabatically extrapolating the observed upperatmosphere temperature down to the deep atmosphere. It has been proposed by Tremblin et al. (2016, 2019) that these theories overestimate the deepatmosphere 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_{2} and NH_{3}, Leggett et al. 2017, 2019), with growing evidence that the increased midinfrared flux cannot be explained by clouds for these objects. However, we emphasise again that the temperaturegradient 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_{2}, and their separation into atmospheric layers with the associated boundary between them, is capable of providing the meanmolecularweight gradient required for the establishment of Ledoux convection. Both CO_{2} 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 largescale 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 rockyexoplanet secondary atmospheres. Second, the chemical timescale needs to be evaluated with chemical networks validated for the conditions of rockyexoplanet 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 welldefined when the system is close to chemical equilibrium. This leads to the following equation for c_{v}: (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 temperaturegradient 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 deepatmosphere 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 T_{0} = 2522 K simulation. The impact of the meanmolecularweight 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_{2} 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 presentday young or irradiated terrestrial exoplanet atmospheres but also paleoclimatology, 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 meanmolecularweight 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 meanmolecularweight gradient, a result that directly impacts how the deepatmosphere 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 JoliotCurie@TGCC under the allocation 2019100925 made by GENCI and of the supercomputer JeanZay@IDRIS under the “grand challenge” 2019100951.
References
 Allard, F., Hauschildt, P. H., Alexander, D. R., Tamanai, A., & Schweitzer, A. 2001, ApJ, 556, 357 [NASA ADS] [CrossRef] [Google Scholar]
 Apai, D., Radigan, J., Buenzli, E., et al. 2013, ApJ, 768, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Arakawa, A., & Jung, J.H. 2011, Atmos. Res., 102, 263 [NASA ADS] [CrossRef] [Google Scholar]
 Artigau, É. 2018, Handbook of Exoplanets (Berlin: Springer), 94 [Google Scholar]
 Biller, B. 2017, Astron. Rev., 13, 1 [Google Scholar]
 Bodenschatz, E., Pesch, W., & Ahlers, G. 2000, Ann. Rev. Fluid Mech., 32, 709 [NASA ADS] [CrossRef] [Google Scholar]
 Brown, J. M., Garaud, P., & Stellmach, S. 2013, ApJ, 768, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Buenzli, E., Apai, D., Morley, C. V., et al. 2012, ApJ, 760, L31 [NASA ADS] [CrossRef] [Google Scholar]
 Busse, F. H. 1978, Rep. Prog. Phys., 41, 1929 [Google Scholar]
 Chabrier, G., Baraffe, I., Allard, F., & Hauschildt, P. 2000, ApJ, 542, 464 [Google Scholar]
 Chalons, C., Girardin, M., & Kokh, S. 2016, Commun. Comput. Phys., 20, 188 [Google Scholar]
 Denissenkov, P. A. 2010, ApJ, 723, 563 [NASA ADS] [CrossRef] [Google Scholar]
 Edwards, H. C., Trott, C. R., & Sunderland, D. 2014, J. Parallel Distrib. Comput., 74, 3202 [CrossRef] [PubMed] [Google Scholar]
 ElkinsTanton, L. 2008, Earth Planet. Sci. Lett., 271, 181 [NASA ADS] [CrossRef] [Google Scholar]
 Forget, F., & Leconte, J. 2014, Phil. Trans. R. Soc. Ser. A, 372, 20130084 [Google Scholar]
 Garaud, P., & Brummell, N. 2015, ApJ, 815, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Garaud, P., Medrano, M., Brown, J. M., Mankovich, C., & Moore, K. 2015, ApJ, 808, 89 [Google Scholar]
 Garaud, P., Kumar, A., & Sridhar, J. 2019, ApJ, 879, 60 [NASA ADS] [CrossRef] [Google Scholar]
 Gregg, M. 1988, SmallScale Turbulence and Mixing in the Ocean, eds. Nihoul J., & Jamart B., (Amsterdam, The Netherlands: Elsevier), 46 [Google Scholar]
 Hurlburt, N. E., Toomre, J., & Massaguer, J. M. 1984, ApJ, 282, 557 [NASA ADS] [CrossRef] [Google Scholar]
 Ito, Y., Ikoma, M., Kawahara, H., et al. 2015, ApJ, 801, 144 [NASA ADS] [CrossRef] [Google Scholar]
 Lebon, G., & Jou, D. 2008, Understanding Nonequilibrium Thermodynamics: Foundations, Applications, Frontiers (Berlin, Heidelberg: Springer) [Google Scholar]
 Leggett, S. K., Tremblin, P., Esplin, T. L., Luhman, K. L., & Morley, C. V. 2017, ApJ, 842, 118 [NASA ADS] [CrossRef] [Google Scholar]
 Leggett, S. K., Dupuy, T. J., Morley, C. V., et al. 2019, ApJ, 882, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Madhusudhan, N. 2019, ARA&A, 57, 617 [NASA ADS] [CrossRef] [Google Scholar]
 Marley, M. S., Saumon, D., & Goldblatt, C. 2010, ApJ, 723, L117 [NASA ADS] [CrossRef] [Google Scholar]
 ModirroustaGalian, D., Ito, Y., & Micela, G. 2021, Icarus, 358, 114175 [CrossRef] [Google Scholar]
 Morley, C. V., Fortney, J. J., Marley, M. S., et al. 2012, ApJ, 756, 172 [NASA ADS] [CrossRef] [Google Scholar]
 Nukiyama, S. 1934, J. Japan Soc. Mech. Eng., 37, 367 [Google Scholar]
 Padioleau, T., Tremblin, P., Audit, E., Kestener, P., & Kokh, S. 2019, ApJ, 875, 128 [CrossRef] [Google Scholar]
 Saumon, D., & Marley, M. S. 2008, ApJ, 689, 1327 [Google Scholar]
 Stern, M. E. 1960, Tellus, 12, 172 [NASA ADS] [CrossRef] [Google Scholar]
 Stevenson, K. B., Lewis, N. K., Bean, J. L., et al. 2016, PASP, 128, 094401 [NASA ADS] [CrossRef] [Google Scholar]
 Tinetti, G., Drossart, P., Eccleston, P., et al. 2018, Exp. Astron., 46, 135 [NASA ADS] [CrossRef] [Google Scholar]
 Tremblin, P., Amundsen, D. S., Mourier, P., et al. 2015, ApJ, 804, L17 [NASA ADS] [CrossRef] [Google Scholar]
 Tremblin, P., Amundsen, D. S., Chabrier, G., et al. 2016, ApJ, 817, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Tremblin, P., Chabrier, G., Baraffe, I., et al. 2017, ApJ, 850, 46 [NASA ADS] [CrossRef] [Google Scholar]
 Tremblin, P., Padioleau, T., Phillips, M. W., et al. 2019, ApJ, 876, 144 [Google Scholar]
 Tremblin, P., Phillips, M. W., Emery, A., et al. 2020, A&A, 643, A23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tsuji, T., Ohnaka, K., Aoki, W., & Nakajima, T. 1996, A&A, 308, L29 [NASA ADS] [Google Scholar]
 Turner, J. 1967, Deep Sea Res. Oceanograp. Abs., 14, 599 [Google Scholar]
 Venot, O., Hébrard, E., Agúndez, M., et al. 2012, A&A, 546, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yano, J.I. 2014, DyAtO, 67, 1 [NASA ADS] [Google Scholar]
 Zahnle, K. J., & Marley, M. S. 2014, ApJ, 797, 41 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1 Diagram illustrating the parametrisation of the massmixing ratio and the subsequent impact on the convectively unstable region of the simulation domain. This parametrisation directly leads to a meanmolecularweight gradient that drives the convection. At the limits of the temperature range explored, there is convective stability in the domain, as ∇μ(X_{eq}) = 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 496^{3} 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 Pseudo2D profiles constructed from 1D profiles of the criterion for Ledoux (left) and Schwarzschild (right) convection as a function of both deepatmosphere temperature and altitude. The profiles are averaged in both the x and y directions. Every deepatmosphere 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 T_{0}, both in the initial conditions and in steady state. This is not the case for Schwarzschild convection, with the initial conditions stable for all T_{0}, 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 and 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 deepatmosphere temperatures (coloured from yellow to deep blue according to decreasing temperature), the average velocity magnitude, total energy, average deepatmosphere 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 Deepatmosphere forcing (equilibrium) temperature and the ratio of timescales explored, illustrating convergence. Left: relationship between the deepatmosphere forcing (equilibrium) temperature and the resulting deepatmosphere 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 T_{0} = 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 smallscale features are evident, including RayleighTaylorlike and KelvinHelmholtzlike instabilities. Despite the high degree of mixing via these turbulent processes, the largescale 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 deepatmosphere 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 largescale 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 temperaturegradient 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 deepatmosphere 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 T_{0} = 2522 K simulation. The impact of the meanmolecularweight gradient is to homogenise entropy. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.