© The Authors 2023

1 Introduction

Mixing-length theory (MLT) for turbulent convection predicts efficient stirring by space-filling homogenous convection (e.g., Chang & Sofia 1987). On this basis, a common assumption in astrophysics and planetary science has been that convec-tive layers in stars and in fluid planets should be well mixed. Condensing species were expected to become uniformly mixed rapidly at depths where they are entirely in vapor form. However, in Jupiter’s deep atmosphere, ammonia appears to be generally depleted and variable in abundance down to at least four scale heights below its condensation level (Li et al. 2017; de Pater et al. 2019; Bolton et al. 2021). Water also has a variable abundance much deeper than its cloud base, as shown by Galileo probe measurements (Wong et al. 2004), infrared spectroscopy (Bjoraker et al. 2015, 2018), and indirect measurements from the Juno Microwave Radiometer (MWR; Li et al. 2020). In the Sun, the convective flow velocities measured by helioseismology appear to be poorly described by MLT (Hanasoge et al. 2012). The possibility that localized downdrafts maintain their coherence without mixing over significant distances has thus been proposed in both the Jovian (Guillot et al. 2020a) and solar (Anders et al. 2019) cases: In Jupiter, storms driven by the condensation of water have long been known to be a powerful source of energy, leading to spatially localized intermittent events that transport large amounts of energy in a short time (e.g., Gierasch et al. 2000). Storms yield precipitation in the form of rain, ice, or, for Jupiter, ammonia-water hail (“mushballs”). Upon evaporation, precipitation leads to a localized increase in mean-molecular weight and decrease in temperature, both effects strongly favoring the formation of a localized downdraft (Guillot et al. 2020b,a).

Regarding the Sun, there is a long history of modeling the effect of compressibility on convection (e.g., Hurlburt et al. 1984; Brummell et al. 2002), focused primarily on large-scale flows. However, more recent helioseismic measurements indicate flow velocities that are orders-of-magnitude smaller than what would be expected from MLT in order to be compatible with the Sun’s luminosity (Hanasoge et al. 2012). It has been suggested that a modification to MLT that includes an additional nonlocal mixing term called “entropy rain” can reconcile these observations with theory (Brandenburg 2016). Entropy rain refers to localized entropy anomalies that can sink through long distances, even through the whole convective region of the Sun, without mixing with their surroundings as MLT would ordinarily posit. Therefore, the total convective heat flux can be greater than the component of heat flux transported by the large-scale flow field to which helioseismic measurements are sensitive. Such a mechanism could reconcile the Hanasoge et al. (2012) observations with the Sun’s observed luminosity. Simulations of laminar downwelling thermals in a highly compressible environment find that thermals form coherent vortex rings that maintain their concentration and do not mix with their environment (Anders et al. 2019), offering a plausible mechanism for entropy rain. This phenomenon motivates more careful analysis of the effect of compressibility on small-scale downdrafts beneath the resolution of larger-scale simulations. We therefore compare our findings here, which consider the influence of turbulence, to this prior work that considers the laminar case.

Global circulation models of Jupiter that include sub-scale moist convection parametrized from the best available models of the Earth’s atmosphere have led to abundance variations that are modest and very rapidly negligible below the water condensation level (DEL Genio & McGrattan 1990). But giant planets bear at least two important differences from Earth: the absence of a surface and the low molecular weight of the dominant gas species. The absence of a surface forces us to seriously grapple with the chaotic dynamics of moist convection operating over multiple scale heights, a situation that never arises on Earth. Furthermore, the low molecular weights of hydrogen and helium imbue condensing species with an additional forcing mechanism, weighing down wet updrafts (see Guillot 1995) and imbuing wet downdrafts with greater power. Also, while the physics of compressible convection has been investigated in some detail in the solar case (e.g., Hurlburt et al. 1984; Brummell et al. 2002; Nordlund et al. 2009), the influence of intermit-tency has not received sufficient attention. We therefore seek to model meso-scale convection from isolated sources in abyssal adiabatic atmospheres to determine whether we should expect the rainfall from a concentrated, energetic storm to efficiently mix, or whether the rainfall can remain coherent and unmixed for some distance.

The topic is important: our Solar System has four planets with abyssal atmospheres, planets with substantial gas envelopes probably constitute the majority of planets in the Galaxy (e.g., Zhu & Dong 2021), and most of them have clouds (e.g., Helling 2019; Loftus et al. 2019). Furthermore, the physics we investigate here should also apply to the atmospheres of brown dwarfs (e.g., Helling & Casewell 2014). Storm activity is also routinely observed in all of the Solar System’s giant planets: thunderstorms have been extensively observed on Jupiter (e.g., Borucki et al. 1982; Little et al. 1999; Kolmasova et al. 2018; Becker et al. 2020); Saturn has an extreme case of decadally intermittent storm activity in its Great White Spot (Heath & McKim 1992; Fischer et al. 2011; Li & Ingersoll 2015); and large methane storms have been observed from Earth on Uranus (de Pater et al. 2015) and Neptune (Molter et al. 2019). The storms we can see from visible clouds may only scratch the surface; abyssal rock storms on Jupiter could be responsible for its seismicity (Markham & Stevenson 2018), and lightning strikes inferred from Voyager’s plasma wave instrument (Gibbard et al. 1999) suggest moist convection may be active at the water cloud layer on Neptune, hundreds of kilometers beneath the visible surface. Such ubiquity calls for a deeper theoretical understanding of abyssal convection, which appears to be key to properly con-textualizing the relationship between stellar and giant planet atmospheres and their deep interiors.

To attack this problem, we begin by introducing some simple intuitive models for negatively buoyant thermals in Sect. 2 that we can use as a heuristic basis, modifying existing self-similar entrainment theories from Turner (1973) to account for the tremendous compressibility of abyssal atmospheres. Armed with some basic intuition, we then use three-dimensional fluid dynamical simulations in Sect. 3 from Simulating Non-hydrostatic Atmospheres in Planets (SNAP; Li & Chen 2019) to account for the effect of self-induced turbulence of strongly negatively buoyant downdrafts. In Sect. 4 we summarize our results, including the computational performance and nondimensional scaling relationships that we then apply to particular celestial objects in Sect. 5. In Sect. 6 we address some of the deficiencies of our simplified model assumption; for example, we inspect the importance of our choice of initial conditions and the effect of environmental turbulence and vertical wind shear. We constrain the degree to which these complications affect our conclusions. Finally, in Sect. 7 we discuss our results in the context of Jupiter and beyond, outline the need for further observations of Uranus, and chart a course for future steps to more fully understand moist convection in abyssal atmospheres.

2 Analytical model

Our analytical model closely follows the classic work Turner (1973), which deftly addresses many of the main problems of the fluid dynamics of buoyancy from isolated sources. This framework provides not only a convenient heuristic, but also experimental corroboration and determination of free parameters that would be prohibitive to constrain on the basis of self-similar theory alone. We modified Turner’s self-similar equations to account for the compressibility of abyssal atmospheres, and used these results as a starting point to compare against more detailed, physics-based simulations. We modeled rainy down-drafts as downward-propagating thermals. Thermals have been observed to be a general feature of convection, constituting the smallest unit of self-organized convection. Following Turner (1963), we characterize the thermal as an entraining vortex ring with a definite size. A vortex ring is a flow structure characterized by a mean propagation velocity and circulation around an azimuthally symmetric, dynamically stable ring. The simplest analytical model for a vortex ring is Hill’s spherical vortex, whose stream function obeys $$\psi =\{\begin{array}{ll}-\frac{3w}{4}\left(1-\frac{4r^2}{b^2}\right)r^2{\sin }^2\theta \hfill & r\le b/2\hfill \\ \frac{w}{2}\left(1-\frac{b^3}{8r^2}\right)r^2\sin \theta \hfill & r\ge b/2\hfill \end{array},$$(1)

where w is the propagation velocity and b is the vortex ring diameter. Thermals have been shown to propagate as vortex rings in laboratory studies (e.g., Turner 1957; Shusser & Gharib 2000; Gao & Yu 2016), and cumulus storms on Earth likewise appear to resemble vortex rings (Wang & Geerts 2013). Figure 1 shows a schematic of the model, with streamlines of Hill’s spherical vortex shown in gray lines with important parameters labeled. We assume the thermal should remain self-similar as it propagates (Turner 1957). The central assumption of the analytical model is that the mean rate of inflow into propagating vortex ring should be proportional to its propagation velocity according to some dimensionless entrainment coefficient α.

We next modeled the evolution of the thermal as it propagates under these assumptions, with the additional assumption that the environment is of uniform density as well as the Boussinesq approximation (we relax these assumptions later). Under our assumptions the entrainment rate is $$\frac{\text{d}m}{\text{d}t}=\frac{\pi }{2}{b}^2\alpha \rho w,$$(2)

where w is the propagation velocity, b is the diameter of the vortex ring, and m is the mass enclosed within the Hill’s spherical vortex. The entrainment coefficient has been empirically determined to match laboratory results for thermals when α = 1/4 (Turner 1973), but it can take other values in different physical systems and in general is an unconstrained scaling coefficient of order unity that must be empirically determined. Under these assumptions, we can write the evolution equation for the vortex ring size as it propagates, $$\frac{\text{d}{b}^3}{\text{d}t}=3\alpha b^{2}w.$$(3)

From Eq. (3), we can immediately see that b ∝ αz, where w = dz/dt. In other words, db/dz = a. So according to this model, in an incompressible medium the thermal diameter increases linearly with propagation distance. Thus α is equivalent to the radian half-angle spread of the cone that encapsulates the entraining spherical vortex as it propagates.

Now we can solve for the equation of motion using Newton’s second law for momentum, so that $\frac{\text{d}}{\text{d}t}\left[mw\right]=g\overline{\delta \rho }\pi b^3/6$, where $m=\left(\rho +\overline{\delta \rho }\right)\pi b^3/6$ is the mass enclosed within the thermal. Then under the Boussinesq approximation, $$\frac{\text{d}\left(b^{3}w\right)}{\text{d}t}=b^3{g}^{\prime },$$(4)

where $g^{\prime }\equiv \frac{\delta \rho }{{\rho }_1}g$. In this model we assume gravity g to be constant, although g′ can vary as the density perturbation in the thermal varies.

Finally, we solved for the evolution of the buoyancy, assuming entrainment according to Eq. (2) and zero detrainment: $\frac{\text{d}{g}^{\prime }}{\text{d}t}=\frac{g}{\rho }\frac{\text{d}\left(\text{Δ}\rho \right)}{\text{d}t}=\frac{g}{\rho }\frac{\text{d}{\rho }_1}{\text{d}m}\frac{\text{d}m}{\text{d}t}$ where ρ₁ = ρ + Δρ. Then if the mass in the thermal is m = ρ₁V, m + dm = (ρ₁ + dρ₁)(V + dV) = (ρ₁ + dρ₁)(V + dm/ρ) if the pressure is the same inside and outside of the thermal. Solving for dρ₁ to first order in dm, we find $\frac{\text{d}{\text{g}}^{\prime }}{\text{d}t}=-\frac{3\alpha }{b}\frac{{\rho }_1}{\rho }{g}^{\prime }w$ from Eq. (2). Then under the Boussinesq approximation $b^3\frac{\text{d}{g}^{\prime }}{\text{d}t}=-g^{\prime }\frac{\text{d}{b}^3}{\text{d}t}$ from Eq. (3). Therefore, from the chain rule, $$\frac{\text{d}\left(b^3{g}^{\prime }\right)}{\text{d}t}=0.$$(5)

Intuitively, this means any volume increases due to entraining surrounding material are compensated by a corresponding decrease in Δρ so that the integrated total buoyancy anomaly is conserved. This result applies to the neutrally stratified medium. In general for a thermal propagating through a stratified medium, one must account for the Brünt-Väisälä frequency of the background material, discussed in Sect. 6.4.

We now modify these equations to account for the compressibility of an ideal gas adiabatic atmosphere. From hydrostatic equilibrium and the ideal gas equation of state, the density profile obeys $$\rho ={\rho }_0{\left(1+\frac{\gamma -1}{\gamma }\frac{z}{H_p}\right)}^{1/\left(\gamma -1\right)},$$(6)

where ρ₀ is the density at z = 0, γ is the Grüneisen parameter, and H_p = k_BT₀/µg is the corresponding pressure scale height. We now solve for the volume of the thermal as it propagates. Two processes affect the volume: entrainment increases the total mass of the thermal as it draws in fluid from its surroundings, while adiabatic compression increase its density. Expressed formally using V = m/ρ, we can write $\frac{\text{d\hspace{0.17em}ln\hspace{0.17em}}V}{\text{d}z}=\frac{\text{d\hspace{0.17em}ln\hspace{0.17em}}m}{\text{d}z}-\frac{\text{d\hspace{0.17em}ln\hspace{0.17em}}\rho }{\text{d}z}$ From Eq. (6), if we assume the pressure is equal inside and outside the thermal and that both the atmosphere and the thermal obey the ideal gas equation of state, then $\frac{\text{d\hspace{0.17em}ln\hspace{0.17em}}\rho }{\text{d}z}=\frac{\text{d}}{\text{d}z}\left[\frac{1}{\gamma -1}\ln \left(1+\frac{\gamma -1}{\gamma }\frac{z}{H_p}\right)\right]=\frac{1}{\gamma H_p+z\left(\gamma -1\right)}$ Therefore, using Eq. (2), $$\frac{\text{d\hspace{0.17em}ln\hspace{0.17em}}V}{\text{d}z}=\frac{3\alpha }{b}-\frac{1}{\gamma H_p+z\left(\gamma -1\right)}.$$(7)

Written another way, $$\frac{\text{d}b}{\text{d}z}=\alpha -\frac{b}{3H},$$(8)

where Η = γH_p + z(γ − 1) is the z-dependent density scale height. Written this way, we can clearly see how the compressibility of the atmosphere affects the thermal’s evolution. While entraining outside gas tends to increase the thermal mass and therefore size, the environment does work on the thermal to compress it and decrease its size. If b ≪ H, then this correction can be negligible. However, if b ~ H, the effect can be substantial and can even overwhelm the effect of entrainment such that db/dz < 0 so that the thermal shrinks rather than grows as it propagates. We plot the thermal size as a function of depth in Fig. 2. We see in the Boussinesq case (black curve) that the thermal grows in size linearly as it sinks. When 0 < b₀ ≪ H₀, this model remains accurate, with a small correction. When b₀ ~ H₀, the thermal can propagate some distance while maintaining its size and concentration; it may even shrink and grow more concentrated when b₀ > 3αH₀. In general for the atmosphere described in Eq. (6), the curves shown in Fig. 2 are concave up because Η increases with depth. So even for initially large b₀, the increase in Η always eventually dominates and leads to a growth and dilution of the thermal.

From Eq. (8) we directly see why atmospheric compression might keep downwelling thermals spatially localized. We can inspect the problem more completely by fully solving for the equations of motion. We seek to modify Eqs. (3)–(5) for a thermal falling through an adiabatic atmosphere described by Eq. (6). We write the dynamic evolution equations for b, w, and g′ expressed as coupled first order ordinary differential equations so that we can straightforwardly solve the system using a standard fourth-order Runge-Kutta integrator. From Eq. (8), $$\frac{\text{d}b}{\text{d}t}=w\left(\alpha -\frac{b}{3H}\right).$$(9)

Next we modify Eq. (4) to obtain an expression for w. We can use the same assumption of conservation of momentum to obtain an expression similar to Eq. (4), but without the Boussinesq approximation that discards terms ∝ g′/g we obtain the more exact expression $$\frac{\text{d}w}{\text{d}t}=g^{\prime }+\frac{\alpha w^2}{b}\left(4\frac{g^{\prime }}{g}-3\right).$$(10)

Equation (10) is equivalent to Eq. (4) when g′ ≪ g.

Finally, we modify Eq. (5) to obtain an expression for g′. One can follow the same derivation up to the final step before invoking the Boussinesq approximation, just allowing ρ to vary with p and assume $\frac{\text{d\hspace{0.17em}}\ln {\rho }_1}{\text{d}p}=\frac{\text{d\hspace{0.17em}}\ln \rho }{\text{d}p}$ as is appropriate for an adiabatic process for an ideal gas. Then using Eq. (9), we obtain $$\frac{\text{d}{g}^{\prime }}{\text{d}t}=-\frac{3w{g}^{\prime }}{b}\left(\alpha -\frac{b}{3H}\right)\left(1+\frac{g^{\prime }}{g}\right),$$(11)

where the factor of (α − b/3H) arises from Eq. (9) and the factor of (1 + g′/g) arises from relaxing the Boussinesq approximation. Equation (11) reduces to Eq. (5) when b ≪ H and g′ ≪ g.

Equations (9)-(11) can be rewritten in their nondimensional forms, $$\frac{\text{d}\beta }{\text{d}\tau }=\omega \left(\alpha -\frac{\beta }{3}\frac{b_0}{H}\right)$$(12) $$\frac{\text{d}w}{\text{d}\tau }={\tilde{g}}^{\prime }+\frac{\alpha {\omega }^2}{\beta }\left(4{\tilde{g}}^{\prime }\frac{{g^{\prime }}_0}{g}-3\right)$$(13) $$\frac{\text{d}{\tilde{g}}^{\prime }}{\text{d}\tau }=-\frac{3\omega {\tilde{g}}^{\prime }}{\beta }\left(\alpha -\frac{\beta }{3}\frac{b_0}{H}\right)\left(1+{\tilde{g}}^{\prime }\frac{{g^{\prime }}_0}{g}\right),$$(14)

by substituting β = b/b₀, ${\tilde{g}}^{\prime }=g^{\prime }/{g^{\prime }}_0$, $\omega =w/{\left(b_0{g^{\prime }}_0\right)}^{1/2}$, where a 0 subscript indicates the parameter value at τ = 0. Assuming ω = 0 at τ = 0, then, the subsequent dynamics are fully specified by the ratios b₀/H and g′₀/g. Of particular interest to note, when g′₀ ≪ g, then Eqs. (12)-(14) become approximately independent of g′₀, because g′₀ only appears as a ratio with g. Thus the nondimensional behavior of these equations of motion are approximately independent of the initial buoyancy perturbation if g′₀ ≪ g. As we will see in Sect. 4, this result is reproduced in our numerical simulations.

We show how the compressive corrections modify the behavior compared to the Boussinesq case in Fig. 3. The figure shows the evolution of non-dimensionalized parameters b, w, and g′ and compares the results from Eqs. (3)–(5) against Eqs. (9)– (11). We see that this model predicts the thermal to remain more spatially localized, thereby enhancing the retained buoyancy density and a faster velocity. This finding is in agreement with Anders et al. (2019), which also modeled downwelling thermals in highly compressible atmospheres. In that work they characterized “stalling” versus “falling” regimes, distinguished by whether entrainment tended to dilute the buoyancy perturbation and slow subsidence (stalling), orwhethercompres-sion dominates so that the thermal maintains its concentration and subsides at larger velocities (falling). According to the model we present here, the distinction between these regimes depends on the relative magnitude of the entrainment coefficient α and the ratio of the thermal size to density scale height b/H.

It is reasonable to estimate at which point we expect this simplified analytic model to break down. To assess this, we consider the relative importance of shear stress and buoyancy. This ratio has been canonically expressed using the nondimensional Richardson number (Ri): $$\text{Ri}\equiv \frac{g}{\rho }\frac{\frac{\partial \rho }{\partial z}}{{\left(\frac{\partial u}{\partial z}\right)}^2}.$$(15)

When Ri is large, buoyant forces are much greater than shear stresses and thus tend to organize the flow in such a way that Reynolds stresses are a small correction. When Ri becomes sufficiently small, shear stresses begin to dominate and the density current becomes susceptible to large-scale instability. In other applications, the flow will become Kelvin-Helmholtz unstable at Ri< 1/4. It is therefore a reasonable first-guess to compute the Richardson number predicted by our analytic model in order to estimate the depth at which the flow may become unstable and thus break our assumption of self-similarity.

Using the stream function from Eq. (1), we can estimate the relationship between the nondimensional propagation velocity ω of the vortex ring and the corresponding shear it creates in its environment. The stream function is related to velocity according to $$v=\frac{1}{r^2\sin \theta }\frac{\partial \psi }{\partial \theta }\widehat{r}-\frac{1}{r\sin \theta }\frac{\partial \psi }{\partial r}\widehat{\theta }.$$(16)

We are interested in the shear flow between the vortex ring and its surroundings. In our simple model for the Hill’s spherical vortex, there is no flow across the spherical boundary of diameter b. Therefore, the flow at the boundary must be entirely in the $\widehat{\theta }$ direction. Using Eqs. (1) and (16), we can write this component of the velocity outside of the vortex ring as $${\upsilon }_{\theta }=-\left(1+\frac{b^3}{16r^3}\right)w\sin \theta .$$(17)

Then at the vortex-environment boundary, the shear flow is $$\frac{\partial {\upsilon }_{\theta }}{\partial r}=\frac{3w}{b}\sin \theta ,$$(18)

which is maximum at the mid-plane θ = π/2 to give us a characteristic shear of about $\frac{3w}{b}$. Now we estimate the corresponding buoyancy gradient. If we take $\frac{\partial \rho }{\partial z}\sim \frac{\text{Δ}\rho }{b/2}$, then $\frac{g\frac{\partial \rho }{\partial z}}{\rho }\sim \frac{2g^{\prime }}{b}$. This gives us an expression for the scaling coefficient of the Richardson number in terms of our analytic parameters, $$\text{Ri}\sim \frac{2b{g}^{\prime }}{9w^2}.$$(19)

As argued above about Eqs. (9)-(11), in the limit where g′ ≪ g the non-dimensionalized equations of motion are approximately independent of our choice of the magnitude of the initial buoyancy perturbation g′₀, implying that the evolution of the nondimensional Richardson number is likewise independent of this choice. This suggests that the penetration depth of down-drafts is actually not sensitive to the magnitude of the buoyancy perturbation. As we will see in Sect. 4.4, this conceptual result is actually reproduced by our numerical simulations.

In Fig. 4 we show how the behavior of the Richardson number evolves with time, and the corresponding depth of the vortex ring. We see that Ri < 1/4 around τ ~ 2.1 corresponding to z/b₀ ~ 1.25. Beyond that depth, the flow may be susceptible to instability. Determining the nature of any such instability, as well as the timescale on which it self-amplifies to substantially disrupt the self-similar flow structure outlined here, is a complex fluid dynamical problem that requires a more sophisticated model than we have presented in this section. As we demonstrate in Sect. 4, the buoyant perturbation actually continues sinking to a much greater depth before finally breaking apart. We therefore compare the results from this simplified framework with fully hydrodynamic simulations in Sect. 4. We also use the analytical framework from this section as a basis to inspect additional confounding factors in Sect. 6.

Fig. 1 Diagram showing the important quantities in our problem overlaid on the stream function of a Hill’s spherical vortex from Eq. (1). A spherical vortex ring of diameter b falls at speed w through a medium of density ρ. Within the vortex, there is a mean density perturbation of δρ, with $g^{\prime }\equiv g\overline{\delta \rho }/\left(\rho +\overline{\delta \rho }\right)$. The vortex entrains external fluid at a rate $\dot{m}=\frac{\pi }{2}{b}^2\alpha \rho w$.

Fig. 2 Thermal diameter, b, as a function of depth for different compressibility regimes. The more compressible, the more spatially coherent this theory predicts the thermal to remain.

Fig. 3 How the three non-dimensionalized variables of interest evolve for the uniform Boussinesq case (dashed, Eqs. (3)–(5)) and the compressible case (solid, Eqs. (9)–(11)). b is the thermal vortex ring diameter, ω = w/(b0g′0)1/2 is its non-dimensionalized mean propagation velocity, $g^{\prime }\equiv \frac{\overline{\delta \rho }}{{\rho }_1}g$ is its mean buoyancy perturbation, and τ = t/(b/g′0)1/2 is the nondimensional time. The subscript ’ denotes initial values. For this figure we use the following parameters: b0/H0 = 1, ω0 = 0, g′0/g = 10%, and α = 1/4.

Fig. 4 Richardson number as a function of τ of the propagating vortex ring using Eq. (19) and the computed parameters from Fig. 3 (red curves). Dashed and solid curves follow the same convention from Fig. 3. The black curves show the corresponding depth, ζ ≡ z/b0. The dashed horizontal line corresponds to Ri= 1/4, at which point the flow becomes susceptible to Kelvin-Helmholtz instabilities.

3 Numerical model

While simple models are useful for building intuition for complex phenomena, their efficacy must in general be compared against experimental results. The framework we followed in Sect. 2 was tested in a laboratory using tanks of water and colored brine (Morton et al. 1956). For the highly compressible environments spanning multiple density scale heights of interest to us, laboratory experiments are impractical. The Earth’s troposphere possesses only a single density scale height when we are interested in semi-infinite atmospheres with relevant dynamics playing out over three or more density scale heights (orders of magnitude different densities). Even for the densest gases, building a column of multiple density scale heights would be prohibitively impractical. Using sulfur-hexafluoride with µ = 146 g mol⁻¹, for example, building a column of three density scale heights would be H = γk_BT/µg ~ 7 km in height using Earth’s gravity and surface temperature, taller than any manmade structure and rivaling the height of Mount Everest. Furthermore, it would require an immense heat source to sustain an adiabatic gradient over that distance because its adiabatic lapse rate would be far steeper than the surrounding atmosphere’s. Given the obvious impracticality of constructing such an apparatus, we instead employ numerical models to perform computational experiments to test our model and account for the highly chaotic nature of turbulence. We seek to do this in a straightforward, easily reproducible way. However, cognizant of the possible biases that result from our assumptions and the limits of fluid dynamical simulations, we discuss confounding factors and additional tests in Sect. 6.

3.1 Settings

We used the flexible SNAP software (Li & Chen 2019) to simulate a localized negative density anomaly in a deep atmosphere by numerically integrating Euler’s equations using a Riemann solver, with total energy as the conserved quantity. Because of the uncertainty in initial conditions for a storm, we assumed all condensible material has been vaporized and is therefore fully miscible in its environment. Because it has been explicitly verified against benchmark results, and to maximize simplicity and reproducibility of our results, we opted to use the initial conditions of the three-dimensional Straka problem (Straka et al. 1993). First, we initiated an isentropic, ideal gas atmosphere defined by a given temperature, T₀, at a given pressure level that obeys Eq. (6) with an ideal gas equation of state. Embedded in the atmosphere, we introduced a buoyancy anomaly centered at z = x = y = 0 that obeys $$\delta T=\{\begin{array}{ll}-T_0\frac{g^{\prime }}{g}\frac{{\pi }^3}{6\left({\pi }^2-8\right)}\cos \left(\frac{\pi r}{b}\right)\hfill & r\le b/2\hfill \\ 0\hfill & r>b/2\hfill \end{array},$$(20)

where r is the distance from the center of the bubble, and b is the size of the bubble using notation consistent with Sect. 2. The temperature perturbed region is isobaric at each vertical level with its environment, such that δΤ/T = δρ/ρ. The choices of coefficients normalizes Eq. (20) so that definition of ${\displaystyle \int \frac{\delta T}{T_0}\text{d}V=\frac{g^{\prime }}{g}}$ holds. Then our non-dimensionalized results are independent of our particular choices for T₀ and g if g′ ≪ g. We calculate the specific entropy perturbations as $$\delta s=c_p\ln \left(\frac{T}{T_{\text{i}}}\right)+R\ln \left(\frac{p}{p_{\text{i}}}\right),$$(21)

where T_i and p_i are the unperturbed initial temperature and pressure profile of the isentropic background atmosphere, and c_p and R are the isobaric specific heat capacity and the specific gas constant, respectively. In Figs. 5 and 6 colors represent entropy perturbation surface densities. These figures show a two-dimensional representation of the three-dimensional simulation by integrating the total entropy perturbation along the third dimension: ∫_Y ρδsdy, where Y is the simulation domain in the y dimension. We discuss the validity of this treatment in more detail in Sect. 6.

We seek to determine the extent to which explicit simulations agree with the idealized model from Sect. 2. The compressibility of the problem can be described by the dimensionless ratio b₀ /H₀. Additionally, the timescale of the problem is set by initial buoyancy perturbation g′ , the potential energy that generates the motion in w₀ as well as circulation and turbulence. For each simulation, we hold the atmospheric structure constant, varying only the initial conditions of the parameters. We therefore test a suite of choices for b₀ and track the subsequent evolution of the down-drafts. Between each simulation, we keep the relative resolution fixed so that if b₀ changes by a factor of two, so too does our grid spacing. In an incompressible environment, then, the simulation output should be scale-independent for flow of constant Reynolds number (which we discuss in more detail in Sects. 3.2 and 6), a general property of the Navier-Stokes equation. We would expect in the incompressible case that the solution should be identical for constant Reynolds number if we use nondimen-sional time τ = t(g′/b₀)^1/2 and nondimensional space ζ = z/b₀. When comparing our results in this way, then, we can isolate the effect of compressibility b₀/H₀ on our simulation outputs. We report our results compare them to the expectations from Sect. 2 in Sect. 4.

Fig. 5 Comparison of a snapshot in time at τ = 8.8 for a sinking bubble simulation with an initial size of b0/H0 = 2/5 for different choices of grid resolution, where ζ ≡ z/b0 is the dimensionless depth. Colors correspond to the integrated line-of-sight specific entropy perturbations, i.e., ∫Y ρδsdy, where δs is defined in Eq. (21), y is the dimension coming out of the page, and Y is its simulated domain. Units on the color bar are arbitrary, linearly proportional to entropy surface density. The dynamical instability observed for Δx ≥ b0/40 is not observed for the lowest-resolution (leftmost) simulation, whose vortex ring remains stable until colliding with the bottom boundary of the simulation.

3.2 Viscosity and turbulence

We now comment on the validity of assuming constant Reynolds number between flows of different scales. As we see from Eq. (10), the evolution of the velocity w0 depends on the thermal size b in a scale-dependent way. The true Reynolds number of our problem is enormous. As we will show in Sect. 4, the velocity scales can be on the order of hundreds of meters per second for length scales on the order of tens of kilometers, while the actual viscosity of gases under the relevant conditions should be tiny, v ~ 10⁻⁵ m² s⁻¹. This characteristic Reynolds number of the problem would be of order Re = w₀b/v ~ 10¹¹. Fully modeling the microscopic dissipation of such vigorous turbulence in non-equilibrium flow is computationally impossible. Prior works (e.g., Anders et al. 2019; Zhou et al. 2020) have addressed this problem by imposing an artificially low Reynolds number (or equivalently, an artificially high viscosity) to force the flow to remain laminar, or for the turbulence to remain sufficiently weak that it can be fully resolved. In this work, we use the standard Roe-Linearized Riemann solver to reconstruct the hydrodynamic fluxes from the volume-averaged quantities (Roe 1981; Toro 2013). To resolve hydrodynamic fluxes between adjacent grid cells, this method transforms the nonlinear Euler equations into quasi-linear equations using the Jacobian matrix of the flux vector fields, then solves the equations as if they were truly linear. This method introduces a resolution-dependent computational viscosity by effectively dissipating sub-grid-scale motions. We therefore find resolution dependence in our outputs that does not converge at the limits of our computational capacities using the Licallo supercomputing cluster at l’Observatoire de la Côte d’Azur. If we use sufficiently coarse resolution, the flow remains laminar. However, at finer turbulence-resolving resolutions we observe new dynamical instabilities not predicted by the analytic model (further discussion in Sect. 4).

3.3 k− ϵ turbulence closure model

We must ensure that our results are reproducible with different choices for turbulent closure. While imposing a global hyper-viscosity stabilizes flow that should be unstable, we wanted to investigate whether a spatially dependent eddy viscosity generated by the flow itself could produce results that are both physically realistic and fully resolvable. We seek a closure model of turbulence that can return resolution-independent results. Because of the nature of our problem, it would be inappropriate to use an imposed global eddy viscosity; we expect the turbulence to be time- and space-dependent because it is generated by the density currents themselves. We therefore seek a closure model that expresses these aspects of our problem. We chose to use a k − ϵ closure model for turbulence, following Launder & Spalding (1972) and Pope (2000). The k − ϵ model introduces two new dynamical variables: k, representing the unresolved kinetic energy of the turbulence, and ϵ representing the microscopic dissipation of turbulent energy. The turbulent kinetic energy is generated by shear in the flow, while the dissipation of turbulent kinetic energy is generated by the turbulence itself. While the details are not fully based on physics, the motivation is consistent with, for example, a Kolmogorov conception of isotropic fully developed turbulence, with turbulent flow generated at larger scales and being dissipated at smaller scales. The dynamical equations are $$\frac{Dk}{Dt}=\frac{1}{\rho }\frac{\partial }{\partial x_j}\left[\frac{{\mu }_t}{{\sigma }_k}\frac{\partial k}{\partial x_j}\right]+\frac{{\mu }_t}{\rho }\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\frac{\partial u_i}{\partial x_j}-ϵ$$(22) $$\frac{Dϵ}{Dt}=\frac{1}{\rho }\left[\frac{{\mu }_t}{{\sigma }_{ϵ}}\frac{\partial ϵ}{\partial x_j}\right]+\frac{C_1{\mu }_t}{\rho }\frac{ϵ}{k}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)\frac{\partial u_i}{\partial x_j}-C_2\frac{{ϵ}^2}{k},$$(23)

where µ_t is the turbulent viscosity, σ_k and σ_ϵare the effective turbulent Prandtl numbers for turbulent energy and dissipation, respectively, C₁ and C₂ are empirical coefficients, and subscripts on u and x are Einstein summation notation. From these quantities, k, which represents the mean kinetic energy of the unresolved turbulent flow, and ϵ, which represents its microscale dissipation, we can dimensionally reconstruct an eddy viscosity, µ_t = k² /ϵ, which originates non-arbitrarily from the density currents of our problem. Importantly, we can see from Eq. (22) that the turbulent kinetic energy is generated by inter-grid point shear – important given our interpretation from Sect. 2 that shear is the source of instability. This same quantity is dimensionally consistent with diffusivity, and can therefore be reasonably invoked to represent diffusion of material as well as momentum. We test the outcomes of two methods: a diffusive case, where the same diffusion coefficient µ_t diffuses both momentum and density differences, and the non-diffusive case, where µ_t acts as a viscosity only. Computationally we introduce buoyancy perturbations as temperature perturbations. Therefore, the diffusion of buoyancy takes the form of a thermal conductivity. However, because of the considerable temperature gradient in the background atmosphere, we cannot conduct heat directly and must instead conduct a temperature analog using dry static energy Ε = c_pT + gz.

Fig. 6 Sample of snapshots in time for entropy perturbations of a downwelling thermal with b0/H0 = 2/5 and Δx = b0/80. The visualization technique is identical to that of Fig. 5. Darker colors correspond to a larger negative entropy perturbation, integrated over the line of sight to visualize the three-dimensional structure.

4 Results

We observe the thermal to form a coherent vortex ring, initially sinking and remaining spatially concentrated. The subsequent behavior depends on the model assumptions.

4.1 Low-resolution validation of the analytic model

At low resolutions, as argued in Sect. 3.2, the effective Reynolds number of the problem is reduced. In this case, the simulation results are in fair agreement with the prediction from Sect. 2. In Fig. 7a we evaluate the vertical distribution of entropy at each time step according to $$\frac{\text{d}S}{\text{d}z}={\displaystyle {\int }_X{\displaystyle {\int }_Y\rho \delta s\text{d}x\text{d}y,}}$$(24)

where X and Y are the horizontal domains of the simulation in the x and y dimension, respectively, and δs is calculated according to Eq. (21). A vortex ring forms as expected and subsides coherently. Contrary to the analytic case, which assumes zero detrainment, our simulation finds that some detrainment does occur. Nevertheless, the motion of the spatially localized central vortex ring remains accurately described by the analytic model from Sect. 2, and the vortex ring remains concentrated down to the simulation’s bottom boundary, traversing about three density scale heights. This result agrees with the findings of Anders et al. (2019), which studied a similar problem with a different method. In the low-resolution case, our model agrees with their assumption of laminar flow, and our results are consistent.

Fig. 7 Evolution of the initial entropy perturbation with time for different resolutions, ζ = z/b0 is the depth, and τ is the dimensionless time. Colors indicate the magnitude of the entropy perturbation at each height $\frac{\text{d}S}{\text{d}\zeta }$ according to Eq. (24). The dashed curve shows the ζ(τ) obtained from integrating Eqs. (9)–(11) using b0/H0 = 2/5, α = 1/4, and g′/g = 7.5%, consistent with simulation parameters. The dashed curves show the Hill’s spherical vortex boundaries, ζ ± b/2b0.

4.2 Dynamical instability observed at higher resolutions

In Sect. 2 we discussed the susceptibility of the flow to a shear instability. At low resolutions, no instability occurs and the flow remains well-described by the simple model from Sect. 2. However, at higher resolutions we observe a dynamical instability that arrests the subsidence of the thermal. We interpret this to be due to the influence of viscosity. At sufficiently low-resolution, small-scale motions are not well-resolved and therefore the influence of turbulence is weak (see Sect. 3.2). However, at higher resolutions the flow is less artificially stabilized and we observe the dynamical instability suspected in Sect. 2.

After the onset of the instability, the flow begins to rapidly detrain most of the buoyancy perturbation contained within the thermal, effectively mixing with the environment around that depth (see Fig. 7b). We refer to this depth as the “penetration depthö, z_p. In Fig. 5 we see the higher-resolution runs exhibit a dynamical instability in which much of the bubble material detrains from the primary vortex ring. In the low-resolution case, this instability is not observed. One therefore must take care to use sufficient resolution (or equivalently low enough viscosity) to capture dynamically important behaviors. Among models with sufficient resolution to capture the instability, the precise details of the instability are somewhat resolution-dependent. However, for sufficiently high resolutions we observe qualitatively similar behavior, albeit with somewhat quantitatively different details. For our purposes in this study we are concerned with constraining the orders of magnitude of rainy downdrafts, and we find convergence in both qualitative behavior and order-of-magnitude penetration depth.

Figure 6 shows the evolution of a sample simulation that initially forms a vortex ring, propagates coherently for some distance, before undergoing dynamical de-coherence due to its own self-induced turbulence. The vortex ring instability is associated with rapid detrainment from the main thermal. After this, the subsidence of the centroid slows considerably.

Fig. 8 Sample output for the k − ϵ model (cf. Fig. 7) using the same parameters and visualization technique.

4.3 Toward resolution independence: The k − ϵ model

We find that the two methods for accounting for turbulence using a k − ϵ model can capture aspects of the dynamical instability for the thermal. The instability is characterized by both a braking in the subsidence rate dense material, and a rapid diffusion of that material. When we use the k − ϵ model as a pure eddy viscosity model, we reproduce the braking effect, but not the rapid diffusion of material. When we include diffusion of dry static energy as described in Sect. 3.3, the diffusion of material is better modeled, but the subsidence rate is faster than using the Roe-linearization method. These results can be seen in Fig. 8. Though imperfect, the k − ϵ method is able to reproduce these behaviors at resolutions as low as Δx = bo/10. Therefore, if one is interested in a particular aspect of the problem, such as finding the depth at which subsidence will appreciably slow, then the k − ϵ model can be useful. However, its output is not fully robust, and for more detailed calculations, high-resolution, turbulence-resolving simulations should be used instead. We plot the comparison between the motion of centroid compared to the pure Riemann solver method in Figs. 9 and 10; we plot the penetration depth z_p (see following subsection) as a function of atmospheric compressibility b₀/H₀ in Fig. 11.

Fig. 9 Tracking the position of the cen-troid for a variety of model types, varying resolution, and computational treatment of turbulence. “High resolution” refers to Δx = b0/80, while “Med resolution” refers to Δx = b0/40. The k − ϵ models were conducted using Δx = b0/40.

Fig. 10 Tracking the nondimensional vertical velocity, ω = w/(b0g′0)1/2, of the centroid as a function of nondimensional time across a variety of parameters using different computational methods (see Fig. 9).

4.4 The penetration depth, z_p

We define a penetration depth that roughly describes the depth to which most of the material comprising the initial buoyancy perturbation reaches before mixing with its surroundings. While there is no perfect way to describe the output of a complex hydro-dynamic simulation with a single number, here we suggest a couple of quantitative methods to do so in order to quickly compare and communicate the output of simulations with different parameters. One method is to quantify for how long the majority of material remains spatially concentrated. Using Eq. (24), we can bin the entropy perturbations into bins of size Δz at any specified time. We can use this to plot a histogram and cumulative distribution function of entropy perturbation as shown in Fig. 12. We note that total entropy is not conserved, but rather always increases according to the second law of thermodynamics. Therefore, ΔS will always be increasing as material mixes, and as kinetic energy is dissipated as heat. We note that while higher-resolution runs have already largely mixed into a quasi-static distribution, the low-resolution case remains strongly peaked in a vortex ring that continues to sink rapidly (see previous subsections). Another way to visualize the simulation is to demarcate the distribution of material using contours of standard deviation, by tracking the centroid. Figure 13 shows how the dense material distributes itself as the downwelling thermal propagates. The dashed curve shows the location of the centroid, defined such that $${\displaystyle {\int }_{z_{\text{btm}}}^{z_{\text{centroid}}}\frac{\text{d}S}{\text{d}z}\text{d}z={\displaystyle {\int }_{z_{\text{centroid}}}^{z_{\text{top}}}\frac{\text{d}S}{\text{d}z}\text{d}z},}$$(25)

where z_btm and z_top are the bottom and top boundaries of the simulation. The dashed curve shows the analytic prediction for the center downwelling thermal, while the solid curve shows the numerically computed simulation centroid. We see the analytic prediction initially agrees closely with the analytic prediction, but abruptly branches off at some characteristic penetration depth characterized by the dynamical instability visualized in Fig. 6 around τ = 7. Despite this, the leading edge of the downdraft appears to remain roughly consistent with analytical predictions even after the disruption. The analytic prediction is shown as a dotted curve in Fig. 13. Evidently, a coherent down-welling thermal continues to propagate despite having detrained a significant fraction of its initial mass. At late times, interaction with the bottom boundary of the simulation box becomes important. Figure 13 also shows the distribution of entropy. If we take the distribution of entropy $\frac{\text{d}S}{\text{d}z}\text{Δ}z$ to be a histogram with bin size Δz according to the simulation resolution, then the shadows of different opacities show the 1, 2, and 3σ distributions of material. Most of the initial entropy perturbation, then, is contained in the darkest shadow of Fig. 13. We see, then, that for this particular simulation, most of the initial buoyancy perturbation penetrates to large ζ; for this example most material ends up at ζ > 6.

So one quantitative option could be to track the dispersion of material as the downdraft subsides, and define z_p according to the point at which the dispersion of material has increased by some factor. While the thermal begins spatially concentrated, it spreads out as it mixes with the environment. For example, we can quantify the vertical distribution of material around the centroid as seen in Fig. 13. We can define the penetration depth to be the point when the standard deviation of material distribution around the centroid increases by a factor of e: $$\sigma \left(z_p\right)=e{\sigma }_0.$$(26)

An alternative option could be to define z_p as the depth at which subsidence appreciably slows, so that the subsidence of the thermal no longer dominates compared to random turbulent motions in the surrounding environment. We can track the cen-troid of entropy perturbations for the simulation, shown in Fig. 9. After the thermal vortex ring undergoes dynamical disruption, the vertical velocity of the centroid branches off of the analytic expectation, plateauing at a much lower value and mixing with the environment (see also Fig. 13). We can then compute dimen-sionless velocity ω = w/(b₀g′₀)^1/2 as the temporal derivative of the cubic spline interpolation of the centroid position shown in Fig. 9. The result is shown in Fig. 10. We can define the penetration depth to be the point at which the subsidence rate reduces to a value comparable to the characteristic velocity environmental turbulence, $$w_0\left(z_p\right)\sim u_0,$$(27)

at which point we say that it can be disrupted and efficiently mixed by the environment. We find that both methods of evaluating z_p yield comparable results when using u₀/(g′b₀)^1/2 ~ 3/10.

Fig. 11 Penetration depth of rainy downdrafts using different computational methods (see Fig. 9), according to our definitions. Red curves use the velocity flattening criterion Eq. (27) with u0/(g′0b0)1/2 = 0.29 (see also Fig. 10), while black curves use the vertical variance criterion from Eq. (26), see also Fig. 13.

Fig. 12 Sample showing a histogram (left) and cumulative distribution (right, integrating from the top) of entropy at a snapshot in time for low-resolution Δx = b0/20 (red), medium-resolution Δx = b0/40 (green), and high-resolution Δx = b0/80. $\text{Δ}S={\displaystyle \int \frac{\text{d}S}{\text{d}z}\text{d}\zeta }$, where $\frac{\text{d}S}{\text{d}z}$ obeys Eq. (24), either integrating between bin boundaries (left panel) or cumulatively from the simulation top boundary (right panel). Simulation input parameters at the same as in Figs. 7–8.

Fig. 13 Tracking the centroid along with the 1σ, 2σ, and 3σ shadows, which represent the distribution of material for b0/H0 = 2/5 and Δx = b0/40. The dashed curve represents the subsidence of the centroid, while the dotted line shows the analytical model prediction for the depth of the leading edge of the vortex ring from Sect. 2, zlead ≡ zcentroid + b(z)/2. The solid line represents the bottom boundary of the simulation. We use dimensionless depth, ζ = z/b0, and dimensionless time, τ = t/(b0/g′)1/2. The dashed-dotted line shows the penetration depth, zp, computed using Eq. (26) for reference.

4.4.1 Dependence on initial buoyancy perturbation, g′₀’

The behavior, both qualitatively and quantitatively, does not appear to be sensitive to the magnitude of the perturbation g′ for a perturbation of constant size. We tested g′₀/g of 1%, 2.5%, and 7.5%. We found the motion of the centroid, when non-dimensionalized, were practically identical so as to be indistinguishable when over-plotted. Likewise, we found the difference in penetration depths between the 1% and 7.5% case to be affected by <1%. This may appear counter-intuitive; one might expect that a denser mass anomaly should be able to sink faster and therefore penetrate deeper than a less dense counterpart. However, if we assume the self-generated shear stresses to be the source of instability, then this result is unsurprising. A smaller buoyancy anomaly creates less rapid fluid motions, and therefore weaker shear. Indeed, in Sect. 2 we argued that our analytic model predicted that the results should be independent of g’ if g′₀ ≪ g. We note however that our model does not consider environmental turbulence; environmental turbulence may make a difference because it does not necessarily scale with the parameters of the problem, b₀ and g′₀ (see Sect. 6.2).

4.4.2 Dependence on compressibility, b₀/H₀

We wish to inspect the behavior while varying the input parameter b₀/H, and determine the extent to which our results depend on different choices of computational technique. To simplify the visualization for different choices of b₀/H₀ and different computational techniques, we show in Fig. 9 the centroid tracking for these different choices, equivalent to the solid curve in Fig. 13. Then in Fig. 10 we show the corresponding vertical velocity, essentially the time derivative of Fig. 9.

The depth at which this occurs depends on the compressibility of the problem. These qualitative results are robust across a range of computational simulations and methods, shown in summary in Fig. 10. We seek to define the penetration depth, the depth at which the velocity slows considerably and a large fraction of the thermal’s buoyancy perturbation detrains from the central vortex ring. This process is continuous and highly complex, so there is no perfect way to quantify where this happens. We chose two methods to define the penetration depth, as follows.

We find that our results depend somewhat on the choice for defining the penetration depth, as well as the choice of computational techniques. However, the top-line order-of-magnitude results are consistent and generally show z_p/b₀ to decrease with increasing b₀/H₀. However, the magnitude of the decrease is less than the corresponding increase in b₀; in particular, z_p/b₀ decreases only by a factor of ~2 while b₀/H₀ increases by an order of magnitude (Fig. 11). The results from the k − ϵ model appear instructive, if imperfect. The advantage of the k − ϵ model is resolution independence. However, the disadvantage is that only part of the dynamics are captured. As discussed in Sect. 4.3, the k − ϵ model that ignores thermal diffusion, the dynamical braking around the correct z_p is observed even at coarse resolutions that do not capture the instability using a standard Roe-linearized Riemann solver. However, the model does not correctly reproduce the profusion of the density perturbation. For this reason, the k − ϵ model yields the correct z_p when using Eq. (27) but not Eq. (26), (Fig. 11). On the other hand when thermal diffusion is introduced, the dynamical braking is suppressed but the profusion of material is captured; in that case, the model yields the correct z_p when using Eq. (26) but not (27). Therefore, if one wishes to use the k − ϵ as a computational expedient for a problem like this, one must ensure it properly captures the relevant physics of interest. For our problem, the main takeaway is that the results are largely in agreement, regardless of these choices when accounting for the above caveats.

5 Applications

We have so far introduced an analytic and numerical framework to model rainy downdrafts in abyssal atmospheres using dimensionless parameters. We now seek to apply these findings to concrete celestial bodies.

5.1 Jupiter and Saturn

Among our primary motivations for this study were Juno’s observations of ammonia depletion below the cloud level on Jupiter (Li et al. 2017). We seek to estimate how our findings for dimen-sionless parameters from Sect. 4 can be dimensionalized in the context of Jupiter. First we seek to estimate δρ/ρ. In Sect. 4 we considered these values between δρ/ρ ∈ [0.5%, 7.5%]. Recall that our choice within this range does not affect the simulated penetration depth z_p. We want to estimate how the presence and vaporization of condensates affects the density of a gas parcel. For an ideal gas, the density is $$\rho =\frac{p\mu }{RT}=\frac{p\left[x{\mu }_c+\left(1-x\right)\mu d\right]}{R{T}_0-x{\mu }_{c}L},$$(28)

where T₀ is the initial temperature of the background atmosphere before adding and vaporizing condensates, µ_d and µ_c are the mean molecular weights of the dry gas and the condensate vapor, respectively, L is the specific latent heat of vaporization, and x is the condensate mole fraction. Notice that there are two ways that increasing the amount of condensate x increases the density: first by changing the mean molecular weight (numerator), and second by reducing the temperature through the latent heat of vaporization (denominator). Then we calculated $$\frac{\delta \rho }{\rho }=\frac{1+x\left(\frac{{\mu }_c}{{\mu }_d}-1\right)}{1-\frac{x{\mu }_{c}L}{R{T}_0}}-1.$$(29)

Next we considered two cases: ammonia condensing in the stratosphere, and water condensing a few bars beneath the photosphere around 300 K. Considering the ammonia level corresponds to ammonia snowing as an isolated species, while the water level would imply that water and ammonia mix, perhaps as a simple fluid or perhaps following the microphysics of mushballs (Guillot et al. 2020b,a). Starting with ammonia, we find between 3-10 × x_⊙ (where x_⊙ is the solar abundance) the δρ/ρ ∈ (1.3%, 4.5%), comfortably within our considered range. For water again considering between 3-10 × x© we find δρ/ρ e (2.5%, 8.5%) again comparable to our considered range.

We then sought to constrain b₀/H₀ for this problem. The density scale height is $H_0=\frac{\gamma R{T}_0}{\mu g}$. For Jupiter this is about 36 km at the ammonia cloud level, and 63 km at the water cloud deck. We can estimate b₀ based on previous simulations that find a characteristic size of ~25 km (Hueso et al. 2002). Fronts can be larger than individual cells, with fronts observed to be of order 50–100 km wide and at times up to 400 km long. Using 20–50 km as a characteristic size range, then based on Fig. 11 we expect the penetration depth of vaporized ammonia to be roughly between 80–150 km. Assuming the ammonia first vaporizes around the 1 bar level, this corresponds to a pressure level of ~ 10–30 bars. These findings suggest that ammonia snow on its own, if it rains from the 1 bar level as a concentrated thermal in a fashion comparable to what we describe in this paper, could plausibly penetrate to a depth of 10 s of bars, as observed by Juno (Li et al. 2017) without invoking any interactions with water. Using water we have H₀ that is smaller by about a factor of two, with an estimate of b₀ that is relatively unconstrained, because JunoCam cannot see down to the water cloud level. Mushball microphysics would further increase the expected penetration depth, as the vaporization of ammonia will be delayed to greater pressures. If this mechanism is indeed at play, it is not unreasonable to suspect that water may behave similarly. If we assume b₀ to be fixed, the penetration depth should be comparable, of order 100–150 km below the water cloud deck, corresponding to a pressure level of 30–100 bar. Based on these arguments, it is plausible that water may likewise be depleted below its cloud deck.

We can use the same argument to predict that Saturn may exhibit a similar depletion in water and ammonia below their respective cloud decks. However, the situation on Saturn may be more complicated; Saturn’s quasi-periodic Great White Spot erupts roughly every 30 Earth-years rather than quasi-continuously as on Jupiter. If a significant fraction of Saturn’s moist convective activity occurs intermittently on multi-decadal timescales, that leaves long quiescent times during which water vapor can diffuse upward again. Convective inhibition may also be at play on Saturn (Li & Ingersoll 2015), further complicating the story.

5.2 Uranus and Neptune

Because Uranus is a high priority target for a future Solar System mission, we consider the implications of our findings there, and on its Solar System cousin Neptune. We consider methane, for which we have some constraints on abundance, estimated to be of order ~40 × x_⊙ on both Uranus (Lindal et al. 1987) and Neptune (Conrath et al. 1991). Using these numbers and taking the methane cloud layer to be at 100 K, we obtain an estimate for δρ/ρ at least 40%, much larger than on Jupiter. Therefore, we can expect that downdrafts associated with such extreme methane enrichment may be denser and stronger than the down-drafts on Jupiter, and additional considerations in violation of our assumption that g′ ≪ g for Eqs. (10)–(11). For the following calculations will will extrapolate our finding that the penetration depth is insensitive to δρ/ρ. Although methane storm activity has been observed on Uranus (de Pater et al. 2015) and Neptune (Molter et al. 2019), there have not been close observations of such events with the resolution and detail of JunoCam. If we assume b₀/H₀ ~ 1, comparable to Jupiter, and calculating the density scale height on Uranus at the methane cloud level to be around 50 km, we calculate the penetration depth to be 150– 200 km for methane storms, corresponding to a pressure level of 20–40 bar. The numbers would be the same on Neptune, owing to Uranus and Neptune’s nearly identical effective temperatures and similar surface gravity. The possibility of a convectively inhibited atmosphere on Uranus and Neptune with a statically stable radiative region (Markham & Stevenson 2021), however, may complicate the dynamics compared to the neutrally stable case investigated in this work. Methane has been unambiguously detected in the atmosphere of Uranus and Neptune, but if methane is depleted below the cloud level like ammonia is on Jupiter, then the abundances that have been measured may not be representative of its bulk abundance. Additionally, methane may not be well-mixed and could include compositional gradients (cf. Sromovsky et al. 2011). Indeed, there is already evidence for latitudinal variation in the distribution of methane on Uranus (Molter et al. 2021) and Neptune (Tollefson et al. 2021). Whether methane is substantially depleted to tens of bars beneath its cloud level may be tested as part of an upcoming mission. If its structure resembles Jupiter, that observation would favor a fluid dynamical explanation for the observed depletion.

5.3 The Sun

Convection on the Sun is somewhat simplified by the absence of condensates undergoing phase transitions. We follow Anders et al. (2019) to estimate the relevant buoyancy and length scales of thermals from the solar photosphere. For a background atmosphere of temperature T₀ ~ 6000 K, the downflow temperature deviation is expected to be δΤ~ − 500 K so that δρ/ρ ~ 8.3%, roughly consistent with our simulations. The thermal perturbation size should be close to the width of the thermal downflow lanes ~ 100 km. For solar gravity at 274 m s⁻² for monatomic hydrogen and helium, this implies b₀/H₀ ~ 1. Therefore, we expect the penetration depth of “entropy rain” to be 300–400 km, or < 0.2% of the total depth of the convective zone. Therefore, in contrast to the findings of Anders et al. (2019), we find that atmospheric compression alone is not a sufficient mechanism to preserve the coherence of a downwelling thermal. The difference between our results stems from our different treatment of turbulence. Anders et al. (2019) considered the laminar case; in that case, our results agree. However, when resolving turbulence using the Roe-Linearized Riemann solver method at high resolution, or using a turbulence closure model like k − ϵ, the density current is arrested at a finite penetration depth z_p that prevents the downwelling thermal from maintaining its coherence through the convective envelope of the Sun. In order for convective heat transport by entropy rain to persist (as described by Brandenburg 2016), some other focusing mechanism, (perhaps magnetic fields; see e.g., Cattaneo et al. 2003; Weiss et al. 2004), or different assumptions about the nature of the problem would be required.

5.4 Exoplanets

Beyond our own Solar System, the basic fluid dynamics we describe here should also apply to exoplanets. The specific length scales of relevance would depend on the details of the system, and we do not attempt to make exhaustive quantitative predictions here. However, our non-dimensionalized results from Sect. 4 can be straightforwardly applied to an exoplanetary system using the methods from this section. Furthermore, we note that our mechanism for depleting volatile vapor below the cloud level should serve as a caution not to mistake measured atmospheric abundances of volatile species for bulk interior abundances in planets with gas envelopes.

6 Confounding factors

Our main analysis has employed a highly simplified picture of rainy downdrafts, aimed at identifying some of the main qualitative features of this type of flow and constraining the order-of-magnitude of downdraft penetration depth. The simplified assumptions employed were used as a starting point to minimize the number of free parameters that could affect our results, as well as to maximize the simplicity of reproducing our results. However, we must inspect whether our idealized results can be reasonably applied to real planetary environments. In this section we address a variety of concerns one might pose.

6.1 Initial conditions

In our computations, we chose somewhat contrived initial conditions based on the three-dimensional Straka problem. One might reasonably ask whether our results would differ if some alternative initial conditions were specified. To inspect this problem, we compare our results for the three-dimensional Straka case to a different arbitrary initialization. As one example, we begin with identical initial conditions to the three-dimensional Straka problem, but use a uniform ΔT in lieu of Eq. (20). As another example, we initialize a cube of material with density perturbations distributed randomly in a cube instead of a centrally condensed bubble as in the Straka case. The cube model in particular is notably different in that its initialization does not involve spherical symmetry. For each grid space, we assign the temperature perturbation to lie randomly between 0 K and 30 K over a cube sized such that the expectation value for the total entropy perturbation is comparable to the 15 K Straka bubble. We find that despite the initialization of the random cube model bearing no spherical symmetry, a vortex ring will still spontaneously form (see Fig. 14). Moreover, we find that the top-line qualitative and order-of-magnitude results hold – the new initial conditions likewise allow the initial perturbation to sink down many multiples of its initial size before substantially mixing with its environment. The penetration depth increases by less than 30%, which is of small consequence in the context of our order-of-magnitude interests. We show how penetration as measured by Eq. (26) changes with different initializations and methods in Table 1. We compare the different numerical methods for the Straka problem medium resolution, high resolution, and the k−ϵ model also plotted on Fig. 11 to results from simulations with different initial conditions. For different initial conditions we investigate the uniform sphere and random cube case described above. While we do not find perfect agreement, we do find consistent orders-of-magnitude, indicating that our overall findings are not catastrophically sensitive to initial conditions. Of course, the choice of initial conditions will affect the subsequent dynamics, but our top-line results are not too sensitive to different arbitrary choices. Modeling a realistic storm, including the development of upwelling convective columns, the microphysics of condensation and coagulation, the precipitation and vaporization, in highly compressible environments is the subject of ongoing and future work. For the moment, our results apply to spatially localized downwelling thermals.

Fig. 14 Integrated line-of-sight entropy perturbations after some initial evolution of the random cube model with a top-down view, using the same plotting technique as, e.g., Fig. 6. We observe that, despite a lack of spherical symmetry upon initialization, a vortex ring as approximately described in Sect. 2 spontaneously forms, albeit retaining some of the stochasticity from the initial conditions.

6.2 Environmental turbulence

We have considered a neutrally stratified adiabatic background atmosphere, but real planetary atmospheric dynamics is driven by cooling from the top. The ensuing turbulent convection that we expect, as well as driven thermal winds, shear between belts and zones, and other dynamic processes in the atmosphere should excite significant turbulence in the environment external to the downwelling thermal. We neglected this effect in our prior analysis. We do not attempt to model such an uncertain and complex phenomenon in our simulations, but we do point to theoretical and experimental work that has been done in the past.

6.2.1 Disruption of a thermal by turbulent surroundings

Following Turner (1962), we can model environmental turbulence in the following way. In the quiescent environment, the flow is driven by the motion of the thermal itself, so that entrain-ment is proportional to its propagation speed w. In a turbulent environment, there exists flow driven by the motion of the thermal, but also environmental flow driven by some outside dynamics. If the turbulence is homogenous and has characteristic length scales that are small compared to the thermal, we can model its effect on the thermal to be detrainment of material from the thermal into the environment. We can specify a characteristic turbulent velocity u₀. In this case, Eqs. (3)–(5) become $$\frac{\text{d}{b}^3}{\text{d}t}=3\alpha b^2{w}_0-3b^2{u}_0$$(30) $$\frac{\text{d}\left(b^3{w}_0\right)}{\text{d}t}=\frac{2}{3}{b}^3{g}^{\prime }-b^2{u}_0{w}_0$$(31) $$\frac{\text{d}\left(b^3{g}^{\prime }\right)}{\text{d}t}=-b^2{u}_0{g}^{\prime }.$$(32)

The above equations modify the evolution equation for the thermal diameter b to be b − b₀ = αz − u₀t, where z is the depth. Introducing the scaling parameter of buoyancy flux F_* = b³g′₀, we can non-dimensionalize the above equations so that $$b_1=\alpha z_1-u_0{t}_1,$$(33)

where $b_1\equiv b{u}_0/F_{\ast }^{1/2}$, $z_1\equiv z{u}_0/F_{\ast }^{1/2}$, and $t_1\equiv t{u}_0^2/F_{\ast }^{1/2}$. We can non-dimensionalize the other quantities using these same scaling parameters: w₁ = w/u₀, and f = b³g′/F_*. The evolution of this system is shown in Fig. 15.

Initializing b₀ = 0 for a point source, the thermal initially grows, before reaching a maximum size and beginning to shrink. At this point, the environmental detrainment begins to overwhelm the thermal entrainment rate (i.e., w ~ u₀). The thermal will continue to propagate as it shrinks, before fully disappearing at some height. We can solve for this numerically, finding z₁ ~ 0.52. Notice that we use b₁(t₁) = 0 here in contrast to b(τ = 0) = b₀ used in the prior sections. This difference arises because the non-dimensionalization in this problem uses u₀ and F_*, in contrast to our use of g′₀ and b₀ in prior sections. This difference arises due to the nature of the problem, because u₀ is a new fundamental scaling relationship for the problem. The statement that b₁ → 0 is essentially saying that (b₀u₀)² ≪ F_*, in the beginning; in other words, the gravitational potential energy of the thermal is much greater than the turbulent kinetic energy of the surroundings.

We now must estimate this quantity for a planet like Jupiter. For a thermal of size b = 40 km, then we expect the penetration depth when accounting for environmental turbulence to be $$z_p\sim 4\times {10}^3{\left(\frac{u_0}{1\text{m\hspace{0.17em}}{\text{s}}^{-1}}\right)}^{-1}\left(\frac{\overline{\text{Δ}T}}{15\text{\hspace{0.17em}K}}\right)\text{\hspace{0.17em}km.}$$(34)

Using u₀ ~ 1 m s⁻¹ is consistent with what we expect for equilibrium turbulent convection. However, one can use a larger value to account for meteorological activity to obtain a shallower penetration depth. Furthermore, if $\overline{\text{Δ}T}$ is smaller, this likewise would decrease the penetration depth. When we compare this value and scaling to our findings using simulations scaled to Jupiter, we find the depth dictated by expected environment turbulence to be one or two orders of magnitude larger that the depth dictated by self-induced turbulence. Evidently, for weak environmental turbulence the effect of instability wrought by the flow itself is more important than the effect of environmental turbulence. Still, for highly turbulent environments or small initial buoyancy perturbations, the effect of environmental turbulence may become important. We therefore add this possibility as a caveat to our findings from Sect. 4.

We have not modified the behavior to account for the compressibility of the environment. We will, however, comment on the behavior we expect qualitatively. Following Fig. 3, we see that in a compressible environment we expect faster propagation velocities. Since the dynamics of the thermal in a turbulent environment are governed by the relative scale of w and u₀, in general for downwelling thermals in compressive environments we expect w to dominate to a deeper depth in the analytic case. We do not endeavor to model compressive convection in a turbulent environment fully in this work, as we have just demonstrated that environmental turbulence should be less important than self-induced turbulence for an environment like Jupiter’s. Atmospheric compressibility does not change this conclusion, and we therefore continue to emphasize our main result that finds a penetration depth caused by a self-induced turbulent instability rather than primarily by mean environmental turbulence.

Fig. 15 Evolution of dimensionless parameters for a thermal propagat ing in a turbulent environment (cf. Turner 1962). $b_1\equiv b{u}_0/F_{\ast }^{1/2}$, $z_1\equiv z{u}_0/F_{\ast }^{1/2}$, $t_1\equiv t{u}_0^2/F_{\ast }^{1/2}$, w1 =w/u0, and f = b3g′/F*, with F*, = b3g′0 The simulation truncates when b1 → 0, indicating that its entire initia buoyancy perturbation has detrained into the environment.

6.2.2 Redistribution of material by eddy diffusion

In this work, we primarily focused on one rather straightforward problem: the dynamics of an initially localized thermal propagating through a neutrally stratified atmosphere at rest. However, in reality, such events would only be one part of the much more extensive dynamics at play in the atmosphere. A fuller picture of the atmosphere would include our results as one component of a larger atmospheric model. As one example, using the eddy-diffusion mass-flux (EDMF) models that have been applied to the Earth (e.g., Siebesma et al. 2007; Suselj et al. 2019), we have investigated plumes in one direction, but have not really discussed eddy diffusion. While a fully robust atmospheric model of this kind is beyond the scope of this work, we can make some comments on the applicability of this work to future endeavors to incorporate our results into more general models of atmospheric circulation. EDMF models employ a modified diffusion equation of the form $$\frac{\partial \overline{\varphi }}{\partial t}=\frac{-\partial }{\partial z}\left[-K\frac{\partial \overline{\varphi }}{\partial z}+M\left({\varphi }_{\text{NL}}-\overline{\varphi }\right)\right],$$(35)

where ϕ refers to the quantity of interest (for example the concentration of a volatile species), $\overline{\varphi }$ refers to the mean abundance at a given vertical coordinate, and ϕ_NL refers to nonlocal mixing from strong updrafts, or downdrafts in our case. Solving this equation in general requires some model for the mass flux M, as well as a relationship between ϕ_NL and the propagation velocity. Our results from Sect. 2 as well as Fig. 10 provide some constraints that can be used for these models. In addition, a current complication when endeavoring to employ EDMF models on giant planets is the lack of an obvious bottom boundary. While the top boundary of the stratosphere holds for giant planets as well, on models of the Earth, the Earth’s surface acts as an obvious bottom boundary. On Jupiter, no such bottom boundary exists. Therefore, when modeling plumes in abyssal environments, one must have some estimate for the distance over which it can be treated according to the analysis in Sect. 2, and at which point the flow may become unstable and simply mix with its environment. Fig. 11 offers reasonable estimates to use when constructing such a model.

We now seek to estimate how efficiently rainy downdrafts can remove volatile species, compared to the efficiency of large-scale turbulent mixing. If we think of the atmosphere under an eddy-diffusion paradigm, then retaining a steady-state compositional gradient requires sources and sinks of composition to accommodate the eddy diffusion flux. According to Fick’s law, the steady-state flux is $J\sim -K\frac{\partial \varphi }{\partial z}$, using the same definition of ϕ from Eq. (35). If we regard the composition sink to be rain-out at the top boundary, and the composition source to be the penetration of rainy downdrafts at the bottom boundary, then we can estimate the efficiency of this exchange using a simple diffusion model.

Consider Jupiter. Juno observed a compositional gradient to vary between about 100 ppm at 1 bar to 373 ppm at 50–60 bar (Li et al. 2017), corresponding to a distance of approximately 60 km. The mean compositional gradient would thus be ~4.5 ppm km⁻¹. If we use a characteristic eddy diffusivity for Jupiter of 2 × 10⁴ m² s⁻¹ by using Jupiter’s scale height of 20 km and vertical convective velocities of order 1 ms⁻¹ (Guillot et al. 2004), then we estimate this eddy diffusion flux to be of order 3 × 10⁻⁷ moles of NH₃ per square meter per second around the 1 bar level. Integrated over the surface of Jupiter, this corresponds to around 10⁸ kg of ammonia per second. One plume as explored in this model, for example the nominal model of ΔT ~ 15 K with b = 40 km corresponds to about 2.4 × 10⁹ kg of ammonia around the 1 bar level. Therefore, to an order of magnitude based on these rough numbers, we would expect to need such a plume to descend about once every four minutes somewhere on Jupiter in order to accommodate the observed compositional gradient. This estimate is consistent with Juno’s observed lightning statistics, detecting large lightning strikes on approximately the same intermittency timescale (Brown et al. 2018).

This work can be regarded as a modest first step toward building a more sophisticated atmospheric model, for example and EDMF model, for Jupiter’s atmosphere.

Fig. 16 Comparing the vertical location of the centroid for different models that address different confounding factors: changing the gravitational field, changing the initial conditions, and introducing vertical wind shear. It is difficult to distinguish by eye between the two assumptions for gravity, because they overlap each other almost perfectly in non-dimensionalized form.

6.3 Vertical wind shear

The background velocity field in our simulations is initialized to be zero. Besides environmental turbulence, this likewise neglects the possibility of vertical wind shear. On the Earth, vertical wind shear is correlated with increased storm activity. If we expect the same principle to apply to Jupiter and other planets, then the regions where stormy rainy downdrafts are occurring probably possess some amount of vertical shear, and we must assess the extent to which this modifies our conclusions. In Fig. 16 we compare zero shear to moderate vertical wind shear (wind velocity changing by 10 m s⁻¹ over one density scale height). The results are practically unaffected. We find that for extreme vertical wind shear (100 m s⁻¹ over one density scale height, not plotted), the perturbation is completely ripped apart and efficiently mixed with the environment. In that case, the extreme vertical wind shear produces Kelvin-Helmholtz instabilities, acoustic radiation, and vigorous dissipative turbulence such that such an extreme shear would require a large amount of energy to sustain. For more moderate vertical wind shear, our results are perturbed but not profoundly affected. On Jupiter, beneath about 4 bars (above the vaporization level for water), the vertical wind shear is very slight (Seiff et al. 1997). Therefore, we can say with some confidence that our results are not too unrealistic for atmospheres that possess realistic vertical wind shear.

6.4 Stratification

Our investigations in this work focus on a neutrally stratified environment. However, some observations suggest Jupiter’s atmosphere may actually be weakly stably stratified (Allison & Atkinson 2001; Wong et al. 2011; Guillot et al. 2020a; Bolton et al. 2021). In this case, the thermal’s acceleration will be modified by the change in potential temperature of its surroundings. If we specify a Brünt-Väisälä frequency, N, associated with the stratification, then Eq. (5) is modified to become (Turner 1973) $$\frac{\text{d}\left(b^3{g}^{\prime }\right)}{\text{d}z}=-b^3{N}^2.$$(36)

For positive real N², the solution to this equation behaves like a converging exponential decay. Therefore, for a constant stratification gradient, there exists some depth at which the thermal will become neutrally buoyant as a combination of its dilution from entrainment and its propagation to a higher density environment. The depth at which this will occur can be determined using dimensional analysis: $$z_{\max }=5{g^{\prime }}_0{}^{1/4}{b}_0^{3/4}{N}^{-1/2}.$$(37)

The coefficient of proportionality can be determined empirically so that this result agrees with observations over a large range of scales (Briggs 1969). If we use characteristic values for Jupiter of b₀ ~ 25 km, g′₀ ~ 2 m s⁻², and N ~ 0.02 s⁻¹ (Wong et al. 2011), we obtain z_max ~ 84 km. This depth should be larger in the case of compression, because the effect of dilution from entrainment will be counteracted by the effect of compression as the thermal falls. Nevertheless, if Jupiter’s atmosphere is really stably stratified, this could substantially affect our results. In general, we expect the penetration depth in a stably stratified environment to be less than for a neutrally stratified environment.

In abyssal atmospheres, some part of the atmosphere may be stably stratified even if it is convective or neutrally stratified above and below that layer. It is then natural to wonder whether a vortex ring could pierce such a layer of stability before continuing to propagate downward. Assuming the layer of static stability is encountered at some depth z < z_p so that the thermal vortex ring structure remains intact, we can ask two types of questions. The first is whether the compositional gradient is extreme enough to decelerate the downwelling thermal. If the magnitude of δρ/ρ within the thermal exceeds the change in buoyancy associated with the compositional gradient. One way to do this is to compare the virtual potential temperature $${\theta }_{\upsilon }\equiv \frac{T}{1-x\left(1-{\mu }_c/{\mu }_d\right)}{\left(\frac{p}{p_{\text{ref}}}\right)}^{-\left(\gamma -1\right)/\gamma },$$(38)

where p_ref is some reference pressure. If θ_υ is the thermal is less than θ_υ at the base of the stable layer, then to first approximation the thermal should pierce through to the other side. On Jupiter where we expect relatively weak stratification, this should occur, although in planets like Uranus with stronger stratification there may be a level of neutral buoyancy within the stable layer. Then one must compare the potential energy barrier from the level of neutral buoyancy to the kinetic energy within the thermal to determine whether it can pierce the stable layer. within the thermal and at the base of the stable layer. On Uranus or Neptune, this would require downdraft velocities in excess of 100 m s⁻¹ (Markham & Stevenson 2021).

6.5 Applicabilityto different environments

We conducted the majority of our tests for Figs. 9, 10, and 11 for a particular set of environmental conditions, namely a hydrogen/helium atmosphere with an Earth-like acceleration of 9.8 m s⁻², releasing the thermal from 4 bars with a temperature of 200 K. However, we claim that our results are robust to different assumptions about the atmosphere. The results are only sensitive to your choice of g and H₀; when properly non-dimensionalized, the results will be equivalent. We demonstrate this in Fig. 16 by changing the gravity of the simulation, finding the non-dimensionalized simulation output to be practically indistinguishable. The results are the same when changing other atmospheric parameters such as the Grüneisen parameter or molecular weight.

6.6 Compositional effects

While the dynamics are robustly unaffected in a composition-ally uniform atmosphere, we have not in this work addressed how compositional differences between the thermal and the environment may affect the dynamics. In this work we modeled density perturbations as temperature perturbations. This is appropriate so long as the Grüneisen parameter is identical inside and outside the thermal, and if the condensing species is fully vaporized. If these assumptions are not met, however, the dynamics are affected. For example, ammonia and water both have more molecular rotational degrees of freedom than hydrogen and helium. They therefore have a smaller Grüneisen parameter, and a less steep adiabatic gradient. If the environmental lapse rate is steeper than the adiabatic lapse rate within the thermal, the dynamics behave like an unstably stratified environment. This may further enhance the downward propagating potential of rainy downdraft laden with vapors like water, ammonia, and methane. The importance of this effect, however, depends on the concentration of vapor in the downdraft, something that is beyond the scope of this work. Additionally, microphysics may allow some fraction of the precipitation to remain condensed, coexisting alongside its dense vapor phase. Follow-up work will investigate the consequences of these dynamics in greater detail, modeling the storms from their initiation through their eventual rainy downdraft phase.

7 Conclusion and discussion

7.1 The penetration depth

Through modeling and simulation, we have arrived at the robust conclusion that spatially localized strong density perturbations can propagate to ~3–8 times their initial size without substantially mixing with their environment (see Figs. 6, 13, and 11). The propagation depth depends on the ratio of the initial perturbation size to the local density scale height, b₀/H₀. We define the penetration depth either according to the vertical spreading of material increasing by a factor of e above its original distribution, or according to velocity flattening where the centroid subsidence velocity falls below some defined value. We find good agreement between the two criteria. However, both criteria are arbitrary; really in neither case does subsidence stop entirely. The material continues to subside beyond the point at which detailed simulations become computationally impractical. Eventually, unmodeled effects, such as environmental turbulence or stratification, are needed to halt the subsidence – but there is nothing fundamental preventing a buoyancy anomaly from subsiding indefinitely. That being said, the downdraft remains quantifiably coherent and concentrated down to a significant depth without significantly mixing with its environment. Ther-mals on the order of the size of a density scale height can propagate about three density scale heights down without mixing. For much smaller perturbations ten times smaller than a density scale height, they can still propagate coherently down to eight times their initial size, or almost a full density scale height.

7.2 Depletion of volatiles in giant planet atmospheres

Applied to Jupiter, this lends credence to the notion suggested in Guillot et al. (2020a) that, even after the vaporization of volatile species, downdrafts can continue to sink coherently without efficiently mixing with their environment. Our results also suggest that the compositional gradients observed in the volatile abundance on Jupiter may exist on other planets, even without invoking the ammonia-water phase mixture relevant on Jupiter. A rainy downdraft that vaporizes quickly upon reaching its boiling point as it descends will still represent a spatially localized density perturbation. The subsequent dynamics would play out along the lines of our simulated results, with the perturbation maintaining some coherence over relatively long length scales -on Jupiter, on the order of a hundred kilometers.

Our results indicate that a pond of rain originating from a source on this scale can be expected to propagate down to tens of bars before efficiently mixing with its environment. It will be of interest to assess the particularity or specificity of the compositional gradients in Jupiter’s volatile distribution observed by Juno. Here we present a plausible mechanism by which a planet with a simpler hydrological cycle could exhibit similar behavior. Uranus, for example, exhibits energetic methane storms (de Pater et al. 2015). Rain from such events could plausibly penetrate to a great depth before mixing with its surroundings. Measuring whether Uranus’s deep methane abundance beneath the clouds is homogenous or variable will be of great interest for an upcoming mission. Comparative planetology between Uranus and Jupiter will be important for interpreting spectra from exoplanets, which may be depleted in volatiles in the observable atmospheres if they are stormy. These results are also of interest to the question of the depth of the weather layer on giant planets. Our results suggest that the weather layer associated with a condensing species should extend substantially below its cloud layer. In the context of Jupiter, our results suggest the weather layer should extend at least to tens of bars.

7.3 On nonlocal convection in the Sun

We now compare our results and conclusions to prior work that has investigated similar phenomena. In particular, Anders et al. (2019) investigated the concept of a localized thermal in a highly compressible environment in the context of the Sun, motivated by the notion of entropy rain. In that work, the authors find results in agreement with our predictions from Sect. 2, but did not produce our observed braking instability from Sect. 4. The primary cause for the differences between our findings was our different treatment of turbulence; their work modeled the laminar case, while this work considers the turbulent case. We discuss these differences in more detail in Sects. 3.2 and 5.3. Concerning entropy rain as described by Brandenburg (2016), our results do not directly corroborate the Anders et al. (2019) finding that the phenomenon can be a simple consequence of atmospheric compression. However, our findings do not rule out the possibility of nonlocal convection if coherence can be enforced by some other mechanism, for example by considering rotation or magnetic field effects.

7.4 Outstanding problems and future work

In this paper we focus particularly on tracking the motion of vaporized rainfall. Future steps for a more complete understanding will involve modeling the complete hydrological cycle, including upwelling columns and homogenous turbulent convection. As discussed in Sect. 6, our results are impacted by more complex but realistic planetary environments. Layers of stable stratification, likely to exist at least intermittently as they do on Earth, as well as environmental turbulence and vertical wind shear likely have an effect on the output. Additionally, the initial conditions we chose for this study were largely ad hoc; a fully self-consistent model that includes storm development and precipitation microphysics would be a reasonable next step. Here we have found, however, that simply invoking a mixing length on the order of a scale height for the downdraft coherence length scale may actually be an underestimate. Indeed, downdrafts can maintain coherence without mixing with their environment down to multiple scale heights, in the context of Jupiter on the order of a hundred kilometers or tens of bars. Planets are natural laboratories of exotic physical processes. We must visit more to learn more about the poorly understood physics, and the range of possible planetary climates.

Acknowledgements

This work has been funded by the CNES postdoctoral fellowship program. Computations were performed using the Licallo supercomputing cluster hosted by Observatoire de la Côte d’Azur. We would like to thank Jeremie Bec and Holger Homann for instructive discussions about resolution dependence and closure models for turbulence. We would also like to thank the anonymous reviewer for constructive feedback. Data availability: Hydrodynamic output data in netCDF format and accompanying analysis software available upon request.

References

All Tables

All Figures

Fig. 1 Diagram showing the important quantities in our problem overlaid on the stream function of a Hill’s spherical vortex from Eq. (1). A spherical vortex ring of diameter b falls at speed w through a medium of density ρ. Within the vortex, there is a mean density perturbation of δρ, with $g^{\prime }\equiv g\overline{\delta \rho }/\left(\rho +\overline{\delta \rho }\right)$. The vortex entrains external fluid at a rate $\dot{m}=\frac{\pi }{2}{b}^2\alpha \rho w$.

In the text

	Fig. 2 Thermal diameter, b, as a function of depth for different compressibility regimes. The more compressible, the more spatially coherent this theory predicts the thermal to remain.
In the text

Fig. 3 How the three non-dimensionalized variables of interest evolve for the uniform Boussinesq case (dashed, Eqs. (3)–(5)) and the compressible case (solid, Eqs. (9)–(11)). b is the thermal vortex ring diameter, ω = w/(b0g′0)1/2 is its non-dimensionalized mean propagation velocity, $g^{\prime }\equiv \frac{\overline{\delta \rho }}{{\rho }_1}g$ is its mean buoyancy perturbation, and τ = t/(b/g′0)1/2 is the nondimensional time. The subscript ’ denotes initial values. For this figure we use the following parameters: b0/H0 = 1, ω0 = 0, g′0/g = 10%, and α = 1/4.

In the text

	Fig. 4 Richardson number as a function of τ of the propagating vortex ring using Eq. (19) and the computed parameters from Fig. 3 (red curves). Dashed and solid curves follow the same convention from Fig. 3. The black curves show the corresponding depth, ζ ≡ z/b0. The dashed horizontal line corresponds to Ri= 1/4, at which point the flow becomes susceptible to Kelvin-Helmholtz instabilities.
In the text

Fig. 5 Comparison of a snapshot in time at τ = 8.8 for a sinking bubble simulation with an initial size of b0/H0 = 2/5 for different choices of grid resolution, where ζ ≡ z/b0 is the dimensionless depth. Colors correspond to the integrated line-of-sight specific entropy perturbations, i.e., ∫Y ρδsdy, where δs is defined in Eq. (21), y is the dimension coming out of the page, and Y is its simulated domain. Units on the color bar are arbitrary, linearly proportional to entropy surface density. The dynamical instability observed for Δx ≥ b0/40 is not observed for the lowest-resolution (leftmost) simulation, whose vortex ring remains stable until colliding with the bottom boundary of the simulation.

In the text

	Fig. 6 Sample of snapshots in time for entropy perturbations of a downwelling thermal with b0/H0 = 2/5 and Δx = b0/80. The visualization technique is identical to that of Fig. 5. Darker colors correspond to a larger negative entropy perturbation, integrated over the line of sight to visualize the three-dimensional structure.
In the text

Fig. 7 Evolution of the initial entropy perturbation with time for different resolutions, ζ = z/b0 is the depth, and τ is the dimensionless time. Colors indicate the magnitude of the entropy perturbation at each height $\frac{\text{d}S}{\text{d}\zeta }$ according to Eq. (24). The dashed curve shows the ζ(τ) obtained from integrating Eqs. (9)–(11) using b0/H0 = 2/5, α = 1/4, and g′/g = 7.5%, consistent with simulation parameters. The dashed curves show the Hill’s spherical vortex boundaries, ζ ± b/2b0.

In the text

	Fig. 8 Sample output for the k − ϵ model (cf. Fig. 7) using the same parameters and visualization technique.
In the text

	Fig. 9 Tracking the position of the cen-troid for a variety of model types, varying resolution, and computational treatment of turbulence. “High resolution” refers to Δx = b0/80, while “Med resolution” refers to Δx = b0/40. The k − ϵ models were conducted using Δx = b0/40.
In the text

	Fig. 10 Tracking the nondimensional vertical velocity, ω = w/(b0g′0)1/2, of the centroid as a function of nondimensional time across a variety of parameters using different computational methods (see Fig. 9).
In the text

	Fig. 11 Penetration depth of rainy downdrafts using different computational methods (see Fig. 9), according to our definitions. Red curves use the velocity flattening criterion Eq. (27) with u0/(g′0b0)1/2 = 0.29 (see also Fig. 10), while black curves use the vertical variance criterion from Eq. (26), see also Fig. 13.
In the text

Fig. 12 Sample showing a histogram (left) and cumulative distribution (right, integrating from the top) of entropy at a snapshot in time for low-resolution Δx = b0/20 (red), medium-resolution Δx = b0/40 (green), and high-resolution Δx = b0/80. $\text{Δ}S={\displaystyle \int \frac{\text{d}S}{\text{d}z}\text{d}\zeta }$, where $\frac{\text{d}S}{\text{d}z}$ obeys Eq. (24), either integrating between bin boundaries (left panel) or cumulatively from the simulation top boundary (right panel). Simulation input parameters at the same as in Figs. 7–8.

In the text

Fig. 13 Tracking the centroid along with the 1σ, 2σ, and 3σ shadows, which represent the distribution of material for b0/H0 = 2/5 and Δx = b0/40. The dashed curve represents the subsidence of the centroid, while the dotted line shows the analytical model prediction for the depth of the leading edge of the vortex ring from Sect. 2, zlead ≡ zcentroid + b(z)/2. The solid line represents the bottom boundary of the simulation. We use dimensionless depth, ζ = z/b0, and dimensionless time, τ = t/(b0/g′)1/2. The dashed-dotted line shows the penetration depth, zp, computed using Eq. (26) for reference.

In the text

Fig. 14 Integrated line-of-sight entropy perturbations after some initial evolution of the random cube model with a top-down view, using the same plotting technique as, e.g., Fig. 6. We observe that, despite a lack of spherical symmetry upon initialization, a vortex ring as approximately described in Sect. 2 spontaneously forms, albeit retaining some of the stochasticity from the initial conditions.

In the text

Fig. 15 Evolution of dimensionless parameters for a thermal propagat ing in a turbulent environment (cf. Turner 1962). $b_1\equiv b{u}_0/F_{\ast }^{1/2}$, $z_1\equiv z{u}_0/F_{\ast }^{1/2}$, $t_1\equiv t{u}_0^2/F_{\ast }^{1/2}$, w1 =w/u0, and f = b3g′/F*, with F*, = b3g′0 The simulation truncates when b1 → 0, indicating that its entire initia buoyancy perturbation has detrained into the environment.

In the text

	Fig. 16 Comparing the vertical location of the centroid for different models that address different confounding factors: changing the gravitational field, changing the initial conditions, and introducing vertical wind shear. It is difficult to distinguish by eye between the two assumptions for gravity, because they overlap each other almost perfectly in non-dimensionalized form.
In the text