Sedimentation Efficiency of Condensation Clouds in Substellar Atmospheres

AM01 calculates the molar mixing ratio of condensed material, q_c, by solving

$$\begin{eqnarray}&&-{K}_{{zz}}\displaystyle \frac{\partial {q}_{t}}{\partial z}-{f}_{\mathrm{sed}}{\omega }_{* }{q}_{c}=0,\end{eqnarray} \tag{ 1 }$$

where q_t is the total condensate molar mixing ratio (including both the condensed and vapor phases), z is altitude, and ω_* is the convective velocity scale. Equation (1) states that the upward mixing of condensate vapor and particles is balanced by the downward sedimentation of particles. Heuristically, under the assumption that q_c/q_t and the mixing length, L, are constant with altitude and given a q_t value below the cloud, q^below_t, the q_t profile above the cloud base is given by

$$\begin{eqnarray}&&{q}_{t}(z)={q}_{t}^{\mathrm{below}}\exp \left(-{f}_{\mathrm{sed}}\displaystyle \frac{{q}_{c}}{{q}_{t}}\displaystyle \frac{z}{L}\right),\end{eqnarray} \tag{ 2 }$$

with all supersaturated condensate vapor condensing and z = 0 at the cloud base. In the full model, L is defined as the ratio of the local lapse rate to the dry adiabatic lapse rate times the scale height H, with a lower limit of 0.1H. The eddy diffusivity K_zz is then related to L by

$$\begin{eqnarray}&&{K}_{{zz}}=\displaystyle \frac{H}{3}{\left(\displaystyle \frac{L}{H}\right)}^{4/3}{\left(\displaystyle \frac{{{RF}}_{h}}{{\mu }_{a}{\rho }_{a}{C}_{p}}\right)}^{1/3},\end{eqnarray} \tag{ 3 }$$

where R is the universal gas constant, μ_a is the atmospheric molecular weight, ρ_a is the atmospheric mass density, and C_p is the specific heat capacity of the atmosphere at constant pressure. F_h is the local heat flux carried by convection, which approaches the interior heat flux σT_eff⁴ in the deep, convective interior of a planet, where σ is the Stefan–Boltzmann constant and T_eff is the effective temperature. In radiative regions, K_zz is set to a specified minimum value. For simplicity, we set K_zz to be a constant throughout the atmosphere in our study (see Section 2.3). The convective velocity scale ω_* is then given by K_zz/L.

The particle size distribution is assumed to be lognormal, with the effective (area-weighted) particle radius given by

$$\begin{eqnarray}&&{r}_{\mathrm{eff}}={r}_{g}\exp \left(\displaystyle \frac{5}{2}{\mathrm{ln}}^{2}{\sigma }_{g}\right),\end{eqnarray} \tag{ 4 }$$

where σ_g is the geometric standard deviation of the lognormal distribution, set to 2 to crudely span the condensation and coagulation modes of cloud particles, and r_g is the geometric mean particle radius, defined through integration of the lognormal distribution by

$$\begin{eqnarray}&&{r}_{g}={r}_{w}{f}_{\mathrm{sed}}^{1/\alpha }\exp \left(-\displaystyle \frac{\alpha +6}{2}{\mathrm{ln}}^{2}{\sigma }_{g}\right).\end{eqnarray} \tag{ 5 }$$

The system is closed analytically by locally fitting the dependence of the particle radius r on the sedimentation velocity of particles v_f using a power law around the radius r_w such that v_f(r_w) = ω_*:

$$\begin{eqnarray}&&{v}_{f}={w}_{* }{\left(\displaystyle \frac{r}{{r}_{w}}\right)}^{\alpha },\end{eqnarray} \tag{ 6 }$$

with v_f given for Stokes flow by

$$\begin{eqnarray}&&{v}_{f}=\displaystyle \frac{2}{9}\displaystyle \frac{{\rm{\Delta }}\rho {{gr}}^{2}\beta }{\eta },\end{eqnarray} \tag{ 7 }$$

where Δρ is the difference in density between the background atmosphere and the density of the particle; g is the gravitational acceleration; β = 1 + 1.26 Kn is the Cunningham slip correction factor, with Kn the Knudsen number of the particle; and η is the atmospheric dynamic viscosity. At Reynolds numbers between 1 and 1000, v_f is augmented by a standard parameterization with respect to the drag coefficient C_d and Reynolds number, while at Reynolds numbers greater than 1000, v_f is given by

$$\begin{eqnarray}&&{v}_{f}=\beta \sqrt{\displaystyle \frac{8{\rm{\Delta }}\rho {gr}}{3{C}_{d}{\rho }_{a}}}.\end{eqnarray} \tag{ 8 }$$

We use the Community Aerosol and Radiation Model for Atmospheres (CARMA) as the "standard" to which we compare AM01. CARMA is a one-dimensional aerosol microphysics model that solves the discretized continuity equation for aerosol particles, subject to vertical transport due to sedimentation and eddy diffusion and production and loss due to particle nucleation (homogeneous and heterogeneous), condensation, evaporation, and coagulation. CARMA was initially used to model Earth's stratospheric sulfate aerosols (Toon et al. 1979; Turco et al. 1979), and has since been generalized to a variety of applications both on Earth and in other planetary atmospheres, including modeling polar stratospheric clouds to inform ozone depletion (Toon et al. 1988), the characteristics of the particles stemming from the eruption of Mount Pinatubo (Zhao et al. 1995), various tropospheric cloud features on Earth (Ackerman et al. 1993, 1995; Jensen et al. 1994), the sulfuric acid clouds of Venus (James et al. 1997; McGouldrick & Toon 2007; Gao et al. 2014), water ice clouds on Mars (Colaprete et al. 1999), and photochemical hazes on Titan (Toon et al. 1992), Pluto (Gao et al. 2017), and ancient Earth (Wolf & Toon 2010). Extension of the model to three-dimensions (Toon et al. 1989) has allowed for the study of Martian dust storms (Murphy et al. 1993) and meteoric dust in Earth's mesosphere (Bardeen et al. 2008), among others. In the rest of this section, we describe qualitatively CARMA's treatment of the various cloud processes and we refer the reader to Appendix for a full description of the physics that is included in CARMA.

CARMA discretizes the vertical extent of the atmosphere into layers and resolves the particle size distribution using mass bins rather than assuming a particular size distribution shape. The sedimentation velocity is computed in the same way as in AM01, but the particle velocities are calculated for each mass bin individually. Similarly, the diffusion velocity associated with eddy diffusion (given a user-defined K_zz) is calculated for particles in each mass bin and for condensate vapors according to Toon et al. (1988). Both the sedimentation and diffusion velocities are then combined to solve for the vertical particle distribution at each time step.

The formation of clouds in CARMA begins with homogeneous or heterogeneous nucleation. Homogeneous nucleation is the formation of stable clusters of condensate molecules directly from the vapor that can then grow to larger cloud particles. The rate of homogeneous nucleation is controlled by (1) the flux of molecules to the cluster, which is dependent on the abundance of condensate vapor, and (2) the material properties of the condensate, such as its surface energy and molecular weight. High surface energy and molecular weight materials tend to nucleate slower than low surface energy and molecular weight materials given the same supersaturation and local temperature. Unlike homogeneous nucleation, heterogeneous nucleation requires foreign surfaces on which stable clusters can form; such foreign surfaces are provided by other aerosol particles in the atmosphere, typically described as condensation nuclei. Thus the size and abundance of these particles strongly impact the rate of heterogeneous nucleation. In addition, the nucleation rate is dependent on the interaction between the condensate and the surface, characterized by the contact angle between the condensate cluster and the surface, the energy needed by a condensate molecule to desorb from said surface, and the oscillation frequency of the condensate molecule on said surface, which is related to the desorption energy (Pruppacher & Klett 1978). The values of the contact angle and the desorption energy are unknown for the substances considered here and must be estimated. Typical contact angles for water range from ∼0° for highly hydrophilic surfaces to >100° for highly hydrophobic surfaces (e.g., Ethington 1990; Yuan & Lee 2013); small contact angles lead to high heterogeneous nucleation rates, and large contact angles lead to low heterogeneous nucleation rates. As our study is focused on trends in cloud distribution instead of simulating any specific known system, we use a small contact angle (01) for all our heterogeneous nucleation cases to increase the efficiency of cloud formation and the optical depth of the resulting cloud. Thermal desorption energies for simple molecules (e.g., H₂O, CO₂, CH₄) over silicate grains in the interstellar medium are typically on the order of 0.1 eV (Seki & Hasegawa 1983; Suhasaria et al. 2015, 2017), though energies >1 eV are also possible for metallic surfaces (e.g., desorption of potassium from a nickel surface; Błaszczyszyn et al. 1995). We adopt a desorption energy of 0.18 eV and corresponding oscillation frequency of 10¹³ Hz for that of water over silicate (Seki & Hasegawa 1983) for all of our simulations. While the desorption energy could potentially be higher, it would not strongly impact our results, as the magnitude of the nucleation rate depends more on the exponential term in the rate equation (Equation (24)).

Following nucleation, cloud particles can grow via condensation or shrink by evaporation. The condensation and evaporation rates in CARMA are partially determined, like nucleation, by the flux of condensate molecules toward or away from the surface of cloud particles. In addition, growth may be limited by the rate with which particles can conduct away the latent heat released upon condensation, while the opposite process can be applied to limit evaporation. Particles can also grow via coagulation, where cloud particles physically stick to each other upon collision and grow larger at the expense of decreasing number densities. For simplicity we do not consider coagulation in this work, as it is not clear whether exoplanet and brown dwarf particles would coagulate to form larger spheres or aggregates (Lavvas & Koskinen 2017). If the former is the case, then the resulting larger particles would increase the corresponding f_sed value of the cloud. In contrast, formation of low density aggregates may result in clouds that are highly extended vertically, decreasing f_sed. We will investigate these effects in a future publication.

We use a set of exoplanet and brown dwarf atmospheres to conduct our model comparisons. These atmospheres originate from a larger model grid detailed in M. S. Marley et al. (2018, in preparation), and are composed of temperature–pressure-composition profiles in radiative–convective-thermochemical equilibrium, assuming no external insolation and H–He-dominated atmospheres with solar metallicity.

We consider objects with an effective temperature of 400 K and log g = 3.25, 4.25, and 5.25 (g in cgs units), corresponding to masses of 0.72M_J, 8.47M_J, and 44.54M_J, respectively. Figure 1 shows the pressure–temperature profiles of our model atmospheres, which cover the pressure range from 64 bars to 0.18 mbar. For the log g = 5.25 case, we also test 3 separate constant-with-altitude K_zz profiles with values of 10⁶, 10⁷, and 10⁸ cm² s⁻¹ to assess the impact of changing K_zz on the value of f_sed. As K_zz is calculated from Equation (3) in AM01 and we are replacing it with a constant profile in our tests, we use Equation (3) to calculate L, given our constant K_zz profiles to ensure self-consistency. The effective temperature of 400 K was chosen to allow for the formation of potassium chloride (KCl) clouds, which is the only type of cloud we consider in this study. KCl is an optimal choice because it undergoes simple phase transition between vapor and solid or liquid upon cloud formation, rather than relying on chemical reactions (Morley et al. 2012) that could complicate the cloud nucleation process. Note that objects with an effective temperature of 400 K and log g = 5.25 do not yet exist, as objects of such high mass takes longer than the age of the universe to cool to such low temperatures. This is acceptable in the present study since we are only interested in determining the general trend in the cloud distribution, rather than predicting the actual cloudiness of real objects.

Zoom In Zoom Out Reset image size

Given the reliance of the nucleation rate on the surface energy of the condensate and, for heterogeneous nucleation, the size and abundance of condensation nuclei, we also conduct two additional sets of tests—one to investigate how f_sed of the KCl cloud changes with different KCl surface energies, and the other to investigate how the size and downward flux of condensation nuclei affect the f_sed of the resulting KCl cloud. For the surface energy tests, we alter the surface energy of KCl as a function of temperature T (Westwood & Hitch 1963; Janz & Dijkhuis 1969),

$$\begin{eqnarray}&&{\sigma }_{s}^{\mathrm{KCl}}(\mathrm{erg}\,{\mathrm{cm}}^{-2})=160.4-0.07{\rm{T}}({\rm{K}}),\end{eqnarray} \tag{ 9 }$$

by decreasing it by factors of 2 and 4, and increasing it by factors of 2, 3, and 4. For the heterogeneous nucleation tests we assume condensation nuclei radii of 0.1, 1, and 10 nm, and downward fluxes of 10, 100, and 1000 cm⁻² s⁻¹, assuming the condensation nuclei to be composed of silicates with density ∼2 g cm⁻³, similar to the meteoritic dust generated by ablating meteorites in Earth's atmosphere. By comparison, Earth receives 100–300 metric tons of interplanetary dust particles per day (Plane 2012), which corresponds to a flux of 27,000–81,000 cm⁻² s⁻¹ of 1 nm particles, though only a fraction of this mass would be capable of acting as condensation nuclei. The full list of test cases are given in Table 1.

Comparison of the two models are accomplished by first running CARMA to equilibrium. CARMA is initialized with a model atmosphere devoid of KCl vapor or cloud particles. KCl vapor is then allowed to diffuse upwards from depth, given a fixed lower boundary mixing ratio of 0.22 ppmv, corresponding to the total abundance of K in a solar metallicity atmosphere (Lodders 2010). Thus we assume that all K is locked in KCl, which is valid given the pressure–temperature space occupied by the KCl clouds in our model upon their formation (see Figure 1 in this work and Figure 2 in Lodders 1999). Upon reaching saturation, with the saturation vapor pressure estimated by Morley et al. (2012) to be

$$\begin{eqnarray}&&\mathrm{log}{p}_{s}^{\mathrm{KCL}}(\mathrm{bars})=7.611-11382/T({\rm{K}}),\end{eqnarray} \tag{ 10 }$$

KCl nucleates homogeneously and forms a cloud deck. Figure 1 compares our model pressure–temperature profiles to KCl's condensation curve; regions of the model atmospheres cooler than the condensation curve are amenable to KCl cloud formation. After nucleation, the KCl cloud particles are free to grow by condensation, shrink by evaporation, and be transported by sedimentation and diffusion.

We set a zero-flux boundary condition at the top of the model domain in the homogeneous nucleation cases, while a finite flux of condensation nuclei is set for the heterogeneous nucleation cases. A total of 65 particle mass bins are used, with the smallest bin corresponding to particles with radius 1 Å, and the mass doubling for successive bins. The time step sizes for the different cases were tuned to avoid numerical errors in the nucleation and condensation processes, which occur when the time step size is too large, and to avoid extremely long run times when the time step size is too small. Model run times are dominated by the mixing timescale, which is related to K_zz and the scale height. Thus cases with high K_zz and/or high gravity were run using short time steps, and vice versa for cases with lower K_zz and/or low gravity. The time step sizes are given in Table 1. Changing the time step size does not lead to any significant changes to the model results. The models were deemed converged once the column mass of particles and vapor became constant with time.

We calculate q_c for the KCl clouds computed by CARMA as a function of pressure level by summing the masses of all particles within a pressure level and dividing it by the local atmospheric mass density, where is the ratio of the condensate molecular weight to the mean atmospheric molecular weight. r_eff at each pressure level is calculated by finding the weighted average particle radii, with the total cross-sectional area of particles at that level as a function of particle radii acting as weights,

$$\begin{eqnarray}&&{r}_{\mathrm{eff}}=\displaystyle \frac{\int {N}_{r}{r}^{3}{dr}}{\int {N}_{r}{r}^{2}{dr}},\end{eqnarray} \tag{ 11 }$$

where N_r is the number density of particles as a function of particle radius. We combine q_c and r_eff to calculate the optical depth Δτ per atmospheric layer Δz for geometric scatterers by using (AM01)

$$\begin{eqnarray}&&{\rm{\Delta }}\tau =\displaystyle \frac{3}{2}\displaystyle \frac{\epsilon {\rho }_{a}{q}_{c}}{{\rho }_{p}{r}_{\mathrm{eff}}}{\rm{\Delta }}z,\end{eqnarray} \tag{ 12 }$$

where ρ_p is the particle mass density. Summing Δτ from the top of the atmosphere downwards gives the cumulative optical depth.

To find the optimal f_sed value that "fit" the cloud distributions computed by CARMA, we minimize the difference in the pressure level P_0.1 where the two cloud distributions from the two cloud models reach cumulative optical depths of 0.1. This is an optimal comparison criterion because the more diffuse, upper portions of the cloud do not strongly impact transmission, reflection, and emission observations, and thus it is irrelevant to assert that the two models are equal there. Though we do not assert that the two cloud distributions are similar at pressures greater than P_0.1, nor do we constrain r_eff in any way, we find that our comparison criterion ensures that these other quantities are very similar between the two models as well, unless the total cumulative optical depth of the cloud distribution computed by CARMA is less than 0.1. The actual "fitting" of f_sed is done by running the AM01 model 10,000 times, with f_sed increasing by steps of 0.001 each time from f_sed = 0.001 to 10, and comparing P_0.1 for the AM01 model with that of CARMA. The f_sed value corresponding to the minimum in difference of P_0.1 between the two models is deemed the f_sed that best "fits" the CARMA cloud distribution.

We find that the best-fit f_sed decreases with increasing K_zz (Figure 2; Table 1). While increasing the mixing strength within the atmosphere should produce more extended clouds, our results show that microphysical processes serve to further enhance vertical extension. The increased mixing leads to higher fluxes of condensate vapor into the cloud forming region, resulting in increased rates of cloud particle nucleation and the formation of numerous small particles. This decreases the cloud's sedimentation efficiency, thereby creating a more vertically extended cloud. Though the cloud particles do grow by condensation from the upward flux of vapor, most of the growth is isolated to near the cloud base, where r_eff reaches a maximum for the CARMA clouds. This is not seen in the AM01 clouds, where r_eff is roughly constant with depth. The increased vertical extent of the cloud with increased K_zz leads to an upward movement of the τ = 0.1 surface in the atmosphere, such that observations probe decreasing column densities and lower temperatures as K_zz increases (assuming no thermal inversions).

Zoom In Zoom Out Reset image size

Unlike with K_zz, the relationship between f_sed and gravity is non-monotonic (Figure 3; Table 1). This can be understood as the action of multiple processes largely canceling each other out. For example, the sedimentation and mixing timescales of cloud particles can be estimated by

$$\begin{eqnarray}&&{t}_{\mathrm{sed}}=\displaystyle \frac{H}{{v}_{f}}\propto \displaystyle \frac{1}{{g}^{2}}\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{t}_{\mathrm{mix}}=\displaystyle \frac{{H}^{2}}{{K}_{{zz}}}\propto \displaystyle \frac{1}{{g}^{2}},\end{eqnarray} \tag{ 14 }$$

and therefore their combined effects cancel out, assuming constant K_zz. However, condensate vapor is not directly affected by sedimentation, and so the decrease in t_mix at the highest log g case could be due to increased vapor fluxes leading to higher rates of nucleation, as with the case with increased K_zz. In addition, due to lower log g atmospheres having higher temperatures despite the same effective temperature (note, for example, the decreasing altitude of the cloud base with increasing log g), part of the (small) difference in f_sed could be due to temperature effects on the nucleation and growth rates of particles. For instance, lower temperatures may lead to lower rates of nucleation (though this could be compensated by higher saturation ratios), and therefore a slight increase in f_sed for the higher log g case (4.25 versus 3.25) may be expected. Curiously, the pressure level at which τ = 0.1 does not change considerably between the three log g cases. Given that greater scale heights lead to higher column densities, observations should probe higher column densities and temperatures with decreasing gravity.

Zoom In Zoom Out Reset image size

The f_sed values for the cases where K_zz and g are varied (∼0.1) are all considerably smaller than those inferred for brown dwarf clouds from observations (∼1–5; e.g., Stephens et al. 2009), despite physical conditions (e.g., high-g) matching that of brown dwarfs. As Figure 4 (and Table 1) shows, this may be due to material properties that control the nucleation rate of particles. The homogeneous nucleation rate (Equation (19)) is heavily dependent on the value of the exponent in its definition, −F/kT, where F is the formation energy of the nucleated critical cluster and k is the Boltzmann constant—the greater the absolute value of this term, the lower the rate. Expanding out the exponent in the relevant physical variables, we find that

$$\begin{eqnarray}&&-\displaystyle \frac{F}{{kT}}=-\displaystyle \frac{16}{3}\pi \displaystyle \frac{{m}^{2}{\sigma }_{s}^{3}}{{\rho }_{p}^{2}{({kT})}^{3}{\mathrm{ln}}^{2}S},\end{eqnarray} \tag{ 15 }$$

where m is the mass of a condensate molecule and S is the condensate vapor saturation ratio. The strong dependence of the exponent on the surface energy translates to a strong dependence of the cloud distribution on condensate material properties. As σ_s increases, the nucleation rate drops precipitously and the peak of the nucleation rate shifts to cooler (lower pressure) regions of the atmosphere where higher saturation ratios somewhat offset the increase in F from the increased σ_s (Figure 5). Beyond a certain threshold, increasing σ_s no longer results in any significant cloud formation, assuming homogeneous nucleation. The lower nucleation rate with increasing σ_s leads to diminutive clouds composed of mostly a small number of large particles and high f_sed values. The increase in the size of cloud particles is due to the combined effect of the decrease in the number density of cloud particles and the increase in the condensate vapor density, leading to increased growth rates per cloud particle (Figure 6). As σ_s can vary by as much as a factor of 10—for example, between KCl (Equation (9)) and ZnS (1672 erg cm⁻²; Celikkaya & Akinc 1990)—the homogeneous nucleation rate and thus the cloud distribution can vary substantially between different condensates despite similar background atmospheric conditions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We find that the cloud distribution arising from heterogeneous nucleation is dependent on the supply rate and size of condensation nuclei (Figure 7; Table 1). The cloud mass density scales with condensation nuclei flux roughly linearly in logspace, though this effect could change at higher condensation nuclei fluxes due to the limited supply of condensate vapor upwelled from depth. The dependence on condensation nuclei radius r_CN arises from the Kelvin curvature effect, where the saturation vapor pressure over a curved surface p_s^cur is greater than that over a flat surface by an amount defined by

$$\begin{eqnarray}&&\mathrm{ln}\displaystyle \frac{{p}_{s}^{\mathrm{cur}}}{{p}_{s}}=\displaystyle \frac{2M{\sigma }_{s}}{{\rho }_{p}{r}_{\mathrm{CN}}{RT}},\end{eqnarray} \tag{ 16 }$$

where M is the molecular weight of the condensate. Thus smaller condensation nuclei lead to greater saturation vapor pressures and lower nucleation rates. However, this effect is only relevant at sufficiently small r, as the right-hand side of Equation (16) approaches 0 at large r. This is reflected in the lower 2 plots of Figure 7, where a much greater difference exists between clouds nucleating on 0.1 and 1 nm condensation nuclei than between clouds nucleating on 1 and 10 nm condensation nuclei.

Zoom In Zoom Out Reset image size

As the particle size distribution is fully resolved in CARMA using mass bins while AM01 relies on a lognormal parameterization of the size distribution, it is useful to assess how much this difference impacts the resulting observable parameters. In order to isolate the effect of the difference in size distribution shape only, we treat only the CARMA results, and compare the binned results with a lognormal distribution computed from them. We also assume that KCl has a constant real refractive index of 1.5, and that it is purely scattering (i.e., its imaginary refractive index is 0 at all wavelengths). The lognormal parameterization of the CARMA results is computed by first calculating q_c and r_eff at every pressure level, according to the procedure given in Section 2.4, and then calculating r_g from r_eff following Equation (4). The total number density N is then computed using Equation (14) from AM01,

$$\begin{eqnarray}&&N=\displaystyle \frac{3\epsilon {\rho }_{a}{q}_{c}}{4\pi {\rho }_{p}{r}_{g}^{3}}\exp \left(-\displaystyle \frac{9}{2}{\mathrm{ln}}^{2}{\sigma }_{g}\right),\end{eqnarray} \tag{ 17 }$$

with σ_g set to 2 as with AM01.

The resulting lognormal distribution is drastically different from the bimodal structure of the actual size distribution as computed by CARMA (Figure 8), where the lognormal parameterization fails to capture the small particle population originating from nucleation. However, this has minimal impact on the optical properties of the cloud, as the "nucleation mode" particles are far too small to significantly affect the oft-observed wavelengths (e.g., visible to thermal-IR). This is shown by the particle cross-sectional area distribution in the top right plot of Figure 8, which ultimately controls the optical depth of the cloud. Nonetheless, the existence of the nucleation mode shifts the lognormal distribution toward smaller radii compared to the "condensation mode" in the CARMA results, which creates differences in optical properties between the binned size distribution and the lognormal parameterization of roughly 100%, or a factor of 2 at most across all wavelengths, and much smaller differences (<20%) at wavelengths corresponding to particle radii where the optical depth is actually significant (≤2 μm).

Zoom In Zoom Out Reset image size

These differences in optical properties are mostly independent of altitude but are dependent on wavelength following the differences in particle number densities and cross-sectional areas between the two size distributions as a function of particle size. Particles with radius r notably impact the extinction of light of wavelength λ when the shape factor 2πr/λ ≥ 1, as shown by the shaded region in the top plots of Figure 8 for λ = 0.35 μm, the lower wavelength bound of the bottom two plots. As such, most of the differences seen in the bottom two plots are due to particles with radii in the shaded region, and the difference with increasing wavelength can be explained by which size distribution dominates as the shaded region moves to larger radii. For example, at 0.35 μm, the lognormal parameterization slightly dominates over the binned distribution, leading to a negative difference in the optical depth and asymmetry factor. However, as the wavelength is increased such that particles smaller than ∼0.4 μm no longer contribute appreciably to extinction, the higher cross-sectional area of the lognormal parameterization versus the binned distribution for particle radii between 0.05 and 0.4 μm becomes irrelevant, and the binned distribution dominates, leading to a positive difference. At longer wavelengths, there are significantly more large particles in the lognormal parameterization than the binned distribution, and so the former dominates again, leading to negative differences. No difference exists for the single scattering albedo, since we have assumed no absorption.

As pointed out in Section 3.2, the f_sed values associated with brown dwarfs (∼1–5) are considerably higher than the values we have found that best fit our KCl clouds (∼0.1). We inferred, then, that this may be due to brown dwarf clouds having higher surface energies. The surface energies of high temperature condensates iron and forsterite, for example, are 1720 erg cm⁻² (Halden & Kingery 1955) and 1280 erg cm⁻² (Miura et al. 2010), respectively. We can compare the value of the factor −F/kT (Equation (15)) between KCl and iron and forsterite to quantify the effects of material properties and atmospheric conditions on the nucleation rate. For KCl, we set ρ_p = 1.99 g cm⁻³; T = 850 K, the condensation temperature of KCl at ∼1 bar in a solar metallicity atmosphere; m = 74.5m_p, where m_p is the proton mass; and use the surface energy expression in Equation (9). For iron and forsterite, we set ρ_p to 7.875 and 3.27 g cm⁻³, respectively; T = 1800 K, roughly the condensation temperature of iron and forsterite at ∼1 bar; and m to 55.845m_p and 140.69m_p, respectively. The resulting −F/kT values for iron and forsterite are factors of 19 and 284 higher than that of KCl, which greatly decrease the homogeneous nucleation rates of these two species, as −F/kT is in the exponential. Iron and forsterite clouds could also form via heterogeneous nucleation on condensation nuclei composed of TiO₂, as has been proposed by, for example, Helling and Woitke (2006).

Morley et al. (2012) showed that Na₂S clouds forming at 600 K with f_sed ∼ 4–5 could explain the color of T dwarfs, indicating that Na₂S's surface energy (currently unknown) may also be large, especially if it is similar to ZnS. Note that Morley et al. (2012) also included KCl in their model, but the high f_sed rendered their KCl clouds optically thin. Comparing the KCl clouds between their model and our results is also difficult, as the K_zz profile they used is not given. In contrast, Morley et al. (2013, 2015) found that KCl and ZnS clouds in the atmosphere of GJ 1214b must have low f_sed values (∼0.1–0.01) to fit the observed transmission spectrum, which is consistent with our results. However, those models had higher metallicities, which may not be completely compatible with the results of our solar metallicity simulations. It remains to be seen what the effect of having low f_sed for KCl—or different f_sed for different clouds—would be on brown dwarf colors. A thicker cloud would lead to redder colors, but the smaller particles could result in the clouds becoming optically thin in the NIR, making the objects appear bluer.

A particular application of our results is on the brown dwarf L/T transition, where the near-IR colors of brown dwarfs switch from a reddening trend with decreasing luminosity to a trend toward the blue (see reviews by Kirkpatrick 2005; Helling & Casewell 2014). Saumon & Marley (2008) and Stephens et al. (2009) found that the AM01 model can reproduce the L/T transition if the f_sed values of clouds therein increased from ∼2 to ≥4–5. Alternatively, Marley et al. (2010) showed that the L/T transition can be caused by the appearance of holes in the clouds, with the f_sed values of the remaining (patchy) clouds staying ∼2. As the silicate and iron clouds sink to higher pressure levels in the atmosphere, the increase in local temperature should lead to an increase in the nucleation rate and thus smaller particles. However, this effect may be countered by the increased abundance of condensable material at depth, thereby growing particles at a faster rate, leading to larger particles and higher f_sed values. This latter process may be more dominant, as the increase in condensation temperature is small compared to the exponential increase in atmospheric density with increasing depth. To test how f_sed evolves with decreasing T_eff for a single cloud species, we examine the behavior of corundum (Al₂O₃), which is an important condensate in L dwarf atmospheres. For our singular purpose of examining trends in f_sed, we assume particle formation via homogeneous nucleation, as we have done for KCl. However, for future, more detailed exploration of corundum cloud formation, we note that equilibrium chemistry arguments suggest that other aluminum-bearing species, such as Al₂O₂ and AlOH, may be more abundant in the gas phase, limiting corundum cloud formation by direct phase change (Helling & Fomins 2013).

For the saturation vapor pressure of corundum, we use

$$\begin{eqnarray}&&\mathrm{log}{p}_{s}^{{\mathrm{Al}}_{2}{{\rm{O}}}_{3}}(\mathrm{bars})=17.7-45892.6/T({\rm{K}}),\end{eqnarray} \tag{ 18 }$$

which is derived from the condensation curve of corundum given in Wakeford et al. (2017) and the solar abundance of Al from Lodders (2010). We use a surface energy of 900 erg cm⁻² (Dobrovinskaya et al. 2009). The model atmospheres used are from the same grid as those used in the other simulations in this work, with log g = 4.75. We assume a constant K_zz value of 10⁸ cm² s⁻¹ for all cases. As T_eff is decreased, the cloud base sinks deeper into the atmosphere, the particles at the cloud base become larger from the increased supply of condensate vapor, and the best-fit f_sed increases (Figure 9). The more massive iron and silicate clouds behaving similarly could thus explain the increase in f_sed with stellar type at the L/T transition.

Zoom In Zoom Out Reset image size

Alternatively, if silicate and iron clouds depend on condensation nuclei to form, then the L/T transition could be explained by their disappearance, such as due to cold-trapping of the material making up the condensation nuclei (e.g., TiO₂; Lee et al. 2015) at depth with decreasing effective temperature, essentially "starving" the clouds.

The reliance on condensation nuclei could explain another anomalous observation: the extreme reddening of young, directly imaged giant exoplanets, which has been attributed to persistent cloudiness even at low luminosity (see, for example, Figure 7 in Bowler 2016). While low gravity objects are intrinsically hotter than high gravity objects at the same effective temperature, the cloudiness can perhaps be enhanced by a high flux of interplanetary dust from their young stellar systems that create large volumes of meteoritic smoke particles within their atmospheres, which can serve as condensation nuclei. If this is the case, then young substellar objects located in young, dusty stellar systems should be cloudier than their isolated counterparts, regardless of mass. We will investigate this mechanism in a future publication.

Classical theories of homogeneous and heterogeneous nucleation are used for computing the rates of new particle generation (Pruppacher & Klett 1978; Lavvas et al. 2011). For homogeneous nucleation, the rate J_hom in units of new particles per volume per unit time is given as

$$\begin{eqnarray}&&{J}_{\hom }=4\pi {a}_{c}^{2}{\rm{\Phi }}{Zn}\exp (-F/{kT}),\end{eqnarray} \tag{ 19 }$$

where n is the number density of condensate vapor molecules, k is the Boltzmann constant, and T is the temperature. a_c is the critical particle radius defined as

$$\begin{eqnarray}&&{a}_{c}=\displaystyle \frac{2M{\sigma }_{s}}{{\rho }_{p}{RT}\mathrm{ln}S},\end{eqnarray} \tag{ 20 }$$

where M, σ_s, ρ_p, and S are the molecular weight, surface tension (energy), mass density, and saturation ratio of the condensate, respectively, and R is the universal gas constant. F is the energy of formation of a particle with radius a_c, given by

$$\begin{eqnarray}&&F=\displaystyle \frac{4}{3}\pi {\sigma }_{s}{a}_{c}^{2}.\end{eqnarray} \tag{ 21 }$$

In other words, classical homogeneous nucleation theory relates the rate of particle formation with the energy of formation. The greater the energy (i.e., the lower the supersaturation), the lower the nucleation rate. The energy of formation is defined as the balance between the increase in energy associated with enlarging a particle's surface to add vapor molecules and the decrease in energy due to the increase in the volume of the particle. Particles with a radius of a_c are at the cusp of this balance, where further growth leads to a net decrease in energy of formation, and therefore continued existence. Φ is the diffusion rate of vapor molecules to the forming particle, defined as

$$\begin{eqnarray}&&{\rm{\Phi }}=\displaystyle \frac{p}{\sqrt{2\pi {mkT}}},\end{eqnarray} \tag{ 22 }$$

where p is the atmospheric pressure and m is the mass of a vapor molecule. Finally, Z is the Zeldovich factor, which takes into account non-equilibrium effects, such as the evaporation of just-formed particles occurring at the same time as particle formation. Z is given by

$$\begin{eqnarray}&&Z=\sqrt{\displaystyle \frac{F}{3\pi {{kTg}}_{m}^{2}}},\end{eqnarray} \tag{ 23 }$$

where g_m is the number of molecules in particles with radius a_c.

The rate of heterogeneous nucleation, J_het, is defined as

$$\begin{eqnarray}&&{J}_{\mathrm{het}}=4{\pi }^{2}{r}_{\mathrm{CN}}^{2}{a}_{c}^{2}{\rm{\Phi }}{c}_{\mathrm{surf}}Z\exp (-{Ff}/{kT}),\end{eqnarray} \tag{ 24 }$$

where a_c is the critical radius of the initial cluster of condensate molecules ("germ") on the nucleating surface; r_CN is the radius of the condensation nuclei; f is a shape factor given by

$$\begin{eqnarray}&&2f=1+{\left(\displaystyle \frac{1-\mu x}{\phi }\right)}^{3}+{x}^{3}(2-3{f}_{0}+{f}_{0}^{3})+3\mu {x}^{2}({f}_{0}-1),\end{eqnarray} \tag{ 25 }$$

where μ is the cosine of the contact angle between the condensate and the nucleating surface and

$$\begin{eqnarray}&&x=r/{a}_{c}\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\phi =\sqrt{1-2\mu x+{x}^{2}}\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{f}_{0}=(x-\mu )/\phi ;\end{eqnarray} \tag{ 28 }$$

and c_surf is the number density of condensate molecules on the nucleating surface, given by

$$\begin{eqnarray}&&{c}_{\mathrm{surf}}=\displaystyle \frac{{\rm{\Phi }}}{\nu }\exp ({F}_{\mathrm{des}}/{kT}),\end{eqnarray} \tag{ 29 }$$

where ν is the oscillation frequency of adsorbed vapor molecules on the nucleating surface and F_des is the desorption energy of that molecule. The unit for J_het is critical germs per condensation nucleus, and thus to get the number of newly nucleated particles per volume per unit time, J_het must be multiplied by the number density of condensation nuclei.

Condensational growth and evaporation in CARMA takes into account the diffusion of condensate molecules to and away from the cloud particle, such that the rate of change in mass m_p of the particle is

$$\begin{eqnarray}&&\displaystyle \frac{{{dm}}_{p}}{{dt}}=4\pi {R}_{d}^{2}D\displaystyle \frac{d{\rho }_{v}}{{{dR}}_{d}},\end{eqnarray} \tag{ 30 }$$

where R_d is the distance away from the center of the particle, D is the molecular diffusion coefficient of the condensate vapor through the background atmosphere, and ρ_v is the mass density of the condensate vapor. Addition of molecules to a particle releases latent heat, and vice versa for the removal of molecules. The rate of change in temperature of the particle, T_p, arising from condensation and evaporation is

$$\begin{eqnarray}&&{m}_{p}{C}_{p}\displaystyle \frac{{{dT}}_{p}}{{dt}}=L\displaystyle \frac{{{dm}}_{p}}{{dt}}-\displaystyle \frac{{dQ}}{{dt}},\end{eqnarray} \tag{ 31 }$$

where C_p is the heat capacity of the particle, L is the latent heat of evaporation of the condensate, and dQ/dt is the conductive cooling rate of the particle given by

$$\begin{eqnarray}&&\displaystyle \frac{{dQ}}{{dt}}=-4\pi {R}_{d}^{2}{\kappa }_{a}\displaystyle \frac{{dT}}{{{dR}}_{d}},\end{eqnarray} \tag{ 32 }$$

where κ_a is the thermal conductivity of the atmosphere. Combining Equations (30)–(32), along with the ideal gas law and the Clausius–Clapeyron equation,

$$\begin{eqnarray}&&\displaystyle \frac{{{dm}}_{p}}{{dt}}=\displaystyle \frac{4\pi {{rDp}}_{s}(S-1)}{\tfrac{{RT}}{M}+\tfrac{{{DLp}}_{s}}{{\kappa }_{a}T}\left(\tfrac{{LM}}{{RT}}-1\right)},\end{eqnarray} \tag{ 33 }$$

where p_sat is the saturation vapor pressure of the condensate (see Jacobson 2005, for a full derivation). CARMA makes the assumption that LM/RT − 1 ∼ LM/RT, and takes into account several additional effects, such that Equation (33) becomes

$$\begin{eqnarray}&&\displaystyle \frac{{{dm}}_{p}}{{dt}}=\displaystyle \frac{4\pi {rD}^{\prime} {p}_{s}(S-{A}_{k})}{\tfrac{{RT}}{{{MF}}_{v}}+\tfrac{D^{\prime} {{ML}}^{2}{p}_{s}}{{\kappa }_{a}^{\prime} {{RT}}^{2}{F}_{t}}},\end{eqnarray} \tag{ 34 }$$

where A_k is the Kelvin factor, which takes into account the curvature of the particle surface, and is given by

$$\begin{eqnarray}&&{A}_{k}=\exp \left(\displaystyle \frac{2M{\sigma }_{s}}{{\rho }_{p}{RTr}}\right).\end{eqnarray} \tag{ 35 }$$

D' and κ_a' are the molecular diffusional coefficient of the condensate vapor through the atmosphere and the thermal conductivity of the atmosphere, respectively, modified to account for gas kinetics effects near the surface of the particle, and are expressed as

$$\begin{eqnarray}&&D^{\prime} =\displaystyle \frac{D}{1+\lambda {\mathrm{Kn}}^{c}}\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}&&{\kappa }_{a}^{{\prime} }=\displaystyle \frac{{\kappa }_{a}}{1+{\lambda }_{t}{\mathrm{Kn}}_{t}^{c}},\end{eqnarray} \tag{ 37 }$$

where λ and λ_t are given by

$$\begin{eqnarray}&&\lambda =\displaystyle \frac{1.33{\mathrm{Kn}}^{c}+0.71}{{\mathrm{Kn}}^{c}+1}+\displaystyle \frac{4(1-{\alpha }_{s})}{3{\alpha }_{s}}\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&{\lambda }_{t}=\displaystyle \frac{1.33{\mathrm{Kn}}_{t}^{c}+0.71}{{\mathrm{Kn}}_{t}^{c}+1}+\displaystyle \frac{4(1-{\alpha }_{t})}{3{\alpha }_{t}},\end{eqnarray} \tag{ 39 }$$

where α_s and α_t are the sticking coefficient and the thermal accommodation coefficient, respectively, which are set to 1, and Kn^c and Kn^c_t are Knudsen numbers of the condensing gas with respect to the particle, defined as

$$\begin{eqnarray}&&{\mathrm{Kn}}^{c}=\displaystyle \frac{3D}{r}\sqrt{\displaystyle \frac{\pi M}{8{RT}}}\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&{\mathrm{Kn}}_{t}^{c}=\displaystyle \frac{{\mathrm{Kn}}^{c}{\kappa }_{a}}{{rD}{\rho }_{a}\left({C}_{p}-\tfrac{R}{2{\mu }_{a}}\right)},\end{eqnarray} \tag{ 41 }$$

where ρ_a is the atmospheric mass density and μ_a is the mean molecular weight of the atmosphere. Finally, F_v and F_t are the ventilation factors that account for the air density variations around a particle as it sediments in an atmosphere, which increases density in front of it and lowers the density behind it (Toon et al. 1989; Lavvas et al. 2011).

CARMA uses the piecewise parabolic method for treating advection of particles (Colella & Woodward 1984) to apply the rate of change in mass dm_p/dt to evolving the binned particle distribution. In other words, CARMA treats the fluxes between mass bins as if they were fluxes between altitude levels.

Coagulation of particles occurs when the particles bump into each other due to Brownian motion and stick. The rate of coagulation of particles from two populations is given by the product of the coagulation kernel, K_b, the number density of one of the groups of particles, and the number density of the other group of particles. Given particle populations 1 and 2, the coagulation kernel is defined by

$$\begin{eqnarray}&&{K}_{b}=\displaystyle \frac{4\pi ({D}_{1}+{D}_{2})({r}_{1}+{r}_{2})}{\tfrac{{r}_{1}+{r}_{2}}{{r}_{1}+{r}_{2}+\sqrt{{\delta }_{1}^{2}+{\delta }_{2}^{2}}}+\tfrac{4({D}_{1}+{D}_{2})}{({r}_{1}+{r}_{2})\sqrt{{V}_{1}^{2}+{V}_{2}^{2}}}},\end{eqnarray} \tag{ 42 }$$

where D_x, r_x, V_x, and δ_x, x = 1, 2 are the molecular diffusion coefficients, radii, thermal velocities, and interpolation factors of particles from populations 1 and 2, respectively. D_x is defined as

$$\begin{eqnarray}&&{D}_{x}=\displaystyle \frac{{kT}\beta }{6\pi \eta {r}_{x}},\end{eqnarray} \tag{ 43 }$$

where β is the Cunningham slip correction factor (see below) and η is the dynamic viscosity. V_x is given by

$$\begin{eqnarray}&&{V}_{x}=\sqrt{\displaystyle \frac{8{kT}}{\pi {m}_{{px}}}},\end{eqnarray} \tag{ 44 }$$

where m_px, x = 1, 2 is the mass of the particle. Finally, δ_x is expressed as

$$\begin{eqnarray}&&{\delta }_{x}=\displaystyle \frac{{(2{r}_{x}+{\lambda }_{x})}^{3}-{(4{r}_{x}^{2}+{\lambda }_{x}^{2})}^{3/2}}{6{r}_{x}{\lambda }_{x}}-2{r}_{x},\end{eqnarray} \tag{ 45 }$$

where λ_x is the particle mean free path and is given by

$$\begin{eqnarray}&&{\lambda }_{x}=\displaystyle \frac{8{D}_{x}}{\pi {V}_{x}}.\end{eqnarray} \tag{ 46 }$$

δ_x interpolates between the continuum and kinetic regimes (Lavvas et al. 2010).

The sedimentation velocity v_f of a spherical particle of radius r is given by Stoke's fall velocity,

$$\begin{eqnarray}&&{v}_{f}=\displaystyle \frac{2}{9}\displaystyle \frac{{\rho }_{p}{{gr}}^{2}\beta }{\eta },\end{eqnarray} \tag{ 47 }$$

where g is the gravitational acceleration and β again is the Cunningham slip correction factor given by

$$\begin{eqnarray}&&\beta =1+1.246\mathrm{Kn}+0.42\mathrm{Kn}\,{e}^{-0.87/\mathrm{Kn}},\end{eqnarray} \tag{ 48 }$$

where Kn is the Knudsen number of the particle defined as the ratio of the atmospheric mean free path l to r, where l can be written as

$$\begin{eqnarray}&&l=\displaystyle \frac{2\eta }{{\rho }_{a}}\sqrt{\displaystyle \frac{\pi {\mu }_{a}}{8{RT}}}.\end{eqnarray} \tag{ 49 }$$

Note that Kn is different from Kn^c and Kn^c_t, defined in Equations (40)–(41). β ∼1 in the continuum regime (Kn ≪ 1). In the kinetics regime, β is large and mostly linear with Kn. Taking this into account, the sedimentation velocity becomes

$$\begin{eqnarray}&&v=A\displaystyle \frac{{\rho }_{p}{gr}}{{\rho }_{a}}\sqrt{\displaystyle \frac{\pi {\mu }_{a}}{2{RT}}},\end{eqnarray} \tag{ 50 }$$

where A is a constant that is ∼0.5.

The velocity associated with advection is calculated from a user-defined wind speed via the piecewise parabolic method (Colella & Woodward 1984). Diffusive velocities are assumed to be dominated by eddy diffusion, and are calculated and combined with advective velocities to solve for new particle and gas distributions using the algorithms outlined in Toon et al. (1988). The eddy diffusion coefficient as a function of altitude is defined by the user.