Convective dynamics and disequilibrium chemistry in the atmospheres of giant planets and brown dwarfs

Disequilibrium chemical processes have a large effect upon the spectra of substellar objects. To study these effects, dynamical disequilibrium has been parameterized using the quench and eddy diffusion approximations, but little work has been done to explore how these approximations perform under realistic planetary conditions in different dynamical regimes. As a first step in addressing this problem, we study the localized, small scale convective dynamics of planetary atmospheres by direct numerical simulation of fully compressible hydrodynamics with reactive tracers using the Dedalus code. Using polytropically-stratified, plane parallel atmospheres in 2- and 3-D, we explore the quenching behavior of different abstract chemical species as a function of the dynamical conditions of the atmosphere as parameterized by the Rayleigh number. We find that in both 2- and 3-D, chemical species quench deeper than would be predicted based on simple mixing length arguments. Instead, it is necessary to employ length scales based on the chemical equilibrium profile of the reacting species in order to predict quench points and perform chemical kinetics modeling in 1-D. Based on the results of our simulations, we provide a new length scale, derived from the chemical scale height, which can be used to perform these calculations. This length scale is simple to calculate from known chemical data and makes reasonable predictions for our dynamical simulations.


INTRODUCTION
As the number of known substellar objects has grown into the thousands, a unique testing ground has developed for our understanding of planetary atmospheres. Substellar objects with masses similar to or in excess of that of Jupiter, subject to little or no irradiation, comprise an especially interesting sample for study, yielding strong observational signals while simultaneously affording relatively simple atmospheric chemistry and dynamics for modeling (Moses 2014;Showman et al. 2010). Jupiter, directly imaged giant exoplanets, and brown dwarfs all share a similar general atmospheric structure, with a lower convective zone driven by cooling of the interior and a cool overlying radiative zone (Burrows et al. 1997;Showman et al. 2010). This structure can facilitate the generation of especially strong molecular features (Madhusudhan et al. 2014), which has led to the (unambiguous) detection of CH 4 , H 2 O, and CO in the atmospheres of directly imaged giant exoplanets (Patience et al. 2010;Barman et al. 2011Barman et al. , 2015Oppenheimer et al. 2013;Konopacky et al. 2013;Janson et al. 2013;Snellen et al. 2014;Chilcote et al. 2015;Macintosh et al. 2015;Samland et al. 2017) and CH 4 , NH 3 , H 2 O, CO and CO 2 in the atmospheres of brown dwarfs (e.g., Geballe et al. 1996Geballe et al. , 2009Oppenheimer et al. 1995Oppenheimer et al. , 1998Schultz et al. 1998;Noll et al. 1997Noll et al. , 2000Leggett et al. 2000;Saumon et al. 2000Saumon et al. , 2007Yamamura et al. 2010). In particular, the abundances of these molecules are sensitive to the strength of vertical mixing in these atmospheres, making it essential to accurately incorporate the effect of dynamics into the interpretation of these observations (Saumon et al. 2003;Zahnle & Marley 2014;Moses et al. 2016).
To explain observations of CO in the atmosphere of Jupiter, Prinn & Barshay (1977) (hereafter PB77) used mixing length theory to describe atmospheric transport, and found that in the limit of small chemical scale height (the scale height of the chemical reaction rate), the predicted abundance will quench to the equilibrium value at the pressure where chemical and dynamical timescales become equal. This approximation provides a framework for tying the strength of atmospheric mixing to observed abundances in the form of the eddy diffusion coefficient, where w is a vertical velocity and L(z) is the characteristic length scale of the atmosphere, which in PB77 is taken to be a density scale height. Further work by Smith (1998) demonstrated that the characteristic length scale for quenching is a function of properties of the equilibrium profile and chemical and dynamical timescales. Modern 1-D atmospheric models utilizing advanced chemical kinetics (e.g., Moses et al. (2011)) build upon this work by directly employing K zz in the evolution equation for each chemical species, where n i is the concentration, or number density, of species i, D is the molecular diffusion coefficient, and P i and L i are production and loss rates . To estimate profiles of K zz for these models, a combination of analytical approximations and simulated velocity profiles from 3-D general circulation models (GCM) are typically used. While these 1-D models have demonstrated the presence of quenching behavior for all of the species discussed above in the atmospheres of Jupiter, directly imaged giant exoplanets, and brown dwarfs (e.g., Moses et al. 2005;Zahnle & Marley 2014;Moses et al. 2016), they are fundamentally limited by their dependence upon accurate prescriptions for the dynamics. Various efforts have been made to explore and describe the transport properties of these atmospheres in 2-and 3-D simulations that employ more realistic dynamics using passive and reactive tracers. One approach has been to include the evolution of a reactive (Cooper & Showman 2006) or settling passive (Parmentier et al. 2013) tracer to a 3-D GCM solving the primitive equations in the radiative zone. While Cooper & Showman (2006) found that quenching occurred as predicted based on Smith (1998), Parmentier et al. (2013) saw that measured eddy diffusion coefficients differed from those estimated conventionally using Equation (1) by orders of magnitude. Another approach, taken by Freytag et al. (2010) has been to consider the evolution of a reactive tracer in a 2-D local box model solving the fully compressible equations of hydrodynamics in a coupled radiative-convective atmosphere. Freytag et al. (2010) observed similar discrepancies between measured and estimated eddy diffusion coefficients. Given these discrepancies, it worth a closer look into the applicability of the eddy diffusion approximation in chemical mixing, and the current models for chemistry in dynamic disequilibrium.
In the present work, we take a different approach from previous dynamics studies. Using 2-D and 3-D local convective box models, we hold the majority of atmospheric parameters constant, and explore the transport properties of a reactive tracer as a function of increasingly realistic dynamical forcing regimes. These simple convection experiments are the regime for which the mixing behavior should be well-described by mixing length theory and eddy diffusion. We compare our results to the predictions of PB77 and Smith (1998), as well as a 1-D model employing a vertical eddy diffusion coefficient profile determined using the results of our simulation. We find an agreement with Smith (1998) that our measured quench points are significantly deeper than would be predicted by PB77, in both 2-and 3-D. We find similarly good agreement with a length scale, H chem,eq , in inferring the quench point from data and predicting it from 1-D models. This length scale is simple to calculate from known chemical data. We conclude with a discussion of our results and their implications for the modeling of chemical disequilibrium.

Equations and assumptions
We solve the fully compressible equations of hydrodynamics, where ρ, u, σ, e, and T are the fluid mass density, velocity, stress tensor, specific internal energy, and temperature respectively, g is the gravitational acceleration, and χ is the thermal diffusivity. We make the assumption of a Newtonian stress tensor, with constant dynamic viscosity µ = νρ, and a constant Prandtl number, The viscosity is defined in terms of the Rayleigh and Prandtl numbers as, where L z is the depth of the convective zone, ∆s is the entropy jump across the domain, and c p is the specific heat capacity at constant pressure. Rayleigh numbers specified in this work are defined at the top of the atmosphere, but increase by a factor of e 2n H (where n H is the number of scale heights across the domain) from top to bottom. The Rayleigh number is left as a free parameter that is explored in this study, as realistic Rayleigh numbers for planetary atmospheres (≈ 10 20 ) are not within reach of current simulations. We hold gravitational acceleration to be constant, but in our formulation of these equations, we do not make the hydrostatic approximation, which can compromise the accuracy of interpretation of phenomena which occur on smaller horizontal scales (Mendonça et al. 2016).
We assume an ideal gas equation of state with a ratio of specific heats γ = 7/5, to emulate a primarily diatomic molecular gas. Our initial conditions and background are of a polytropic atmosphere. For a convective region, this has a polytropic index of, where is the superadiabatic excess, which we take to be 10 −4 based on work by Anders & Brown (2017) demonstrating that ∝ Ma 2 and based on MLT calculations finding that typical Ma for these objects should be ∼ 10 −2 . Polytropic atmospheres have background temperature profiles (T 0 ) that are linear with altitude, which leads to the profile shown in the left panel of Figure 1. For an ideal gas, as we work with here, with P = ρT , the initial density profile is ρ = T m and the initial pressure profile is P = T m+1 . Our initial conditions are in hydrostatic and thermal equilibrium, but are superadiabaticallystratified and unstable to convection. As shown in Figure 1, the polytropic temperature-pressure profile does a good job of approximating the realistic profile for a directly imaged giant planet used in Moses et al. (2016).

Reactive tracer field
To this system, we add a reactive tracer c that reacts according to the reaction, where X 2 is presumed to be the dominant component of the atmosphere (i.e. of density approximately equal to ρ 0 , the background density), and k is the forward reaction rate constant, constructed as an unmodified Arrhenius rate law, where τ (z = 0) is a fixed value for the timescale of the reaction at the bottom of the atmosphere, with T act the fixed activation temperature for the reaction. The choice of these two values sets the ratio of chemical to density scale heights and the vertical profile of the chemical timescale. The reactive tracer is evolved along with the dynamics, where c is the mole fraction of our reactive tracer, ρ 0 is the background density profile, and c eq is the equilibrium profile of c. D is the molecular diffusivity, set by the Schmidt number, which for simplicity's sake we set to one (i.e., ν = D). Background profiles are used for the concentration of X 2 , and the Arrhenius rate law due to the fact that fluctuations in the density variable will be of order (Anders & Brown 2017). We chose to work with two different equilibrium profiles in this work, as shown in the right panel of Figure  1, aimed at bracketing the two extremes of equilibrium profiles relevant to real chemical species. Specifically, those that vary slowly with height (i.e. approximately linear or constant profiles) and those that vary rapidly with height (i.e. profiles that can be approximated using tanh functions). By including two different choices of equilibrium profiles, we are able to explore speciesdependent effects on quenching behavior.

Numerics
In this work, we utilize the open-source framework Dedalus 1 (Burns et al. 2017), which has been successfully used to study atmospheric dynamics (e.g., Lecoanet et al. 2014Lecoanet et al. , 2015Lecoanet et al. , 2016bAnders & Brown 2017). Dedalus is a pseudo-spectral code that employs implicitexplicit timestepping and solves linear terms implicitly in spectral space, and nonlinear terms explicitly in physical space (Burns et al. 2016(Burns et al. , 2017. Our simulations resolve the physical viscosity of the simulation on the grid, obviating the need for artificial hyperdiffusive terms, or filters that suppress artificial modes in the numerical solution (e.g., Cooper & Showman 2006). In our simulations, all field variables are represented with a Fourier 1 Dedalus is available at http://dedalus-project.org basis in the horizontal, and a Chebyshev basis in the vertical (Boyd 2001). The code evaluates nonlinear terms on a grid with a factor 3/2 more points than Fourier coefficients. Resolutions are reported in Table 1.
In our simulations, the horizontal boundaries are periodic. On the upper and lower boundaries we employ stress-free, impenetrable, fixed thermal flux, and zero tracer flux boundary conditions. We initiate our simulations by setting the fluctuating temperature component field to be a noise field of small perturbations (of order 0.01 or smaller). Initially, the fluctuating component of density is set equal to zero, and the tracer fields begin in equilibrium. We have performed experiments using different sets of initial conditions for the tracer fields and found no dependence of the tracers' evolution on their initial conditions. 3. RESULTS

Dynamic evolution
All simulations are dynamically evolved to a steady state characterized by a relatively constant rms vertical velocity, with evolution times described in Table  1. With increasing Rayleigh number, our 2-D simulations transitioned from laminar flows which settled on a constant rms vertical velocity quickly and showed no long term evolution (10 4 ≤ Ra < 10 6 ) to flows initially dominated by multiple large-scale circulating structures which eventually underwent merger events prior to very slowly evolving towards a constant rms vertical velocity (10 6 < Ra ≤ 10 7 ). In contrast, in 3-D simulations, plume structures at the top and bottom of the domain dominate the transport, and large-scale circulating structures never develop due to the lack of a pre- ferred axis. These three dynamic regimes are showcased in Figure 2, which compares two different 2-D Ra cases captured at the same time in their dynamical evolution, and a 3-D case that has reached a similarly evolved state after a shorter evolutionary period. The demand for higher resolution and substantial long term evolution of the higher Ra 2-D cases only further highlights the significant computational challenges of accessing realistic regimes for planetary atmospheres (Showman et al. 2011). This is due to the establishment of the large-scale circulating structures mentioned above that span the domain in 2-D (which are not present in 3-D). In 3-D, shorter evolution times, due to the lack of these slowly evolving, circulating structures make runs less challenging. The increased resolution requirements due to the extra dimension, however, adds significantly to computational cost. Given that realistic Ra are inaccessible to atmospheric modeling, we explore the rela- tionship between Ra and the mean transport properties of the atmosphere while acknowledging the caveat that further dynamic regimes may exist that do not scale with what we observe and infer here.

Disequilibrium chemistry
When each simulation was evolved to the dynamically steady state described in Section 3.1, reactive tracers were injected, and the simulations were evolved forward a further five hundred buoyancy times in the 2-D cases, and roughly one hundred buoyancy times in the 3-D cases, where, τ B = L z /g . In all 2-D cases, these reactive tracers were observed to reach a steady state, as measured by a lack of evolution in various measures of the quench point, within approximately fifty buoyancy times. In the 3-D cases, this convergence happened even more rapidly, taking approximately thirty buoyancy times.
Two different metrics are used to measure the quenching behavior of the reactive tracers in our simulations. The first metric is an "observed" quench point, where the mole fraction of the reactive tracer at the top of the atmosphere is used as a proxy for what would be observed in an atmosphere remotely. The corresponding pressure for where the equilibrium profile of the reactive species takes on that value is taken to be the quench point, making the standard assumption that the quenching species immediately takes on the equilibrium value where it begins to quench. The second metric is a measured departure point, which is the location at which the reactive species profile moves above the equilibrium profile. This metric allows for comparison of the value at which the tracer is observed to quench to the location at which the profile actually begins to diverge from the equilibrium profile. Figures 3 and 4 show the 2-D and 3-D median vertical profiles of vertical velocity and the mole fraction of the reactive tracer, averaged horizontally and medianed in time over the end portions of their evolution. In 2-D, with increasing Ra, the vertical velocity profile shifts to slightly higher values (left panel), leading to an observed quench point deeper in the atmosphere (right panel). As the vertical velocity profiles increase in magnitude, the quench point inferred from the top of the atmosphere moves deeper into the atmosphere. In contrast, in 3-D, the vertical velocity profiles decrease in magnitude with increasing Ra, and the observed quench point rises in the atmosphere. In the 2-D case, there is  Figure 3 for the 3-D cases. The tracer profiles are medianed over between 25 and 150 buoyancy times after they reach a steady state. These cases were not evolved for as long as the 2-D cases, as they reached a steady state more quickly and with significantly less variation. In the 3-D cases, the median vertical velocity of the dynamics decreases with increasing Rayleigh number, and the quench point inferred from the top of the simulation moves higher in the atmosphere. also a pronounced transition region between the departure point and the height in the atmosphere at which the reactive tracer profile takes on a constant value. This behavior, which is reduced in the 3-D case, is likely an effect of the constrained 2-D dynamics described in Section 3.1. This is demonstrated in Figure 5, which shows the two-dimensional distribution of the reactive tracer and the dynamical flow patterns. In the 2-D case (top), the large-scale circulating flow structures trap fluid in their relatively stagnant centers (bright regions), leading to regions without significant quenching and therefore lower concentrations of the reactive tracer. Above these stagnant centers, horizontal flows homogenize the upper levels of the profile to the value determined by the region of strongest vertical flow. In the 3-D case (bottom), these large-scale circulating structures are weaker and more local, and so strong vertical flows mix the fluid to different quench values that are not smoothed out across the entire domain as is true in the 2-D case. The region of strongest vertical flow, however, still dominates in setting the quench point determined by the mean profile. Figures 6 and 7 shows various measured and predicted quench points relative to the classical quench point predicted using the theory described by PB77. Three salient trends emerge in Figure 6 for both equilibrium profile cases in 2-D: 1) significantly deeper quenching than predicted by PB77, 2) the steady approach of the departure point towards some depth intermediate between the classic PB77 prediction and the departure measured at the lowest Ra probed, and 3) the slow convergence between the departure point and the measured quench point.
Based on these trends, some results of the quench approximation may reasonably be expected to hold at realistic Ra. Namely, that the quenched profile will be rapidly achieved, such that the measured quench point will accurately track the point at which the equilibrium and actual chemical profiles diverge. Comparing the two cases, the slowly varying linear profile clearly has a wider spread of predicted and measured quench points than the rapidly varying tanh profile. This speaks to the fact that for species where quenching occurs in a relatively constant region of their equilibrium profile, measuring and predicting quench points is necessarily more difficult due to the fact that a very small difference in mole fraction may reflect a large change in pressure.
In our 3-D simulations, trends are present, but not as clearly as in our 2-D cases. Both the measured quench point and the departure point continue the trend shown in 2-D of quenching deeper in the atmosphere than predicted by PB77, but whether or not the departure point is converging towards a point intermediate between that measured at the lowest Ra and the prediction of PB77, and towards the measured quench point is less clear. Although the measured quench point and the departure point may converge at much higher Ra, it is not clear from the Ra probed that this is the case, making it difficult to predict from the quench point where the chemical profile actually diverges from the equilibrium profile. Similarly to the 2-D cases, however, we see that the tanh profile has a narrower range of measured quench depths and departure points, indicating that our interpretation from our 2-D results, that it is easier to measure quenching in a rapidly varying region of a chemical profile, holds in 3-D.
As was described in Smith (1998) and confirmed in Cooper & Showman (2006), a significantly smaller characteristic length scale, H Smith , than the density scale height is required to match our results to predictions made using mixing length theory. We find, however, that the Smith length scale is still too large to accurately predict the behavior of our reactive tracer in 2-D.
A slightly better metric is the pure chemical length scale described in PB77, which in our formulation takes the form, where H k and H ρ are the scale heights of the rate constant and density respectively. A general derivation of this length scale is given in the Appendix. This length scale is simpler to derive than the Smith length scale in the limit that the departure point accurately predicts the observed quench point (i.e. if c varies very slowly with height after quenching), such that the third term can be taken to be zero, An even better estimation, as shown in Figure 6 as H chem,eq , involves approximating the third term by setting c equal to its equilibrium profile (as should be valid at the quench point), which is only slightly more difficult to calculate from known chemical data. In 3-D, it is not clear from our results that one length scale does significantly better than another in predicting the measured quench point, as it crosses the predictions of the various length scales with increasing Ra. In this case, H chem,eq seems to still do the best job of predicting the departure point, but it is difficult to infer if this will continue to be the case with increased Ra.

1-D Chemical Modeling
In addition to our 2-and 3-D models, we created a 1-D eddy diffusion model that solves for the equilibrium profile of the reactive tracer according to, Figure 7. The same measurements made in Figure 6 for the 3-D simulations. where w is the median profile of the vertical velocity calculated over the course of the last five hundred buoy- Figure 9. The same measurements made in Figure 8 for the 3-D simulations.
ancy times of our simulations' evolution, and the eddy diffusion is calculated using that same velocity profile and the various length scales probed in Figures 6 and  7 following Equation 1. These models are solved as a nonlinear boundary value problem with ∂/∂t = 0. In 2-D and 3-D, as shown in Figures 8 and 9, 1-D models using eddy diffusion coefficients involving the chemical equilibrium profile and rate law (e.g., H Smith or H chem ) do a relatively good job of modeling the quenching behavior we measured from our simulations. No one eddy diffusion model (based on the different length scales) stands out strongly as the preferred model to use for 1-D chemical modeling. In 2-D, there is some spread in the results of the different models, though the Smith length scale seems to do marginally better at predicting the observed quench point and departure point. In 3-D, there is also some spread to the results, but the best model to use is not clear. The fact that in the case of the 1-D models no one length scale wins out speaks to the difference between modeling small scale dynamics where the disequilibrium chemistry is occurring at, and the dynamics driving transport at large scales.

CONCLUSION
In this work we take a different tack from previous dynamical studies exploring the transport properties of the atmospheres of weakly irradiated giant plan-ets and brown dwarfs, and examine the relationship between chemical disequilibrium behavior and the dynamic regime being simulated. We parameterize the dynamic driving by the Rayleigh number, Ra. Using 2-and 3-D local box models of vertically stratified convective regions in a parameter space relevant to these objects, we have confirmed that the simple assumption of a chemically-independent characteristic length scale is not sufficient to explain the quenching behavior of a reactive tracer. Instead the length scale depends on the chemical properties of the reacting species (as was demonstrated for the 1-D case by Smith (1998)). Further, we have shown that there exists another length scale, the chemical scale height, H chem,eq (Equation 12), which more accurately predicts the measured quench point, can be incorporated effectively into an eddy diffusion coefficient for 1-D chemical kinetics modeling, and is simpler to calculate than the Smith length scale. Finally, we have shown that for the Ra that are currently accessible to numerical simulation, the transport properties of the atmosphere depend on the dynamic regime being simulated. This is an important consideration as the field moves forward towards coupling complex chemistry more directly with more realistic atmospheric dynamics (Heng & Showman 2015;Drummond et al. 2016).
In this paper we have focused on validating the quench approximation in the convective regions of substellar atmospheres. Overlying these convective regions, however, there are the stably stratified radiative zones. In these zones, the basis of the quench approximation, the eddy diffusion approximation, does not necessarily hold, and the effect of the transport properties of these zones upon quenching species may be very different (Parmentier et al. 2013). Additionally, even if the eddy diffusion approximation does reasonably describe the transport properties of the radiative zone, for species that quench within these regions, differences in the structure of the temperature-pressure profile may lead to quenching regimes where current models for estimating the eddy diffusivity no longer apply (such as regions where the chemical timescale increases very slowly). We look forward to further investigation into the dynamics of these regions, and their interaction with the chemistry of the atmospheres of giant planets and brown dwarfs.
Computations were conducted, with support from the NASA High End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center, on Pleiades with allocations GID s1647 (PI: Brown) and GID s1419 (PI: Toomre). Financial support was provided by University of Colorado Boulder and Bates College startup funding. Additionally, this work benefited from the Exoplanet Summer Program in the Other Worlds Laboratory (OWL) at the University of California, Santa Cruz, a program funded by the Heising-Simons Foundation. The authors would like to thank Mark Rast, Julianne Moses, and Vivien Parmentier for useful discussions. We thank the anonymous referee for useful comments that improved the quality of the paper.
the concentration of the reactant A will evolve (in the absence of dynamics), according to, where P i and L i are the production and loss rates for other reactions in the chemical network. If the forward reaction given by Equation 1 is the rate-limiting step for the destruction or production of a species X, then the timescale for a system to relax to equilibrium from some initial concentration of X, [X] 0 , can be estimated as, where the assumption is made that the reverse reaction is very slow compared to the forward reaction and that the forward reaction is part of the only significant mechanism for X production/destruction. For a species to quench, an atmospheric parcel dynamically perturbed from equilibrium must be replenished by the dynamics on a timescale,