© P. Tremblin et al. 2021

1 Introduction

With an increasing number of observations of exoplanet atmospheres taking place using transmission spectroscopy (e.g., Sing et al. 2016), a growing concern regarding cloud cover has been abound among the exoplanet community. Clouds may lead to “flat” spectra of the atmosphere of super-Earths and mini-Neptunes (Kreidberg et al. 2014), potentially reducing our ability to obtain information about the gas phase but also providing clues of where and how condensible materials are forming and evolving.

Cloud modeling remains a very important challenge in exoplanetology given the complexity of the phenomena that link chemistry with complex microphysics, radiative transfer, and hydrodynamics. Several complementary approaches have been developed in the past few years in response to this challenge: global 1D models with either simplified (e.g. Fortney et al. 2008; Tan & Showman 2019) or complex (e.g. Helling et al. 2019) microphysics; and 3D global circulation models (GCMs) with passive (e.g. Parmentier et al. 2016) or active (e.g. Lines et al. 2018) clouds. Early studies such as that of Gierasch et al. (1973) and more recent ones (e.g. Tan & Showman 2021a,b) have also identified the key role of radiative instabilities at a global scale – in planetary atmospheres (Venus, Jupiter) and in the atmosphere of exoplanets and brown dwarfs.

However, all these approaches use a simplified approach to hydrodynamics, either because they are 1D or because they are global and cannot adequately capture small-scale, non-hydrostatic buoyancy effects. We therefore propose in this paper another complementary approach based on theory with the inclusion of source terms in the Rayleigh-Taylor stability analysis in a similar way to the recent theoretical development of diabatic convection proposed in Tremblin et al. (2019). The classical (incompressible) Rayleigh-Taylor instability is triggered by a jump in density with a heavy fluid on top of a light fluid (see Chandrasekhar 1961; Zhou 2017a,b), but cannot account for the stabilizing-destabilizing effect of source terms similarly to Schwarzschild convection. We then use local small-scale simulations in order to properly study the interplay between buoyancy and radiative transfer with opacity discontinuities and its impact on the structure of clouds.

Tremblin et al. (2019) recently proposed a new paradigm for the development of convective motions in the presence of compositional and thermal source terms, namely, via diabatic convection. This paradigm can describe many convective systems, such as moist convection in Earth atmosphere, thermohaline convection in Earth oceans, and CO/CH₄ radiative convection in the atmosphere of brown dwarfs and giant exoplanets. In this paper, we study in a similar way the dynamical behavior of an opacity discontinuity between a “higher-opacity” and “lower-opacity” medium subject to radiative heating and cooling. Such a discontinuity could be unstable or stable due to the impact of radiative transfer on buoyancy similarly to the impact of thermal source terms in diabatic convection. This radiative Rayleigh-Taylor instability (RRTI, hereafter) is, however, an interface instability that is different from diabatic convection, akin to the difference between the standard Rayleigh-Taylor instability and Schwarzschild convection in the adiabatic case.

In Sect. 2 and Appendix A, we present the linear stability analysis predicting the RRTI growth rate and we present our numerical setup in Sect. 3. In Sect. 4, we propose a first idealized study in the brown-dwarf regime and provide a simplified understanding of the RRTI mechanism. In Sect. 5, we study the impact of the temperature gradient of the atmosphere on this instability for an Earth-like idealized setup. In Sect. 6, we discuss the possible implications of this instability on the cloud structure of different atmospheres. We then give our conclusions in Sect. 7 and discuss the possible expectations for future exoplanet observations.

2 Linear stability analysis for the radiative Rayleigh-Taylor instability

We start from the equations of hydrodynamics with gravity and a simple heating source term H(X, T) in the total energy equation that depends on temperature, T, and composition, X, to model the radiation: $$\begin{array}{lllllllll}\frac{\partial \rho }{\partial t}+\nabla \left(\rho u\right)\hfill & =\hfill & 0\text{cr}\frac{\partial \rho u}{\partial t}+\nabla \left(\rho u\otimes u+P\right)\hfill & =\hfill & \rho g\text{ cr}\frac{\partial \rho X}{\partial t}+\nabla \left(\rho Xu\right)\hfill & =\hfill & 0\text{ cr}\frac{\partial \rho \mathcal{E}}{\partial t}+\nabla \left(u\left(\rho \mathcal{E}+P\right)\right)\hfill & =\hfill & \rho c_{p}H(X,T),\hfill \end{array}$$(1)

with the total energy as $\mathcal{E}=e+u^2/2+\varphi$, the internal energy, e, the velocity u, the gravitational potential ϕ, along with the equation of state (EOS) of an ideal gas as ρe(γ − 1) = P, and the gravity as g = −∇ϕ (aligned with the y axis in the rest of the paper). With the ideal gas law, the temperature and the (constant) mean molecular weight are linked via P = ρk_bT∕μ. We assume an hydrostatic background with: $$\frac{\partial P_0}{\partial y}=-{\rho }_{0}g,H(X_0,T_0)=0,$$(2)

assuming a discontinuity of composition in two parts of the domain with $X_0^{+}$ in the upper half and $X_0^{-}$ in the lower half and with a continuous density at the interface ρ₀. We define H_T,X the partial derivative of the source term with respect to temperature and composition, respectively, and we assume for simplicity (and for the linear stability analysis in Appendix A) that H_T,X is constant at the interface.

The classic Rayleigh-Taylor analysis predicts no instability in the absence of density discontinuities. However, following Tremblin et al. (2019), we show in Appendix A a linear stability analysis of the diabatic Rayleigh-Taylor instability in the Boussinesq regime, which gives a radiative Rayleigh-Taylor instability for the system (1): $${\omega }^2=\frac{gk}{2}\frac{H_X}{T_0{H}_T}\left(X_0^{-}-X_0^{+}\right),$$(3)

with k as the horizontal wavelength of the perturbation and ω the growth rate of the instability. Equation (3) is a simplified version of a more general diabatic Rayleigh-Taylor growth rate given in appendix (see Eq. (A.30)). The effect of a stable stratification in both side of the interface is neglected here but can be taken into account in the more general growth rate. We refer to the appendix for more details.

We highlight that this growth rate is the discontinuous version of the continuous case studied in Tremblin et al. (2019), with the criterion for diabatic convection without mean-molecular-weight gradient nor compositional source terms: ∇_X H_X > 0 with ∇_X = ∂logX∕∂logP. Qualitatively, the instability in the continuous and discontinuous cases are easy to understand and rely on buoyancy: when a parcel of fluid moves upward, if a compositional change during the motion induces heating, the temperature increases and the density decreases, which is destabilizing relative to the buoyancy. On the contrary, if the compositional change induces cooling, the temperature decreases and the density increases, which is stabilizing for an upward motion. To trigger the instability in a stratified medium, the impact of the source term has to induce heating for upward motions and cooling for downward motions and also be sufficiently large to overcome the stabilizing effect of a stable stratification.

Remarkably, we get an instability without a discontinuity of density at the interface, but with a heating source term that depends on composition with a discontinuity of composition. This situation is typical of cloud interfaces in planetary atmospheres and, therefore, we explore this instability with numerical simulations in brown dwarf and rocky exoplanet regimes.

3 Numerical method and setup

We use the code ARK (Padioleau et al. 2019) to solve the hydrodynamics equations in two dimensions (x and y) with a radiative source term as presented in Tremblin et al. (2019). In this paper, we ignore changes in the mean molecular weight and we keep it constant. Then we trace the opacity of the medium with a scalar field, X. The equations solved in this setup are the system Eq. (1) with $$H(X,T)=\frac{4\pi \kappa (X)}{c_p}(J-\sigma T^4/\pi ),$$(4)

which is equivalent to the divergence of the radiative flux (see Mihalas & Mihalas 1984). The opacity is κ(X) = Xκ₁ + (1 − X)κ₂, tracing gases with two different opacities, κ₁ and κ₂. The mean gray intensity is J = (I^↑ + I^↓)∕2 with I^↑ and I^↓ the upward and downward intensities computed using the radiative transfer equation in a standard two-stream approximation (see Tremblin et al. 2019, for details). Here, we only treat absorption opacities and scattering is ignored.

The hydrodynamics solver is a well-balanced and all-regime solver extensively described and tested in Padioleau et al. (2019). These properties of the numerical method are essential for this study. The well-balanced property means that the method is able to maintain the hydrostatic balance at machine precision. This gives us a very precise control of the instability: even when it is initialized on an unstable equilibrium, the simulation does not develop velocities unless we explicitly add a perturbation. The all-regime property means that the solver has a low-Mach correction to reach a high accuracy in the regime of low velocities. This low-Mach correction is activated in all the simulations presented in this paper and is also essential for capturing the instability in the regime we explore here.

There is no diffusivity, viscosity, or sub-grid turbulence in the model. The dissipation relies only on the numerical diffusion at the grid scale of the all-regime numerical method. Thanks to the low-Mach correction, that dissipation is significantly reduced, as shown in Padioleau et al. (2019).

The initial conditions are all chosen in a similar way: we initialized the scalar tracing opacities X = X₀ to 0 and 1 or vice versa in the upper and lower half of the domain. We then pre-compute a pressure-temperature profile that satisfy the discrete version of the hydrostatic balance and energy conservation in Eq. (2), starting from an imposed base pressure P_bot. As explained above with the well-balanced property, the simulation will remain static in the absence of perturbations even if these initial conditions are unstable.

The boundary conditions for the density and pressure are imposed by a linear extrapolation of the temperature and an extrapolation of the hydrostatic balance at the top and bottom boundaries of the domain in the y direction. In addition, we impose zero velocities. For the radiative variables, we impose the downward radiative flux at the top $\pi I_{\text{top}}^{\downarrow }$ and the net radiative flux at the bottom $\pi (I_{\text{bot}}^{\uparrow }-I_{\text{bot}}^{\downarrow })$. The downward radiative flux at the top represents the emission of the layers of the atmosphere that are above our computational domain. Imposing the net radiative flux at the bottom implies that the downward intensity at the bottom is re-emitted upward by deeper layers or a planetary surface. The magnitude of these fluxes are chosen in conjunction with the gray opacities, κ₁ and κ₂, to ensure a realistic temperature profile, furthermore the adiabatic index, γ, is adjusted such that this idealized setup is stable with regard to Schwarzschild convection, and allows for the study of other instabilities. The boundary conditions in the x direction are periodic. The initial perturbation is adapted differently in the brown dwarf and Earth-like regimes; thus, we provide details on this choice in Sects. 4 and 5.

4 RRTI in the brown dwarf regime

We first choose our parameters to model a stably stratified part of the upper atmosphere of a L dwarf with effective temperature T_eff ~ 1600 K and surface gravity log(g) = 5. All the parameters including the opacities and radiative boundary conditions are given in Table 1 and are set to reproduce a temperature gradient decreasing with height stable to Schwarzschild convection. We assume the presence of an opacity interface e.g. coming from condensation of silicates or iron, with an opacity ratio of a factor of 10. We use a perturbation in the initial vertical velocity: $$u_{0,y}(x,y)=A{c}_{\text{s}}\sin (m\pi x/x_{\text{max}})e^{-{(y-0.5y_{\text{max}})}^2/{(w{y}_{\text{max}})}^2},$$(5)

with c_s the local sound speed, x_max, y_max the horizontal and vertical extent of the box, respectively.

Figure 1 shows the maps of the opacity tracer at t = t_end when we start from the higher-opacity medium on top of the lower-opacity one (left panel) and vice versa (right panel), clearly demonstrating that higher-opacity over lower-opacity is stable while lower-opacity over higher-opacity is unstable. We highlight that higher-opacity versus lower-opacity only refers to κ₂ versus κ₁ and not to an optically thick versus thin transition: in this setup, both medium are optically thin and the total optical depth after crossing both media is 0.18. This instability can be directly linked to the dependence of the opacity to a local tracer given in Sect. 2.

We show in the top panel of Fig. 2, the initial and final profile of potential temperature $\theta =T{(P_{\text{ref}}/P)}^{(\gamma -1)/\gamma }$ in the unstable simulation. It shows that the potential temperature is always increasing with height, hence, the simulations are always stably stratified, namely, they are stable with regard to thermal convection. The bottom panel shows the evolution of the averaged internal, gravitational, and kinetic energy in the simulation. We can see that the kinetic energy is equal to the opposite of the gravitational energy: this is expected for a buoyancy instability that converts gravitational potential energy into kinetic energy. The evolution of the internal energy is not conservative, which is expected because of the radiative transfer source term. Its evolution leads to conditions prone to the RRTI instability, which we will explore in more detail below.

In order to explore the link to RRTI, we show in Fig. 3 the heating rate computed from our initial conditions when the interface is artificially displaced by 10% up or down in the box. Since the heating rate is exactly zero when the interface is at the middle of the box, this artificial displacement allows us to probe numerically the evolution of H with composition and get an estimate of H_X in our simulations. As explained above, we can expect the flow to be unstable to buoyancy when an upward displacement leads to heating and a downward displacement to cooling, while upward displacements inducing cooling would correspond to a stable situation. Indeed, Fig. 3 shows that vertical velocities and the instability appear when the displacement of the interface is unstable to buoyancy. We can provide a theoretical estimate of the growth rate using Eq. (3) and by numerically computing H_X,T at the interface from the initial conditions. Such an estimate gives a growth rate of 0.04 s⁻¹ while the measured growth rate in the simulation is 0.07 s⁻¹ (70% higher, measured between t = 10 and t = 50 s before reaching saturation). The theoretical estimate provides the correct order of magnitude and the agreement can be considered as relatively good given the approximations in the linear stability analysis that partly ignore the complexity and non-locality of radiative transfer.

Figure 3 also shows that the displacement of the interface leads to a non-local heating and cooling preferentially in the higher-opacity medium. Such a non-locality is not surprising: when the medium is optically thin, a local modificationof its opacity has an impact in the entire domain. Heating and cooling happens preferentially in the higher-opacity domain simply because the heating rate is proportional to the opacity (see Eq. (4)). The mean velocity profile shows that motions are indeed appearing preferentially in the entire higher-opacity medium and not only at the interface. This non-locality and asymmetry between the two media can explain the rapidly non-linear and strongly asymmetric development of the perturbations in Fig. 1: the downward column is narrow and fast while the upward motions are spatially large and slow.

In order to study this radiative heating and cooling more in depth, we need to explore the evolution of the upward and downward intensities as a function of the interface displacement. We simplify the problem by assuming that the radiative transfer equation leads to a simplified solution of the type $I(y)=e^{-\kappa \rho \Delta y}I(y_0)+(1-e^{-\kappa \rho \Delta y})\sigma T{(y)}^4/\pi ,$ with a vanishing opacity in the lower-opacity medium, namely, I(y) = I(y₀) and a small opacity in the higher-opacity medium, such that I(y) is in the linear regime as a function of y. At equilibrium in the initial condition, the radiative flux is constant and the heating rate is zero, this implies that I^↑ − I^↓ is constant and (I^↑ + I^↓)∕2 = σT⁴∕π, that is, the intensity profiles increase or decrease with altitude similarly to temperature. Under these assumptions, Fig. 4 shows the typical equilibrium downward and upward intensity profiles in gray when the temperature is decreasing with altitude. We detail the two different cases.

Higher-opacity over lower-opacity: starting from the top boundary, the downward intensity is increasing when going down, which corresponds to emission in the higher-opacity medium; it is then constant when propagating through the lower-opacity medium. Then, starting from the bottom boundary the upward intensity is constant through the lower-opacity medium, its value being fixed by the net radiative flux boundary condition, and then decrease in the higher-opacity medium which corresponds to absorption. When the interface is moving up (red curves), the optical path in the higher-opacity medium is decreasing, corresponding to a decrease in emission for the downward photons which decreases first the downward intensity and then the upward intensity because of the imposed net radiative flux at the bottom boundary. The decrease in both intensities hence results in a decrease in the mean gray intensity J. Using Eq. (4) together with the fact that T is constant indicates that there is a net cooling of the atmosphere when the interface is moving up, which has a stabilizing effect. Lower-opacity over higher-opacity: starting from the top boundary, the downward intensity is constant when going down in the lower-opacity medium, then increases in the higher-opacity medium, which corresponds to emission in the higher-opacity medium. Then starting from the bottom boundary the upward intensity is decreasing in the higher-opacity medium which corresponds to absorption and is constant through the lower-opacity medium. When the interface is moving up, the optical path in the higher-opacity medium is increasing, corresponding to an increase in emission for the downward photons which increases first the downward intensity and then the upward intensity because of the imposed net radiative flux at the bottom boundary. The increase in both intensities hence results in an increase in the mean gray intensity J. Using Eq. (4) together with the fact that T is constant, means there is a net heating of the atmosphere when the interface is moving up, which is destabilizing.

In summary, when the temperature is decreasing with altitude, a decrease in the optical path in the higher-opacity medium results in a decrease in emission hence radiative cooling and an increase in the optical path results in an increase in emission hence radiative heating. Interestingly, it can be expected from Fig. 4 that this behavior is inverted in the presence of a temperature inversion. Figure 5 shows the expected intensity profiles if the temperature is increasing with height. The same analysis can be repeated but with the higher-opacity medium absorbing the downward photons (since the intensity is decreasing along the path). Consequently, a decrease in the optical path in the higher-opacity medium results in a decrease in absorption (hence, radiative heating), while an increase in the optical path in the higher-opacity medium results in an increase in absorption (hence, radiative cooling). It is then clearthat we should expect the instability to be inverted in that case, namely: when the temperature is increasingwith altitude, a higher-opacity medium on top of a lower-opacity medium should be unstable and lower-opacity over higher-opacity should be stable. We explore this possibility in an Earth-like regime in the next section.

Fig. 1 Final two-dimensional (2D) maps of the opacity tracer in a brown-dwarf regime with a negative vertical temperature gradient. Simulations are started from: left, a higher-opacity medium (yellow) on top of a lower-opacity medium (dark blue). Right, a lower-opacity medium on top of an higher-opacity medium.

Fig. 2 Global properties of the simulations in the brown-dwarf regime. Top: initial and final potential temperature profiles for the simulation with a lower-opacity medium on top of a higher-opacity one. Bottom: evolution of the averaged kinetic, gravitational and internal energies during the simulation. We plot the evolution of the differences from their initial values.

Fig. 3 Brown-dwarf regime with a negative vertical temperature gradient. Top: higher-opacity over lower-opacity. Bottom: lower-opacity over higher-opacity. Left: heating rate profiles when the position of the interface is displaced up or down by 10% from the initial condition. Right: mean vertical velocity profile in the simulation at t = tend ∕2.

Fig. 4 Simplified downward and upward intensity profiles in an atmosphere with a negative vertical gradient of temperature. Grey profiles are those set for an initial condition at radiative equilibrium with the interface at the middle of the box. Red and blue profiles corresponds to intensities when the interface has been moved up and down.

Fig. 5 Simplified downward and upward intensity profiles in an atmosphere with a positive vertical gradient of temperature. Grey profiles are those set for an initial condition at radiative equilibrium with the interface at the middle of the box. Red and blue profiles corresponds to intensities when the interface has been moved up and down.

5 Earth-like regime: dependence on the temperature gradient

In order to explore the behavior of the RRTI with a temperature inversion, we performed simulations in an Earth-like regime for which a positive vertical gradients of temperature could arise either because of irradiation condensation and evaporation or hydrodynamical effects. We therefore need to impose a temperature inversion by adding a forcing term in the heating source term $$H(X,T)=\frac{4\pi \kappa (X)}{c_p}(J-\sigma T^4/\pi )-\frac{T-T_{\text{forcing}}}{{\tau }_{\text{forcing}}},$$(6)

with T_forcing an imposed profile with a positive vertical gradient and τ_forcing the timescale for the forcing. Prior to this, we first considered a purely radiative setup (i.e., τ_forcing →∞) and adjusted the parameters (see Table 2) to reproduce the unstable behavior in a negative vertical gradient of temperature (see Fig. 6). The radiative balance in Earth atmosphere is typically in the infrared at wavelength around 10 μm. We assume a ratio of opacity of 1000 for Earth clouds (Kokhanovsky 2004) and adjust the radiative boundary conditions to get a temperature-pressure profile that matches the international standard atmosphere (ISA) values at 2 km altitude and is stable to Schwarzschild convection. The instability appears to be weaker in the Earth-like regime than in the brown-dwarf regime, a velocity perturbation in that context tends to trigger sound waves that are reflected on the boundaries of the domain and strongly interfere with the interface. To overcome this numerical limitation, we have used an interface perturbation rather than a velocity perturbation: $$y_{\text{int}}(x)=0.5y_{\text{max}}(1+A\sin (m\pi x/x_{\text{max}})).$$(7)

This type of perturbation appears to have a much weaker interaction with the boundaries of the simulation. Figure 6 shows that we recover the behavior expected from the brown-dwarf setup: higher-opacity over lower-opacity is stable and lower-opacity over higher-opacity is unstable. Figure 7 shows the behavior of the heating rate when we displace the interface in the initial conditions and the velocity profile at the middle of the simulation. The heating rate behaves similarly to the brown-dwarf regime: in the lower-opacity over higher-opacity case, an upward displacement of the interface results in radiative heating and a downward displacement to radiative cooling, hence, the interface is unstable to buoyancy. The difference with the brown-dwarf regime is that the opacity jump is larger and the higher-opacity medium is optically thick. Indeed we can also see in Fig. 7 (left) that heating and cooling are localized close to the interface which is a consequence of optically-thick radiative transfer that degenerates to local thermal diffusion. The right panel also shows that the motions tend to be localized closer to the interface. This difference does not impact the general behavior of the instability.

In order to explore the behavior of the instability with a temperature inversion, we used a simple Newtonian forcing as proposed in Eq. (6). This Newtonian cooling is a very simplified model to force a temperature inversion that could arise from irradiation, phase change, or dynamics. The initial profile is computed to be at equilibrium H = 0 for Eq. (6), that is, including Newtonian cooling. We point out that the argument (Sect. 4) meant to explain the instability in the presence of a temperature inversion still applies here because Newtonian cooling depends only on temperature. We use a linear forcing profile with altitude and parameterize the top and bottom forcing temperature (see Table 3). We also choose a forcing timescale of the form τ_forcing = τ_ref × (κ_ref∕κ(X)). This form allows us to invert the higher-opacity and lower-opacity medium in the simulation and keep a continuous temperature profile stable to Schwarzschild convection. It is also reasonable to assume that irradiating and dynamical forcing will be more efficient in the higher-opacity medium. All the parameters for the forcing are listed in Table 3 and we show in Fig. 8 the opacity tracer at the end of the simulation. As expected from Fig. 5, the instability is inverted: higher-opacity medium over lower-opacity medium is unstable and lower over higher-opacity is stable when the temperatureis increasing with altitude. We also show in Fig. 9 the profiles of the heating rate when the interface is displaced upward and downward. It shows that in the higher-opacity over lower-opacity case, an upward displacement results in radiative heating and a downward displacement in radiative cooling and confirm the instability regime. As in Fig. 7, we can also see that heating and cooling are localized close to the interface because of the optically thick regime in the higher-opacity medium similarly to the motions triggered by the instability.

We schematically summarize in Fig. 10 all the different possibilities of the stable and unstable regimes assuming the presence of a higher-opacity cloud layer in an atmosphere with different temperature structures. Essentially, the top part of the cloud layer is unstable if the temperature decreases with altitude, stable with a temperature inversion; whereas the base is stable if the temperature decreases with height and unstable with a temperature inversion. This might lead to the cloud cover being patchy. Interestingly, there is only one case in which a cloud cover would be stable: negative temperature gradient at the base and a temperature inversion at the top. We discuss in the next section the possible implications for cloud covers in different objects.

Fig. 6 Final 2D maps of the opacity tracer in an Earth-like regime with a negative vertical temperature gradient. Simulations are started from: (left) a higher-opacity medium (yellow) on top of a lower-opacity medium (dark blue); (right) a lower-opacity medium on top of a higher-opacity medium.

Fig. 7 Earth-like regime with a negative vertical temperature gradient. Top: higher-opacity over lower-opacity. Bottom: lower-opacity over higher-opacity. Left: heating rate profiles when the position of the interface is displaced up or down by 10% from the initial condition. Right: mean vertical velocity profile in the simulation at t = tend ∕2.

Fig. 8 Final 2D maps of the opacity tracer in an Earth-like regime with a positive vertical temperature gradient. Simulations at left: higher-opacity medium (yellow) on top of a lower-opacity medium (dark blue). Right: lower-opacity medium on top of a higher-opacity medium.

Fig. 9 Earth-like regime with a positive vertical temperature gradient. Top: higher-opacity over lower-opacity. Bottom: lower-opacity over higher-opacity. Left: heating rate profiles when the position of the interface is displaced up or down by 10% from the initial condition. Right: mean vertical velocity profile in the simulation at t = tend ∕2.

Fig. 10 Schematics of the possible stable and unstable situations for a higher-opacity cloud layer depending on the temperaturestructure of the atmosphere.

6 Interpretations of the structure of cloud covers

6.1 Cold ice- and gas-giant planets

The first important distinction here is made by checking if clouds are expected to be radiatively active or passive, meaning that the radiative timescale is smaller or larger than the advective timescale. If the radiative timescale is small, the RRTI will develop quickly compared to other dynamical effects and will strongly impact the shape of the clouds. If this timescale is long, wind and convection will impact the higher-opacity layer before the radiative instability develops and clouds can be seen as radiatively passive scalar following the dynamics (ignoring for the moment condensation or evaporation; see Sect. 7). We estimate in Fig. 11 the radiative timescale and advective timescale for the cloudy Solar System planets, the extrasolar giant planet HD 209458b, and the T dwarf 2M 1047. We estimate the radiative timescale with ${\tau }_{\text{rad}}=c_{\text{p}}/\kappa \sigma T_{\text{1bar}}^3$, assuming a cloud opacity at κ = 1 cm² g⁻¹ and a temperature from the pressure/temperature structure at 1 bar. The global advective timescale is estimated with τ_adv = R_p∕u_wind with R_p the radius of the object and u_wind the typical wind velocities in the atmosphere. The wind velocity for 2M 1047 is taken from the recent wind measurement in Allers et al. (2020). All the values used to get these estimations are given in Table 4. We emphasize that this is an order-of-magnitude estimation to assess roughly the relative importance between radiative and dynamical effects.

We can see that the cold ice- and gas-giant planets in the Solar System are essentially in the regime in which the advective timescale is smaller than the radiative timescale which implies that clouds are essentially radiatively passive and follow the flow as a passive scalar. The radiative timescale is much smaller for the giants than for Earth and Venus because they are colder, but also because the specific heat capacity in these objects is around one order of magnitude larger. As a result, the radiative timescale is almost two orders of magnitude larger. Surprisingly, HD 209b and the T dwarf 2M 1047 fall in a similar regime as Earth and Venus, mainly because the radiative timescale is small as a result of the high temperature that characterizes these objects. Although the Solar System’s giant planets are often seen as our closest analog to hot giant exoplanets and brown dwarfs, they seem to be bad proxies with regard to understanding the cloud dynamics on these objects compared to Earth and Venus, at least as far as the RRTI is concerned. Similarly to the Solar System giants, we can expect that cold exoplanets and cold Y dwarfs, such as Wise 0855, will also be shown to have relatively radiatively passive clouds.

Fig. 11 Comparison between the radiative and advective timescale for HD 209458b, the T dwarf 2M 1047, Venus, Earth, Jupiter, Saturn, Uranus, and Neptune.

6.2 Earth and Venus

The dynamics of clouds in the atmosphere of Earth and Venus provides a possible test bench for a mechanism such as the RRTI by comparing the correlation between stability or instability and temperature gradients to Fig. 10.

Mammatus clouds in Earth’s atmosphere are opaque lobes hanging at the base of a cloud, typically beneath the anvil part of cumulonimbus when they overshoot the tropopause or also at the base of thunderstorm clouds (see Fig. 12 and Winstead et al. 2001). The understanding of the formation mechanism of these unusual cloud structures remains a challenge and many processes have been proposed so far based on evaporative cooling, cloud-base radiative heating, dynamical instabilities (Kelvin-Helmholtz, “standard” Rayleigh-Taylor) and many others (Schultz et al. 2006, 2007; Garrett et al. 2010). However, none of these processes seem to provide a definitive answer since they do not seem to always lead to the formation of mammated clouds. However, among all the observed properties of the mammatus, the correlation between their formation and the presence of a temperature inversion is certainly the one that is shared most conspicuously. This correlation was reported as early as 1911 by Clayton, using kite soundings at the Harvard College Astronomical Observatory.

Another interesting case is the main cloud cover in the atmosphere of Venus. Since a stable cloud deck is present, preventing us from seeing the surface at visible wavelength, the RRTI would predict this cloud cover to be stable only if a temperature inversion is present at the top of the cloud deck. We provide in Fig. 13, the pressure/temperature profiles measured by two of the four Pioneer probes (Seiff 1983) and the location of the main cloud deck (Formisano et al. 2006). It does appear that a temperature inversion is present at the location of the top boundary of the clouds explaining the stability to the RRTI and the formation of the cloud cover on the dayside of the planet. Yet such temperatureinversions are not possible without irradiation on the nightside of the planet. We would therefore expect the top cloud cover to be destabilized there. This is also what has been observed by Venus Express with relatively recent analysis of the cloud dynamics on the nightside of the planet (Peralta et al. 2017). While the cloud cover is smooth and follows super-rotation on the dayside (which is well reproduced by global circulation models (GCM) Lebonnois et al. 2016; Garate-Lopez & Lebonnois 2018), more complex dynamics and cloud structure with stationary wave patterns seems to be present on the nightside and is not well reproduced by GCM simulations.

Interestingly, both types of phenomena seem to be qualitatively in agreement with our predictions for the behavior of the RRTI as a function of the temperature gradient (see Fig. 10): the base of a cloud cover would be unstable to the RRTI when a temperature inversion is present (similar to the case of Mammatus) and the top of the cloud cover would be stable with a temperature inversion (dayside of Venus), whereas it would be unstable with temperature decreasing with height (nightside of Venus). More detailed studies are needed to confirm this link, such as with local convective simulations, as in Lefèvre et al. (2017, 2018), where the RRTI mechanism could be identified and characterized.

Fig. 12 Appearance of mammatus clouds over Montparnasse tower in Paris on the 21st of November 2016.

Fig. 13 Pressure and temperature profiles of the atmosphere of Venus measured by the Pioneer probes (Seiff 1983) and location of the main cloud deck (Formisano et al. 2006).

7 Discussion and conclusions

7.1 Speculations for exoplanets

Based on the insights provided by the RRTI and its possible role in the dynamics of Earth and Venus clouds, we can provide expectations for the cloud cover of irradiated rocky and giant exoplanets. Similarly to the dayside of Venus, we could expect a stable cloud cover to be possible only when irradiation create a temperature inversion at the top of the clouds (this inversion can be caused by the absorption of stellar light by the clouds themselves). If the cloud cover can grow to a large vertical extension on the dayside, it can then be advected on the nightside and survive the instabilities created at the top by the disappearance of the temperature inversion and consequently cover the entire planet. This possibility will strongly depend on the rotation period and wind velocity: if locally the clouds are submitted to a rapid day-night forcing they could be rapidly destabilized before growing significantly at the planetary scale, while for slow rotation or wind speed they could grow significantly on the dayside and survive on the nightside. Our mechanism may contribute to shaping certain cloud properties on Earth and Venus, and on irradiated exoplanets in general, although multiple mechanisms are crucial in determining the weather system on Earth and Venus.

For isolated objects such as brown dwarfs or low-irradiated hot exoplanets as the type usually observed by direct imaging, the formation oftemperature inversions appears difficult and is not really expected to take place. We therefore expect the top of a cloud cover to be unstable to the RRTI everywhere in the atmosphere of these objects. As a consequence of the RRTI, patchy cloud coversmay be ubiquitous in both L and T dwarfs and not only at the L/T transition. To illustrate this point, we show in Fig. 14 a simulation with an unstable cloud layer leading to patches of opacity. We perform this simulation on longer time-scale to probe the saturated steady state, the solution does not evolve much between 150 and 300 s, which is, therefore, much longer than the turnover timescale associated with the initial transitional phase. In that context, 1D atmospheric models with homogeneous cloud cover may not be a sufficiently realistic approximation.

Fig. 14 Opacity tracer in a simulation with an unstable cloud cover. Top panel: initial state. Middle: at t = 150 s. Bottom: at t = 300 s.

7.2 Conclusions

By using an analytical stability analysis and 2D radiative-hydrodynamical simulations, we have shown the existence of an instability of an opacity discontinuity between a lower-opacity and a higher-opacity medium: the Radiative Rayleigh-Taylor instability, namely, a particular case of the general diabatic Rayleigh-Taylor instability.

• In an atmosphere with a negative vertical gradient of temperature, a lower-opacity medium on top of a higher-opacity medium is unstable and a higher-opacity medium on top of a lower-opacity medium is stable.

• In an atmosphere with a positive vertical gradient of temperature, the behavior is inverted: a lower-opacity medium on top of a higher-opacity medium is stable and an higher-opacity medium on top of a lower-opacity medium is unstable.

• Applied to a higher-opacity cloud layer, this mechanism predicts that the base of a cloud can be unstable in the presenceof a temperature inversion and a cloud cover can only be stable if a temperature inversion is present at the top.

This mechanism could shed some light on the interpretation of several cloud structures in different objects:

• The radiative timescale of the ice- and gas-giant Solar System planets and cold exoplanets and Y dwarfs is long. In that context, clouds can be expected to be radiatively passive and they essentially follow the dynamics as long as only radiative transfer is concerned.

• The RRTI predicts that the base of a cloud is unstable with a temperature inversion, this could offer a possible mechanism to explain the formation of mammatus clouds in Earth’s and Earth-like atmospheres. Furthermore, a cloud cover can only be stable with a temperature inversion at the top, which seems to be the case for Venus dayside temperature profiles. By generalizing this finding with regard to irradiated exoplanets, we may expect stable large-scale cloud covers only in irradiated planets with slow rotation or wind speed, so that the cloud cover can grow on large scales on the dayside (and possibly be advected on the nightside).

• Isolated and low-irradiated objects like brown dwarfs and hot exoplanets observed by direct imaging are not expected to have temperature inversions. In that context, patchy cloud covers may be ubiquitous in their atmospheres.

We have, of course, used a very simplified setup in this paper. Some limitations regarding our current approach are:

• We partly neglected condensation or evaporation leading to the formation or destruction of the higher-opacity material and irradiation. We only explored with a simplified thermal forcing the possible formation of a temperature inversion because of evaporative cooling or irradiation. As a consequence, we have ignored the feedback of the instability on the forcing causing the temperature inversion.

• We used a simplified approach for radiative transfer: a gray model neglecting scattering. We also used the two-stream approximation in column, which neglects lateral radiative coupling.

• We used 2D simulations: further studies ought to consider this instability in 3D and assess its impact on a opaque layer of finite thickness. Longer timescales are also needed to probe the turbulent steady state.

Some of these limitations might not be important, that is, condensation for the formation of the higher-opacity material might happen in a first phase and be relatively negligible in the following evolution governed by radiation hydrodynamics. Scattering will probably not have much impact when the radiative balance is mainly in the infrared, and in any case, it does not impact radiative heating and cooling. However, a better modeling of evaporation, irradiation, and 3D radiative transfer (e.g., with the M1 model González et al. 2007) are certainly the next steps toward a better understanding of the RRTI and cloud dynamics in the context of exoplanets and brown dwarfs. With the upcoming arrival of the James Webb Space Telescope and observational data with a large spectral coverage, these types of models might be essential for the understanding of clouds in the atmosphere of giant and rocky exoplanets.

Acknowledgements

The authors are thankful to the anonymous referee for his/her valuable comments that help improving this manuscript. P.T. acknowledges supports by the European Research Council under Grant Agreement ATMO 757858. P.T. also thanks A. Barazzutti for providing Fig. 12, for her natural curiosity about natural phenomena without which this work would not have been possible, and her unconditional support over the years.

Appendix A Linear stability analysis of the diabatic Rayleigh-Taylor instability

A.1 General equations

For the linear stability analysis, we will use the Boussinesq approximation, following Tremblin et al. (2019). We decompose all the fields as e.g. ρ(t, x, y, z) = ρ₀(z) + δρ(t, x, y, z) and the background state is given by the static equations: $$\begin{array}{lllllllll}\frac{\partial P_0}{\partial z}\hfill & =\hfill & -{\rho }_{0}g\text{ cr}R(X_0,T_0)\hfill & =\hfill & 0\text{ cr}u\hfill & =\hfill & 0\text{ cr}H(X_0,T_0)\hfill & =\hfill & 0,\hfill \end{array}$$(A.1)

with R(X, T) the source term in the advection/reaction equation of the fluid concentration, X, and ρC_pH(X, T) the source term in the energy equation i.e. thermal conduction or radiative transfer heating rate.

The linearized Boussinesq approximation leads to the following system: $$\begin{array}{lllllllll}\nabla \left(\delta u\right)\hfill & =\hfill & 0\text{ cr}\frac{\partial \delta X}{\partial t}+\delta u\cdot \nabla \left(X_0\right)\hfill & =\hfill & \frac{\partial R}{\partial X}\delta X+\frac{\partial R}{\partial T}\delta T\text{ cr}\frac{\partial {\rho }_0\delta u}{\partial t}+\nabla \left(\delta P\right)-\delta \rho g\hfill & =\hfill & 0\text{cr}\frac{\partial \delta T}{\partial t}+T_0\delta u\cdot \nabla \log {\theta }_0\hfill & =\hfill & \frac{\partial H}{\partial X}\delta X+\frac{\partial H}{\partial T}\delta T,\hfill \end{array}$$(A.2)

with θ₀ the potential temperature $T_0{(P_{\text{ref}}/P_0)}^{(\gamma -1)/\gamma }$ (P_ref is a reference pressure and γ the adiabatic index). The system is closed with the equation of state for an ideal gas (by removing δP in the Boussinesq limit): $$0=\frac{\delta \rho }{{\rho }_0}+\frac{\delta T}{T_0}-\frac{\partial \log {\mu }_0}{\partial X}\delta X$$(A.3)

and we consider in the rest of the analysis an interface between two sub-domains located in z = 0. Before going to the Boussinesq regime, we recall first the demonstration in the incompressible regime, without composition nor source terms in order to highlight the modification of the demonstration for the Boussinesq limit.

A.2 Incompressible limit

In the incompressible case, the system is reduced to: $$\begin{array}{l}\nabla \left(\delta u\right)=0\\ \frac{\partial {\rho }_0\delta u}{\partial t}+\nabla \left(\delta P\right)-\delta \rho g=0.\end{array}$$

And we can introduce an extra equation on the advection of density: $$\frac{\partial \delta \rho }{\partial t}+\delta u\nabla {\rho }_0=0.$$(A.6)

Then we assume the form exp(ωt + i(k_xx + k_yy)) for the disturbance: $$i{k}_x\delta u+i{k}_y\delta v+{\partial }_z\delta w=0\text{ cr}\omega {\rho }_0\delta u+i{k}_x\delta P=0\text{cr}\omega {\rho }_0\delta v+i{k}_y\delta P=0\text{cr}\omega {\rho }_0\delta w+{\partial }_z\delta P+\delta \rho g=0\text{cr}\omega \delta \rho +\delta w{\partial }_z{\rho }_0=0,$$(A.7)

which is reduced to (defining $k^2=k_x^2+k_y^2)$: $$k^2\delta P+\omega {\rho }_0{\partial }_z\delta w=0\text{ cr}\omega {\rho }_0\delta w+{\partial }_z\delta P+\delta \rho g=0\text{cr}\omega \delta \rho +\delta w{\partial }_z{\rho }_0=0.$$(A.8)

We then need to find the differential equation for δw: $$\omega {\rho }_0\delta w-\omega {\partial }_z\left({\rho }_0{\partial }_z\delta w\right)/k^2+\delta \rho g=0,$$(A.9)

with the relation that links δρ and δw in the incompressible limit: $$\omega \delta \rho +\delta w{\partial }_z{\rho }_0=0.$$(A.10)

We explicitly keep the relation between δρ and δw here, because this is what is going to change in the Boussinesq limit with source terms. If we assume that in each subdomain, the density is constant, the differential equation is then $${\partial }_z^2\delta w-k^2\delta w=0,$$(A.11)

which has the solution δw(0^±)e^∓kz and assuming the continuity of velocity at the interface leads to δw(0^±) = δw(0).

We now replace δρ in the differential equation Eq. (A.9) and integrate around the interface between z = ±ϵ to get the jump relation: $${\omega }^2{k}^2{\displaystyle {\int }_{-ϵ}^{+ϵ}{\rho }_0}\delta w\text{d}z-{\omega }^2{[{\rho }_0{\partial }_z\delta w]}_{-ϵ}^{+ϵ}-g{k}^2{\displaystyle {\int }_{-ϵ}^{+ϵ}\delta }w{\partial }_z{\rho }_0\text{d}z=0.$$(A.12)

The first term can be decomposed as follow: $$\begin{array}{lllll}{\displaystyle {\int }_{-ϵ}^{+ϵ}{\rho }_0}\delta w\text{d}z\hfill & =\hfill & {\displaystyle {\int }_{-ϵ}^0{\rho }_0}\delta w\text{d}z+{\displaystyle {\int }_0^{+ϵ}{\rho }_0}\delta w\text{d}z\text{ cr}\hfill & \approx \hfill & -ϵ{\rho }_{-}\delta w(0)+ϵ{\rho }_{+}\delta w(0),\hfill \end{array}$$(A.13)

which tends to zero as ϵ tends to zero. This result holds for any term that cannot be written as a derivative of a discontinuous function: only those terms (i.e., Dirac functions) contribute to the jump relation. The second term can be written as follows: $$\begin{array}{lllll}{[{\rho }_0{\partial }_z\delta w]}_{-ϵ}^{+ϵ}\hfill & =\hfill & -k{\rho }_0(ϵ)\delta w(ϵ)e^{-kϵ}-k{\rho }_0(-ϵ)\delta w(-ϵ)e^{+k(-ϵ)}\text{ cr}\hfill & \to \hfill & -k\delta w(0)({\rho }_{-}+{\rho }_{+}),\hfill \end{array}$$(A.14)

δw can be removed from the integral in the third term since it is continuous at the interface and we can perfom the integration similarly to the second term to obtain the classical result: $${\omega }^2=gk\frac{{\rho }_{+}-{\rho }_{-}}{{\rho }_{+}+{\rho }_{-}}.$$(A.15)

We now introduce the Boussinesq limit, without composition nor source terms.

A.3 Boussinesq limit

In the Boussinesq limit, the differential equation in Eq. (A.9) is the same but the link between δρ and δw is different. Assuming no composition nor source terms, the Boussinesq limit adds the following relations: $$\begin{array}{lllll}\omega \delta T+\delta w{T}_0{\partial }_z(\log {\theta }_0)\hfill & =\hfill & 0\text{cr}\delta \rho /{\rho }_0\hfill & =\hfill & -\delta T/T_0,\hfill \end{array}$$(A.16)

which gives the following relation between δρ and δw: $$\omega \delta \rho /{\rho }_0=\delta w{\partial }_z(\log {\theta }_0),$$(A.17)

We note that the relation between δρ and δw is proportional to the instability criterion in the continuous case (∂_z logθ₀ < 0). We assume the log-density gradient and the log-potential temperature gradient constant in each subdomain. This assumption might not be true depending on the equation of state but we will assume that the wavelength of the perturbation is sufficiently small so that both gradients can be considered constant. The differential equation Eq. (A.9) with Eq. (A.17) then become: $$\begin{array}{lllll}{\omega }^2{\partial }_z({\rho }_0{\partial }_z\delta w)-{\omega }^2{k}^2{\rho }_0\delta w-g{k}^2{\rho }_0{\partial }_z(\log {\theta }_0)\delta w\hfill & =\hfill & 0\text{cr}{\partial }_z^2\delta w+{\partial }_z(\log {\rho }_0){\partial }_z\delta w-k^2(1+g{\partial }_z(\log {\theta }_0)/{\omega }^2)\delta w\hfill & =\hfill & 0,\hfill \end{array}$$(A.18)

whereby the solution is $\delta w(0^{\pm })e^{q^{\pm }z}$ with $$\begin{array}{ll}{q}^{\pm }\hfill & =-{({\partial }_z(\log {\rho }_0))}^{\pm }/2\hfill \\ \hfill & \mp \sqrt{{({\partial }_z(\log {\rho }_0))}^{\pm ,2}/4+k^2\left(1+g{({\partial }_z(\log {\theta }_0))}^{\pm }/{\omega }^2\right)}\hfill \end{array}.$$(A.19)

We note that q^± is independent of ω if the sub-domains are neutral to buoyancy (∂_z logθ₀ = 0).

To write the jump conditions we need to write the non-conservative product ρ₀ ∂_z logθ₀, with discontinuous variables only in the derivative (using the equation of hydrostatic equilibrium): $${\rho }_0{\partial }_z\log {\theta }_0=-{\partial }_z{\rho }_0-{\rho }_0^{2}g/(P_0\gamma ).$$(A.20)

The second term can be discontinuous but involves no derivatives, therefore, it does not contribute to jump relations (by integration of this term on each side of the interface). It is only the first term (equivalent to a dirac function) that contributes. The jump condition by integration of Eq. (A.18) between ± ϵ is then given by: $${\omega }^2{[{\rho }_0{\partial }_z\delta w]}_{-ϵ}^{+ϵ}+{[{\rho }_0]}_{-ϵ}^{+ϵ}\delta w(0)g{k}^2=0;$$(A.21)

hence, the growth rate in the Boussinesq limit is given by: $${\omega }^2=g{k}^2\frac{{\rho }_{+}-{\rho }_{-}}{{\rho }_{-}{q}^{-}-{\rho }_{+}{q}^{+}},$$(A.22)

We note that we recover the incompressible result if the log-density and log-potential temperature gradients are negligible in Eq. (A.19), that is, q^± = ∓k.

A.4 Diabatic Rayleigh-Taylor instability

Now we go back to the full system with composition and source terms. We go on to add to Eq. (A.9) the following equations: $$\begin{array}{lllllll}\omega \delta X+\delta w{\partial }_z{X}_0\hfill & =\hfill & R_X\delta X+R_T\delta T\text{cr}\omega \delta T+\delta w{T}_0{\partial }_z(\log {\theta }_0)\hfill & =\hfill & H_T\delta T+H_X\delta X\text{ cr}\delta \rho /{\rho }_0\hfill & =\hfill & -\delta T/T_0+{\partial }_X(\log {\mu }_0)\delta X.\hfill \end{array}$$(A.23)

In that case, $$\begin{array}{lllllllllll}({\omega }_T^{\prime }-\omega )\delta w{\partial }_z{X}_0\hfill & =\hfill & (H_T-\omega +H_X/(T_0{\partial }_X\log {\mu }_0))\text{cr}\hfill & \hfill & \times ((R_X-\omega )\delta X+R_T\delta T)\text{cr}({\omega }_X^{\prime }-\omega )\delta w{\partial }_z(\log {\theta }_0)\hfill & =\hfill & ((R_X-\omega )/T_0+R_T{\partial }_X(\log {\mu }_0))\text{cr}\hfill & \hfill & \times ((H_T-\omega )\delta T+H_X\delta X)\text{ cr}\delta \rho /{\rho }_0\hfill & =\hfill & -\delta T/T_0+{\partial }_X(\log {\mu }_0)\delta X,\hfill \end{array}$$(A.24)

with ${\omega }_X^{\prime }=R_X+R_T{T}_0{\partial }_X\log {\mu }_0$ and ${\omega }_T^{\prime }=H_T+H_X/(T_0{\partial }_X\log {\mu }_0)$. This implies the relation between δρ and δw: $$\begin{array}{lllllll}\delta \rho /{\rho }_0\hfill & =\hfill & \delta w({\alpha }_T{\partial }_z\log {\mu }_0-{\alpha }_X{\partial }_z(\log {\theta }_0))\text{cr}{\alpha }_T\hfill & =\hfill & ({\omega }_T^{\prime }-\omega )/((H_T-\omega )(R_X-\omega )-H_X{R}_T)\text{ cr}{\alpha }_X\hfill & =\hfill & ({\omega }_X^{\prime }-\omega )/((H_T-\omega )(R_X-\omega )-H_X{R}_T).\hfill \end{array}$$(A.25)

We then assume α_T and α_X constant in the whole domain. We can then carry out the following redefinition: $$\begin{array}{lllll}\log \psi \hfill & =\hfill & -{\alpha }_X\log {\theta }_0+{\alpha }_T\log {\mu }_0\text{ cr}\delta \rho /{\rho }_0\hfill & =\hfill & \delta w{\partial }_z(\log \psi ).\hfill \end{array}$$(A.26)

The differential equation Eq. (A.9) with Eq. (A.26) then becomes: $$\begin{array}{lllll}{\omega }^2{\partial }_z({\rho }_0{\partial }_z\delta w)-{\omega }^2{k}^2{\rho }_0\delta w-g{k}^2{\rho }_0\omega {\partial }_z(\log \psi )\delta w\hfill & =\hfill & 0\text{cr}{\partial }_z^2\delta w+{\partial }_z(\log {\rho }_0){\partial }_z\delta w-k^2(1+g{\partial }_z(\log \psi )/\omega )\delta w\hfill & =\hfill & 0.\hfill \end{array}$$(A.27)

In each sub-domain we can solve the differential equation in δw to get $\delta w(0^{\pm })e^{q^{\pm }z}$, with $$\begin{array}{lllll}{q}^{\pm }\hfill & =\hfill & -{({\partial }_z(\log {\rho }_0))}^{\pm }/2\text{cr}\hfill & \mp \hfill & \sqrt{{({\partial }_z(\log {\rho }_0))}^{\pm ,2}/4+k^2\left(1+g{({\partial }_z(\log \psi ))}^{\pm }/\omega \right)}.\hfill \end{array}$$(A.28)

Remarkably, we can get a discontinuity (and instability) in log ψ with a continuous density at the interface. We therefore assume ρ₀ to be continuous at the interface. We also need to employ this hypothesis in order to have a well-defined jump relation for the term ρ₀ ∂_z(logψ), with a discontinuous density as this is a non-conservative product that is a-priori not well defined. We highlight that this hypothesis asserts that the interface is neutral to the classic Rayleigh-Taylor instability, but does not prevent the medium from becoming stratified on each side of the interface (with non-zero ${\partial }_z{(\log {\theta }_0)}^{\pm }$ in Eqs. (A.26) and (A.28).

The jump condition, via the integration of Eq. (A.27) between ± ϵ, becomes: $$\omega {\rho }_0{[{\partial }_z\delta w]}_{-ϵ}^{+ϵ}-{[\log \psi ]}_{-ϵ}^{+ϵ}{\rho }_0\delta w(0)g{k}^2=0,$$(A.29)

which gives the diabatic Rayleigh-Taylor growth rate for a continuous density: $$\omega =g{k}^2\frac{\log {\psi }^{+}-\log {\psi }^{-}}{q^{+}-q^{-}}.$$(A.30)

We note that this growth rate takes into account the stabilizing effect of stratification in the medium ${\partial }_z{(\log {\theta }_0)}^{\pm }\ne 0$. In the limit of very large stratification $q^{\pm }\to \mp k\sqrt{g{({\partial }_z(\log \psi ))}^{\pm }/\omega }$, which gives the following limit of the growth rate: $${\omega }^{3/2}=\sqrt{g}k\frac{\log {\psi }^{+}-\log {\psi }^{-}}{\sqrt{{({\partial }_z(\log \psi ))}^{+}}+\sqrt{{({\partial }_z(\log \psi ))}^{-}}}.$$(A.31)

Within the limit of small source terms H_X,T, R_X,T → 0, we have ${({\partial }_z(\log \psi ))}^{\pm }\to ({({\partial }_z(\log {\theta }_0))}^{\pm }-{({\partial }_z(\log {\mu }_0))}^{\pm })/\omega$, which shows that ω → 0 when ∂_z(logθ₀)^)±→ +∞.

A.5 Radiative Rayleigh-Taylor instability

Assuming R = 0 and no mean molecular weight gradient, we get: $$\begin{array}{lllllll}\omega \delta X+\delta w{\partial }_z{X}_0\hfill & =\hfill & 0\text{cr}\omega \delta T+\delta w{T}_0{\partial }_z(\log {\theta }_0)\hfill & =\hfill & H_T\delta T+H_X\delta X\text{ cr}\delta \rho /{\rho }_0\hfill & =\hfill & -\delta T/T_0,\hfill \end{array}$$(A.32)

assuming only a discontinuous composition X, we get $$\begin{array}{lllll}\log \psi \hfill & =\hfill & -H_X{X}_0/(\omega T_0(H_T-\omega ))\text{cr}\delta \rho /{\rho }_0\hfill & =\hfill & \delta w{\partial }_z(\log \psi ),\hfill \end{array}$$(A.33)

which gives the growth rate, assuming continuous ρ₀: $${\omega }^2=g{k}^2\frac{H_X}{T_0(H_T-\omega )}\frac{X_0^{-}-X_0^{+}}{q^{+}-q^{-}}.$$(A.34)

We can simplify this expression assuming ω ≪ H_T and neglecting the log-density and logψ gradients such that q^± = ∓ k. In this limit, the radiative Rayleigh-Taylor growth rate is simply given by: $${\omega }^2=\frac{gk}{2}\frac{H_X}{T_0{H}_T}\left(X_0^{-}-X_0^{+}\right).$$(A.35)

A discontinuity of composition associated with a source term that depends on composition (i.e., a discontinuity in XH_X) can therefore lead to an instability similar to the adiabatic Rayleigh-Taylor instability even if the density is continuous at the interface.

References

All Tables

All Figures

	Fig. 1 Final two-dimensional (2D) maps of the opacity tracer in a brown-dwarf regime with a negative vertical temperature gradient. Simulations are started from: left, a higher-opacity medium (yellow) on top of a lower-opacity medium (dark blue). Right, a lower-opacity medium on top of an higher-opacity medium.
In the text

	Fig. 2 Global properties of the simulations in the brown-dwarf regime. Top: initial and final potential temperature profiles for the simulation with a lower-opacity medium on top of a higher-opacity one. Bottom: evolution of the averaged kinetic, gravitational and internal energies during the simulation. We plot the evolution of the differences from their initial values.
In the text

	Fig. 3 Brown-dwarf regime with a negative vertical temperature gradient. Top: higher-opacity over lower-opacity. Bottom: lower-opacity over higher-opacity. Left: heating rate profiles when the position of the interface is displaced up or down by 10% from the initial condition. Right: mean vertical velocity profile in the simulation at t = tend ∕2.
In the text

	Fig. 4 Simplified downward and upward intensity profiles in an atmosphere with a negative vertical gradient of temperature. Grey profiles are those set for an initial condition at radiative equilibrium with the interface at the middle of the box. Red and blue profiles corresponds to intensities when the interface has been moved up and down.
In the text

	Fig. 5 Simplified downward and upward intensity profiles in an atmosphere with a positive vertical gradient of temperature. Grey profiles are those set for an initial condition at radiative equilibrium with the interface at the middle of the box. Red and blue profiles corresponds to intensities when the interface has been moved up and down.
In the text

	Fig. 6 Final 2D maps of the opacity tracer in an Earth-like regime with a negative vertical temperature gradient. Simulations are started from: (left) a higher-opacity medium (yellow) on top of a lower-opacity medium (dark blue); (right) a lower-opacity medium on top of a higher-opacity medium.
In the text

	Fig. 7 Earth-like regime with a negative vertical temperature gradient. Top: higher-opacity over lower-opacity. Bottom: lower-opacity over higher-opacity. Left: heating rate profiles when the position of the interface is displaced up or down by 10% from the initial condition. Right: mean vertical velocity profile in the simulation at t = tend ∕2.
In the text

	Fig. 8 Final 2D maps of the opacity tracer in an Earth-like regime with a positive vertical temperature gradient. Simulations at left: higher-opacity medium (yellow) on top of a lower-opacity medium (dark blue). Right: lower-opacity medium on top of a higher-opacity medium.
In the text

	Fig. 9 Earth-like regime with a positive vertical temperature gradient. Top: higher-opacity over lower-opacity. Bottom: lower-opacity over higher-opacity. Left: heating rate profiles when the position of the interface is displaced up or down by 10% from the initial condition. Right: mean vertical velocity profile in the simulation at t = tend ∕2.
In the text

	Fig. 10 Schematics of the possible stable and unstable situations for a higher-opacity cloud layer depending on the temperaturestructure of the atmosphere.
In the text

	Fig. 11 Comparison between the radiative and advective timescale for HD 209458b, the T dwarf 2M 1047, Venus, Earth, Jupiter, Saturn, Uranus, and Neptune.
In the text

	Fig. 12 Appearance of mammatus clouds over Montparnasse tower in Paris on the 21st of November 2016.
In the text

	Fig. 13 Pressure and temperature profiles of the atmosphere of Venus measured by the Pioneer probes (Seiff 1983) and location of the main cloud deck (Formisano et al. 2006).
In the text

	Fig. 14 Opacity tracer in a simulation with an unstable cloud cover. Top panel: initial state. Middle: at t = 150 s. Bottom: at t = 300 s.
In the text