Formation and dynamics of water clouds on temperate sub-Neptunes: the example of K2-18b

Hubble (HST) spectroscopic transit observations of the temperate sub-Neptune K2-18b were interpreted as the presence of water vapour with potential water clouds. 1D modelling studies also predict the formation of water clouds at some conditions. However, such models cannot predict the cloud cover, driven by atmospheric dynamics and thermal contrasts, and thus their real impact on spectra. The main goal of this study is to understand the formation, distribution and observational consequences of water clouds on K2-18b and other temperate sub-Neptunes. We simulated the atmospheric dynamics, water cloud formation and spectra of K2-18b for H2-dominated atmosphere using a 3D GCM. We analysed the impact of atmospheric composition (with metallicity from 1*solar to 1000*solar), concentration of cloud condensation nuclei and planetary rotation rate. Assuming that K2-18b has a synchronous rotation, we show that the atmospheric circulation in the upper atmosphere essentially corresponds to a symmetric day-to-night circulation. This regime preferentially leads to cloud formation at the substellar point or at the terminator. Clouds form for metallicity>100*solar with relatively large particles. For 100-300*solar metallicity, the cloud fraction at the terminators is small with a limited impact on transit spectra. For 1000*solar metallicity, very thick clouds form at the terminator. The cloud distribution appears very sensitive to the concentration of CCN and to the planetary rotation rate. Fitting HST transit data with our simulated spectra suggests a metallicity of ~100-300*solar. In addition, we found that the cloud fraction at the terminator can be highly variable, leading to a potential variability in transit spectra. This effect could be common on cloudy exoplanets and could be detectable with multiple transit observations.


Introduction
Detection surveys revealed a high abundance of exoplanets with masses intermediate between the Earth and Neptune, called super-Earths and sub-Neptunes. Statistical analysis of Kepler data suggests a transition between these two populations at ∼1.8 R Earth (Fulton et al. 2017), compatible with models of photoevaporation of H 2 -dominated atmospheres surrounding rocky cores (Owen & Wu 2017;Lehmer & Catling 2017) or formation model of water world Zeng et al. (2019). Following the observed trend in giant planets of the Solar System and prediction of planetary formation model, one would expect that the fraction of heavy elements (called the metallicity) in primary atmospheres decreases with the planetary mass (Kreidberg et al. 2014a;Fortney et al. 2013). Sub-Neptunes are thus expected to be enriched in heavy elements reaching typically 100-1000× solar metallicity. Measuring the atmospheric composition in particular the water abundance of sub-Neptunes would provide major constraints on planetary formation and evolution. Unfortunately, past transit spectroscopic observations of warm sub-Neptunes have been unsuccessful for measuring molecular abundances because of the presence of high and thick clouds or hazes (Kreidberg et al. 2014b;Knutson et al. 2014;Benneke et al. 2019a).
K2-18b is a moderately irradiated (1.06±0.06 × that of the insolation on Earth) sub-Neptune planet Montet et al. 2015;Benneke et al. 2017). It is on an orbit slightly inside the inner limit of the classical Habitable Zone (Kopparapu et al. 2014). It is however probably in the Habitable Zone for slow synchronous rotating terrestrial planets around Mstars 1 (Yang et al. 2013). Its measured mass and radius (8.6±1.4 M Earth and 2.71±0.07 R Earth from Cloutier et al. 2019) are quite similar to GJ 1214b, and suggest it is surrounded by a low molec-ular weight (H 2 /He-dominated) atmosphere with potentially a fraction of water vapour. Recent analysis of HST/WFC3 transit observations revealed an absorption at 1.4 µm, interpreted as water vapour absorption (Tsiaras et al. 2019;Benneke et al. 2019b) or methane absorption (Bézard et al. 2020) in a low molecular weight (H 2 /He-dominated) atmosphere. Benneke et al. (2019b) also suggested that water clouds (most likely icy, but possibly liquid for a very narrow range of parameters) may form in the atmosphere of K2-18b. K2-18 b appears as a prime target to characterize the atmospheric composition of super-earth/mini-Neptune and the climate dynamics of moderately irradiated worlds. Previous atmospheric modelling studies of K2-18b are based on 1-dimensional models (Benneke et al. 2019b;Madhusudhan et al. 2020;Scheucher et al. 2020;Bézard et al. 2020;Blain et al. 2020). Such models cannot predict the cloud cover, driven by atmospheric dynamics and thermal contrasts, and thus the real impact of clouds on spectra and on the planetary albedo. Previous 3D climate studies for terrestrial planets showed how the atmospheric circulation and cloud radiative effects strongly affect the cloud distribution and the habitability (Leconte et al. 2013b,a;Yang et al. 2013;Wolf & Toon 2014). In particular, 3D models predict a strong cloud formation on the dayside of slow synchronous rotating terrestrial planets around M-stars. For such conditions, the associated high planetary albedo could place the inner edge of the Habitable Zone closer to the host star (Yang et al. 2013). This intense cloud formation could also limit the detectability of water vapour and other species from transmission spectroscopy (Fauchez et al. 2019;Komacek et al. 2020;Suissa et al. 2020). These studies illustrate the importance of 3D modelling to assess the effect of clouds on observational spectra.
Here, we use a three-dimensional General Circulation Model (GCM) to simulate the atmosphere of K2-18b for H 2 -dominated atmosphere with water clouds. In Section 2, we describe the 3D model and the conditions used in this study. In Section 3, we present our results on the thermal structure, the atmospheric dynamics, the cloud distribution and variability. We analyse the impact of atmospheric metallicity (from 1×solar to 1000×solar), concentration of cloud condensation nuclei and rotation rate. We discuss the impact of clouds on transit spectra in Section 4. We finish with a summary and conclusions in Section 5.

The LMD Generic GCM
We simulated the atmosphere of K2-18b using the LMD Generic GCM. This model has been specifically developed for exoplanet and paleoclimate studies. In particular, it has been used for studying the atmospheres of moderately irradiated terrestrial exoplanets (Wordsworth et al. 2011;Leconte et al. 2013a,b;Turbet et al. 2016Turbet et al. , 2018 and warm sub-Neptunes (Charnay et al. 2015a,b). The model is derived from the LMDZ Earth (Hourdin et al. 2006) and Mars (Forget et al. 1999) GCMs. It solves the primitive hydrostatic equations of meteorology using a finite difference dynamical core on an Arakawa C grid. The model includes schemes for adiabatic convective adjustment, large-scale water cloud condensation, moist convection, precipitation and evaporation (see the description of the water cycle in Charnay et al. (2013) and in the next subsection). Mayne et al. (2019) showed that the full Navier Stokes equations can produce some differences compared to the primitive hydrostatic equations for simulations of the warm sub-Neptune GJ1214b. These differences are significant for cases with strong day-night thermal contrasts. However, the thermal contrast should be very small on K2-18b (see section 3.1.) and its atmospheric scale height is approximately half that of GJ1214b for a similar planetary radius. We expect the shallow fluid approximation to remain fully acceptable under K2-18b's regime.
In this paper, simulations were performed with a horizontal resolution of 64×48 (corresponding to resolutions of 3.75 • latitude by 5.625 • longitude). We also did a test with a 128×96 resolution but we did not notice significant differences. For the vertical discretization, the model uses pressure coordinates. In this work, we used 40 layers equally spaced in log pressure, with first level at 80 bars and top level at 0.2 mbar.
The radiative scheme is based on the correlated-k method with the absorption data calculated directly from high resolution spectra. We used the k-coefficients computed by Blain et al. (2020) with bins of 200 cm −1 (spectral resolution R∼50 at 1 µm). We included absorption by H 2 O, CH 4 , CO, CO 2 , NH 3 , PH 3 , H 2 S, HCN, K, Na, FeH, TiO and VO. We used HITEMP 2010 for H 2 O, CO and CO 2 (Rothman et al. 2010), TheoReTS for CH 4 (Rey et al. 2018), and ExoMol for NH 3 and PH 3 (Coles et al. 2019;Yurchenko 2015;Sousa-Silva et al. 2015).
The H 2 -H 2 and H 2 -He collision-induced absorptions (CIA) from the HITRAN database (Karman et al. 2019) were included, as well as water vapour continuum from Clough et al. (1989). Rayleigh scattering by H 2 , H 2 O and He was included, based on the method described in Hansen & Travis (1974). Radiative transfer is computed with the Toon et al. (1989) scheme including water cloud optical properties.

Water clouds
The model simulates water cloud formation including condensation, evaporation, coalescence and sedimentation. We used the routines for the water cycle described in Charnay et al. (2013). The mixing ratio of water vapour at the first pressure level is fixed (see values in Table 1). Water vapour is advected by the atmospheric circulation and condenses where its partial pressure exceeds the water vapour saturation pressure. We fixed the density number of cloud condensation nuclei (CCN), exploring values from 10 4 to 10 7 CCN/kg. We based this range on cloud particle concentrations for terrestrial water clouds, with typically 10 7 − 10 8 droplets/kg for marine cumulus and 10 4 − 10 5 ice particles/kg for cirrus (Wallace & Hobbs 2006). The CCN concentration is a key parameter that controls the properties of the water clouds. Exploring the sensitivity of the results to this poorly known parameter allows us to account for most of the uncertainties related to the water clouds microphysics and particles size distribution. For reference simulations, we chose a concentration of 10 5 CCN/kg corresponding to conditions for cirrus formation. For a sub-Neptune like K2-18b, CCN sources for water clouds could be micrometeorites, photochemical hazes or NH 4 Cl cloud particles. The latter are salts and could be efficient CCN. NH 4 Cl clouds should form at ∼0.1 bar for K2-18b's conditions (Blain et al. 2020).
In all simulations, water clouds form for water vapour pressures always lower than 6.11 mbar (i.e. below the triple point of water). They are thus expected to be composed of ice particles only. We assumed spherical particles or crystal shapes of rimed dendrites (geometric parameters from Heymsfield (1977)). For K2-18b's conditions, the air flow around falling cloud particles can have a high Reynolds number. We modified our sedimenta-10 0 10 1 10 2 10 3 10 4 Particle radius (melted water, in m) Sedimentation velocity (m/s) Sedimentation velocity for ice particles on K2-18b Spherical particles (VStokes) Spherical particles (Vdrag) Non-spherical particles (Vdrag for rimed dendrites) Fig. 1. Sedimentation velocity as a function of particle radius (equivalent radius of melted water) computed for K2-18b at 10 3 Pa. The blue line follows the Stokes law for spherical particles. The yellow line shows the terminal velocity for the general drag coefficient. The red line is computed with the general drag coefficient and for rimed dentrites (parametrization from Heymsfield (1977)).
tion routine (based on the Stokes law) to take into account these regimes. The terminal velocity of cloud particles is given by: where r c is the particle radius, ρ atm is the atmosphere density, ρ particle is the cloud particle density and C D is the drag coefficient, expressed as (Clift & Gauvin 1971): Re 1 + 0.15Re 0.687 + 0.42 1 + 42500 is the Reynolds number and η is the dynamic viscosity. Fig.1 shows the deviations between the Stokes law and the general terminal velocity for spherical particles (blue and yellow lines). Differences are significant for particles larger than 200 µm. For rimed dendrites, the drag is increased compared to spherical particles because of a higher drag surface for the same mass. The terminal velocity is reduced by a factor of ∼4 for rimed dendrites compared to spherical particles with a 100 µm radius.
The model takes into account the conversion of cloud ice particles to snowflakes by coalescence with other particles following the parameterization from Boucher et al. (1995). For K2-18b, coalescence generally dominates for water ice precipitation. Condensation and evaporation of water clouds also include latent heat release, which can potentially trigger moist convection. However, moist convection rarely occurs in our simulations of K2-18b because water condenses at high altitude, where the atmosphere is almost isothermal. In addition, moist convection could be inhibited for sub-Neptunes with high fractions of water vapour due to the vertical gradient of mean molecular weight appearing with water condensation (Leconte et al. 2017). Table 1 indicates the critical water vapour mixing ratios above which moist condensation should be inhibited. The water vapour mixing ratios exceed these limits for all cases where water condenses in K2-18b's atmosphere (i.e. 100×, 300× and Composition r vap (mol/mol) q vap (kg/kg) q crit (kg/kg) 100×solar 0.057 0.22 0.06 300×solar 0.088 0.2 0.08 1000×solar 0.14 0.2 0.144 Table 1. Water vapour volume mixing ratio (r vap ), mass mixing ratio (q vap ) and maximal mass mixing ratio for moist convection (q crit ) for the different atmospheric compositions.

Model parameters and configurations
where Q is the planet's tidal dissipation quality factor, ω p is the planet's primordial rotation rate (ω p ∼ 1.7 × 10 −4 for Jupiter) and D is the semi-major axis. From formula (3), τ ∼ 17 Ma for Q = 100 (typical value for terrestrial planets) and τ ∼ 17 Ga for Q = 10 5 (typical value for gas giant planets). We can expect that sub-Neptunes have intermediate values for Q, likely lower than 10 4 , giving τ 1.7 Ga. With an estimated age of 2.4 ±0.6 Ga (Guinan & Engle 2019), K2-18b would likely be tidally locked. Note that with its potential non-zero eccentricity (e ∼ 0.09 in Cloutier et al. (2019)), K2-18b is a good candidate for a spinorbit resonance as Mercury. For most of our simulations, we assumed that the planet is tidally-locked with a synchronous rotation around its host star. We also performed tests with different rotation rates (see Section 3.4). The simulations were performed for H 2 -rich atmospheres with 1×, 10×, 100×, 300× and 1000×solar metallicity. The atmospheric composition is based on calculations with the 1D model Exo-REM, which has already been applied to K2-18b (Bézard et al. 2020;Blain et al. 2020). The 1D model was run with non-equilibrium chemistry and an eddy mixing coefficient K zz =10 6 cm 2 /s. Figure A.1 shows the atmospheric composition profiles used for the different cases. We assumed no longitudinal variation of the atmospheric composition except for water. Table  2.3 shows the specific heat capacity (c p ), the atmospheric scale height at the 300 K level (H) and the mean molecular weight for the different compositions. All simulations were initialised with a 1D thermal profile (computed with the 1D version of the model) and were run for more than one thousand K2-18b orbits (∼100 Earth years). Simulations with 300×solar metallicity provide the best fits to transit observations (discussed in Section 4). We used this atmospheric composition as reference in this study to explore sensitivity to cloud particle size and rotation rate. Fig. 2 shows temperature profiles at the pole, at the substellar point and at the equatorial morning terminator for different atmospheric metallicities. We notice that the isothermal region appears at a lower pressure for the high metallicity cases (e.g. ∼0.01 bar for 1000×solar compared to ∼0.1 bar for 1×solar). This is due to the enhanced infrared opacities and greenhouse effect for high metallicity cases. In addition, the thermal gradient and the temperature at 1 bar are reduced for the 1000×solar case. This is due to the specific heat capacity per molecule which increases with the metallicity and the abundances of H 2 O, CO 2 , CH 4 , NH 3 , which have more degrees of freedom (i.e. vibration modes) than H 2 and He.

Thermal structure and atmospheric dynamics
There is little longitudinal/latitudinal temperature variations, apart from the upper atmosphere (at pressures lower than 10 mbar) for high metallicity (≥100×solar). CH 4 is the main absorber of stellar flux in the upper atmosphere. It produces radiative heating and a stratospheric thermal inversion at the substellar point for high metallicity cases. The weak horizontal temperature variations are due to the long radiative timescale compared to the advection timescale. 1D modelling is therefore an excellent approach for computing the thermal structure of such a temperate planetary atmosphere.
Pressure-latitude cross-sections of zonally-averaged zonal wind and mean equatorial zonal winds (Fig. 3) show the presence of a weak tropospheric equatorial superrotation jet (i.e. westerly winds) at pressures between 0.1 and 10 bars. Showman & Polvani (2011) showed that an equatorial superrotation jet develops on tidally-locked planets by the formation of standing planetary-scale equatorial and Rossby waves. A condition for this mechanism to occur is that the equatorial Rossby deformation radius be smaller than the planetary radius. We can define the dimensionless equatorial Rossby deformation length (Leconte et al. 2013b): where H is the atmospheric scale height and N = g T g c p + dT dz is the Brunt-Vaisala frequency. In the troposphere, the vertical temperature gradient deviates from the dry adiabatic gradient (Γ dry = − g c p ) by around 1-10% and L ∼ 0.7. An equatorial superrotation jet forced by stationary planetary waves can develop there. In contrast, the stratosphere is almost isothermal with L ∼ 2. An equatorial jet cannot develop from stationary planetary waves at pressures lower than than ∼0.1 bar. We note that for the case with 10×solar metallicity, superrotation develops at pressures lower than 0.1 bar and at all latitudes. Superrotation is likely triggered by a different process here (e.g. by barotropic waves as on Titan). In any case, the atmospheric circulation for pressures lower than 0.1 bar is dominated by a day-night circulation, with upwelling air on the dayside and downwelling air on the nightside (see Fig. 4a for illustration). Winds are relatively axisymmetric around the substellar-anti-stellar axis. For 300×solar metallicity, horizontal winds are maximum at the terminator and reach 200 m/s at top model level (at 0.2 mbar). The vertical wind is around 0.2 m/s at the substellar point for pressures lower than 0.1 bar.
The dayside is heated by stellar radiation and cooled by ascending air producing adiabatic cooling. The temperature of the nightside is controlled by radiative cooling and by adiabatic warming produced by downwelling air. Counter-intuitively and because of this nightside adiabatic warming, the coldest point in the atmosphere between 0.5 and 5 mbar is not at the anti-substellar point but at the terminator (see Fig. 5). For the 300×solar metallicity, the temperature at longitude ±90 • and at 1 mbar is ∼20 K lower than at the substellar point (see Fig. 5).

Cloud formation
As illustrated in Fig. 2, water clouds can only form on K2-18b for high atmospheric metallicity (100-1000×solar) if the atmosphere is H 2 -dominated with a solar C/O. In this case, water condensation occurs between 2 and 10 mbar, where the atmosphere is almost isothermal. Cloud formation is ruled by the day-night atmospheric circulation described in the previous subsection. For such low pressures, clouds should be composed of water ice. They preferentially form at the substellar point at the tropopause by the vertical advection of water vapour, or at the terminator due to the strong radiative cooling occurring when the air moves towards the nightside (see Fig. 4b). Cloud precipitations evaporate below the cloud layer in undersaturated air (typically in one atmospheric scale height), as virgae on Earth. This increases the water vapour mixing ratio below the cloud layer as illustrated in the bottom left panel of Fig. 5. Fig. 6 shows the equatorial cloud mass mixing ratio and the vertical integrated mass of condensed water for 100×, 300× and 1000×solar metallicity. For 100× and 300×solar, clouds essentially form at the substellar point. They appear at similar pressures at other locations but intermittently and with much smaller mixing ratio. For 1000×solar metallicity, the thermal contrast between the dayside and the nightside is enhanced. This is due to the shorter atmospheric radiative timescale for a higher metallicity, which implies a higher mean molecular weight, a lower specific heat capacity per mass and a photosphere at a lower pressure (see Menou (2012)). The dayside is warmer and less favourable to cloud formation. The shorter radiative timescale makes the terminator colder and very favourable to cloud formation. Optically thick clouds form there circling the whole planet. Despite the strong condensation, the abundance of water vapour above the cloud layer at the terminator is almost unchanged because it is advected by horizontal winds. The adiabatic warming caused by subsiding air combined with the humidity reduced by precipitation at the terminator limits cloud formation on the nightside.
We explored the sensitivity to cloud particle size by changing the CCN concentration for the 300×solar metallicity. Fig. 7 shows the equatorial cloud mass mixing ratio for spherical cloud particles with CCN concentrations of 10 5 , 10 6 and 10 7 nuclei/kg of air. The cloud distribution is significantly changed with a regime of preferential cloud formation at the substellar point for low CCN concentrations (10 4 -10 5 /kg of air) and a regime of preferential cloud formation at the terminator and on the nightside for high CCN concentrations (10 6 -10 7 /kg of air). By fixing the CCN concentration n CCN , the mean cloud particle size is , where q cloud is the mass mixing ratio of condensed water. The particle size is typically r c = 450 µm, 200 µm, 60 µm and 30 µm for n CCN = 10 4 , 10 5 , 10 6 and 10 7 kg −1 respec-Atmospheric composition c p (J kg −1 K −1 ) H(km) Mean molecular weight (  tively. For n CCN = 10 4 and 10 5 kg −1 , the sedimentation velocity of spherical particles is ∼10-30 m/s at 10 mbar (see Fig. 1), much higher than the vertical wind speed (w ∼ 0.2 m/s). The sedimentation timescale (τ sed = H/v sed ) is also much shorter than the advection timescale (τ adv = R p /v). Clouds are not efficiently transported and are limited to their initial formation location (essentially the substellar point). For n CCN = 10 6 and 10 7 , the sedimentation velocity of spherical is ∼0.3-1 m/s and the sedimentation timescale is relatively close to the advection timescale (τ sed /τ adv = 0.1 − 0.3). This ratio is enhanced by upward vertical winds on the dayside and reach values close to 1 for n CCN = 10 7 kg −1 taking into account vertical winds. Clouds are then quite efficiently advected horizontally to the terminator, reducing the efficiency of the dayside cold trapping. This enhanced global water transport favours cloud formation at terminator and on the nightside. Surprisingly, cloud formation is reduced at the dayside. This is due to a cloud radiative feedback. Clouds formed at terminator or on the nightside warm the atmosphere below the cloud deck by greenhouse effect. The greenhouse effect is illustrated in Fig. 6 with a reduction of the outgoing longwave radiation where clouds are present. Because of the efficient heat redistribution, the dayside becomes warmer and potentially too warm for cloud formation, suppressing the associated cold trapping and the horizontal/vertical gradient of humidity. This reinforces cloud formation at terminator and on the nightside. In contrast, we found that the cloud albedo effect is limited. For K2-18's spectrum (M2.8 stellar spectral type) and for K2-18b's CH 4 -and H 2 O-rich atmosphere, a large part of the stellar radi-Article number, page 5 of 16 A&A proofs: manuscript no. article_K2-18b_Charnay_arxiv ation is absorbed by CH 4 and H 2 O above clouds. The planetary albedo is always lower than 0.1 and the cloud distribution only marginally changes this value.
We also performed a simulation with n CCN = 10 5 and with non-spherical ice particles (rimed dendrites, see Section 2). Sedimentation speed is reduced by a factor of ∼4 and the results are intermediate between the case with n CCN = 10 5 and the case with n CCN = 10 6 . The morning (west) terminator appears more cloudy than the evening (east) terminator because it is slightly colder.
Exploring a large range of CCN concentration is a simple way to cover the possible impact of microphysics and particle shape.
From these simulations, we conclude that the cloud distribution on K2-18b is controlled by the global dayside-to-nightside circulation, particle size and cloud radiative effects. Radiative feedbacks can profoundly alter cloud distribution, making it very sensitive to cloud microphysics. In any case, the terminator is at least partially cloudy for metallicity 100×solar with clouds confined between 2 and 10 mbar.

Ascending region
Stellar heating Adiabatic cooling

Dayside-nightside circulation for tidally-locked planet
Radiative cooling Radiative cooling b) Fig. 4. Illustration of the day-night atmospheric circulation around K2-18b for a synchronous rotation. Panel a shows the circulation and warming/cooling zones. Panel b shows the preferential location of cloud formation.

Cloud variability
Cloud infrared opacity produces a significant greenhouse effect warming the atmosphere below the cloud layer. As discussed in the previous paragraph, this cloud radiative effect has a strong impact on the spatial cloud distribution. We also found that it produces time variability in the cloud fraction. We computed the global cloud fraction and the cloud fraction at the terminator for the different atmospheric compositions assuming that the atmosphere is cloudy when the cloud optical depth in the visible range (τ cloud ) is greater than 1. The cloud optical depth in transit (τ cloud H ) is related to the normal cloud optical depth (τ cloud V ) by (Fortney 2005): η∼65 for 300×solar metallicity. In the optical regime (r λ), the particle extinction cross-section is equal to 2πr 2 and the normal cloud optical depth can be expressed as: where m cloud is the vertical integrated mass of condensed water. Clouds are optically thick (τ cloud H 1) for m cloud 2r c ρ ice 3η . For r c = 200 µm, this corresponds to m cloud 1.9×10 −3 kg/m 2 . From   Fig. 6, the vertical integrated mass of condensed water is generally much larger than 0.01 when clouds form, meaning that water clouds on K2-18b are optically thick at the limbs. Note that what matters is the cloud optical depth in the visible part of the atmosphere. Our simplified calculation is valid because clouds form high enough in the atmosphere with little latitudinal variations of the vertical cloud extend at the terminator. Using the time-averaged GCM outputs, clouds would be optically thick at the terminator at all latitudes for all cases. Note that atmospheric columns close to the terminator can have a significant contribution to the limb optical depth. We computed the mean optical depth over an opening angle of 15 to 30 • longitude, depending on the atmospheric scale height and based on the work by Caldas et al. (2019) (see Fig. 2 in their paper). If we use the instantaneous GCM outputs, the mean cloud fraction is around 99% for 1000×solar and around 13% for 100×solar and 300×solar metallicity (see Fig. 8). For the last two cases, the cloud fraction varies between 0% and 50%. This difference between the instantaneous and time-averaged atmospheric state can have important implications for observations as discussed in Section 4. Figure 9 is an illustration of the atmospheric transmission at 0.4 micron for the 300×solar metallicity with 10 6 nuclei/kg (relatively cloudy case), computed with Pytmosph3R (Caldas et al. 2019). It shows the altitude were the atmosphere becomes opaque. One limb (east limb, on the right in the plot) is more opaque (difference of 50 km) because of the presence of clouds.

Asynchronous rotation
As discussed in Section 2.3, K2-18b may have a nonsynchronous rotation, since it is close to the limit for tidallocking. We tested the impact of rotation rate on the atmospheric dynamics and cloud distribution for 2:1, 4:1 and 10:1 spin-orbit resonance. Fig. 10 shows zonally-averaged zonal wind and maps of the vertical integrated mass of condensed water. For the 2:1 resonance, two zonal jets appear at high latitudes. The cloud distribution is mostly unchanged, with a preferential formation at the substellar point but the amount of cloud at the substellar point is reduced by a factor of ∼2, due to a weaker day-night circulation. The cloud distribution for a 3:2 spin-orbit resonance would therefore be similar to our reference 1:1 case, with a reduction of the amount of cloud at the substellar point. For a 4 times faster rotation, the high latitude jet are reinforced. The cloud distribution is significantly changed with preferential formation at the terminator. Cloud maximum is shifted eastward due to the easterlies at pressures lower than 10 mbar, where clouds form. For a 10 times faster rotation, the superrotation is well developed at mid and high latitudes, and easterlies appear at the equator in the troposphere. The circulation is closer to an Earth-like circulation. Again, clouds preferentially form at the terminators but their maximum is shifted eastward due to the superrotation, which is developed everywhere in the stratosphere.

Observational transit spectra
Using the outputs of the 3D GCM, we computed transit spectra using Exo-REM (Baudino et al. 2015(Baudino et al. , 2017Charnay et al. 2018;Blain et al. 2020) and we compared them to HST data from Tsiaras et al. (2019) and Benneke et al. (2019b) in order to constrain the atmospheric metallicity and to highlight implications for future observations. We computed spectra with k-coefficient bands of 20 cm −1 (R∼500 at 1 µm) using the temperature and cloud profiles from the GCM at the terminators. Fig. 11 shows the transit spectra with no cloud for 10×, 100× and 1000×solar Article number, page 7 of 16 A&A proofs: manuscript no. article_K2-18b_Charnay_arxiv Mean equatorial (averaged between latitude -30 • and +30 • ) zonal winds, vertical wind, temperature, radiative warming/cooling, water vapour mixing ratio and water cloud mixing ratio as a function of longitude and pressure for 300×solar metallicity. metallicity compared to HST/WFC3 data from Tsiaras et al. (2019). As shown by Bézard et al. (2020) and also discussed in Blain et al. (2020), the transit spectrum in HST/WFC3 band and in particular the 1.4 µm band is dominated by CH 4 absorption for a H 2 -dominated atmosphere with solar C/O. Spectral features are too deep for the 1×solar and 10×solar metallicity compared to observations. These two cases are ruled out if they are cloudfree as predicted by 3D model.
We computed the transit spectrum with partial cloud cover as the sum of a cloudy spectrum and a clear spectrum pondered by the cloud fraction. This approximation is acceptable because the temperature profile and the altitude of clouds do not significantly change with latitude. The transit depth is given by (see also Line & Parmentier (2016)): where D clear and D cloudy are the clear and fully cloudy transit depth, and f c is the cloud fraction at the terminator. Fig. 12 shows simulated transit spectra from 0.4 to 2 µm with f c = 0, f c = 1 and with the cloud fraction computed in Section 3.3, including the 1 sigma error (blue areas). The cases with 100×solar and 300×solar metallicity are weakly impacted by clouds. Note that for 100×solar, clouds are not optically very thick even for f c = 1. For the case with 300×solar metallicity and 10 6 CCN/kg, the cloud top is a slightly higher and the cloud fraction is much higher and cloud, flattening the transit spectrum in the visible range and reducing the depth in spectral windows. Absorption bands at 1.15, 1.4 and 1.8 µm are still strong. For the case with 1000×solar, the almost complete cloud fraction at the terminator strongly flattens the transit spectrum masking almost completely the molecular absorptions.
For each case, we computed the reduced chi-squared for HST data from Tsiaras et al. (2019), HST data from Benneke et al. (2019b), HST+K2+Spitzer data HST data from Benneke et al. (2019b). Table 4 summarizes these results with the planetary bond albedo, the global cloud fraction and the cloud fraction at the terminator for the different atmospheric compositions. For all datasets, the minimum of chi-squared is obtained for the 300×solar metallicity. A higher concentration of CCN weakly changes the chi-squared. Using HST from Tsiaras et al. (2019), the cases with 300×solar and 100×solar metallicity are within 1 sigma. The case with 1000×solar metallicity is within 2 sigma while the cases with 1×solar and 10×solar metallicity are ruled out. The chi-squared is increased for all cases using   Atm. Transmittance Fig. 9. Transmittance map at a given time for the 300×solar metallicity (CCN concentration of 10 6 nuclei/kg) as seen for a primary transit. The difference between the west terminator (left in the plot) and the east terminator (right in the plot) is due to clouds. metallicity of K2-18 b is likely ∼100-300×solar if it has a solar C/O ratio.
As discussed in Section 3.3, the cloud fraction at the terminator of K2-18b could be highly variable. This would produce variability in the transit depth, correlated to spectral windows. It would be maximum at 0.77 µm. For the most variable case from our simulations (i.e. 300×solar with CCN=10 6 /kg), the transit depth varies with a standard deviation of 14 ppm at 0.77 µm and 7 ppm at 1.07 µm. The transit depth uncertainty of individual HST-WFC3 transits is ∼80 ppm, too large to search for spectral variability even with 9 transits. However, this variability could be detectable with multiple transits observed by JWST-NIRISS. The variability of transit depth due to a variable cloud fraction can be expressed as: D cloudy = D clear at wavelengths probing above the cloud layer and D cloudy is constant in spectral windows (in the optical regime with λ r c ). δD=0 when the stellar light is absorbed above the cloud layer. Measurements with high SNR of transit variability could thus infer D cloudy and provide constraint on the altitude/pressure of the cloud top. If the atmospheric composition can be retrieved from transit depth above the cloud layer, then the cloud fraction at the terminator, its variability and the whole clear transit spectrum could be derived too. Fig. 13 shows the thermal emission flux from 1D simulations with Exo-REM with no cloud. We can notice a reduction of the thermal flux at 15 µm for high atmospheric metallicity due to the absorption by CO 2 . At 15 µm with high metallicity, the secondary eclipse depth is around 40 ppm. Using the JWST Exposure Time Calculator, we found that the photometric uncertainty for one secondary eclipse with MIRI-Imaging is around 70 ppm with F1280W filter (centered at 12.8 µm) and 80 ppm with F1500W filter (centered at 15 µm). Several eclipses would be required to detect K2-18b, which is thus too cold for observations in thermal emission.
Finally complementary information on K2-18b's atmosphere could be provided by high-resolution Doppler spectroscopy from ground-based instruments. Such observations can be efficient to probe cloudy atmospheres and to constraint abundances and cloud-top pressure (Gandhi et al. 2020). As shown by Blain et al. (2020), one full transit of K2-18b observed with VLT-CRIRES+ might be sufficient to detect CH 4 .

Summary and conclusions
In this study, we analysed the atmospheric dynamics, cloud formation and observational implications for K2-18b with a H 2dominated atmosphere with solar C/O. We explored the effects of atmospheric metallicity, CCN concentration and rotation rate. Assuming a synchronous rotation, we found that the atmospheric circulation in the upper atmosphere (above 0.1 bar) corresponds essentially to a symmetric day-to-night circulation. The heat transport is very efficient and there are only modest horizontal temperature changes. This simple regime leads to preferential cloud formation between 2 and 10 mbar at the substellar point or at the terminator, which is the coldest region. We found that water clouds never form for low atmospheric metallicity (i.e. 1×solar and 10×solar). Clouds always form for metallicity 100×solar but the cloud cover is never total. For 100-300×solar metallicity, they preferentially form at the substellar point and the cloud fraction at the terminators is small. For 1000×solar metallicity, the weaker heat redistribution leads to an absence of cloud on the dayside but very thick clouds at the terminator. Due to the high fraction of water vapour, we predict large ice cloud particles with radii of 30-450 µm for a realistic range of CCN concentration. We found that the cloud particle size can strongly impact the cloud distribution because of cloud radiative feedbacks. Increasing the rotation rate to that for a 2:1 resonance does not significantly impact the cloud distribution, although it tends to decrease cloud formation at the substellar point and to enhanced it at the terminators. Because of the inhomogeneous cloud cover and the absorption of stellar radiation by CH 4 and H 2 O high in the atmosphere, the planetary albedo is very low (lower than 0.1).
Comparing transit spectra simulated from the outputs of the 3D model to HST observations, we found that data are compatible with a 100-300×solar metallicity, similar to the conclusions from Bézard et al. (2020) and Blain et al. (2020). A 100-300×solar metallicity would be consistent with the massmetallicity trend from the Solar System (Kreidberg et al. 2014a;Blain et al. 2020), as shown in Fig. 14. For such a composition, clouds have a relatively small impact on transit spectra in the near-infrared, even with high CCN concentrations. An important implication of the day-night circulation is that the cloud formation at the terminator does not affect the abundance of water vapour above clouds. The atmospheric metallicity and the C/O  ratio could be well retrieved with future observations (e.g. JWST, ARIEL, VLT-CRIRES+, ELTs). Unfortunately, it would be very difficult to distinguish cases with the same metallicity but different CCN concentrations or different rotation rates since the cloud distribution is very sensitive to the thermal structure and to the different parameters, leading to degenerated solutions. In contrast to a fast rotation rate, a slow synchronous rotation should exhibit a transmission spectrum blueshifted by the day-night circulation, as observed on hot Jupiters (Snellen et al. 2010;Brogi et al. 2016). But it would only be by ∼100 m/s (10 times lower than the precision obtained with VLT-CRIRES on HD 209458b by Snellen et al. (2010)). Finally, we found that for some parameters, the cloud fraction at the terminator is highly variable. This produces variability in transit spectra correlated with spectral windows. Similar spectral variability is observed on brown dwarfs and attributed to inhomogeneous cloud cover. To our knowledge, transit spectral variability due to clouds has never been studied. We can expect that future JWST observations of cloudy exoplanets could reveal such a variability. Its detection could be used to distinguish condensate clouds from photochemical hazes, which should not present strong horizontal or temporal variability. Another way to distinguish water clouds from photochemical hazes would be transmission spectroscopy in the visible range. Water clouds should produce a flat transit spectrum for the cloudy part while sub-micrometric haze particles should produce a slope due to Mie-scattering Lavvas et al. (2019).
To conclude, K2-18b is a unique target for studying the composition and formation of sub-Neptunes thanks to its relatively clear atmosphere. Interestingly, laboratory work by Hörst et al. (2018) suggest that the photochemical haze production rate is relatively low for 100×solar metallicity gas mixture at 300 K, much lower than for gas mixtures at 400-600 K. Temperate sub-Neptunes might thus be more promising than warm sub-Neptunes for transit spectroscopy. The atmospheric circulation and the cloud formation on temperate sub-Neptunes should also have many similarities with those on rocky planets in the habitable zone of low-mass stars. With the major role played by water clouds on the climate, the characterisation of temperate sub-Neptunes may lead to major advances in the understanding of the habitability of exoplanets. Transit spectra for 100×, 300× and 1000×solar metallicity, for cloud-free cases (solid red lines), fully cloudy cases (dotted red lines), and with partial cloud cover (blue regions). The simulations were performed with 10 6 CCN/kg for the bottom left panel and with 10 5 CCN/kg for the other panels. K2 and HST data from Benneke et al. (2019b) are indicated with black dots. Atmospheric composition profiles used for the radiative transfer in the LMDG GCM and for 1×, 10×, 100×, 300× and 1000×solar metallicity. Chemical abundances are computed with Exo-REM with non-equilibrium chemistry for an eddy-mixing coefficient K zz =10 6 cm 2 /s and for solar C/O ratio. Note that the mixing ratio of water is computed separately in the GCM with coud condensation.