The impact of phase equilibrium cloud models on GCM simulations of GJ 1214b

ABSTRACT

We investigate the impact of clouds on the atmosphere of GJ 1214b using the radiatively coupled, phase-equilibrium cloud model EddySed coupled to the Unified Model general circulation model. We find that, consistent with previous investigations, high metallicity (100 × solar) and clouds with large vertical extents (a sedimentation factor of f_sed = 0.1) are required to best match the observations, although metallicities even higher than those investigated here may be required to improve agreement further. We additionally find that in our case which best matches the observations (f_sed = 0.1), the velocity structures change relative to the clear sky case with the formation of a superrotating jet being suppressed, although further investigation is required to understand the cause of the suppression. The increase in cloud extent with f_sed results in a cooler planet due to a higher albedo, causing the atmosphere to contract. This also results in a reduced day–night contrast seen in the phase curves, although the introduction of cloud still results in a reduction of the phase offset. We additionally investigate the impact the Unified Model’s pseudo-spherical irradiation scheme on the calculation of heating rates, finding that the introduction of nightside shortwave heating results in slower mid-latitude jets compared to the plane-parallel irradiation scheme used in previous works. We also consider the impact of a gamma distribution, as opposed to a lognormal distribution, for the distribution of cloud particle radii and find the impact to be relatively minor.

scattering, planets and satellites: atmospheres, planets and satellites: gaseous planets; GJ 1214b

1 INTRODUCTION

GJ 1214b, the archetypal warm Neptune, was one of the first mini-Neptunes/super-Earths to be discovered (Charbonneau et al. 2009) and is notable for its featureless transmission spectrum (Gillon et al. 2014; Kreidberg et al. 2014). It has been the focus of general circulation model (GCM) simulations for a number of years (e.g. Menou 2012; Kataria et al. 2014; Charnay, Meadows & Leconte 2015a; Charnay et al. 2015b; Drummond et al. 2018), and has served as a test case for the limits of applicability of the primitive equations in mini-Neptunes (Mayne et al. 2019) as well as tests of model convergence (Wang & Wordsworth 2020). Despite the inference of clouds in its atmosphere (Gillon et al. 2014; Kreidberg et al. 2014), GCM simulations of GJ 1214b have mostly ignored clouds, leaving their investigation to one-dimensional models. The exception to this is Charnay et al. (2015b) who modelled the atmosphere of GJ 1214b including clouds consisting of radiatively active cloud droplets of fixed radii. They found that particles of |$0.5\, \mathrm{\mu m}$| in size could be lofted high in the atmosphere and their impact on the transmission spectrum improved agreement with observations.

While not capturing the complete dynamics, one-dimensional models allow for improved modelling of the chemistry and cloud physics. Morley et al. (2013) modelled GJ 1214b with a photochemical kinetics model with a parametrized vertical mixing set through the eddy mixing rate, K_zz, including both KCl and ZnS clouds as well as hazes, and found that enhanced metallicity could explain the observations if the sedimentation efficiency was sufficiently low to result in large cloud scale heights. They also found that photochemical haze, with |$1{{\ \rm per\ cent}}$| to |$5{{\ \rm per\ cent}}$| of haze precursors (e.g. C₂H₂, C₂H₄, etc.) forming haze particles, could explain the observations. The latter mechanism does not require atmospheric mixing to carry the particles high into the atmosphere as hazes are expected to form at the low pressures needed to explain observations, and the small particle sizes and associated sedimentation time-scales allow hazes to persist in the upper atmosphere. More sophisticated modelling of the cloud formation process has been done using the carma cloud microphysics code by Gao & Benneke (2018) who modelled KCl and ZnS clouds on GJ 1214b. They found that, in the absence of haze, metallicities of 1000 × solar and mixing rates of |$K_{zz}=10^{6}\, \mathrm{m^2s^{-1}}$| are required to adequately explain the flat transmission spectrum.

The modelling of clouds in three dimensions often takes the form of post-processing GCM models due to the simplicity and relative computational ease (e.g. Robbins-Blanch et al. 2022), allowing for the generation of synthetic observables; however, without the self-consistent coupling of the cloud layer to the atmosphere, the full impact of clouds may be missed. GCM simulations of radiatively coupled clouds on gas giants have primarily focused on hot Jupiters and have been done adopting a range of complexities, from modelling the cloud formation microphysics (e.g. Lee et al. 2016; Lines et al. 2018b) to the use of parametrized cloud models (e.g. Lines et al. 2019; Christie et al. 2021) or fixed particle sizes (e.g. Charnay et al. 2015b; Parmentier, Showman & Fortney 2020; Roman et al. 2021). Although the increased complexity comes at computational expense compared to post-processing, these models find that radiatively coupled clouds can have an impact on the thermal structure of the atmosphere, motivating their necessity.

In this paper, we study clouds in the atmosphere of GJ 1214b using the phase-equilibrium eddysed cloud model radiatively coupled to the Met Office’s unified model (UM) GCM. The use of a one dimensional model coupled to a GCM provides a complementary view of cloud formation to Charnay et al. (2015b) in that while it does not couple the advection and sedimentation directly to the atmosphere via tracers, it allows for a size distribution to be modelled with the particle sizes varying throughout the atmosphere. We use this setup to investigate the impact of clouds on both the atmospheric dynamics as well as the synthetic observations. Overall, we find that supersolar metallicities and large cloud scale heights are required to significantly impact the dynamics and the observables, consistent with previous studies, and that this increased cloud also results in a cooling and contraction of the atmosphere due to the increased albedo.

The structure of the paper is as follows: In Section 2, we outline the components of the simulations. In Section 3, we present the result from the simulations with conclusions and a summary found in Section 4. We also include a number of appendices. In Appendix  A, we demonstrate the impact of switching from a plane-parallel zenith-angle approximation in attenuating incoming stellar radiation to a more physically motivated pseudo-spherical geometry method. In Appendix  B, we explore the impact of particle sizes being distributed in a gamma distribution instead of a lognormal distribution. In Appendix  C, we investigate the impact of KOH on KCl cloud formation. In Appendix  D, we discuss the convergence properties of the simulations presented in the paper. In Appendix  E, we present a number of additional plots.

2 THE MODEL

To understand the impact of clouds on the atmosphere of GJ 1214b, we simulate the atmosphere using the um coupled to the EddySed phase-equilibrium cloud model. We introduce these models and present the specifics of the simulations below.

2.1 The Unified Model GCM

The um solves the full, deep-atmosphere, non-hydrostatic Navier–Stokes equations (Mayne et al. 2014; Wood et al. 2014) and has been adapted to model hot Jupiters (Mayne et al. 2014, 2017) and mini-Neptunes (Drummond et al. 2018; Mayne et al. 2019). Radiative transfer is done using the socrates radiative transfer code based on Edwards & Slingo (1996) which has been adapted and benchmarked in Amundsen et al. (2014, 2017). We assume an atmosphere dominated by H₂ and He with the gas-phase opacity sources being H₂O, CO, CO₂, CH₄, NH₃, Li, Na, K, Rb, Cs, and H₂-H₂ and H₂-He collision-induced absorption (CIA). Opacities are computed using the correlated-k method and ExoMol line lists (Tennyson & Yurchenko 2012; Tennyson et al. 2016). Unless discussed below, the setup is the same as our previous simulations using a coupled version of eddysed (Christie et al. 2021) except with the planetary parameters appropriate for GJ 1214b. The common parameters used across all simulations are found in Table 1.

Table 1.

Common parameters.

	Value
Grid and timestepping
Longitude cells	144
Latitude cells	90
Hydrodynamic time-step	30 s
Radiative transfer
Wavelength bins	32
Wavelength minimum	0.2 μm
Wavelength maximum	322 μm
Damping and diffusion
Damping profile	Polar
Damping coefficient	0.2
Damping depth (ηs)	0.7
Diffusion coefficient	0.158
Planet
Intrinsic temperature	100 K
Initial inner boundary pressure	200 bar
Semimajor axis a	1.23 × 10−2 au
Stellar constant (at 1 au)	3.9996 |$\mathrm{W\, m^{-2}}$|

Previous studies which coupled EddySed to the um (Lines et al. 2018b; Christie et al. 2021) did not include CO₂ as a source for gas-phase opacity as it has a relatively low abundance in atmospheres with solar metallicity; however, as mini-Neptune atmospheres may have enhanced metallicities, meaning CO₂ abundances may become significant, we follow Charnay et al. (2015a) in taking the abundances of CO, CO₂, H₂O, and CH₄ to be in chemical equilibrium through the reactions

$$\begin{eqnarray*} &\mathrm{CO} + 3\mathrm{H_2} \rightleftharpoons \mathrm{H_2O} + \mathrm{CH_4} \nonumber \, \, ,\\ &\mathrm{CO_2} + \mathrm{H_2} \rightleftharpoons \mathrm{H_2O} + \mathrm{CO}\, \, , \end{eqnarray*}$$

(1)

with the solution found using a modification of the analytic chemistry method presented in the appendix of Burrows & Sharp (1999).

Due to the uncertainty in atmospheric composition of GJ 1214b, we consider two different metallicities, solar and 100 × solar, with a focus on the 100 × solar case. As differing atmospheric metallicities can result in very different atmospheric scale heights, thus requiring different computational domain heights, we situate each computational domain such that the mid-point of the vertical axis is located at the observed planetary radius (⁠|$1.7\times 10^7\, \mathrm{m}$|⁠). While this does result in differing gravitational accelerations g at the inner boundary, as the um accounts for spatial variation in g, the gravitational accelerations in regions where the two domains overlap agree. Gao & Benneke (2018) found that a metallicity of 1000 × solar best fit the observations. For such a high metallicity, the assumption that the atmosphere is dominated by hydrogen and helium breaks down, and properly accounting for these atmospheric compositions requires modifications to the socrates configuration. These higher metallicites may be investigated in a future work, but are beyond the scope of the investigation presented here.

To compute the gas constant R, heat capacity c_p, and mean molecular weight, we compute equilibrium chemistry profiles for GJ 1214b for solar and 100 × solar metallicities using the 1D radiative-convective code atmo (Tremblin et al. 2015; Drummond et al. 2016) and average the properties over the model atmosphere (see Table 2). These 1D profiles also serve as the initial conditions for our simulations.

Table 2.

Metallicity dependent parameters.

	Solar	100 × Solar
Gas parameters
Specific gas constant R (J kg−1 K−1)	3.513 × 103	1.565 × 103
Specific heat capacity cp (J kg−1 K−1)	1.200 × 104	5.632 × 103
Mean molecular weight (g mol−1)	2.36	5.31
Simulation parameters
Inner boundary (m)	1.46 × 107	1.58 × 107
Domain height (m)	4.8 × 106	2.4 × 106
g at inner boundary (m s−2)	12.2	10.4
Kzz, 0 (m2 s−1)	7 × 102	3 × 103
Radiation time-step (s)	150	60/150

We also make use of the um’s new pseudo-spherical irradiation scheme (Jackson et al. 2020). In our previous works (e.g. Christie et al. 2021), a plane-parallel scheme was used for attenuating the incoming shortwave radiation. When using the old scheme, the incoming shortwave radiation goes to zero at the terminator, precluding any possibility of nightside heating. The new pseudo-spherical scheme remedies this by computing the attenuation using spherical shells. We find that in allowing for nightside heating there is a reduction of wind speed near the poles which may be physically relevant in scenarios where there is a mid-latitude or polar jet. We present the results from our tests of the new scheme in Appendix  A.

The radiative time-step for most simulations is |$150\, \mathrm{s}$|⁠; however, in the case 100 × solar metallicity simulations with f_sed = 0.5 and 0.1 the increased opacity in the upper atmosphere necessitates a shorter radiative time-step of |$60\, \mathrm{s}$|⁠.

2.2 The eddysed cloud model

To include clouds in our simulations, we couple to the um the one-dimensional, phase-equilibrium, parametrized cloud model eddysed (Ackerman & Marley 2001), treating each vertical column within the GCM independently. The benefit of this approach is that it allows for a level of sophistication beyond the assumption of fixed cloud decks and particle sizes, while still being fast enough to allow long simulation times. The underlying assumption is that clouds exist within an equilibrium where vertical mixing is balanced by sedimentation,

$$\begin{eqnarray*} \frac{\partial q_\mathrm{t}}{\partial z} = -\frac{\left\langle v_\mathrm{sed} \right\rangle }{K_{zz}}q_\mathrm{c} , \end{eqnarray*}$$

(2)

where q_c is the condensate mass mixing ratio, q_v is the vapour mass mixing ratio, and q_t = q_c + q_v. eddysed further assumes that the mass-averaged sedimentation velocity <v_sed > is proportional to a characteristic mixing velocity w_⋆ = K_zz/L_mix, where L_mix is the characteristic mixing scale, which we take to be equal to the pressure scale height H and the constant of proportionality is f_sed.¹ Combined with equation (2), this yields the governing equation

$$\begin{eqnarray*} \frac{\partial q_\mathrm{t}}{\partial z} = -f_\mathrm{sed}\frac{w_\star }{K_{zz}} q_\mathrm{c} = -\frac{f_\mathrm{sed}}{L_\mathrm{mix}}q_\mathrm{c}\, \, \end{eqnarray*}$$

(3)

where z is the vertical coordinate. The cloud distribution is thus determined by a cloud scale height |$L_\mathrm{cloud}=f_\mathrm{sed}^{-1}L_\mathrm{mix}$| and is independent of K_zz, except through the impact of K_zz on the particle sizes and thus the heating rates. It is assumed that any vapour in excess of the local vapour saturation mass mixing ratio q_versus immediately becomes condensate,

$$\begin{eqnarray*} q_\mathrm{c} = \max (0,q_\mathrm{t} - q_\mathrm{vs})\, \, . \end{eqnarray*}$$

(4)

The sedimentation factor f_sed is taken to be constant, as in Lines et al. (2018b) and Christie et al. (2021); however, we note that Rooney et al. (2022) have investigated variations of the eddysed model with an altitude dependent f_sed. Unfortunately, as the model adopted in Rooney et al. (2022) require parameters that are poorly constrained for GJ 1214b, we opt to retain f_sed as a fixed parameter subject to a parameter study.

Particle sizes are assumed to be lognormal with the peak of the distribution r_g given by

$$\begin{eqnarray*} r_\mathrm{g} = r_w f_\mathrm{sed}^{1/\alpha }\exp \left(-\frac{\alpha + 6}{2}\ln ^2\sigma \right), \end{eqnarray*}$$

(5)

where σ parametrizes the distribution width, r_w is the particle radius at which the sedimentation speed is equal to w_⋆ (v_sed(r_w) = w_⋆ = K_zz/H), and α is the scaling of v_sed with particle radius, v_sed ∝ r^α, estimated at r = r_w. While r_g corresponds to the peak of the distribution, the area-weighted effective radius r_eff, more relevant to radiative properties, is given by

$$\begin{eqnarray*} r_\mathrm{eff} = r_w f_\mathrm{sed}^{1/\alpha }\exp \left(-\frac{\alpha +1}{2}\ln ^2\sigma \right). \end{eqnarray*}$$

(6)

The assumption of lognormality of the radius distribution results in the exponential dependence of r_g and r_eff on model parameters. As a point of comparison, we derive the eddysed particle sizes for a gamma distribution and investigate the impact of such a choice in Appendix  B.

The particle sedimentation speed v_sed, used in the computation of r_w, is computed in the terminal velocity limit using

$$\begin{eqnarray*} v_\mathrm{sed}(r) = \frac{2}{9}\frac{\beta g r^2\Delta \rho }{\eta }\, \, , \end{eqnarray*}$$

(7)

where β = 1 + K_n(1.256 + 0.4exp (− 1.1/K_n)) is the Cunningham slip factor, K_n = λ/r is the Knudsen number, λ is the mean free path, g is the gravitational acceleration, Δρ = ρ_c − ρ_a is the relative density of the cloud particles, and

$$\begin{eqnarray*} \eta = \frac{5}{16}\frac{\sqrt{\pi m k_\mathrm{B} T}\left(k_\mathrm{B}T/\epsilon \right)^{0.16}}{1.22\pi d^2}\, \, \end{eqnarray*}$$

(8)

is the dynamic viscosity within the atmosphere (Rosner 2000). The molecular diameter of H₂ is |$d = 2.827\times 10^{-8}\, \mathrm{cm}$|⁠, |$\epsilon =59.7 k_\mathrm{B}\, \mathrm{K}$| is the depth of the Lennard–Jones potential well, k_B is the Boltzmann constant, and m is the mean mass of a particle in the gas phase within the atmosphere.

As in Christie et al. (2021), we assume the mixing rate K_zz can adequately be described by the global average mixing. For GJ 1214b, we adopt mixing rates from Charnay et al. (2015a) who performed a GCM tracer study and fit the resulting distribution to a one-dimensional analytic model to estimate K_zz,

$$\begin{eqnarray*} K_{zz} = \frac{K_{zz,0}}{P^{0.4}_\mathrm{bar}}, \end{eqnarray*}$$

(9)

where P is the pressure and the specific value of K_{zz, 0} is metallicity dependent (see Table 2). The introduction of clouds may alter the structure of the atmosphere and thus the mixing; however, in our previous study of HD 209458b, we did not find the mixing to be significantly impacted (Christie et al. 2021).

The critical vapour saturation mixing ratios q_vs are adopted from Morley et al. (2012) and Visscher, Lodders & Fegley (2006). We only consider KCl and ZnS as most other condensable species are expected to have condensed near the base of or below our computational boundary. KCl clouds are expected to condense directly from gas phase KCl, which at solar metallicities, is the dominant K-bearing species in the gas phase. The critical saturation vapour pressure is given by

$$\begin{eqnarray*} P_\mathrm{vs,KCl} = 10^{7.6106 - 11382\, \mathrm{K}/T}\, \, \mathrm{bar}, \end{eqnarray*}$$

(10)

following Morley et al. (2012). At metallicities of 100 × solar, however, KOH becomes more abundant in the gas phase than KCl in the limit of chemical equilibrium. To ascertain the of KOH impact on cloud formation, we examine KCl condensation in the presence of KOH using Atmo in Appendix  C. We find that the inclusion of KOH does shift the cloud deck slightly; however, once cloud begins to form, both gas-phase KCl and KOH are quickly depleted and the cloud mixing ratio is insensitive to the inclusion of KOH.

Unlike KCl, ZnS is not expected to exist in the gas phase and as a result forms via the chemical reaction

$$\begin{eqnarray*} {\rm Zn} + {\rm H}_2{\rm S} \longrightarrow {\rm H}_2 + {\rm ZnS(s)} \end{eqnarray*}$$

(11)

directly into the solid state. Due to high surface energies, ZnS clouds likely require nucleation sites as precursors to cloud formation (Gao & Benneke 2018). Within our models we assume that such sites exist and ultimately do not impact the final distribution of cloud particles. We adopt a condensation curve from Visscher et al. (2006) and Morley et al. (2012), including a correction for the atmospheric metallicity, namely,

$$\begin{eqnarray*} P_\mathrm{vs,ZnS} = 10^{12.8117 - 15873\, \mathrm{K}/T - \left[\mathrm{Fe/H}\right]} \, \, \mathrm{bar}\, \, . \end{eqnarray*}$$

(12)

As KCl cloud particles may serve as condensation nuclei for ZnS clouds, considering them as separate cloud particles as done here may not be appropriate. This is not accounted for in the eddysed model as implemented but has been investigated by Gao & Benneke (2018) using the carma microphysics code.

To couple eddysed radiatively to the atmosphere, at each point in the atmosphere and for each cloud species the particle size distribution is resolved into 54 size bins between the sizes of |$10^{-9}\, \mathrm{m}$| to |$1.1\times 10^{-3}\, \mathrm{m}$|⁠, allowing for the radiative properties (opacity, single-scattering albedo, and asymmetry parameter) to be computed using pre-calculated tables. The cloud distribution and radiative properties are computed every radiative time-step (see Table 2). As the coupling of eddysed to the atmosphere is done in such a way that the advection of condensate is not considered, the thermal impact of evaporation and condensation is neglected. Latent heat release was included in the model of Charnay et al. (2015b) and was found not to impact the dynamics. A more general discussion of the impact of latent heating in brown dwarf and exoplanetary atmospheres can be found in Tan & Showman (2017).

2.3 Initial conditions and parameter study

To initialize our simulations we use one-dimensional pressure-temperature profiles for GJ 1214b generated using atmo to construct a three-dimensional atmosphere initially in hydrostatic equilibrium and without any winds. Simulations are then run for 600 d without clouds in order to spin-up the atmosphere at which point clouds are enabled. Simulations are then run for an additional 400 d with clouds enabled, except for the 100 × solar metallicity f_sed = 0.5 and 0.1 cases. These two cases are run for 800 d with clouds enabled due to the greater impact of clouds on the atmosphere. Clear sky simulations are run for 1000 d to allow for an appropriate comparison.

For our main parameter study, we investigate two atmospheric metallicities, solar and 100 × solar, and three values of the sedimentation factor f_sed for the eddysed cloud models: f_sed = 1.0, 0.5, and 0.1. These choices of sedimentation factor are informed by previous one-dimensional models of clouds on GJ 1214b which found that clouds need significant vertical extent to match observations (Morley et al. 2013).

To avoid issues associated with the sudden introduction of a new source of opacity, cloud opacities are linearly increased from zero to their full values over the course of the first 100 d after clouds are enabled (i.e. between days 600 and 700).

3 RESULTS

In this section, we present the results of our simulations, first focusing on the impact of cloud on the atmospheric structure followed by an investigation of the impact of clouds on the synthetic observations. A discussion of the convergence of the simulations is found in Appendix  D.

3.1 Atmospheric structure

3.1.1 Thermal structure and the distribution of clouds

The equatorial temperature profiles for the solar and 100 × solar metallicities are shown in Figs 1 and 2, respectively. In each panel, the solid lines correspond to cloudy simulations and the dashed line is the clear sky case for the appropriate metallicity. In the case of solar metallicity (Fig. 1), the introduction of cloud results in only a small (⁠|$\sim 50 \, \mathrm{K}$|⁠) shift to cooler temperatures between 0.1 and 1.0 bar. Due to the relatively small changes in the temperature structure, in the phase equilibrium limit the distribution of clouds resemble the case of post-processing a clear-sky simulation (not shown) with ZnS clouds forming at |$1\, \mathrm{bar}$| and KCl forming at |$0.1\, \mathrm{bar}$|⁠. The existence of ZnS clouds below the KCl cloud deck depends strongly on the existence of nucleation sites, either of KCl, haze, or another form of dust grain (Gao & Benneke 2018), and as such the presence of ZnS clouds at these pressures should be viewed with an increased level of uncertainty. In the case of 100 × solar metallicity (Fig. 2), the addition of clouds has an increased impact over the lower metallicity case due to the scattering of stellar radiation back into space by KCl clouds, reducing the net heating. While in the clear sky case only 0.8 per cent of the incoming stellar radiation is scattered back through the top of the atmosphere, the cloudy simulations show 7.8 per cent, 15.9 per cent and 44.9 per cent of incoming stellar radiation scattered back for f_sed = 1.0, 0.5, and 0.1, respectively. In the f_sed = 1.0 case where clouds are relatively shallow, this results in a shift of |$\sim 100\, \mathrm{K}$| near the ZnS cloud deck at |$1\, \mathrm{bar}$|⁠. As the clouds increase in vertical extent (i.e. with decreasing f_sed), the impact of the reduced heating extends throughout the atmosphere at pressures below |$1\, \mathrm{bar}$|⁠, with up to a |$\sim 300\, \mathrm{K}$| shift around |$P\sim 0.1\, \mathrm{bar}$| for the f_sed = 0.1 case. This shift in temperature results in the base of the clouds appearing at higher pressures (see Fig. 3).

Figure 1.

Pressure-temperature profiles for the case of a solar metallicity atmosphere. The cloudy cases are shown as solid lines while the clear-sky case is shown as a dashed line, for comparison. Dotted grey lines indicate the condensation curves for KCl and ZnS.

Figure 2.

Pressure-temperature profiles for the case of a 100 × solar metallicity atmosphere. The cloudy cases are shown as solid lines while the clear-sky case is shown as a dashed line, for comparison. Dotted grey lines indicate the condensation curves for KCl and ZnS.

Figure 3.

Cloud mass mixing ratio profiles for the case of a 100 × solar metallicity atmosphere. KCl clouds are shown as solid lines while ZnS clouds are dashed lines.

Above the cloud deck, the vertical distribution of cloud exhibits the pressure dependence |$q_\mathrm{c}\propto p^{f_\mathrm{sed}}$| characteristic of the eddysed cloud model. Horizontally, however, there is not significant variation in the cloud mixing ratio. This is in part due to the fact that the although horizontal variation in temperatures are seen the cloudy simulations (see Fig. 2), the temperatures never approach or cross the condensation curve. The lack of observed horizontal variation in cloud is likely also due, in part, to limitations of the modelling approach. A setup that more directly couples the vapour abundance to the local transport processes, instead of assuming a global average mixing rate, has the potential to show more horizontal variation (see e.g. fig. 2 in Charnay et al. 2015b).

The effective particle sizes within the clouds (see equation 6 and Fig. 4) decrease with increased cloud extent as, by assumption, the reduced rainout necessitates smaller radii for the particles to remain suspended. The variation of particle sizes with pressure takes on a form characteristic of EddySed where initially particle sizes increase with decreasing pressure due to the increase in vertical mixing but eventually turn over and begin to decrease with pressure as the atmosphere becomes increasingly more inefficient at supporting cloud particles. The effective particle sizes peak around pressures of 10⁻² bar with sizes of the order of 1–10 microns. As the particle sizes in eddysed are determined solely by the local conditions, the particle sizes become small as atmospheric pressures decrease. These small particles in the upper atmosphere are responsible for the scattering of the incident stellar radiation.

Figure 4.

Profiles of the effective radius r_eff for the case of a 100 × solar metallicity atmosphere. KCl clouds are shown as solid lines while ZnS clouds are dashed lines.

3.1.2 Velocity structures

There has been a lot of investigation of the zonal velocities in GCM simulations of GJ 1214b and the inconsistencies across simulations, although the differences in heating profiles between temperature forcing (Mayne et al. 2019), dual band grey radiative transfer (Menou 2012; Wang & Wordsworth 2020), and multiband correlated-k radiative transfer (Kataria et al. 2014; Charnay et al. 2015a, b) at least in part explains the differences.² In the simulations presented here, we find an equatorial jet forms in all cases except the 100 × solar metallicity, f_sed = 0.1 case (see Fig. 5). For the solar metallicity cases, we find two additional polar jets with amplitudes ∼1.2 to |$1.3\, \mathrm{km\, s^{-1}}$|⁠, qualitatively consistent with the results of Charnay et al. (2015a), with limited impact from the introduction of clouds. The speed of the of the equatorial jet also remains consistent at |$\sim 1.6\, \mathrm{km\, s^{-1}}$|⁠.

Figure 5.

Zonal velocities for the solar metallicity case (top row) and the 100 × solar metallicity case (bottom row). Positive values indicate superrotation. The black dashed lines indicate the contour of zero zonal velocity.

For atmospheres with 100 × solar metallicity, we find superrotation in the mid-latitude and polar regions in the clear sky, f_sed = 1.0, and f_sed = 0.5 cases; however, there does not exist a jet separate from the equatorial jet. Since the models do not include any drag at the lower boundary, a counterrotating flow also forms below the jets, roughly at pressures of 0.1–1 bar as angular momentum is carried upwards into the jet (see Fig. 5 and Showman & Polvani 2011). The speed of the equatorial jet does appear to increase with cloud scale height; however, in the f_sed = 0.1 case, the superrotating jet eventually decays, forming two dayside gyres and a counterrotating wind (see Fig. 6 as well as additional plots in Appendix  E). The mechanism causing the decay of the jet is unclear, although it does appear to be related to the increased opacity and the contracted temperature profile. To explore this, we ran two additional simulations (not shown) with 100 × solar metallicity and f_sed = 0.1 with the cloud opacity introduced at the start and no opacity ramp included. In the first test, we used the same initial temperature pressure profile outlined in Section 2.3 while in the second we used a cooler profile approximating the pressure-temperature profile in the end state of our f_sed = 0.1 case shown in Fig. 2. In the former case, the atmosphere spins up, forming an equatorial jet, before undergoing the same decay seen in our main simulations. In the latter case, no transient equatorial jet forms, with the final state being qualitatively the same as in Fig. 6, right-hand panel. We leave a more detailed investigation of the underlying mechanism and whether it is numerical in origin for future work.

Figure 6.

The horizontal wind speed at |$0.001\, \mathrm{bar}$| for each of the 100 × solar metallicity cases. The horizontal winds show limited variation between cases except for the f_sed = 0.1 case where the nightside gyres move closer to the equator resulting in flow at the mid-latitudes on the nightside moving in a counter-rotating direction.

3.2 Observational diagnostics

To allow for comparison with observations, we have generated synthetic observations using the socrates radiative transfer routines within the um with the spectral resolution increased to |$10\, \mathrm{cm^{-1}}$|⁠. This approach allows for the synthetic observations to be generated self-consistently using the same routines and opacity data used in the heating and cooling calculations. These routines have previously been used in Lines et al. (2018a, 2019) and Christie et al. (2021) and the specifics of the transmission spectrum calculation are presented in Lines et al. (2018a).

3.2.1 Transmission spectra

We calculate transmission spectra for each of the simulations as in Christie et al. (2021), with the exception that the model results are not scaled to agree with the observed value at 1.4 microns, as was the case in Christie et al. (2021). As the results from the solar metallicity simulations show little impact from the introduction of clouds, we opt to focus on the results from the 100 × solar metallicity simulations here.

The left-hand panel of Fig. 7 shows the synthetic observations for the HST WFC3 bandpass covering 1.1–|$1.7\, \mathrm{\mu m}$|⁠. The observations of Kreidberg et al. (2014) are included, with the data shifted up by 0.002 to facilitate comparison.³ A shallow cloud layer, as in the f_sed = 1.0 and f_sed = 0.5 cases, has only a minor impact on the transmission spectrum. The f_sed = 0.1 case shows the best agreement with the data, in so far as the spectrum is relatively flat, however, the 1.4 micron CH₄ feature is still visible. The sharp decrease in the transmission spectrum of the f_sed = 0.1 case near 1.6 microns is due to insufficient precision in the fixed-point recording of the refraction index data for ZnS and is unlikely to be physical in origin (see Querry 1987).

Figure 7.

The simulated transmission spectra for the HST WFC3 (left) and JWST MIRI/LRS (right) bandpasses for the 100 × solar metallicity simulations. The observational data from Kreidberg et al. (2014) are also included in green, and the data have been shifted upwards by 0.002 to facilitate comparison.

To better understand how the specific cloud species are impacting the transmission spectrum in the f_sed = 0.1 case, we breakdown the spectrum further. Fig. 8 shows the transmission spectrum along with post-processed versions of the same run where different combinations of cloud species contribute to the opacity in the transmission spectrum calculation. The transmission spectrum for the clear sky case is shown in grey for comparison. The impact of the clouds on the thermal profile and the extent of the atmosphere as it contracts thermally can be seen in a relatively uniform shift to smaller values in the transmission spectrum from the true clear sky spectrum to the cloudy spectrum where clouds contribute to the heating and cooling but are transparent in transit. Looking instead at the individual contributions to the transmission spectrum, the largest contribution can be seen to be ZnS. As we have made assumptions favourable to the formation of ZnS clouds in that there are always sufficient condensation nuclei and there are no energy barriers limiting condensation (see Gao & Benneke 2018 for a full discussion), the models may be overestimating the impact of ZnS.

Figure 8.

The simulated transmission spectra for the HST WFC3 (left) and JWST MIRI/LRS (right) bandpasses for the 100 × solar metallicity clear sky case and the cloudy 100 × solar metallicity f_sed = 0.1 case allowing for differing contributions to the transmission spectrum for each of the cloud species. The grey lines at the bottom show the major gas-phase contributors in the indicated parts of the spectrum.

We additionally look at the transmission spectra in the JWST MIRI/LRS bandpass (see Fig. 7, right-hand panel), motivated by the possibility of future observations (Bean et al. 2021). As in the previous analysis of the WFC3 bandpass, clouds show the largest impact in the f_sed = 0.1 case. The cooler atmosphere results in a smaller transit across all wavelengths in the bandpass (see Fig. 8, right-hand panel). In the region between |$6$| and |$9\, \mathrm{\mu m}$|⁠, the spectrum is dominated by CH₄ without any direct impact from clouds seen in our synthetic spectrum. At longer wavelengths, however, KCl clouds, and to a lesser extent ZnS clouds, flatten the spectrum, obscuring spectral features (see e.g. the NH₃ feature between |$10$| and |$11\, \mathrm{\mu m}$| in Fig. 8, right-hand panel). The limited impact of ZnS in this bandpass may be in part to the limitations in the Querry (1987) refraction index data, discussed above.

3.2.2 Phase curves

Although there do not currently exist phase curve observations of GJ 1214b, there are planned observations using JWST MIRI/LRS (Bean et al. 2021). These observations cover the wavelength range between 5 and 12 microns and are expected to be an excellent probe of the thermal structure of the planet due to the minimal impact of scattering by aerosols. Fig. 9 (right-hand panel) shows the synthetic phase curves for each simulation with 100 × solar metallicity. As expected given the small impact of the cloud on the thermal structure in the f_sed = 1.0 case, there is minimal difference between the clear sky and f_sed = 1.0 phase curves. The simulation with the most vertically extended clouds – the f_sed = 0.1 case – shows a decrease in the peak of the phase curve associated with the cooler atmospheres which can be contrasted with previous cloudy models of the hot Jupiter HD 209458b where dayside temperatures increased with cloud scale height (i.e. with decreasing f_sed) resulting in higher peaks in phase curves (Christie et al. 2021). Despite these differences, the introduction of cloud still results in an increase in contrast in the phase curves between the minima and maxima (a factor of 2.1 for the clear case compared to 4.0 for the f_sed = 0.1 case) and a shift of the peak back towards 180°, as discussed in Parmentier et al. (2020) and seen in our previous work (Christie et al. 2021). The lack of offset is also likely due to the lack of equatorial jet pushing the hotspot westward. We also note that, as expected for the given wavelength range, the phase curve is probing thermal emission with only |$\sim 10^{-9}{{\ \rm per\ cent}}$| of the emission coming from reflected starlight in the clear sky case and increasing to |$\sim 2{{\ \rm per\ cent}}$| in the f_sed = 0.1 case.

Figure 9.

Simulated phase curves for each of the 100 × solar metallicity simulations for the HST WFC3 bandpass (left, with a wavelength range between 1.1 and |$1.7\, \mathrm{\mu m}$|⁠) and the JWST MIRI/LRS bandpass (right, with a wavelength range between 5 microns and |$12\, \mathrm{\mu m}$|⁠).

As the effective temperature of GJ 1214 is |$3026\, \mathrm{K}$| (Charbonneau et al. 2009), the stellar blackbody peaks at |$\sim 1\, \mathrm{\mu m}$| allowing for the possibility of scattered starlight by clouds contributing significantly to the HST WFC3 phase curve (see Fig. 9, left-hand panel). For the clear sky, 100 × solar metallicity case where there is limited scattering of starlight, the planetary thermal emission dominates the phase curve, with |$F_\mathrm{p}/F_\mathrm{s}\, \sim 2.3$| to 2.7 × 10⁻⁷. With increased cloud extent, however, the scattered stellar component begins to dominate the phase curve, and for f_sed = 0.1, the phase curve peaks with F_p/F_s = 2.6 × 10⁻⁵, with 99.8 per cent of the planetary flux coming from the scattered stellar component.

4 DISCUSSION AND CONCLUSIONS

In this paper, we have investigated the impact of clouds on the dynamics and observables of the warm Neptune GJ 1214b, specifically through the coupling the one-dimensional phase equilibrium eddysed cloud model to the UM GCM. Consistent with previous investigations, we find that increased metallicity and increased cloud vertical extent (i.e. decreased f_sed) are necessary to impact the pressure temperature profiles and synthetic observables. The most significant impact was for our 100 × solar metallicity case with f_sed = 0.1 where |$\sim 45{{\ \rm per\ cent}}$| of the incident stellar radiation was scattered back into space resulting in cooling and contraction of the atmosphere. These clouds, primarily the ZnS component, also increased the atmospheric opacity within the spectral windows but did not entirely obscure the spectral features (see e.g. 1.4 microns in Fig. 7). To reproduce the flat spectrum of Kreidberg et al. (2014), the cloud content either needs to be further increased through alterations to the cloud model or the metallicity further increased, consistent with the conclusions of Gao & Benneke (2018).

Dynamically, the inclusion of clouds results in an increase in the speed of the equatorial jet, except in the 100 × solar metallicity case with f_sed = 0.1, where we see a suppression of the equatorial jet and the formation of two dayside gyres. It is unclear whether the latter configuration is stable or whether it will continue to evolve, or to what extent it is numerical in origin. Understanding this better likely requires further improving our cloud models, and potentially moving to tracer-based models that better capture the localized dynamics and improved microphysics to properly model the particle size distribution. The move to a tracer-based model would also facilitate the inclusion of photochemical hazes, which do not fit naturally in the eddysed framework.

ACKNOWLEDGEMENTS

We thank the anonymous referee for comments that greatly improved the quality of this paper. The observational data were retrieved from Dr. Hannah Wakeford’s online archive at www.stellarplanet.org. Material produced using Met Office Software. This research made use of the ISCA High Performance Computing Service at the University of Exeter. This work was performed using the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the Science and Technology Facilities Council DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via Science and Technology Facilities Council capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. DiRAC is part of the National e-Infrastructure. This work was partly funded by the Leverhulme Trust through a research project grant [RPG-2020-82], a Science and Technology Facilities Council consolidated grant [ST/R000395/1] and a a UK Research and Innovation Future Leaders Fellowship [grant number MR/T040866/1]. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising.

DATA AVAILABILITY

The research data supporting this publication are archived on the Harvard Dataverse and are available at doi.org/10.7910/DVN/8XAEPY.

Footnotes

1

As discussed in Ackerman & Marley (2001) where the EddySed model is applied to brown dwarfs, in convective atmospheres the characteristic scale L_mix is taken to be the convective mixing scale. As we are applying the model to radiation-dominated regions of the atmosphere, we opt for the simpler assumption of L_mix = H.

2

Also relevant is the investigation of the impact of choice of radiative transfer scheme on models of the hot Jupiter HD 209458b by Lee et al. (2021). Relevant to the discussion here, they found differing temperature structures in the upper atmosphere and differing widths of the central jet between their dual grey and correlated-k models.

3

As the appropriate radius for the |$200\, \mathrm{bar}$| inner boundary is unknown, a uniform shift in the transmission spectrum approximates a move in the location of the inner boundary. To properly address this issue, the appropriate value of the inner boundary radius could be found by searching the parameter space of possible radii; however, this would be computationally expensive for what is likely little actual gain.

4

The gamma distribution is someitmes presented as the Hansen distribution (e.g. Hansen 1971; Burningham et al. 2021) or the potential exponential distribution (e.g. Helling et al. 2008) in different forms in the literature.

APPENDIX A: PSEUDO-SPHERICAL IRRADIATION IN THE UM

In this work, we take advantage of the UM’s implementation of pseudo-spherical irradiation (Jackson et al. 2020) which calculates the attenuation of the incoming stellar radiation using spherical shells. This is in contrast to the previous plane-parallel implementation where the attenuating column is approximated using the vertical column N, corrected using the zenith angle α, yielding an effective attenuating column N/cos α. As a result, in the plane-parallel implementation, shortwave heating goes to zero at the terminator, without any possibility of nightside heating, an issue remedied by the pseudo-spherical treatment. This treatment has also been applied to the study of the impact of flares on ‘Earth–like’ terrestrial exoplanets in Ridgway et al. (submitted to MNRAS). In this section, we compare GJ 1214b clear sky results using the two different irradiation schemes to ascertain the impact of the change, and we demonstrate that for our test case the switch to pseudo-spherical irradiation results in slower jets near the poles.

For this test, we ran each of the two cases for 1000 d using the clear-sky solar-metallicity configuration outlined in the paper. Qualitatively, we see the same velocity structures forming in both cases (see Fig. A1), with a fast equatorial jet and two polar jets with the equatorial jet having a maximum speed of |$\sim 1.5\, \mathrm{km\, s^{-1}}$|⁠. The impact of the new irradiation scheme is primarily felt at the poles where shortwave radiation heating the nightside reduces the temperature contrast around the poles resulting in a slower polar jet (⁠|$1.2\, \mathrm{km\, s^{-1}}$| in the plane-parallel case versus |$1.0 \, \mathrm{km\, s^{-1}}$| in the spherical case). This can be seen in Figs A2 and A3 where the plane-parallel case shows a larger decrease in the temperature approaching the pole.

Figure A1.

The zonally averaged azimuthal velocity for a simulation with plane-parallel irradiation (left) and pseudo-spherical irradiation (right). Both simulations are run without clouds for 1000 d. The black dashed line indicates the contour of zero zonal velocity.

Figure A2.

The zonally averaged temperature for a simulation with plane-parallel irradiation (left) and pseudo-spherical irradiation (right). Both simulations are run without clouds for 1000 d.

Figure A3.

The temperature at |$10^{-4}\, \mathrm{bar}$| for a simulation with plane-parallel irradiation (left) and pseudo-spherical irradiation (right). Both simulations are run without clouds for 1000 d.

We estimate the energy deposited in the atmosphere based on the heating rates and find that the spherical irradiation scheme results in |$\sim 6{{\ \rm per\ cent}}$| more atmospheric heating than the plane-parallel case. While in theory the atmospheric heating should be insensitive to the radiation scheme, losses do occur due to various assumptions being made, and we will investigate this further in a follow-up paper.

We note that while using spherical irradiation makes modest changes to the structure of the atmosphere in the regions of interest for studies like this one, it does better characterize the underlying radiative transfer within the atmosphere, and it may offer some improvement in stability as we see a reduction in both horizontal and vertical velocities near the poles.

APPENDIX B: eddysed WITH A GAMMA DISTRIBUTION

Without direct observations of clouds on exoplanets, it is difficult to constrain the appropriate particle size distribution to use in the modelling of clouds on these planets; however, in water clouds on Earth, where in situ observations of raindrop and snowflake sizes are possible, the size distributions are often fit to mono- or multimodal lognormal (e.g. Ackerman et al. 2000), exponential (e.g. Marshall & Palmer 1948), or gamma distributions (e.g. Ulbrich 1983). Informed by this, models of exoplanetary clouds that do not explicitly track the particle sizes often assume one of these distributions (e.g. Ackerman & Marley 2001; Helling, Woitke & Thi 2008; Charnay 2018).

In this appendix, we look at the impact of switching from a lognormal particle size distribution to a gamma particle size distribution within the eddysed framework. While the assumption of a lognormal particle size distribution is often the simplest from a computational and analytical standpoint, eddysed can be formulated using a gamma distribution⁴ which is presented below.

B1 Derivation

Considering a layer of cloud within the atmosphere, we take the number distribution of cloud particles to be dn/dr = Nf(r; A, B), where N is the total number of particles and f(r; A, B) is the gamma distribution,

$$\begin{eqnarray*} f(r; A,B) = \frac{B^Ar^{A-1}}{\Gamma (A)}\exp \left(-Br\right). \end{eqnarray*}$$

(B1)

The two parameters A and B are more often written as α and β; however, as those are used elsewhere in paper, we opt for this simple substitution in nomenclature. In Ackerman & Marley (2001), the variance of the distribution is |$\operatorname{\operatorname{Var}}\left[\log X\right] = \log ^2\sigma$| where σ parametrizes the width of the distribution. For a gamma distribution, the variance of a gamma-distributed variable X depends only on the parameter A,

$$\begin{eqnarray*} \operatorname{\operatorname{Var}}\left[\log X\right] = \psi ^{(1)}(A), \end{eqnarray*}$$

(B2)

where ψ⁽¹⁾(x) is the trigamma function. Following the method of Ackerman & Marley (2001), we can then fix the value of A based on a prescribed σ,

$$\begin{eqnarray*} \psi ^{(1)}(A) = \log ^2\sigma \, \, . \end{eqnarray*}$$

(B3)

As σ is taken to be constant throughout the atmosphere in the eddysed model, A only needs to be calculated at the beginning of the simulation, limiting computational overhead. For the default value of σ = 2 used here, A ≈ 2.54278. With A fixed, we can now calculate B based on the sedimentation arguments of Ackerman & Marley (2001). The sedimentation factor is defined as

$$\begin{eqnarray*} f_\mathrm{sed} = \frac{\int _{0}^{\infty } r^{3+\alpha }(\mathrm{ d}n/\mathrm{ d}r)\mathrm{ d}r}{r_w^\alpha \int _{0}^{\infty } r^3 (\mathrm{ d}n/\mathrm{ d}r)\mathrm{ d}r}. \end{eqnarray*}$$

(B4)

This can be rewritten as

$$\begin{eqnarray*} f_\mathrm{sed} = \frac{\operatorname{\operatorname{E}}\left[r^{3+\alpha }\right]}{r_w^\alpha \operatorname{\operatorname{E}}\left[r^3\right]} = \frac{1}{(Br_w)^\alpha }\frac{\Gamma (A+3+\alpha)}{\Gamma (A+3)}, \end{eqnarray*}$$

(B5)

where E[x] is the expected value of x. We can then solve for B,

$$\begin{eqnarray*} B = \frac{1}{r_w}\left[\frac{\Gamma (A+3+\alpha)}{f_\mathrm{sed}\Gamma (A+3)}\right]^{1/\alpha }. \end{eqnarray*}$$

(B6)

In the rare cases where α is an integer, the identity Γ(z + 1) = zΓ(z) can be used to eliminate the gamma function entirely and reduce the equation to

$$\begin{eqnarray*} B = \frac{1}{r_w}\left[\frac{(A+2+\alpha)\cdots (A+3)}{f_\mathrm{sed}}\right]^{1/\alpha }\, \, . \end{eqnarray*}$$

(B7)

While this is worth noting, we will not use this simplification for the rest of the derivation as α is generally not an integer.

We can then compute the total number of particles based on mass conservation

$$\begin{eqnarray*} q_\mathrm{c}\rho _a = \frac{4}{3}\pi N\rho _\mathrm{c}\int _{0}^{\infty }r^3 f(r;A,B)\mathrm{ d}r, \end{eqnarray*}$$

(B8)

resulting in

$$\begin{eqnarray*} N = \frac{3q_\mathrm{c}\rho _\mathrm{a}B^3}{4\pi \rho _\mathrm{c}(A+2)(A+1)A}. \end{eqnarray*}$$

(B9)

The area-weighted effective radius r_eff in a cloud layer is given by

$$\begin{eqnarray*} r_\mathrm{eff} = \sqrt{\operatorname{\operatorname{E}}[r^2]} = \left(\frac{A(A+1)}{B^2}\right)^{1/2}\, \, . \end{eqnarray*}$$

(B10)

The scaling of r_eff with the size parameter σ is very different in the case of a gamma distribution compared to the lognormal distribution. For small values of σ, the two distributions scale similarly with the distributions becoming sharply peaked with |$r_\mathrm{eff} \sim r_wf_\mathrm{sed}^{1/\alpha }$|⁠; however, for wider distributions (e.g. for σ greater than 2, see Fig. B1), r_eff decreases slower with σ in the case of the gamma distribution case compared to the lognormal distribution case. For wider distributions, we may, as a result, expect larger differences in particle sizes and optical properties between cloud layers with a lognormal distribution than with cloud layers with a gamma distribution. Coincidentally, we find that for reasonable values of α the effective radii r_eff for the two distributions are similar for σ ∼ 2, the value used here and widely in the literature.

Figure B1.

The effective particle radii r_eff for eddysed formulated with a lognormal distribution (solid lines) and a gamma distribution (dashed lines). For the purposes of this plot, α = 2 has been assumed.

To illustrate how using eddysed with a gamma distribution differs from the default lognormal distribution case, the lognormal and gamma distributions are plotted in Fig. B2 for f_sed = 1.0 and α = 2 and three values of σ. Compared to the lognormal distribution case, both location of the peak of the gamma distribution and the number of particles show a much weaker dependence on σ with the peak of the gamma distribution staying close to r_w and the width being modulated by the low-mass tail. This weaker dependence and the skewness of the gamma distribution means that while varying σ does change the particle numbers, the condensate mass remains more concentrated in the high-mass end of the distribution compared to the lognormal distribution where a wider distribution results in more massive particles being traded for exponentially more particles at the peak of the distribution.

Figure B2.

A comparison of the lognormal (red) and gamma (blue) distributions for three values of σ. In each case, f_sed = 1.0 and α = 2 have been assumed.

B2 Results

To examine the impact of the change in distribution, we focus on the 100 × solar metallicity f_sed = 0.1 case as our analysis has already shown cloud in this case to have a noticeable impact, although we only compare the results after 1000 d. Although larger values of σ may result in larger differences between the two cases, we opt to keep σ = 2 as it represents the standard value used here and in other studies. Comparing the two cases, we find negligible differences in the equatorial temperature structure; however, at mid-latitudes we find cooler temperatures on the night side and along the morning terminator in the case of a gamma distribution (see Figs B3 and B4). The differences occur around |$P\sim 10^{-3}\, \mathrm{bar}$|⁠, with the lognormal distribution case being |$\sim 100\, \mathrm{K}$| hotter than the gamma distribution case. We note that the flow at mid-latitudes is somewhat different between the lognormal and gamma distribution cases, with the equatorial jet decaying at a slower rate in the case of the gamma distribution; however, longer run times are required to understand to what extent these differences are transient. This does not have a noticeable impact on the transmission spectrum, either through the change in atmospheric temperature or through any change in cloud opacity used in the transmission spectrum calculation.

Figure B3.

The temperature profiles at a latitude of 45° for the 100 × metallicity, f_sed = 0.1 case for a lognormal distribution (solid lines) and a gamma distribution (dashed lines) assuming σ = 2. Equatorial profiles are not included as they did not differ noticeably between the two distributions.

Figure B4.

The temperatures at |$1\, \mathrm{mbar}$| for the lognormal distribution case (left) and the gamma distribution case (right).

Switching to a gamma distribution may result in larger differences from the lognormal distribution case in parts of the parameter space not investigated here (e.g. the case of very wide distributions). We have opted to not investigate these cases as wide distributions quickly run into the issue of significant fractions of the particle distribution falling outside of the permitted or physically plausible range of particle sizes at low pressures. Relaxing or modifying the eddysed assumption of a fixed distribution width may allow this issue to be avoided, but also it may provide an avenue for moving beyond eddysed more generally.

APPENDIX C: COMPARISON OF CONDENSATION CURVES WITH atmo

KCl clouds form directly from the condensation of gaseous KCl which, at solar metallicity, is the dominant K-bearing species in the atmosphere. As metallicity increases, however, the abundance of KOH increases, eventually exceeding KCl in abundance. To ascertain what impact this may have, we use the atmo chemistry code to compute the chemical equilibria for the initial 100 × solar temperature profile for two cases. In the first, we allow KOH to form, but in the second, we artificially exclude KOH. When KOH is allowed to form, there is a reduction in gaseous KCl, resulting in it becoming supersaturated at a slightly lower pressure (see Fig. C1). In either case, once KCl clouds begin to form, abundances of both KCl and KOH quickly drop above the cloud deck, resulting in the KCl cloud profile being the same in either case, except for the slight shift in the cloud deck.

Figure C1.

The chemical equilibrium for the major K-bearing species for the 100 × solar case, generated using the atmo code. In one case, KOH has been included (solid line) while in the other KOH has been removed from the possible species (dashed line). Although when KOH is allowed to form it becomes the dominant K-bearing species below the KCl cloud deck, it has a negligible impact on the formation of KCl(s).

For the purposes of this paper where we only consider chemical and phase equilibrium cases, we take this to be sufficient indication that there is only a minor impact from ignoring KOH. Non-equilibrium models may be impacted, but that is beyond the scope of this paper.

APPENDIX D: CONVERGENCE AND CONSERVATION

In this section, we look at the convergence of the simulations. Examining first at the evolution of the peak zonal velocities (Fig. D1), the clear sky simulations quickly approach a near-steady state within the first few hundred days. There remains a gradual continued spin-up of the atmosphere as the models lack any explicit physically motivated dissipation mechanism. While it is possible that the simulations may eventually reach a state in which the numerical dissipation balances the physical forcing, previous experiments (not shown) have found that the simulations become unstable before this occurs.

Figure D1.

The evolution of the zonal velocity in each of the simulations. For each cloudy simulation, the clouds were enabled after 600 d. The transition in the flow in the 100 × solar metallicity f_sed = 0.1 case can be clearly seen in the right plot.

For the solar metallicity case, the introduction of cloud at 600 d results in a shift in the zonal velocities of |$\sim 0.2 \, \mathrm{km\, s^{-1}}$| with the simulations quickly reaching a new near-steady state. The 100 × solar metallicity case, however, sees fluctuations in zonal velocity with time after the introduction of clouds at 600 days, although the zonal velocities do not exhibit a long term diverging trend.

Although the um does not explicitly conserve axial angular momentum (to machine accuracy), we see it remains relatively well conserved for the solar metallicity simulations, losing |$\sim 0.05{{\ \rm per\ cent}}$| over the 1000 d of simulation (see Fig. D2, left-hand panel). The 100 × solar metallicity simulations, on the other hand, show a larger loss of angular momentum after clouds are enabled, with the f_sed = 0.1 case losing |$\sim 1.4{{\ \rm per\ cent}}$| of the initial angular momentum and the majority of that loss occurring during the final 800 d (see Fig. D2, right-hand panel). This is possibly, in part, due to the poor conservation during the contraction of the atmosphere as well as the polar diffusion scheme dissipating angular momentum from the atmosphere. This will be investigated in a future work.

Figure D2.

The evolution of the axial angular momentum within the computational domain, normalized to the initial value, for each simulation presented in the paper.

The transition between the flow structures in the 100 × solar metallicity, f_sed = 0.1 case can be seen in Fig. D1 (right-hand panel), with the transition in the flow strucutre being illustrated in Fig. E2.

Figure E1.

The evolution of the 100 × solar metallicity f_sed = 0.1 simulation at 0.001 bar over the length of the period in which clouds are included.

Figure E2.

The same as in Fig. E1, except for the 60 d around the transition in the flow structure.

APPENDIX E: ADDITIONAL PLOTS

In this appendix, we include a number of additional plots illustrating the transition in flow structure seen in the 100 × solar metallicity f_sed = 0.1 simulation. Fig. E1 shows the evolution of the flow over the length of the simulation beginning at when clouds are enabled. Fig. E2 shows evolution over the period between 1220 and 1280 d when the flow undergoes the largest change in structure.