The Impact of Phase Equilibrium Cloud Models on GCM Simulations of GJ~1214b

We investigate the impact of clouds on the atmosphere of GJ~1214b using the radiatively-coupled, phase-equilibrium cloud model {\sc EddySed} coupled to the {\sc Unified Model} general circulation model. We find that, consistent with previous investigations, high metallicity ($100\times$ solar) and clouds with large vertical extents (a sedimentation factor of $f_\mathrm{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_\mathrm{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_\mathrm{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 the {\sc 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 log-normal distribution, for the distribution of cloud particle radii and find the impact to be relatively minor.


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 (Kreidberg et al. 2014;Gillon 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 et al. 2015a,b;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 (Kreidberg et al. 2014;Gillon 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) which modelled the atmosphere of GJ 1214b including clouds consisting of radiatively-active cloud droplets of fixed radii.They found that particles of 0.5 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 parametrised vertical mixing set through the eddy mixing rate,   , including both KCl and ZnS clouds as ★ E-mail: d.christie@exeter.ac.uk 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% to 5% of haze precursors (e.g., C 2 H 2 , C 2 H 4 , 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 timescales allow hazes to persist in the upper atmosphere.More sophisticated modeling 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   = 10 6 m 2 s −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 selfconsistent coupling of the cloud layer to the atmosphere, the full impact of clouds may be missed.GCM simulations of radiativelycoupled 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 parametrised cloud models (e.g., Lines et al. 2019;Christie et al. 2021) or fixed particle sizes (e.g., Charnay et al. 2015b;Parmentier et al. 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 E S cloud model radiatively coupled to the Met Office's U M (UM) GCM.The use of a one dimensional could 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 log-normal 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.

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

The U M 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(Mayne et al. , 2017) ) and mini-Neptunes (Drummond et al. 2018;Mayne et al. 2019).Radiative transfer is done using the S radiative transfer code based on Edwards & Slingo (1996) which has been adapted and benchmarked in Amundsen et al. (2014Amundsen et al. ( , 2017)).We assume an atmosphere dominated by H 2 and He with the gas-phase opacity sources being H 2 O, CO, CO 2 , CH 4 , NH 3 , Li, Na, K, Rb, Cs and H 2 -H 2 and H 2 -He collision-induced absorption (CIA).Opacities are computed using the correlated-k method and E M 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 E S (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.
Previous studies which coupled E S to the UM (Lines et al. 2018b;Christie et al. 2021) did not include CO 2 as a source for gasphase opacity as it has a relatively low abundance in atmospheres with solar metallicity; however, as mini-Neptune atmospheres may have enhanced metallicities, meaning CO 2 abundances may become significant, we follow Charnay et al. (2015a) in taking the abundances of CO, CO 2 , H 2 O, and CH 4 to be in chemical equilibrium through the reactions 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 midpoint of the vertical axis is located at the observed planetary radius (1.7 × 10 7 m).While this does result in differing gravitational accelerations  at the inner boundary, as the UM accounts for spatial variation in , 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 S 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 , heat capacity   , and mean molecular weight, we compute equilibrium chemistry profiles for GJ 1214b for solar and 100× solar metallicities using the 1D radiative-convective code A (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.
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 timestep for most simulations is 150 s; however, in the case 100× solar metallicity simulations with  sed = 0.5 and 0.1 the increased opacity in the upper atmosphere necessitates a shorter radiative timestep of 60 s.

The E S Cloud Model
To include clouds in our simulations, we couple to the UM the one-dimensional, phase-equilibrium, parametrised cloud model E -S (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, where  c is the condensate mass mixing ratio,  v is the vapour mass mixing ratio, and  t =  c +  v .E S further assumes that the mass-averaged sedimentation velocity  sed is proportional to a characteristic mixing velocity  ★ =   / mix , where  mix is the characteristic mixing scale, which we take to be equal to the pressure scale height  and the constant of proportionality is  sed . 1 Combined with Equation 2, this yields the governing equation 1 As discussed in Ackerman & Marley (2001) where the E S model is applied to brown dwarfs, in convective atmospheres the characteristic scale  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  mix =  .
where  is the vertical coordinate.The cloud distribution is thus determined by a cloud scale height  cloud =  −1 sed  mix and is independent of   , except through the impact of   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  vs immediately becomes condensate,  c = max(0,  t −  vs ) . (4) The sedimentation factor  sed is taken to be constant, as in Lines et al. (2018b); Christie et al. (2021); however, we note that Rooney et al. (2021) have investigated variations of the E S model with an altitude dependent  sed .Unfortunately, as the model adopted in Rooney et al. (2021) requires parameters that are poorly constrained for GJ 1214b, we opt to retain  sed as a fixed parameter subject to a parameter study.
Particle sizes are assumed to be log-normal with the peak of the distribution  g given by where  parametrises the distribution width,   is the particle radius at which the sedimentation speed is equal to  ★ ( sed (  ) =  ★ =   /), and  is the scaling of  sed with particle radius,  sed ∝   , estimated at  =   .While  g corresponds to the peak of the distribution, the area-weighted effective radius  eff , more relevant to radiative properties, is given by The assumption of log-normality of the radius distribution results in the exponential dependence of  g and  eff on model parameters.As a point of comparison, we derive the E S particle sizes for a gamma distribution and investigate the impact of such a choice in Appendix B.
The particle sedimentation speed  sed , used in the computation of   , is computed in the terminal velocity limit using where  = 1 +  n (1.256 + 0.4 exp(−1.1/n )) is the Cunningham slip factor,  n = / is the Knudsen number,  is the mean free path,  is the gravitational acceleration, Δ =  c −  a is the relative density of the cloud particles, and is the dynamic viscosity within the atmosphere (Rosner 2000).The molecular diameter of H 2 is  = 2.827 × 10 −8 cm,  = 59.7B K is the depth of the Lennard-Jones potential well,  B is the Boltzmann constant, and  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   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 onedimensional analytic model to estimate   , where  is the pressure and the specific value of  ,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  vs are adopted from Morley et al. (2012) and Visscher et al. (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  vs,KCl = 10 7.6106−11382 K/ bar, (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 A 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 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); Morley et al. (2012), including a correction for the atmospheric metallicity, namely,  vs,ZnS = 10 12.8117−15873 K/ −[Fe/H] bar .
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 E S model as implemented but has been investigated by Gao & Benneke (2018) using the CARMA microphysics code.
To couple E S 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 m to 1.1 × 10 −3 m , allowing for the radiative properties (opacity, singlescattering albedo, and asymmetry parameter) to be computed using pre-calculated tables.The cloud distribution and radiative properties are computed every radiative timestep (see Table 2).As the coupling of E S 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).

Initial Conditions and Parameter Study
To initialise our simulations we use one-dimensional pressuretemperature profiles for GJ 1214b generated using A to construct a three-dimensional atmosphere initially in hydrostatic equilibrium and without any winds.Simulations are then run for 600 days without clouds in order to spin-up the atmosphere at which point clouds are enabled.Simulations are then run for an additional 400 days with clouds enabled, except for the 100× solar metallicity  sed = 0.5 and 0.1 cases.These two cases are run for 800 days with clouds enabled due to the greater impact of clouds on the atmosphere.Clear sky simulations are run for 1000 days 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  sed for the E S cloud models:  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 days after clouds are enabled (i.e., between days 600 and 700).

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.

Thermal Structure and The Distribution of Clouds
The equatorial temperature profiles for the solar and 100× solar metallicities are shown in Figures 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 (Figure 1), the introduction of cloud results in only a small (∼ 50 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 bar and KCl forming at 0.1 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 (Figure 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% of the incoming stellar radiation is scattered back through the top of the atmosphere, the cloudy simulations show 7.8%, 15.9% and 44.9% of incoming stellar radiation scattered back for  sed = 1.0, 0.5 and 0.1, respectively.In the  sed = 1.0 case where clouds are relatively shallow, this results in a shift of ∼ 100 K near the ZnS cloud deck at 1 bar.As the clouds increase in vertical extent (i.e., with decreasing  sed ), the impact of the reduced heating extends throughout the atmosphere at pressures below 1 bar, with up to a ∼ 300 K shift around  ∼ 0.1 bar for the  sed = 0.1 case.This shift in temperature results in the base of the clouds appearing at higher pressures (see Figure 3).
Above the cloud deck, the vertical distribution of cloud exhibits the pressure dependence  c ∝   sed characteristic of the E S 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 modeling 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, for example, Figure 2 in Charnay et al. 2015b).
The effective particle sizes within the clouds (see Eqn. 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 E S 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 −2 bar with sizes on the order of 1 to 10 microns.As the particle sizes in E S 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.

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 forc-ing (Mayne et al. 2019), dual band grey radiative transfer (Menou 2012;Wang & Wordsworth 2020), and multi-band correlated-k radiative transfer (Kataria et al. 2014;Charnay et al. 2015a,b) at least in part explains the differences 2 .In the simulations presented here, we find an equatorial jet forms in all cases except the 100× solar metallicity,  sed = 0.1 case (see Figure 5).For the solar metallicity cases, we find two additional polar jets with amplitudes ∼ 1.2 to 1.3 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 ∼ 1.6 km s −1 .
For atmospheres with 100× solar metallicity, we find superrotation in the mid-latitude and polar regions in the clear sky,  sed = 1.0, and  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 to 1 bar as angular momentum is carried upwards into the jet (see Figure 5 and Showman & Polvani 2011).
The speed of the equatorial jet does appear to increase with cloud scale height; however, in the  sed = 0.1 case, the superrotating jet eventually decays, forming two dayside gyres and a counter-rotating wind (see Figure 6 as well as additional plots in Appendix E).The mechanism causing the decay of the jet is unclear, although it does 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.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  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  sed = 0.1 case shown in Figure 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 6, right panel.We leave a more detailed investigation of the underlying mechanism and whether it is numerical in origin for future work.

Observational Diagnostics
To allow for comparison with observations, we have generated synthetic observations using the S radiative transfer routines within the UM with the spectral resolution increased to 10 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. (2018aLines et al. ( , 2019)); Christie et al. (2021) and the specifics of the transmission spectrum calculation are presented in Lines et al. (2018a).

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 panel of Figure 7 shows the synthetic observations for the HST WFC3 bandpass covering 1.1 to 1.7 m.The observations of Kreidberg et al. (2014) are included, with the data shifted up by 0.002 to facilitate comparison. 3A shallow cloud layer, as in the  sed = 1.0 and  sed = 0.5 cases, has only a minor impact on the transmission spectrum.The  sed = 0.1 case shows the best agreement with the data, insofar as the spectrum is relatively flat, however, the 1.4 micron CH 4 feature is still visible.The sharp decrease in the transmission spectrum of the  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 et al. 1987).
To better understand how the specific cloud species are impacting the transmission spectrum in the  sed = 0.1 case, we breakdown the spectrum further.Figure 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.
We additionally look at the transmission spectra in the JWST MIRI/LRS bandpass (see Figure 7, right 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  sed = 0.1 case.The cooler atmosphere results in a smaller transit across all wavelengths in the bandpass (see Figure 8, right panel).In the region between 6 m and 9 m, the spectrum is dominated by CH 4 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, for example, the NH 3 feature between 10 m and 11 m in Figure 8, right panel).The limited impact of ZnS in this bandpass may be in part to the limitations in the Querry et al. (1987) refraction index data, discussed above.

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.Figure 9 (right 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  sed = 1.0 case, there is minimal difference between the  clear sky and  sed = 1.0 phase curves.The simulation with the most vertically extended clouds -the  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  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  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 hot spot westward.We also note that, as expected for the given wavelength range, the phase curve is probing thermal emission with only ∼ 10 −9 % of the emission coming from reflected starlight in the clear sky case and increasing to ∼ 2% in the  sed = 0.1 case.
As the effective temperature of GJ 1214 is 3026 K (Charbonneau et al. 2009), the stellar blackbody peaks at ∼ 1 m allowing for the possibility of scattered starlight by clouds contributing significantly to the HST WFC3 phase curve (see Figure 9, left 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  p / s ∼ 2.3 to 2.7 × 10 −7 .With increased cloud extent, however, the scattered stellar component begins to dominate the phase curve, and for  sed = 0.1, the phase curve peaks with  p / s = 2.6×10 −5 , with 99.8% of the planetary flux coming from the scattered stellar component.

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 E S cloud model to the UM GCM.Consistent with previous investigations, we find that increased metallicity and increased cloud vertical extent (i.e., decreased  sed ) are necessary to impact the pressure temperature profiles and synthetic observables.The most significant impact was for our 100× solar metallicity case with  sed = 0.1 where ∼ 45% 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 Figure 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  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 a 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 don't fit naturally in the E S framework.

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 , corrected using the zenith angle , yielding an effective attenuating column /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 days using the clear-sky solar-metallicity configuration outlined in the paper.Qualitatively, we see the same velocity structures forming in both cases (see Figure A1), with a fast equatorial jet and two polar jets with the equatorial jet having a maximum speed of ∼ 1.5 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 km s −1 in the plane-parallel case versus 1.0 km s −1 in the spherical case).
This can be seen in Figures A2 and A3 where the plane-parallel case shows a larger decrease in the temperature approaching the pole.
We estimate the energy deposited in the atmosphere based on the heating rates and find that the spherical irradiation scheme results in ∼ 6% 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 multi-modal log-normal (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 et al. 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 E S framework.While the assumption of a log-normal particle size distribution is often the simplest from a computational and analytical standpoint, E S can be formulated using a gamma distribution4 which is presented below.

B1 Derivation
Considering a layer of cloud within the atmosphere, we take the number distribution of cloud particles to be / =   (; , ) where N is the total number of particles and  (; , ) is the gamma distribution, 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 Var [log ] = log 2  where  parametrises the width of the distribution.For a gamma distribution, the variance of a gamma-distributed variable  depends only on the parameter , where  (1) () is the trigamma function.Following the method of Ackerman & Marley (2001), we can then fix the value of  based on a prescribed ,  As  is taken to be constant throughout the atmosphere in the Emodel,  only needs to be calculated at the beginning of the simulation, limiting computational overhead.For the default value of  = 2 used here,  ≈ 2.54278.With  fixed, we can now calculate  based on the sedimentation arguments of Ackerman & Marley (2001).The sedimentation factor is defined as This can be rewritten as where  [] is the expected value of .We can then solve for , In the rare cases where  is an integer, the identity Γ( + 1) = Γ() can be used to eliminate the gamma function entirely and reduce the equation to 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 resulting in The area-weighted effective radius  eff in a cloud layer is given by The scaling of  eff with the size parameter  is very different in the case of a gamma distribution compared to the log-normal distribution.For small values of , the two distributions scale similarly with the distributions becoming sharply peaked with  eff ∼    1/ sed ; however, for wider distributions (e.g., for  greater than 2, see Figure B1),  eff decreases slower with  in the case of the gamma distribution case compared to the log-normal distribution case.For wider distributions, we may, as a result, expect larger differences in particle sizes and optical properties between cloud layers with a log-normal distribution than with cloud layers with a gamma distribution.Coincidentally, we find that for reasonable values of  the effective radii  eff for the two distributions are similar for  ∼ 2, the value used here and widely in the literature.
To illustrate how using E S with a gamma distribution differs from the default log-normal distribution case, the log-normal and gamma distributions are plotted in Figure B2 for  sed = 1.0 and  = 2 and three values of .Compared to the log-normal distribution case, both location of the 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   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 log-normal distribution where a wider distribution results in more massive particles being traded for exponentially more particles at the peak of the distribution.

B2 Results
To examine the impact of the change in distribution, we focus on the 100× solar metallicity  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 days.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 midlatitudes we find cooler temperatures on the night side and along the morning terminator in the case of a gamma distribution (see Figures B3 and B4).The differences occur around  ∼ 10 −3 bar, with the log-normal distribution case being ∼ 100 K hotter than the gamma distribution case.We note that the flow at mid-latitudes is somewhat different between the log-normal 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.
Switching to a gamma distribution may result in larger differences from the log-normal 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 E S assumption of a fixed distribution width may allow this issue to be avoided, but also it may provide an avenue for moving beyond E S 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 A 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 Figure 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.
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 (Figure 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.
For the solar metallicity case, the introduction of cloud at 600 days results in a shift in the zonal velocities of ∼ 0.2 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 ∼ 0.05% over the 1000 days of simulation (see Figure D2, left panel).The 100× solar metallicity simulations, on the other hand, show a larger loss of angular momentum after clouds are enabled, with the  sed = 0.1 case losing ∼ 1.4% of the initial angular momentum and the majority of that loss occurring during the final 800 days (see Figure D2, right 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.
The transition between the flow structures in the 100× solar metallicity,  sed = 0.1 case can be seen in Figure D1 (right panel), with the transition in the flow strucutre being illustrated in Figure E2.

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  sed = 0.1 simulation.Figure E1 shows the evolution of the flow over the length of the simulation beginning at when clouds are enabled.
Figure E2 shows evolution over the period between 1220 days and 1280 days when the flow undergoes the largest change in structure.E1, except for the 60 days around the transition in the flow structure.

Figure 1 .Figure 2 .
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 3 .Figure 4 .
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.

Figure 5 .Figure 6 .
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.

Figure 8 .
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  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.

Figure 9 .
Figure9.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 m) and the JWST MIRI/LRS bandpass (right, with a wavelength range between 5 microns and 12 m).

Figure A1 .Figure A2 .Figure A3 .
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 days.The black dashed line indicates the contour of zero zonal velocity.

Figure B1 .Figure B2 .
Figure B1.The effective particle radii  eff for E S formulated with a log-normal distribution (solid lines) and a gamma distribution (dashed lines).For the purposes of this plot,  = 2 has been assumed.

Figure B3 .
Figure B3.The temperature profiles at a latitude of 45 • for the 100× metallicity,  sed = 0.1 case for a log-normal 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 D1 .Figure D2 .Figure E1 .
Figure D1.The evolution of the zonal velocity in each of the simulations.For each cloudy simulation, the clouds were enabled after 600 days.The transition in the flow in the 100× solar metallicity  sed = 0.1 case can be clearly seen in the right plot.

Table 1 .
Common Parameters 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.
The same as in Figure