© ESO 2018

1 Introduction

Two decades of exoplanet research has yielded over 3600 confirmed detections¹ and has instilled the knowledge that with great numbers comes great diversity. Planets have been discovered having a wide range of orbital parameters, inferred bulk compositions and, most importantly for the question of habitability, atmospheric properties. With technological constraints limiting our ability to probe small terrestrial planets like our own, some of the best observational data regarding atmospheric composition and dynamics comes from the well studied planetary class of hot-Jupiters. Their large planetary radii coupled with their short-period, and often transiting nature leads to the retrieval of high-resolution spectral datasets. Such data has revealed the presence of a number of chemical species including water, sodium and potassium (e.g. Charbonneau et al. 2002; Sing et al. 2011; Birkby et al. 2017). One of the most important results is that recent observations show, even within this limited sub-set of planets, that there exists a large range in atmosphere types which encompass clear, hazy and cloudy skies (Sing et al. 2016).

Clouds can be challenging because they obscure our view of the underlying atmospheric composition. Their wavelength dependent opacity can act as a blanket that mutes water features (Deming et al. 2013; Helling et al. 2016c) and other key species in the infra-red (IR) region of the transmission spectrum (Pont et al. 2013; Iyer et al. 2016). Helling et al. (2016c) show that obscuration by high altitude cloud in their simulations of HD 209458 b (as opposed to deeper clouds in HD 189733 b) can explain the reduced H₂O amplitude in atmosphere observations of HD 209458 b, and underscores the impact of the vertical cloud structure on direct observables. Additionally, photochemically produced complex hydrocarbons (Zahnle et al. 2016) or small cloud particles lofted into the upper atmosphere both contribute to an enhanced Rayleigh scattering slope in the visual that can, via an increased continuum opacity, mask, sometimes completely, the alkali metal lines (Kirk et al. 2017). The masking of absorption lines in transmission spectra makes placing constraints on absolute and even relative atmosphericchemical abundances difficult. Spectral lines are also key in measuring the doppler shift caused by the fast global winds detected on some of these planets (Snellen et al. 2010; Birkby et al. 2013, 2017; Louden & Wheatley 2015).

Cushing & Roellig (2006) obtain spectra of the L-type brown dwarf DENIS 0255 and infer a silicate feature from a dip in the emission around 10 μm which is not seen in clear skies atmosphere models. No such detection has been made in an exoplanet atmosphere; Richardson et al. (2007) claimed to have detected a broad emission silicate feature at 9.65 μm in HD 209458 b, but it was not found in subsequent studies (Swain et al. 2008). The James Webb Space Telescope (JWST) however will enhance available constraints on models, with the 10 μm wavelength coverage potentially allowing for the first direct detections of clouds via their vibrational mode absorption features (Wakeford & Sing 2015). Sing et al. (2016) suggest that silicate detection requires significant vertical mixing to sustain condensates at observable pressures. However, the primary difficulty in obtaining direct signatures of clouds in the spectra, currently, is due to the quality of the data obtained from exoplanet observations.

Aside from the aforementioned modifications to both the visual and near-IR transmission spectra from cloud scattering and absorption, clouds have also been inferred by other means. One of the more compelling observational indications of inhomogeneous cloud coverage is the westwards offset in the peak optical phase curve, first detected for Kepler-7b by Demory et al. (2013). A westwards shift in the visual curve suggests more efficient backscattering of optical photons near the night-day terminator. Oreshenko et al. (2016) investigates the sensitivity of this offset in Kepler-7b to aerosol and condensate composition. Such a phase curve offset has also been observed in thermal emission (although the presence of clouds is not required to produce this effect); an eastwards offset in the IR phase curve which shows the hot-spot shift from the sub-stellar point (Knutson et al. 2007; Crossfield et al. 2010). Evans et al. (2013) find a decreasing albedo towards longer wavelengths of their HST STIS observations; evidence for day-side optically thick, reflective clouds. A large horizontal temperature gradient, expected for hot-Jupiters, can lead to a cyclic mechanism of night-side condensation to day-side evaporation (Lee et al. 2016). Clouds advected onto the day-side could possibly introduce a higher visual albedo compared to the day-night transition region for a clear skies atmosphere. Further studies have shown that offsets in the visual for Kepler planets are common (Angerhausen et al. 2015; Esteves et al. 2015; Shporer & Hu 2015). Armstrong et al. (2016) observed the first temporal variability of the phase curve from HAT-P-7b which is indicative of variable cloud coverage potentially driven and moderated by global winds. Additionally, Spitzer data from the emission of WASP-43b shows a distinct lack of night-side flux, consistent with increased cloud opacity in this region (Stevenson et al. 2017). For the hottest exoplanets that are expected to exhibit such a hemispherical dichotomy in cloud coverage, ionising processes, for example by the stellar irradiation and cosmic ray impact (Helling et al. 2016a) may open new pathways to form cloud particles even higher in the atmosphere. We would consequently expect the formation of an ionosphere also on exoplanets (Helling & Vorgul 2017). This may be the case for the hot super-earth CoRoT-7b (Helling et al. 2016b; Mahapatra et al. 2017).

Louden & Wheatley (2015) and Brogi et al. (2016) inferred wind speeds of up to 5.3 km s⁻¹ on HD 189733 b, and models have shown that high-velocity, super-rotating jets should be common for tidally locked hot-Jupiters (Showman & Guillot 2002; Mayne et al. 2014b). Zonally asymmetric radiative forcing is responsible for the generation of strong eastwards equatorial flow on these hot-Jupiter exoplanets; a combination of eddy and mean flow momentum transport specifically involved in the jet pumping (Showman & Polvani 2011; Tsai et al. 2014; Mayne et al. 2017). For the formation of clouds this is particularly relevant as the advection of cloud particles from the night-side to the day-side could result in their subsequent evaporation. Likewise, the advection of metal enriched gas can provide condensible material to the night-side. Moreover, the vertical transport of cloud from its coupling to the vertical advection of the gas (from mean flows) and the gravitational settling of cloud particles due to precipitation, can play a large role in the distribution of cloud and gases. As cloud particles drift into the deeper, hotter regions of the atmosphere they can evaporate and enrich specific element abundances (Woitke & Helling 2004; Witte et al. 2009).

We choose to focus on two well studied hot-Jupiters; HD 209458 b and HD 189733 b. Their similar orbital motions, physical size, mass, equilibrium temperatures, and the discovery of water vapour, neutral oxygen, carbon monoxide, and methane, suggest that models indicating they share similar cloud properties are unsurprising (Benneke 2015; Helling et al. 2016c). Given these similarities, the observational differences such as the weak H₂ O absorption depth in HD 209458 b compared to HD 189733 b is surprising (Sing et al. 2016). Helling et al. (2016c) suggest that the vertical distribution of cloud, and hence modification of the atmosphere optical depth, could be responsible for these contrasting observables.

The formation, structure, composition, and evolution of clouds are all sensitive to the local thermal, chemical and dynamical conditions of the atmosphere which in turn respond to the radiative influence of such clouds on their host atmospheres. This complex interplay can be best understood by consistently treating cloud formation as part of a 3D atmosphere or global circulation model (GCM). Much progress has been made recently in modelling cloud-free planetary atmospheres, particularly hot-Jupiters. Simulation results (e.g. Dobbs-Dixon & Lin 2008; Dobbs-Dixon & Agol 2013; Carone et al. 2014; Mayne et al. 2014b; Heng & Showman 2015; Amundsen et al. 2016; Boutle et al. 2017) show a variety of extremely complex atmospheres that can be dominated by radiative flux, momentum-coupled horizontal, and vertical advection and driving from the deep atmospheric thermal flux (in turn driven by overshoot from the convective interior).

Recently, Lee et al. (2016) combined the kinetic cloud formation model developed by Woitke & Helling (2003, 2004); Helling & Woitke (2006); Helling et al. (2008; for a summary see Helling & Fomins) with a 3D global circulation model (Dobbs-Dixon & Agol 2013). The cloud model describes the micro-physics of seed formation (nucleation), surface growth, and evaporation consistently coupled with element conservation, gas phase chemistry, hydrodynamical mixing, and gravitational settling. We utilise the time-dependent representation of the cloud formation model as derived in Helling & Woitke (2006), and couple it to the Met Office GCM, the “Unified Model” (Wood et al. 2014; Mayne et al. 2014b), otherwise known as the UM. We present our model results for the atmospheres of the hot Jupiters HD 209458 b and HD 189733 b. We use this advanced, cloud-coupled UM version in order to present a comparative study between the two similar giant gas planets which only differ significantly in their gravitational strength. Helling et al. (2016c) concentrated on the cloud structures for 1D atmospheres, Lee et al. (2016) presented the first consistent 3D atmosphere solution with cloud formation, and our paper presents a continuation of these modelling efforts, including full radiative treatment (absorption and scattering) of the cloud.

In Sect. 2 we discuss the hydrodynamics, thermodynamics, radiative transfer and cloud formation formalism, including initial and boundary conditions. Our results are presented in Sect. 3 and discussed in the context of previous and complementary studies in Sect. 4. We summarise our results and discuss the model’s limitation and potential improvements in Sect. 5.

2 Numerical methods

We coupled our state-of-the-art numerical 3D GCM with a state-of-the-art microphysical cloud model in order to simulate the dynamics of a cloudy hot-Jupiter atmosphere. To solve for the thermodynamical properties and hydrodynamical motions of the atmosphere, we used the Met Office UM which has been well tested, not only for our own atmosphere (Brown et al. 2012), but also for hot-Jupiter atmospheres in Mayne et al. (2014b, 2017) and Amundsen et al. (2014, 2016).

The microphysical cloud model that we utilised is summarised in Helling & Fomins (2013), and has already been tested in its stationary form for brown dwarfs and giant gas planets by Helling et al. (2008a,b) and Witte et al. (2009, 2011), and by Juncher et al. (2017), in addition to the time-dependent version as part of 2D turbulence modelling in Helling et al. (2001, 2004) and the already mentioned 3D GCM for HD 189733 b in Lee et al. (2016). Figure 1 is a flow schematic of the physical processes of our coupled model which are described in the following sections.

2.1 Hydro- and thermodynamical modelling

The UM uses the dynamical core and finite-difference code ENDGame (Even Newer Dynamics for General Atmospheric Modelling of the Environment, Wood et al. 2014). For each vertical layer the atmosphere is discretised onto a latitude-longitude grid by using a staggered Arakawa-C grid (Arakawa & Lamb 1977). The vertical structure uses a geometric height (as opposed to pressure height) vertical grid using a vertically staggered Charney-Phillips grid (Charney & Phillips 1953). ENDGame is semi-Lagrangian and time stepping uses a semi-implicit scheme (Crank & Nicolson 1996).

The UM’s dynamical core solves five equations that correspond to the solution for a full, deep and non-hydrostatic atmosphere (see Mayne et al. 2014b, for more information): $$\frac{Du}{Dt}=\frac{uv\text{\hspace{0.17em}tan}\varphi }{r}-\frac{uw}{r}+fv-f^{\prime }w-\frac{c_{\text{p}}\theta }{r\text{\hspace{0.17em}cos}\varphi }\frac{\partial \Pi }{\partial \lambda }+D(u),$$(1) $$\frac{Dv}{Dt}=\frac{u^2\text{\hspace{0.17em}tan}\varphi }{r}-\frac{vw}{r}-uf-\frac{c_{\text{p}}\theta }{r}\frac{\partial \Pi }{\partial \varphi }+D(v),$$(2) $$\frac{Dw}{Dt}=\frac{u^2+v^2}{r}+uf-g(r)-c_{\text{p}}\theta \frac{\partial \Pi }{\partial r},$$(3) $$\delta \frac{D{\rho }_{\text{gas}}}{Dt}=-{\rho }_{\text{gas}}\left[\frac{1}{r\text{\hspace{0.17em}cos}\varphi }\frac{\partial u}{\partial \lambda }+\frac{1}{r\text{\hspace{0.17em}cos}\varphi }\frac{\partial (v\text{\hspace{0.17em}cos}\varphi )}{\partial \varphi }+\frac{1}{r^2}\frac{\partial (r^{2}w)}{\partial r}\right],$$(4) $$\frac{D\theta }{Dt}=\frac{Q}{\Pi }+D(\theta ),$$(5) $${\Pi }^{\frac{1-\kappa }{\kappa }}=\frac{R{\rho }_{\text{gas}}\theta }{p_0}.$$(6)

We define the coordinates λ as the longitude, ϕ as the latitude, r as radial distance from the planetary centre, and t as time. The equations are solved for the zonal (u), meridional (v) and vertical (w) wind velocity components as well as the fluid density, ρ_gas, and the potential temperature (temperature at a fixed reference pressure), θ. c_p is the fluid specific heat capacity, R is gas constant, Q is the heating rate, D is the diffusionoperator, κ is the ratio $\frac{R}{c_{\text{p}}}$, Π is the Exner pressure function and the Coriolis parameters f = 2Ωsinϕ and f′ = 2Ωcosϕ, with Ω the rotation frequency of the planet. g(r) is the height-varying gravitational acceleration at radial position r calculated via $$g(r)=g_{\text{p}}{\left(\frac{R_{\text{p}}}{r}\right)}^2,$$(7)

where g_p is the gravitational acceleration, R_p is the radial value for the planetary inner boundary and g(r) is always directed towards the lower boundary in each cell.

To smooth out grid-scale noise and to account for unresolved sub-grid processes such as eddies and turbulence we used the diffusion scheme² in Mayne et al. (2017). Additionally, at the upper boundary we implemented a sponge layer as described in Mayne et al. (2014b).

2.2 The cloud formation model

Cloud formation is a phase transition process from a gas to a solid (or liquid) phase and progresses through two major steps: seed formation (nucleation), surface growth, and evaporation. These processes are complemented by element conservation, gravitational settling, and advection. Helling & Fomins (2013) present the kinetic cloud formation approach derived by Woitke & Helling (2003, 2004), Helling & Woitke (2006), and Helling et al. (2008). We therefore only summarise the main equation body, applicable to this work, that we solved in combination with chemical equilibrium in our 3D radiation-hydrodynamics code. In the following work, “cloud” is used to refer to an ensemble of particles that have a local and a global size distribution, that are made of a mix of materials which can change in time and in space, and which has a certain geometrical extension.

The cloud particle properties in each cell are described in terms of a set of dust moments L_j (r, t) [cm^j g⁻¹] (j = 0, 1, 2, 3). Dust moments are integrated volume-weighted particle size distributions (Woitke & Helling 2003, 2004; Helling & Woitke 2006; Helling et al. 2008): $${\rho }_{\text{gas}}(r,t)L_j(r,t)={\displaystyle {\int }_{V_l}^{\infty }f}(V,r,t)V^{j/3}\text{d}V,$$(8)

where j is the moment index, ρ_gas(r, t) is the gas density, V_l is the seed particle volume [cm³] and f(V, r) is the cloud particle size distribution function [cm⁻⁶].

The dust moments L_j(r, t) are calculated from the dust moment equation (Woitke & Helling 2003): $$\frac{D({\rho }_{\text{gas}}{L}_j)}{Dt}=V_l^{j/3}{J}_{*}+\frac{j}{3}{\chi }^{\text{net}}{\rho }_{\text{gas}}{L}_{j-1}.$$(9)

The evolution of the dust moments in time (left hand side) is dependent on the nucleation of seed particles (first term right hand side) and the net surface growth velocity (second term right hand side). J_* is the nucleation rate [cm⁻³ s⁻¹], and χ^net is the net surface growth velocity [cm s⁻¹] which has a negative value for evaporation (≤0 cm s⁻¹) and positive for growth (≥0 cm s⁻¹). Applying the operator splitting method, Eq. (9) is integrated in time during each hydrodynamic timestep. For this, we applied the implicit ordinary differential equation solver LIMEX (Deuflhard et al. 1987).

The cloud model describes the nucleation process by applying the modified classical nucleation theory to calculate the nucleation rate (e.g. Helling & Fomins 2013; Gail & Sedlmayr 2014; Lee et al. 2015b). The formation of seed particles during the nucleation process provides a surface for further chemical reactions. The cloud particle growth or evaporation velocity, χ^net (r, t) [cm s⁻¹], is calculated considering the summation of the growth of each dust species, s from each contributingreaction r (Gail & Sedlmayr 1986; Helling & Woitke 2006; Helling & Fomins 2013).

The local number density of cloud particles, n_d, [cm⁻³] is calculated by $$n_{\text{d}}(r)={\rho }_{\text{gas}}{L}_0(r),$$(10)

and the mean volume of a cloud particle (which we refer to in the rest of the manuscript as “mean cloud particle volume”) by $$\langle V\rangle =\frac{L_3}{L_0}.$$(11)

Elements are depleted or enriched as a response to the formation or evaporation of cloud. The change in elemental abundance is accountedfor via the following equation, $$\begin{array}{l}\frac{D(n_{\langle \text{H}\rangle }{ϵ}_i)}{Dt}=\sqrt[3]{36\pi }{\rho }_{\text{gas}}{L}_2{\sum }_{r=1}^R\frac{\Delta V_r^s{n}_r^{\text{key}}{v}_r^{\text{rel}}{\alpha }_r}{v_r^{\text{key}}}\left(1-\frac{1}{S_r}\frac{1}{b_{\text{surf}}^s}\right),\hfill \end{array}$$(12)

where n_⟨H⟩ is the total number density of hydrogen molecules and ϵ_i is an element’s relative abundance to hydrogen. For elements where the abundance can be enriched or depleted, the value of ϵ_i is also advected with atmospheric flow, as solved in the UM by Eqs. (1)–(7).

Cloud particles may settle or “precipitate” through the atmosphere with a velocity dependent on the local gas density and cloud particle size. To account for this precipitation, a drift velocity is calculated and used to advect the cloud particles downwards in a column. For the case of free molecular flow (Knudsen number ≥ 1) (Woitke & Helling 2003): $$v_{\text{dr}}=-\frac{\sqrt{\pi }}{2}\frac{g{\rho }_{\text{d}}\langle a\rangle }{{\rho }_{\text{gas}}{c}_s}.$$(13)

Here g is again the gravitational acceleration,ρ_d is the cloud particle mass density, c_s is the speed of sound, and ⟨a⟩ is the mean particle size, $$\langle a\rangle =\sqrt[3]{\frac{3}{4\pi }}\frac{L_1}{L_0}.$$(14)

The sedimentation routine uses a first order Eulerian advection scheme, working downwards in each column from the model top. To maintain stability, an exponential limiter (Rotstayn 1997) is used to prevent violation of the Courant-Friedrichs-Lewy (CFL) condition (Courant et al. 1928).

Fig. 1 The eight-stage DIHRT-UM cloud formation and evolution model flow. Heating rate solutions to the gas and cloud radiative transfer are applied to the UM and feedback onto the dynamics. Absorption and scattering coefficients are found independently for the cloud and gas phase but the radiative transfer equations are only solved once for the combined coefficients.

Fig. 2 Horizontally averaged pressure–temperature profiles of our three atmosphere models, prior to modification by the cloud radiative feedback, at tcloud = 0 days (see Table 1 and Sect. 2.6). Dashed lines show condensation curve data for enstatite (MgSiO3) and forsterite (Mg2SiO4).

Fig. 3 Zonally averaged temperature [K] (left) and zonal wind speed [m s−1] (right) of both the hot HD 209458 b (upper) and standard HD 209458 b (lower) atmospheres at tcloud = 0 days (see Table 1 and Sect. 2.6).

2.3 Cloud specifics

In our model, we consider the seed formation via the homogeneous nucleation of TiO₂. Cloud particle growth from five dust species (TiO₂[s], SiO[s], SiO₂[s], MgSiO₃[s], and Mg₂ SiO₄) is modelled to occur via 22 surface reactions which involves 13 gas phase species forming the solid materials. The elements involved in the cloud formation processes (nucleation, growth and evaporation) are Si, Mg, Ti, and O and these elements (with the exception of oxygen) are depleted or enriched according to Eq. (12) as cloud formation progresses in time and in space. The relevant material data (surface tension, vapour pressures) are the same as in Lee et al. (2016). The initial, unaltered element abundances were chosen to be solar and we adopted the values from Asplund et al. (2009) in order to remain comparable to the work by Lee et al. (2016). The element abundances are input material properties for our gas phase chemistry calculation from which we derive the gas phase composition in local thermodynamic equilibrium (Helling et al. 2001), hence, the number densities of the 13 gas species that are directly involved in the cloud formation processes.

2.4 Gas phase radiative transfer

A chemical equilibrium model solves for the abundances of key species involved in cloud formation. For computational efficiency however, we used the fast analytical formulation of Burrows & Sharp (1999; used in Amundsen et al. 2016) that calculates abundances of key absorbing species involved in our radiative hydrodynamics. This parameterised gas phase chemistry assumes solar metallicity. Analytical formulae for exoplanet chemistry has also been published by Heng & Tsai (2016).In addition, the abundances of the alkali species Na, K, Li, Cs, and Rb are approximated using the same parameterisations described in Amundsen et al. (2016). Since the chemistry scheme used in the radiative transfer stage does not use the elemental abundances depleted by cloud formation, the abundances of key absorbers may well be overestimated.

We used the Suite of Community Radiative Transfer codes based on Edwards and Slingo (SOCRATES) code which is coupled to the UM (Edwards & Slingo 1996). SOCRATES implements the plane-parallel two-stream approximation and uses the sophisticated correlated-k method to avoid computationally demanding line-by-line calculations (Lacis & Oinas 1991; Amundsen et al. 2014, 2016). By solving the radiative transfer equations, we improve upon the Newtonian Cooling approximation used in many previous models (Cooper & Showman 2005; Menou & Rauscher 2009; Heng et al. 2011; Mayne et al. 2014b).

Rayleigh scattering is included for H₂ and He withtheir numerical factors described in Amundsen et al. (2016). Opacities are determined, where possible, from the latest ExoMol line lists (Tennyson & Yurchenko 2012; Tennyson et al. 2016) and the model includes line absorption from H₂ O, CO, CH₄, NH₃, Li, Na, K, Rb, Cs, and H₂ -H₂ and H₂-He collision induced absorption (CIA). Further details regarding the opacity implementation can be found in Amundsen et al. (2014, 2016). Additionally, the treatment of overlapping absorption in the gas phase is described in Amundsen et al. (2017).

A ghostlayer above the numerically modelled domain is used to account for the radiative processes taking place in the atmosphere extending from the lowest pressures of our simulated domain to P ≈ 0 Pa. Thereby theflux at the top boundary is not the exact stellar flux. For more details on this layer, see Amundsen (2015).

2.5 Cloud radiative transfer

Clouds have been shown to have a dramatic effect on atmospheric opacities; their presence can alter the local thermal conditions by way of back-warming (Lee et al. 2016). Mie Theory (Mie 1908) and Effective Medium Theory (EMT) are used to obtain the cloud absorption and scattering coefficients, over the 32 bands described in Amundsen et al. (2014), which are appended to the coefficients for the gas phase and used in the calculation for the full radiative transfer performed by SOCRATES. All cloud properties relevant for the opacity calculations (cloud particle size, cloud particle number density, material composition) are assumed to be homogeneous across a computational grid cell, hence, they are represented by one value. This represents the effect of the grid resolution on our radiative transfer calculation.

The simultaneous condensation of multiple minerals onto the nucleated seed particles leads to a bulk composition that is the sum of the different volume fractions from the individually condensed materials. Since a single particle is composed of multiple materials it contributes mixed optical properties. A quantification of the bulk optical property of the mixed particle is needed, and EMT is used to account for this. We take our material real and imaginary refractive indices and extinction coefficients from Lee et al. (2015a, 2016), using the same extrapolation methods to cover the necessary wavelength range that our radiative transfer scheme extends over. Following the Bruggeman method (Bruggeman 1935), we solve for the average dielectric function of the particle, ϵ_av, using a Newton-Raphson minimisation scheme. For cases of non-convergence, the solution is obtained via the analytic Landau-Lifshitz-Looyenga method (Looyenga 1965). This methodology is applied inboth Lee et al. (2016) and Juncher et al. (2017).

The dielectric function is a necessary input to the Mie Theory solution. Bohren et al. (1983) BHMie routines are used to solve for the absorption, Q_abs(λ,⟨a⟩), and scattering, Q_sca(λ,⟨a⟩), efficiencies. Here, λ is our considered wavelength. Using the expression of size parameter, $$x=\frac{2\pi \langle a\rangle }{\lambda },$$(15)

for particles of similar size to the photon wavelength (10⁻⁶ ≤ x ≤ 10³) where the interaction between the light and particle is non-trivial, the full Mie theory calculation is required. In Lee et al. (2016), analytical expressions are used for the limiting cases of Mie theory; the large particle, hard sphere scattering case (where x ≥ 1000) and the small metallic sphere case (where x ≤ 10⁻⁶). Improved stability in our Mie code means we can instead approximate these regimes using the full Mie code solution but fixing size parameters that exceed these limiting values (for x ≥ 1000, x = 1000 and x ≤ 10⁻⁶, x = 10⁻⁶.)

The resulting absorption and scattering coefficients are scaled by the cloud cross section which is defined as, $$\frac{n_{\text{d}}}{{\rho }_{\text{gas}}}\pi {\langle a\rangle }^2.$$(16)

These cloud scattering and absorption coefficients are then appended to the same coefficients from the gas phase for a total (cloud and gas) opacity.

Using SOCRATES, radiative fluxes and heating rates are calculated using a two-stream solver with a full treatment of scattering for both the stellar and thermal fluxes using the Practical Improved Flux Method (Zdunkowski et al. 1980). Delta-rescaling is performed to reduce errors due to strong forward scattering by cloud particles. This configuration of SOCRATES is also used by the Met Office climate model (Walters et al. 2017) and has undergone extensive validation for clouds on Earth.

2.6 Numerical approach

We performed three simulations that encompass two well-known transiting hot-Jupiters, HD 209458 b and HD 189733 b. Each of our simulations is composed of three numerical phases or “stages”, primarily to isolate the cloud formation dependence on the atmosphere’s thermal profile, and the radiative effect of cloud back onto the atmosphere. These stages are described in detail in Sects. 2.6.2–2.6.4. The first considers a clear-skies, cloudless atmosphere which we term the “cloud free” stage. A further period of time follows in which clouds can form and evolve in response to the atmospheric thermal, chemical and dynamical conditions, but do not radiatively interact with (feedback onto) the atmosphere; this stage we call “transparent cloud”. The final and third stage allows the clouds to both scatter and absorb radiation and hence modify the atmospheric conditions. This stage we term “radiatively active cloud”. To evolve through transient instabilities arising from large heating rates³, caused by steep gradients in the cloud opacity, we introduced a sub-period termed “transient radiative cloud”, which includes a number of numerical aids to maintain model stability during this period. The total combined simulation length extends from time t = 0 to t = t_end days (throughout this paper, “day” refers to an Earth day unless stated otherwise), with t_end, alongside full model parameters, defined in Table 1. Since the length of the cloud free stage differsbetween our three atmosphere models, for ease of comparison we define a new time, t_cloud, which quantifies the length of time since clouds are included.

Fig. 4 Horizontal temperature [K] distribution of both the hot HD 209458 b (upper), standard HD 209458 b (middle), and HD 189733 b (lower) atmospheres at tcloud = 0 days (see Table 1 and Sect. 2.6). Plots at constant pressures of 102 Pa and 2 ×102 Pa for HD 209458 b and HD 189733 b respectively.

2.6.1 Initial conditions

Our models started from rest with an initial, spherically symmetric, temperature–pressure (TP) profile. For HD 209458 b, these profiles represent both non-adjusted and hotter interior atmospheres which we subsequently refer to as the “standard” and “hot” cases respectively (see Amundsen et al. 2016, for details). The TP profile of HD 189733 b was calculated using a similar method to that of Amundsen et al. (2016), but only a standard profile was considered. For P ≤ 10⁵ Pa, and by t = 400 days (under cloud-free conditions), non-evolving wind speed maxima and well resolved flows indicated the atmosphere had reached a quasi-steady-state. For higher pressures the thermal timescale increases significantly Mayne et al. (2014b). The physical motivation for using a hotter “deep” interior is raised by the works of Cooper & Showman (2005); Mayne et al. (2014b); Amundsen et al. (2016) and Mayne et al. (2017) which have shown the deep atmosphere continues to evolve beyond our longest current simulation time of 1200 days, indicating that the thermal structure of a non-adjusted standard atmosphere is likely to be too cool. Further studies have shown that dynamical processes can be responsible for the slow heating of the radiatively inactive regions of the atmosphere, including the downwards flux of kinetic energy (Showman & Guillot 2002), and transportation of potential temperature (Tremblin et al. 2017). Simulations in 1D are unable to capture these physical processes due to the omission of horizontal advection. In 3D, computational feasibility limits the ability to capture such processes acting over long timescales. Therefore, we used a TP profile with +800 K to mimic the increased deep atmosphere temperature effect of Tremblin et al. (2017); this is the method used in Amundsen et al. (2016).

In all models, the zeroth vertical level (z = 0) was set to an initial pressure of P = 2 × 10⁷ Pa. Since the surface gravity of HD 189733 b is larger than that of HD 209458 b, the atmosphere is more compressed and hence the vertical geometric extent of the atmosphere is reduced from ~10 000 km (HD 209458 b) to 3200 km. Correspondingly, the outer boundary pressures vary between our simulations; in HD 189733 b, we simulated up to P ~ 10 Pa, whereas much lower pressures of P ~ 10⁻¹ Pa were probed for HD 209458 b.

The hydrodynamical time step was set to t = 30 s, shorter than that of Mayne et al. (2014b), but consistent with Amundsen et al. (2016). In our case, a shorter time step was required to account for cloud thermal adjustment from advection in high velocity (e.g. jet) regions. Taking a maximum wind speed of 5400 m s⁻¹ which is obtained numerically for HD 209458 b in Mayne et al. (2014b) and observationally for HD 189733 b in Louden & Wheatley (2015), the dynamical timescale of the jet was found to be ~1 day. Hence, the dynamics of the atmosphere are temporally highly-resolved with almost 10⁵ steps per the advective timescale in the most dynamically active regions. Our radiative time step is consistent with the optimised value of t = 150 s in Amundsen et al. (2016), such that radiative transfer is solved every five hydrodynamical steps.

Verticalvelocities were damped using the standard “sponge” technique (see Mayne et al. 2014b; Amundsen et al. 2016). The damping magnitude was set via the damping coefficient, initially set to 0.15. A linear profile was used whereby the damping coefficient is constant, for a given height, with latitude and is implemented from 0.75× z_top (where z_top is the upper boundary height) increasing steeply to z_top. We term this damping geometry, defined for each simulation in Table 1, “linear” (see Mayne et al. 2014a,b, for more information on damping.)

2.6.2 Cloud free

Our cloud free simulations represent the prior state-of-the-art, atmospheric circulation modelling, and included full radiative transfer for gaseous absorbers and scatterers. In this cloud free stage we evolved our initial profiles from rest (v = 0), without clouds (J_* = 0 and χ^net = 0), until the atmosphere has reached a quasi-steady-state, defined by the formation of planetary scale flows (super-rotating equatorial jet) with no significant changes in the maximum zonal wind speed over time. While such conditions are typically met by 400 days, we chose to evolve the cloud-free stage for longer (see Table 1 for details).

2.6.3 Transparent cloud

Subsequent to the cloud free stage, each model ran for a further period of time totalling 50 days (t_cloud = 0–50 days) in which clouds could form and evolve in response to the local atmospheric properties, but were transparent to radiation; no scattering or absorption from cloud particles was considered. This method allowed us to determine the dynamical evolution of the cloud in response to the cloud free temperature profile and to more clearly identify the effects of cloud absorption and scattering on the final atmospheric profile.

One of the difficulties of this self-consistent coupled-cloud model, identified in Lee et al. (2016) is the integration time of the cloud chemistry step. For the first 100 hydrodynamical time steps, we solved for the cloud chemistry at each step. After this period, we then solved the cloud chemistry equations every ten hydrodynamical time steps. This may slightly increase the number of cloud particles that are transported across the night-day terminator, in response to the increased day-side temperature (i.e. cloud that may have evaporated on the next step due to the temperature contrast will remain in the condensed phase for up to nine further hydrodynamical time steps until it is able to respond to the new thermal conditions). This decreased computation frequency means the whole cloud formation complex, including the drift velocity and cloud opacity (both of which are updated after the cloud formation chemistry integration), was also updated less frequently. While this advection of cloud onto the day-side may lead us to slightly overestimate the day-side cloud coverage, we expect this to only apply to the most dynamically active regions (super-rotating equatorial jet) where wind velocities can carry the cloud particles further over the response window of 5 min (10 × 30 s time step). This interval updating was required to obtain the necessary accuracy to resolve the cloud dynamics while maintaining computational feasibility.

Where evaporation occurs quickly (i.e. in response to a sudden thermal change), the solution to the evaporation velocity can become numerically stiff and hence we employed the same technique used in Lee et al. (2016) to avoid this issue. The methodinvolves removing a large amount of the evaporating, volatile material at the start of the step, if the time step required is small. The return of elements to the gas phase, from each contribution dust species s is (Woitke 2006): $${ϵ}_i={ϵ}_i^b+\frac{{\nu }_{i,s}1.427m_u}{V_{0,s}}\Delta L_{3,s}.$$(17)

where ${ϵ}_i^b$ is the element abundance prior to the instantaneous evaporation step, ν_i,s is the stoichiometric coefficient of element i in species s and 1.427m_u (where m_u is the atomic mass unit) is the factor for converting between gaseous mass density ρ_gas and hydrogen nuclei density n_⟨H⟩.

In cells which report a negative growth velocity (χ < 0 [cm s⁻¹]), but for which the solution is not numerically stiff, standard evaporation, described in Sect. 2.2, applies. Material is removed from the particle down to the evaporation limit which is set to the seed particle size. When the evaporation limit is reached, the seed particle still remains. At this point a separate test of the thermal stability of the seed particle (χ of TiO₂) is used. If ${\chi }_{{\text{TiO}}_2}<0$ [cm s⁻¹] then seed particles are removed from the reservoir.

Cloud chemistry was not integrated for cells which contained little to no cloud providing the nucleation rates were small (n_d < 10⁻¹⁰ and J_* < 10⁻¹⁰ [cm⁻³ s⁻¹]). This is computationally efficient and removes problems where if χ^net is large but J_* is small, very large particles occur which have drift velocities that become difficult to handle due to CFL condition concerns. Similarly, if χ^net < 10⁻²⁰ [cm s⁻¹] then there is no particle growth; only significant departures from equilibrium were solved in time. This value has been tested and optimised in Lee et al. (2016).

Fig. 5 Zonally averaged atmospheric profiles of HD 189733 b at tcloud = 0 days (see Table 1 and Sect. 2.6).

2.6.4 Radiatively active cloud

In the radiatively active cloud stage, which immediately follows the transparent cloud stage, we took into account the influence of radiative energy transfer by absorption and scattering from the cloud. Hence, cloud particles could effect the local temperature (hence pressure, density, and velocity) only during this stage. For each of our three gas giant atmosphere models, we simulated this stage for a further period of time totalling 50 days (t_cloud = 50–100 days) which is split into two stages: transient radiative cloud and full radiative cloud.

Transient radiative cloud: during the transient radiative cloud stage, cloud particles were able to scatter and absorb radiation and thus change the thermal properties of the atmosphere by way of heating rates calculated by the radiative transfer scheme. Due to the large cloud opacity gradients, particularly at the cloud “top”, the heating rates can in turn become large. Given our time step, limited by computational feasibility and selected for accuracy and stability, in some circumstances (e.g. regions of low pressure), heating rates in cloud regions can modify the thermal structure of the atmosphere so quickly that transient instabilities appear.

To deal with these instabilities, the following conditions were necessarily applied to HD 209458 b and HD 189733 b for 14 days and 20 days respectively (the time required before the model could be run without the following measures). We introduced a heating rate cap or limit (Q_max = 1.3 × 10⁻³ K s⁻¹) which slows the thermal evolution of the atmosphere. An additional stability enhancing measure included prescribing a number of upper atmosphere layers in which clouds that extend into these regions become transparent to radiation. These transparent layers were set to have a zero cloud opacity (no contribution from scattering and absorption). This stopped the low pressure upper atmosphere from evolving too quickly. To assist with stability we also increased the damping coefficient. Additionally, the sponge layer geometry was switched from a linear profile, to one where the damping effect at polar regions was enhanced; we term this damping geometry “Parametric” (see Drummond et al. 2018a, and Table 1 for more details). These conditions, listed in Table 1, are removed at the end of this transient radiative cloud stage to allow for a further period of time with fully active clouds where no restrictions apply; the heating rate limit and transparent layers were switched off, and the damping magnitude and geometry were returned to their original values.

Conditions were also applied to the cloud model during the radiatively active cloud stage to aid numerical stability and enhance computational efficiency. To avoid expensive Mie calculations for cloud depleted cells we enforced a limit of n_d < 10⁻² cm⁻³ for the opacity calculations. For values less than this the cloud opacity was set to zero. Should the cloud opacity be calculated, a minimum value of κ_cloud = 10⁻⁷ cm² g⁻¹ applied. Drift velocities were calculated only for regions with n_d > 10⁻¹⁰ cm⁻³. Rarely, the cloud chemistry integration can partially fail causing the generation of erroneous cells that contain tiny but non-zero values for situations where no cloud should form. Floating point precision in these almost-cloud-free cells could therefore produce fictitious particles with extremely large sizes. For both these reasons, we limited the radiation transfer coupling to cells containing particles with ⟨a⟩ < 50 μm (since this value vastly exceeded the sizes of particles expected to form) otherwise the opacity was set to zero. These conditions were not removed at the end of the transient radiative cloud stage and continued for the rest of the simulation.

Full radiative cloud: clouds were fully coupled to the radiation scheme during this stage, which continued for a further period of time after the transient radiative stage, until t_end. For all simulations, this stage ended at t_cloud = 100 days, which is the summation of all the cloud-enabled stages; transparent, transient radiative and full radiative. It is worth noting that while the UM does not inherently conserve axial angular momentum (AAM), we tracked its value which was conserved to better than 0.01% in all simulations, by the end of the full radiative cloud stage.

Fig. 6 Time evolution of the horizontally averaged cloud particle number density, nd (left) and mean cloud particle volume, ⟨V ⟩ (right) during the transparent cloud stage of the hot HD 209458 b model at between tcloud = 0 and tcloud = 50 days (see Table 1 and Sect. 2.6).

2.6.5 Extended period

At t_cloud = 100 days, for the hot HD 209458 b model only, we continued the run for a further period of time (six planetary periods) but used higher resolution (500 band) spectral files for a second diagnostic call to the radiation scheme. This allowed us to calculate the emission spectra and phase curves to a much higher accuracy than the standard 32 bands used for the heating rates (see Boutle et al. 2017). This process was additionally carried out for a cloud free atmosphere so that comparisons can be made between the fluxes from the clear and cloudy atmospheres.

3 Results

We begin in Sect. 3.1 by examining the first stage of our simulations; a clear-skies, cloud-free atmosphere. By doing so we can identify the key thermal and dynamical features, as well as the differences between our chosen profiles, which can then be used to explain the transparent cloud results in Sect. 3.2. Finally in Sect. 3.3, we focus on the results of our radiatively active, cloudy atmospheres. As described in Sect. 2.6, we simulate two planetary atmospheres; HD 209458 b and HD 189733 b. We remind the reader that for HD 209458 b we use two temperature profiles, a standard atmosphere (solved for radiative-convective balance in 1D models) and a hot atmosphere with a hotter deep atmosphere to account for predicted dynamical processes, unable to be captured by 1D solutions and unable to be retrieved in practical timeframes in 3D simulations, that drive the deep atmosphere to much higher temperatures (+800 K at P = 10 bar) than the standard solution.

3.1 Cloud free

As discussed in Sect. 2.6, we evolve our three atmospheric profiles (HD 209458 b (Standard + Hot) and HD 189733 b) without clouds with the UM for at least 800 days. Figure 2, shows the horizontally averaged TP profiles for each atmosphere at the end of the cloud free stage.

3.1.1 HD 209458 b

Alongside Fig. 3 which shows the temperature and zonal windspeed of HD 209458 b, we can see that for P < 10⁶ Pa, our hotter HD 209458 b is thermally and dynamically similar to that of the standard model which shows only minor differences compared with the results of Mayne et al. (2014b) and Amundsen et al. (2016). Since the irradiation is the dominant heating source for these tidally locked planets, the energy deposition is largest at the equator, with a corresponding increase in temperature near the equator with respect to higher latitudes. This trend reverses for P > 10⁶ Pa where the polar regions become warmer than the equator (see discussion in Mayne et al. 2014b, 2017). The key identifying feature of the cloud free standard model is the temperature inversion, where similar temperatures are predicted at both P = 10⁵ Pa and P = 10⁷ Pa. This local temperature inversion does not exist for the hot model. At the millibar pressure level there is no change between simulations in the horizontal temperature contrast, shown in Fig. 4. The fast super-rotating equatorial jet is clearly resolved in Fig. 3 (right) with minimal differences between the two thermal profiles. Both jets possess an almost identical shape and velocity (≈6.8 km s⁻¹), as well as the depth, approximately P = 10⁶ Pa, at which thejet core extends into the atmosphere (the penetration depth). Cooler temperatures are found at high latitudes away from the maximum stellar energy deposition at ϕ = 0°, λ = 0°.

Fig. 7 Hot HD 209458 b slices at P = 102 Pa (upper), P = 103 Pa (middle), and P = 105 Pa (lower) for the cloud particle number density, nd [cm−3] and horizontal wind velocity [m s−1] (left) and temperature [K] (right). Sampled at tcloud = 50 days (see Table 1 and Sect. 2.6), prior to cloud radiation coupling.

Fig. 8 Zonal mean of the cloud particle number density, log(nd) [cm−3] for hot HD 209458 b (upper), standard HD 209458 b (middle) and HD 189733 b (lower) atmospheres at tcloud = 50 days (see Table 1 and Sect. 2.6).

3.1.2 HD 189733 b

From the horizontally averaged pressure–temperature profiles in Fig. 2, the HD 189733 b profile is shown to be cooler than the standard HD 209458 b atmosphere, with a maximum temperature difference of Δ T_max ≈−200 K. However, similar to the standard HD 209458 b simulation, for a cloud free atmosphere, a temperature inversion occurs at P = 10⁶ Pa. In Fig. 4 minor differences, compared to HD 209458 b, in the horizontal temperature contrast are observed, with decreased temperatures on the night-side. The vertical distributions of temperature and zonal wind velocity for the clear skies stage are shown in Fig. 5. A comparison between the HD 189733 b and the standard HD 209458 b show very few differences; the most noticeable change is that for HD 189733 b, the jet appears to advect heat more efficiently than HD 209458 b due to the higher equatorial night-side temperatures.

3.2 Transparent cloud

In this section we consider the cloud evolution and structure prior to including cloud scattering and absorption. Cloud formation (seed nucleation and surface growth) and evolution processes (advection, precipitation and compositional changes) are reacting to the dynamical, thermal and chemical properties of our clear skies, cloud-free atmospheres. At this stage, clouds cannot modify the atmospheric thermal and dynamical properties, as this situation applies in the following radiatively active clouds stage. We stress that although transparent to radiation, clouds are not post-processed but form and evolve with the atmosphere’s local thermochemical properties.

3.2.1 Hot HD 209458 b

The convergence tests of Lee et al. (2016) indicate that cloud is in a relatively steady-state throughout their simulation (from 0 to 65 days) with no significant change in the cloud particle number density from 20 days onwards. In Fig. 6, analysis of the cloud particle number density over the first 50 days of our simulation indicate that the cloud is in a quasi-steady-state, with only small changes in particle number density in the region of P ~ 10² Pa. Cloud particles form in dense decks (n_d > 5000 cm⁻³) between P = 8 × 10⁵ and P = 10 Pa, with cloud particle number density tapering off significantly to the upper simulation boundary. The slow evolution of the cloud particles in regions accessible toobservations is most likely attributed to low nucleation rates from the low pressure environment. The evolution of the cloud particle volume, ⟨V ⟩, which reveals the surface growth efficiency after nucleation, is seen in Fig. 6 to increase from P = 10³ Pa to P = 10 Pa. Both evolutionary plots show a swift decline of cloud particles and volume, at the base around 5 × 10⁵ Pa. This is due to the settling of cloud particles into the deeper atmosphere which is too hot for seed nucleation and particle surface growth. We find it important to note that while we may obtain a quasi-steady-state over the fast cloud formation timescale, the slow processes of vertical mixing and settling mean we are unable to probe the long term evolutionary effects of the cloud. Tackling this problem in 3D GCMs for a fully consistent, active cloud model is beyond the scope of this work.

The temperature and cloud particle number density at various pressures within the atmosphere are shown in Fig. 7. In all regions where nucleation and surface growth is possible, bands of cloud particles form outside of the equatorial jet. The distribution of particles in the jet at P = 10² Pa (see Fig. 7a) is time dependent, with the less thermally stable components of the mixed-composition particles evaporating in response to advection over the high temperature, “hot-spot” region (centred at λ = 200°) and re-condensing on the night-side. There are small longitudinal asymmetries in the particle distribution with enhanced numbers on the cooler night side. This is particularly clear in Fig. 7c where at the deeper P = 10³ Pa level the hot-spot offset (advected heat from the maximum irradiation at the sub-stellar point, e.g. Showman & Guillot 2002; Mayne et al. 2014b, 2017) is more clearly defined, leading to a more pronounced eastwards shift in cloud particle depletion. At P = 10⁵ Pa vortices that form across the night-day terminator lead to the trapping of cloud that is seen as enhanced particle number densities in these regions.

The zonally averaged cloud particle number density (Fig. 8) reveals a more complete picture of the cloud structure. The cloud bands are part of two out-of-equator decks that extend from the top of the simulation domain (P ≈ 10 Pa) at the poles, down into the atmosphere to P ≈ 10⁶ Pa for all latitudes. Diverging winds from the jet cause seed particles that form there to advect to higher latitudes, although a decreased number remain trapped at the equator. These winds also constantly advect condensible gas to these higher latitudes where the cooler temperatures on the night side increase the saturation and lead to more efficient cloud formation. This can be seen clearly in Fig. 7 where for all pressures plotted, there are higher cloud particle abundances found on the night-side. This is particularly true for P = 10³ Pa where the hot-spot etches out a contrasting area of low cloud particle number density.

The cloud top has a wing-like shape where the cloud particle number density drops with decreasing latitude due to the increasing temperatures towards the hotter jet region. The cloud base for the hot HD 209458 b is well defined since the temperature increase is linear at pressures of around P = 10⁶ Pa with no strong latitudinal temperature gradient (see Fig. 3). This means the transition at which nucleation is no longer possible (J_* = 0) due to an inhibiting temperature, occurs at the same pressure level. There are two regions of enhanced particle number density; one extending from P = 10³ Pa to P = 5 × 10⁴ Pa and another for a narrow pressure range at the cloud base. The former is likely due to the balance between settling cloud particles from the upper atmosphere and preferential thermal conditions for nucleation. At the cloud base, while the temperature is low enough for seed particle nucleation, it is too hot for the surface growth of silicates. Therefore, all available TiO₂ material either goes into nucleating seed particles (and increasing seed number densities) or is the only species available for surface growth on these seeds. The latter effect is shown to be true in Lee et al. (2015b).

The depletion and enrichment of elemental abundances is a consequence of the cloud formation and evaporation processes. In Fig. 9a the density scaled abundances of our three tracked elements are shown as a function of atmospheric pressure. Titanium is depleted down towards the cloud base at P = 8 × 10⁵ Pa. At higher pressures the cloud particles are unable to form and therefore the elemental abundances remain at their initial solar metallicity values. Very small enhancements are made to the metallicity below the cloud base due to the evaporation of particles settling into the hot-zone. Silicon and magnesium are depleted for P ≤ 4 × 10⁵ Pa leading to the narrow TiO₂[s] only layer at the very base of the cloud between P = 4 × 10⁵ Pa and P = 8 × 10⁶ Pa. The growth velocities between these pressures for our magnesium and silicon based condensates is zero since the temperature in this region is too high for thermal stability, whereas TiO₂[s] is able to exist at the higher temperatures. Figure 9b shows the contribution, by volume, of each surface growth species. This TiO₂ layer is clearly identified, with the volume fraction increasing for pressures lower than P = 4 × 10⁵. A preference for the formation ofSiO₂[s] means that SiO[s] composes only a small fraction of the cloud particles.

The depletion of cloud particles at the equator observed from the horizontal slices of the atmosphere shown in Fig. 7 is seen clearly in Fig. 8. Number densities five times smaller than those seen in the areas of highest concentration, exist in the jet between P = 10² Pa and P = 10⁵ Pa, with these numbers decreasing with time. While a smaller number of particles persist in the jet, many of the particles have advected out of the equatorial region and into the cloud decks. Since the jet region is hotter than at higher latitudes, it does not allow for large surface growth velocities. This means that the particles near the equator are typically smaller in size than those in the higher latitude bands.

The mean particle size (Fig. 10, left) shows a clear correlation with the cloud particle number density; where cloud particles exist in large number densities the cloud particle size is small and where the particle number density is low, the cloud particle sizes are large. In areas with small numbers of particles, more condensible material grows on a single particle, increasing the size. This effect is similarly seen in the mean cloud particle volume (Fig. 10, right) where the volumes are largest in regions with the lowest particle number densities. With the exception of the vortex regions near the cloud base, cloud particles are consistently sub-micron sized. In HD 209458 b, since clouds scatter strongly as a function of the number and size of their particles, a combination of the cloud particle number density and the mean particle size is a more useful metric than the cloud particle volume. In further analysis, we will consider the cloud distribution to be defined primarily by the number of cloud particles and not total volume.

Fig. 9 Chemistry in the hot HD 209458 b atmosphere: depletion of elemental abundances and volume contribution of dust species at tcloud = 50 days (see Table 1 and Sect. 2.6). Panel a: Zonal and meridional mean, density scaled, elemental abundances of Ti, Si, and Mg. The initial, undepleted values are the solar element abundances. Panel b: Zonal and meridional mean dust volume fractions.

Fig. 10 Zonal average of the (left) mean cloud particle size, ⟨a⟩ [μm] and (right) mean cloud particle volume, ⟨V ⟩ [cm−3], for the hot HD 209458 b atmosphere at the end of the transparent cloud stage at tcloud = 50 days (see Table 1 and Sect. 2.6).

3.2.2 Standard HD 209458 b

Since the standard HD 209458 b temperature profile is almost identical to the hotter one for P ≤ 5 × 10⁵ Pa it is not surprising that differences in cloud structure between these simulations are minimal. The zonally averaged particle number density (Fig. 8) shows a cloud deck that extends from the top of the simulation domain to P ~ 10⁶ Pa. The cloud base matches the hottest part of the atmosphere, with the temperature decreasing again at higher pressures. The temperature inversion seen in Fig. 3 (lower, left) is enough to allow for nucleation and surface growth, with high pressures supporting cloud formation. This leads to the formation of two cloud decks, the deep atmosphere one with vastly higher cloud particle number densities.

There is a similar increase in cloud particles at the cloud base, and areas of depletion between P = 10⁴ Pa and the cloud base at polar regions. Unlike for the hot profile, enhanced cloud wings just below P = 10³ Pa do not exist, but we instead see an increase in cloud particles from P = 10³ Pa and upwards. Again, particles persist in the jet and more so at high latitudes in the upper atmosphere, suggesting a less efficient settling process.

3.2.3 HD 189733 b

Unlike the standard HD 209458 b model, the temperature does not increase, prior to the inversion, to levels that inhibit cloud formation. Hence, the large number densities from the pressure supported high-growth region at pressure higher than P = 10⁶ Pa are seen at the base of Fig. 8.

3.3 Radiatively active cloud

After a period of running the simulations with transparent clouds, we couple particles to the radiation scheme. The following results consider the atmosphere at the most evolved point, after the combination of the transient radiative and full radiative cloud stages.

Fig. 11 Thermal, dynamical and cloud properties of the hot HD 209458 b atmosphere after tcloud = 100 days (see Table 1 and Sect. 2.6) with 50 days of radiatively active clouds.

Fig. 12 Longitudinally averaged gas temperature [K] (left) and particle number density (right) for the hot HD 209458 b model at the end of the transparent cloud stage at tcloud = 50 days (see Table 1 and Sect. 2.6) and the end of the radiatively active stage at tcloud = 100 days. For each time, three profiles are plotted that cover latitudes of ϕ = 0° (black), ϕ = 45° (blue) and ϕ = 60° (red).

3.3.1 Hot HD 209458 b

From t_cloud = 50 days clouds begin to scatter and absorb radiation. 14 days of a transient radiative cloud stage are required to maintain simulation stability as the atmosphere approaches a steady state, defined such that the temperatures are no longer significantly changing between successive time steps. A further 36 days of fully radiatively active clouds are then simulated. The final atmospherictemperature is shown in Fig. 11 (upper, left) and the zonal wind speed in Fig. 11 (upper, right). The difference in the zonal mean temperature shown in Fig. 12 shows that the temperature of the atmosphere decreases from the upper simulation boundary to pressures of approximately P = 10⁶ Pa and all latitudes see a temperature decrease of a similar magnitude. This cooling effect is due to the large proportion of the incoming stellar flux that is back-scattered instead of being absorbed (as per cloud free simulations). The largest cooling is found at the top of the atmosphere, in the equatorial region, with Δ T reaching up to − 250 K. At ϕ = 45° and ϕ = 60°, the zonally averaged profiles show small heating around 2 × 10⁵ Pa. This limited heating region traces the cloud base where cloud particles are better absorbers due to their larger (~1 μm) size, and due to their interaction with the planet’s long-wave thermal flux.

At the equator, a large pressure range of more efficient heating (Δ T ≈ +100 K) occurs near 10⁴ Pa, matching the base of the small particles trapped in the jet. Figure 12 also reveals the change in particle number density from cloud radiation coupling. At the equator, almost every pressure level shows the particle number density decreasing when the cloud becomes radiatively active, but this decrease halts towards the end of the 50 days, potentially as a result muted dynamics (reduced vertical and meridional wind speeds) that occur during this radiative stage. This indicates that despite the cooler temperatures providing enhanced conditions for particle growth, the transfer of particles out of the jet from lower to higher latitudes is still the dominating factor. Additionally, the jet becomes depleted of metal rich gas which inhibits further nucleation and growth of cloud particles. Conversely, for ϕ = 45°, probing the cloud decks, the cooler temperatures and advection of particles from the jet cause the particle number density to rise for P ≤ 2 × 10³ Pa. At ϕ = 60°, there is an increase in cloud particle number density above P = 10² Pa, but this value decreases at all other altitudes. The net effect is to increase the contrast between the vertically extensive and particle number dense mid to high latitude cloud decks, and the less cloudy equatorial jet.

While valuable information regarding the atmosphere’s thermal change is obtained from the pressure–temperature profiles, it is useful to understand how the clouds may impact the direct observables. In Fig. 13, the cloud optical depth as a function of atmospheric pressure, for four different wavelength ranges is shown. The efficiency of the Mie scattering from sub-micron particles is reflected in the scattering optical depth increasing with decreasing wavelength, since the Mie scattering strength increases when the particle size matches the interacting radiation. The role of cloud, in the absorption of stellar flux, is limited to the longest wavelengths. For λ < μm, the cloud scatters all incoming stellar irradiation, highlighting the importance of this process in cloudy atmospheres. In Fig. 14, the reflected stellar flux per layer is shown for varying pressure depth. The difficulty of penetration by the stellar irradiation, due to cloud scattering, is clearly seen with the bulk of the returned stellar flux tracing the cloud top (only minimal stellar radiation is absorbed here, as per Fig. 13). Even the most transparent cloud in the equatorial jet allows minimal stellar flux to travel below the millibar pressure level.

Figure 15 plots the spectral energy distributions for the cloudy and clear planetary flux, in addition to the received stellar flux, both normalised to W m⁻² for an observer at a distance of 1 au. The planetary flux shows the total of the reflected stellar flux and the emitted thermal flux from the planet’s atmosphere. We can begin to appreciate the dominance of the scattering processes from the cloud due to the large fractional return of the incoming stellar flux from the emission and reflectance in the SED, from the cloudy atmosphere. In the cloud-free simulation, strong absorption from sodium and potassium is seen in the visual to near-infrared. This absorption is lost in the cloud simulation which sees a proportional return of the incoming stellar flux, and thus highlights the impact of cloud on obscuring chemical species. The reduced stellar heating from the high albedo of the atmosphere leads to the natural consequence of a cooling atmosphere, due to the lowering of the atmosphere’s total energy budget. A combination of the cooler atmosphere, and the increased optical depth to thermal radiation due to the cloud, leads to reduced outgoing flux at thermal infra-red wavelengths compared to the cloud free simulations. Indeed, the flux in the Spitzer bandpass between 3.6 μm and 4.5 μm is reduced by an order of magnitude when clouds are active.

With radiatively active clouds, the highly efficient Mie scattering and hence reduced absorption in the visual leads to an increase in flux in the Kepler bandpass, seen in Fig. 16. We find no westwards offset in the visual phase curve, but this is expected due to cloud particles persisting across all planetary longitudes; a situation made possible by the high condensation temperature of our silicate particles. While the increased albedo from cloud coverage is detectable in the Kepler bandpass, variability is not. This is not true however for the thermal emission, Fig. 17, where we not only retrieve the eastwards offset due to the advected hot-spot, but find a variability in the flux. This flux variability is a consequence of the time-dependent nature of the cloud distribution (opacity) in the equatorial jet region. In Fig. 18 we show the evolution of the cloud particle number density and upwards thermal flux for the upper atmosphere pressure level of P = 1 Pa. By analysing these values in the jet at two times, separated by approximately half the jet advective timescale, features in the inhomogeneity of the cloud coverage and emergent flux are offset by roughly λ = 180°. Large, non-transient deviations from the mean cloud particle distribution, opacity and hence thermal flux therefore evolve with the advective timescale, driving the variability seen in the thermal phase curves. In addition to the variability, the cooling atmosphere and increased IR opacity results in an order of magnitude reduction, compared to clear-skies, in the Spitzer-band flux. This worsens the agreement with observations of emission, which for the day-side is closely matching the clear-skies predictive flux, suggesting that we are missing some physics from our cloudy model.

Figure 11 (upper, right) shows the maximum wind velocity in the equatorial jet has increased by over 100 m s⁻¹, but more noticeable is the change of the jet shape, broadening at higher pressures. The out-of-equator westwards, retrograde flow is replaced with eastwards flow such that the majority of the atmosphere is now in prograde rotation (we note that AAM is conserved to better than 0.01%).

In Fig. 11 (middle, left) the cloud structure remains in a two deck, out-of-equator configuration. The width of the particle depleted region at the equator increases in response to the widening of the jet itself. The top of atmosphere pressure increases due to the changing thermal profile, but the efficient sedimentation of particles produces a well defined cloud top. The depleted cloud region from P = 10⁴ Pa still exists, along with the TiO₂ layer at the base. The deck loses its layered structure however, with a more stochastic or disorganised nature to the particle number density distribution. The lofted particles in the upper atmosphere in the equatorial region have drifted down to higher pressures, possibly in response to less efficient vertical mixing from the dampened winds. The mean cloud particle volume, shown in Fig. 11 (middle, right) indicates that the large growth velocities at higher pressures continues to result in more voluminous cloud particles near the base, and this is reflected by the increased mean particle size in Fig. 11 (lower, left).

Figure 11 (lower, right) shows the cloud particle number density for a horizontal slice at P = 10² Pa. There is little zonal asymmetry in the cloud coverage, clearly indicating the ability for silicate cloud particles to persist across the night and day side of the planet. Since the day-side temperature at this pressure level does not exceed the condensation temperature of the silicate materials used in our simulation, cloud particles can grow and persist at all longitudes. There is a stark contrast with the distribution of cloud particles at the same pressure level for the end of the transparent cloud stage (cf. Fig. 7a). Between ϕ = 30°–60° there still remains a slight bias towards cloud on the cooler night side, but it is more subtle. For ϕ >60° the reverse is true, with a large contrast between the day and night and regions of sparse cloud coverage on the night side. Interestingly, the low cloud particle number density regions correspond to cold spots that would usually indicate favourable conditions for cloud formation. Particles originally trapped at the equator at t_cloud = 50 days still exist but are smoothed out over the full longitude of the atmosphere. On the night side at higher latitudes there are regions of extremely small cloud particle number densities which can be explained by the low night side pressures and hence fast drift velocities that settle the cloud particles more efficiently to higher pressure levels.

The surface growth velocity is displayed in Fig. 19. Cloud particle growth is still occurring at the end of the simulation, although the velocities are low due to the depletion of elemental abundances. Particles continue to grow between P = 0.01 Pa and P = 5 × 10³ Pa, as well as near the cloud base between P = 10⁵ Pa and P = 10⁶. The higher pressures at the base result in more efficient cloud particle growth but this region is also supported by the upwards mixing of elements in the evaporative zone. There is a distinct lack of cloud particle growth between P = 2 × 10³ Pa and P = 10⁵ that could be attributed to higher, unfavourable temperatures that would inhibit cloud formation without the additional support of high pressures.

Where cloud exists, the drift or settling velocity of the particles (Fig. 20a) never drops below 2 cm s⁻¹, and is mostly sustained above 10 cm s⁻¹. The lowest velocities are found at P = 10³ Pa which corresponds to a region of low particle number density. The highest drift velocities (v_dr > 1 m s⁻¹) are found inthe very upper atmosphere due to the low atmospheric densities. To investigate the balance of cloud particle settling and vertical mixing, the difference between the settling velocity and the vertical wind speed is plotted in Fig. 20b, such that the cloud mixing parameter defined as, ζ = v_dr − w, (where w is the vertical wind speed) is positive for net upwards transport of cloud particles. On the planet night-side, vertical equilibrium is achieved for P > 10⁴ Pa, but rarely exceeds this balance. This means that cloud particles around the base are suspended on the night-side, but do not see a net upwards advection. The day-side is far more conducive to vertical mixing due to the faster vertical wind speeds. For 10 Pa < P < 10⁴ Pa, cloud particles are mixed upwards efficiently, despite precipitation velocities reaching 1 m s⁻¹. The transition to net upwards advection correlates with the well defined cloud top at P = 10 Pa. The ability for cloud particles to be lofted into the upper atmosphere, or at least suspended at higher altitudes (we do not resolve turbulent mixing in our simulations), is supported by Parmentier et al. (2013) who perform simple analysis of TiO tracer particles to find that lofting can occur for particles of 5 μm or less. Such a result appears to be necessary to explain the “hazy” transmission spectra of some planets where condensation curves are not crossed at observable pressures; the dredging of particulates from the deeper atmosphere is required to contaminate the observationally accessible upper atmosphere (Sing et al. 2016).

Fig. 13 Meridional mean at λ = 180° of the cloudunitless optical depth. The scattering component of the optical depth are shown by the solid lines and the total opticaldepth (scattering and absorption) from the cloud is shown by the dashed lines.

Fig. 14 Upwards (returned) stellar flux (W m−2) for λ = 180°.

Fig. 15 Spectral energy distributions of the active, cloudy, hot HD 209458 b (black line) and cloud free (blue line) simulations. The spectrum is taken at secondary eclipse (ϕ = 180°) and includes the reflected stellar and emitted thermal flux. Cloud free absorption is from Na and K.

Fig. 16 Visual phase curve of the active, cloudy, hot HD 209458 b simulation for flux at 1 au received in the Kepler bandpass between 400 and 800 nm. Black points correspond to six complete orbits after the simulation end at tcloud = 100 days (see Table 1 and Sect. 2.6). Blue points follow the flux from a cloud free simulation.

Fig. 17 Phase curves of the active, cloudy, hot HD 209458 b simulation for flux at 1 au received between 3.5 and 4.5 μm. Black dots correspond to six complete orbits after the simulation end at tcloud = 100 days. The red line is the flux from the thermal emission of the atmosphere only and is to be compared with the first orbit (uppermost black dots), and does not included reflected stellar flux.

Fig. 18 Upwards thermal flux (upper) and cloud particle number density (lower) at P = 1 Pa in the hot HD 209458 b atmosphere after tcloud = 106 (left) and 106.5 (right) days. Latitude has been confined to the equatorial jet (− 25° < ϕ < 25°).

3.3.2 Standard HD 209458 b

The initial atmospheric temperature profiles, shown in Fig. 2, reveal a similar or in some regions identical thermal structure above P = 10⁶ Pa between the hot and standard HD 209458 b models. Maybe surprisingly then, the radiatively active and transparent temperatures and cloud particle number densities which are plotted in Fig. 21 show for the standard HD 209458 b, a vastly differing picture to the hotter interior simulation. For the standard simulation, a stronger cooling is found for all latitudes, with the cloud deck at the 1 mbar pressure level leading to a shift in temperature by up to 250 K. For higher latitudes, the thermal properties of the atmosphere do not change below 5 × 10⁴ Pa, even though the temperature inversion leads to the formation (in the transparent cloud stage) of a expansive cloud deck (n_d ≫ 10⁶ cm⁻³). At the equator however, atmospheric cooling takes place for P ≤ 10⁶ Pa. Only the equatorial profile sees a small heating, around P = 2 × 10⁴ Pa. The largest difference in particle number density occurs for the equator where cloud particles decrease for P ≤ 10³. In the higher latitude cloudy regions there is a net tendency to decrease the number of cloud particles, but both latitudes increase in number for the lowest pressures.

Fig. 19 Horizontally averaged growth velocity, χnet for the active hot and standard HD 209458 b atmospheres at tcloud = 100 days (see Table 1 and Sect. 2.6).

3.3.3 HD 189733 b

Due to the larger cloud particle number density a 20-day period of transient radiative cloud was necessary to achieve model stability before fully active clouds were enabled for the final 30 days. Changes to the thermal properties of HD 189733 b during the active cloudy stage are shown in Fig. 22 (upper, left) where the zonal mean of the atmospheric temperatureis shown. As per both variants of HD 209458 b, highly scattering cloud returns a large proportion of the incoming stellar flux, reducing the net atmospheric energy budget and leading to cooling. This effect is better seen in Fig. 23 where zonal means of both temperature and cloud particle number density are plotted for three latitudes. Unlike HD 209458 b, the temperature does not vary much with latitude. In the transparent cloud stage at the milli-bar pressure level the temperature contrast between the equator and high latitudes is only 200 K (compared to 500 K for HD 209458 b). Cooling of up to 200 K occurs at the atmosphere top (P = 1 Pa) for all latitudes. Temperatures diverge from the transparent cloud profiles for P ≤ 10⁴ Pa for ϕ = 45° and only P = 10³ Pa for ϕ = 60°, but the equator sees a temperaturechange for P ≤ 10⁶ Pa, showing that the equatorial region is the most effected by the scattering clouds. A maximum atmosphere temperature change of Δ T_max = −200 K is found in the jet for a wide range of pressures around P = 10⁵ Pa. Like the standard HD 209458 b simulation, the presence of the thermal inversion means cloud extends all the way down to the lower simulation boundary at2 × 10⁷ Pa. The lack of a defined cloud base means no regions of heating, with the exception of a small 50 K rise at the equator which coincides with the penetration depth of the jet structure.

Figure 22 (upper, right) shows the zonal mean of the longitudinal or zonal velocity, which when compared to the values at the end of the transparent cloud stage in Fig. 5 (right) reveals significant changes to the dynamics of the atmosphere. In response to the aforementioned thermal changes, at the lowest pressure the super-rotating equatorial jet core appears to narrow from 30° to 20°. However, lower velocities outside of the jet core lead to prograde flow that increases to, at its widest point, 40°. Similar to what is seen in HD 209458 b, the reverse, westwards flow at higher latitudes decreases in velocity and meridional extent, meaning more of the atmosphere experiences a prograde flow. The jet itself has a lower velocity, likely due to the reduced energy pumping from decreased stellar heating.

The zonal mean of the cloud particle number density and volume are shown in Figs. 22 (middle, left) and 22 (middle, right), respectively. As per the standard HD 209458 b atmosphere, the temperature inversion at the deepest pressures leads to cloud forming down to the lower simulation boundary. Since the temperatures in the deep atmosphere are still sufficiently low to enable condensation, the high pressures are supportive of and increase the surface growth (and nucleation) rates leading to a deck with two orders of magnitude more particles than the upper atmosphere, above the inversion. For P ≤ 10⁶ Pa the equator is again less voluminous and numerous in cloud particles than higher latitudes, although there is less of a contrast than for HD 209458 b. The only regions devoid of cloud particles are the lowest pressures at the jet boundary where settling has precipitated the cloud particles faster than it could be mixed upwards. At higher latitudes the cloud extends to the very upper boundary, meaning the cloud top is at P ≈ 1 Pa. A horizontal slice of the cloud particle number density at the P = 100 Pa pressure level (Fig. 22 – lower, right) reveals a more disorganised structure, with a less cloudy jet but high cloud particle number density at the jet boundary and higher latitudes. As expected from the condensation temperature of silicate particles at this pressure, there is little asymmetry in the distribution of cloud particles as a function of longitude. The effect of changing thermal properties of the atmosphere due to active clouds is best viewed via the radiatively transparent and active profiles of the particle number density shown in Fig. 23. Between the atmosphere top and P = 10⁴ Pa, the cooler temperatures lead to an increase in particles at the equator but a small decrease at higher latitudes. The moderately higher temperatures at the equator due to the geometry of irradiation means that, prior to active clouds, the nucleation rate is slower than the surrounding higher-latitude and more saturated areas. Hence the decrease in temperature from the scattering clouds means the remaining condensible gases trapped in the jet at the equator can now form seed particles. From P > 5 × 10⁵ Pa at the equator and for P > 5 × 10³ at higher latitudes, there is no divergence from the pre-active stage seed particle number density. This is because, particularly below the temperature inversion at P ≈ 10⁶ Pa, all condensible gases have formed the maximum number of seed particles. Lower temperatures that could enhance seed nucleation or particle growth are of no use if the elemental reservoir is dry.

The mean particle size is shown in Fig. 22 (lower, left). The largest particles are found around the deep atmosphere polar hot-spots and correspond with a decrease in the number of cloud particles. In these regions the low nucleation rate has produced cells that contain a small number of particles that provided limited sites for further surface growth and as a result exceed 1 μm in size. However, for the majority of the atmosphere the particles are sub-micron in size, consistent with our HD 209458 b atmospheres. The sensitivity of particle or droplet size with the ratio of nucleation to surface growth, in addition to changing reaction rates across both a vertically or horizontally varying temperature profile reveals one of the significant limitations of post-processing cloud (Parmentier et al. 2016).

The lower drift velocities for HD 189733 b compared to HD 209458 b, shown in Fig. 20a indicate that gravitational settling is less efficient in this atmosphere. Indeed the mixing parameter in Fig. 20b reveal this to be a partially correct assumption since unlike HD 209458 b, there is a fast net upwards advection of cloud on the planet nightside for P > 1000 Pa. This is not true for HD 209458 b which does not see upwards advection of cloud particles for any pressure on the nightside. However, there is only lofting of cloud particles on the dayside for P > 100 Pa.

Due to cloud structure similarities (in particular, the lack of zonal asymmetry) between HD 209458 b and HD 189733 b, we do not include in this work, an analysis of the observation properties (emission spectra and phase curves).

Fig. 20 Vertical transport of cloud, considering the settling velocities (left) and the cloud mixing parameter (right). Left panel: Horizontally averaged drift velocities, vdr [m s−1], for both the active hot and standard HD 209458 b atmospheres at tcloud = 100 days (see Table 1 and Sect. 2.6). Right panel: Meridional mean of the cloud mixing parameter (ratio of vertical windspeed to drift velocity) as a function of pressure. The red line follows the night-side at λ = 0° and the blue line follows the day-side at λ = 180° with solid as HD 209458 b and dashed as HD 189733 b.

Fig. 21 Longitudinally averaged gas temperature [K] (left) and particle number density (right) for the standard HD 209458 b model at the end of the transparent cloud stage at tcloud = 50 days (see Table 1 and Sect. 2.6) and the radiatively active stage at tcloud = 100 days. For each time, three profiles are plotted that cover latitudes of ϕ = 0° (black), ϕ = 45° (blue), and ϕ = 90° (red).

Fig. 22 Thermal, dynamical and cloud properties of HD 189733 b after tcloud = 100 days (see Table 1 and Sect. 2.6) with 50 days of active clouds. Black regions correspond to areas with values below the colour-bar minimum.

Fig. 23 Longitudinally averaged gas temperature [K] (left) and particle number density (right) for thestandard HD 189733 b model at the end of the transparent cloud stage at tcloud = 50 days (see Table 1 and Sect. 2.6) and the final active stage at tcloud = 100 days. For each time, three profiles are plotted that cover latitudes of ϕ = 0° (black), ϕ = 45° (blue) and ϕ = 90° (red).

4 Discussion

All of our simulations show a tendency for clouds to exist in thick decks, at high latitudes, outside of the equatorial jet. It would be easy to assume that the higher temperatures in the jet region contribute to conditions that inhibit the formation of cloud. However, at the millibar pressure level, the maximum temperature for the hot HD 209458 b does not exceed 1100 K, well below the T ~1500 K condensation temperature of our silicate condensates (with the exception of the less thermally stable SiO[s]). At the start of the simulation, both the nucleation and growth velocities around the equator are positive, indicating that cloud formation is active in this region. We find that winds diverging from the equatorial jet lead to the meridional transport of cloud particles to higher latitudes. This advective process also applies to gas which is constantly fed to higher latitudes where seed particle nucleation and surface growth is enhanced from the higher supersaturations due to the lower temperatures. This strong dynamical and thermal interplay does not remove all material from the equator, but leaves a smaller number of cloud particles, particularly for HD 209458 b. The role of meridional (and zonal) winds in this process accentuate the limitation of using cloud that is post-processed from 3D atmospheric profiles; coupled models are necessary to handle the atmosphere dynamics and their influence on cloud particles.

Since, except for the high temperature deep atmosphere, the atmospheric temperature does not generally exceed the condensation temperature of our silicate materials, cloud particles can form and persist on both the night and day sides of both planets investigated. A possible explanation for the discrepancy between our results and Lee et al. (2016), who find that the day-side temperatures evaporate the cloud formed on the cooler night-side leading to a cyclic cloud formation process, is that their day-side temperatures are artificially high due to the non-scattering model where the calculated scattering cross section is added entirely into absorption. Our work indicates that the cloud scattering far outweighs absorption and hence this approximation is unlikely to capture the true radiative balance of the atmosphere. Absorption dominates only at the cloud base where the larger particles are more efficient at absorbing the outgoing thermal radiation, leading to small regions of heating beyond the initial temperature. Our distribution of cloud particles actually shows a slight bias towards higher cloud particle number densities on the day-side, with a depleted region at the southern polar region. Such an effect could be due to the pressure difference of the atmosphere between the day and night side, with cloud particles at a set physical height advecting between higher, cloud formation supporting pressures on the day side to lower pressures on the night side. The effect demonstrates the importance of capturing the horizontal dynamics of the cloud.

A result consistent with both HD 209458 b and HD 189733 b is that after clouds are coupled to the radiative transfer scheme and thus begin altering the thermal profile of the atmosphere by scattering and absorption, the distribution of cloud particles becomes more disorganised. This is particularly true when looking at the cloud particle number density at P = 100 Pa. Prior to the active stage, the cloud coverage is largely inhomogeneous with latitude due to the jet dynamics and thermal effects, as well as showing slight bias to night-side formation from cooler temperatures. After clouds are allowed to scatter and absorb for 50 days, a much more stochastic distribution of cloud appears, particularly for HD 189733 b. Variability in the cloud opacity causes fluctuating heating rates that lead to changes in atmospheric density. In HD 189733 b, density enhancements running parallel and just outside of the equatorial jet lead to sharp increases in cloud particle numbers, since the number density is a linear function of the gas density.

Due to the symmetric distribution of cloud particle number density and volume between the east and west limbs of the atmosphere, the cloud albedo is close to constant across the stellar-illuminated atmosphere. The consistency in the reflected flux with longitude means that we find no westwards offset in the visual phase curve. The work of Parmentier et al. (2016) suggests that a planetary equilibrium temperature exceeding 1500 K would be necessary to introduce evaporative regions for MgSiO₃. For temperatures in our HD 209458 b atmosphere, we would need to include the more volatile metal sulphides as a condensate species to have a similar day-night cloud cycle as found in Lee et al. (2016). However, since our particles are mixed composition it is likely the evaporation process would happen on and off the more stable silicate particles we explore here, meaning rather than day-side evaporation of single-specied cloud decks, we would be seeing large bulk changes in the particle volume due to the removal of the volatile components of that particle. Again, the value of capturing the horizontal advection of cloud particles is demonstrated by the advection-driven temporal variability seen in the thermal phase curves, and will be critical in modelling observations of atmospheric variability in planets such as HAT-P-7b (Armstrong et al. 2016).

Despite the large contrast in cloud particle number density between the equator and higher latitudes, Fig. 14 which showsthe outgoing stellar irradiation, reveals the equatorial region is almost as “active” as the more optically thick cloud at higher latitudes, suggesting that only a small covering of tiny particles are necessary to effectively Mie scatter a significant proportion of the irradiation. Our primary result is that our sub-micron sized, silicate-based and hence highly scattering cloud particles result in a large fraction of the stellar irradiation being blocked from penetrating deep into the atmosphere. This stellar irradiation, which would otherwise be absorbed by gas species in our cloud free simulations, no longer contributes to heating. The net energy budget of the atmosphere reduces and the overall effect is for the atmosphere to cool. As a result of the cooler conditions and the increased optical depth due to the cloud, the outgoing thermal radiation is reduced and contributes to a much lower thermal emission flux. The implications for a cooler atmosphere are significant. Firstly, they provide us with what may be the most self-consistent pressure–temperature profiles of HD 209458 b to date. Secondly, the change in conditions invoke us to consider condensates that would not be undergoing active condensation at our pre-active cloud temperatures. It it worth cautioning that our current non-inclusion of condensates with a reduced diaphaneity will likely overestimate the scattering intensity of our cloud. We endeavour to explore a wide range of potential condensates, including iron and corundum, in future work and discuss the implications of this in Sect. 5. A secondary effect to the cooling atmosphere in response to highly scattering, radiatively active clouds is a weakened meridional and vertical flow; the maximum vertical windspeed for HD 209458 b is reduced typically by an order of 5. Interestingly, zonal flows are maintained, with the jet strength stable across the full length of the simulation. The mechanism behind this dynamical change is beyond the scope of this work and will be explored in future studies.

It is a valuable exercise to consider how the results may change when relaxing some of our necessary model limitations. One of these is the simplification that all seed particles are formed via the homo-molecular nucleation of TiO₂. In theory, there would be a wide chemical range of additional cloud condensation nuclei, from additional oxides to complex hydrocarbons (aerosols). Despite the low abundance of titanium, seed particles of this composition form freely, in large quantities throughout the atmosphere, in our simulations. The addition of further seed species is therefore unnecessary, for at least providing condensation sites, as there are plentiful cloud condensation nuclei available for surface growth. For this reason, including more nucleation species is unlikely to change the overall horizontal distribution of the cloud. However, should seed particles form in increased abundance due to multiple nucleation pathways, the mean particle size will likely be affected; increased seed particles will result in surface growth occurring over a larger number of particles, diluting the volume per seed. As cloud particles become smaller in size, their vertical dynamics will change; precipitation will be reduced as a result of lower settling velocities, potentially leading to increased cloud coverage in the upper atmosphere. Additionally the size of the particles is a strong function of their radiative properties. The smaller particles will contribute a large scattering (as shown in Fig. 13), potentially altering the thermal balance of the atmosphere. It is also worth considering the effect additional seed species may have on the overall chemical composition of the atmosphere. Even if additional nucleation does not lead to seed that play a large role in subsequent condensation, such seed particles would remove their component elements from the gas phase.

One of the key problems concerning cloud evolution is understanding the interplay between vertical mixing and gravitational settling or precipitation. Our simulations investigate the short term dynamics of the cloud particles, mostly dominated by the fast horizontal advection. On much longer timescales, if clouds are not supported by upwards advection they may settle out of the observable atmosphere and no longer contribute the same optical properties that we explore here. With our simulation times unable to capture the long timescales of the deep atmosphere precipitation, we must consider the balancing of vertical wind and settling velocities as a way of quantifying the net direction of cloud particles. The problem is more complex than this however as once cloud has settled into the evaporative region and returned to its gas phase constituents, there is no longer a particle to settle and the condensible gases are more easily mixed upwards (since they have no associated settling velocity) into cooler regions where they can re-nucleate. Providing there exists even the most minimal updraft near the cloud base, this should be enough to sustain a cloud deck. The deck thickness is then essentially defined by the “convective mixing” and settling timescales much like in the Ackerman & Marley (2001) model. Again, post-processing procedures for diagnostic cloud severely restricts the ability to capture this cyclic process of settling, evaporation and re-nucleation. In our work the vertical mixing effect can partly be seen via the large growth velocities near the cloud base in Fig. 19. We have also shown via the mixing parameter in Fig. 20b, the balance between settling and updraft velocities, the net direction or mixing efficiency of the cloud. For HD 209458 b in particular, we pay attention to the large net upwards advection experienced by cloud particles on the planet day-side, suspending and in some regions lofting micron to sub-micron particles. This supports the conclusion of Sing et al. (2016) who find that particulate lofting is necessary for some planets to make the upper, observable atmosphere hazy or cloudy. Additionally, the suspension of cloud particles in our simulations is consistent with the results of Parmentier et al. (2013) who find that global advection patterns lead to a well-mixed upper atmosphere (P ≳ 10⁵ Pa). For HD 189733 b this particle lofting also extends to the planet nightside, helping to keep cloud particles at high altitudes. A surprising outcome, perhaps, is that despite the diminished vertical updrafts during the radiatively active cloud stage, many regions of the atmosphere can still either suspend (no net vertical motion) or loft (net upwards advection) cloud particles. For the planetary atmospheres we study here at least, this might attenuate the importance of the turbulently transported cloud particles, a process we are unable to model at these numerical resolutions, but is parametrised in many 1D cloud formation models as a vertical eddy diffusion parameter, K_zz, which includes both bulk and turbulent flow.

In our introduction we discussed the lack of a trend in the level of observable cloudiness with planetary equilibrium temperature. This is particularly true for both HD 209458 b and HD 189733 b where despite very similar atmospheric temperatures, Sing et al. (2016) show that while HD 189733 b shows evidence of cloudiness and haziness from the strong Rayleigh scattering slope, HD 209458 b appears to have a less (but still) cloudy atmosphere with more defined alkali metal absorption, reduced Rayleigh scattering and a stronger H₂ O absorption. Our simulations indicate that both planets would have cloudy atmospheres, although HD 189733 b has much higher cloud particle number densities. The higher surface gravity of HD 189733 b would suggest that precipitation is a much more relevant dynamical component to the evolution of the cloud than for HD 209458 b and therefore it should be expected that cloud particles are more likely to settle out of the observable atmosphere and not contribute towards Rayleigh scattering. However we find that drift velocities in the upper atmosphere of HD 189733 b are comparable to those of HD 209458 b and in fact, coupled with the strong mean vertical atmospheric flows, lead to enhanced lofting of cloud particles in HD 189733 b. Therefore lofting in HD 189733 b, along with the large cloud particle number densities from efficient cloud formation, could lead to the enhanced cloud opacity inferred from the HD 189733 b transmission spectra. It is worth noting that although higher drift velocities may be expected from increased values of the gravitational acceleration, stronger gravity can lead to atmospheric compression, larger local gas densities and hence stifled drift velocities. Therefore it is wise to consider the balance between drift velocities and the vertical wind speeds.

The temperature of the deep atmosphere, below the cloud base, remains constant despite the overall cooling of the atmosphere.This could partly be a result of the long radiative timescales at these pressures which have not yet adjusted to the thermal changes in the upper atmosphere over a short (50 day) period. It is also possible that the cloud base, shown to absorb the intrinsic planetary thermal emission, keeps the temperatures here high by trapping the outgoing convective heat flux. A much longer simulation time, beyond the scope of this work, would be needed to disentangle these effects.

An interesting observation is that more cloud forms in the upper atmosphere of the hot HD 209458 b atmosphere compared to the standard profile. Why should this occur when the thermal structure of these atmospheres is consistent above P = 10⁵ Pa? The thermochemical conditions should be identical between the two. The answer likely lies in the deep atmosphere where the two differ in their thermal structure. For the standard atmosphere, the temperature inversion allows a thick cloud deck to form, exhausting the elemental abundances. For the hotter atmosphere, precipitating particles evaporate at the condensation temperature transition and enrich the gas phase. These elements may be rapidly advected upwards to lower pressure regions that lead to increased surface growth (and nucleation) velocities that produce either more seed particles or larger grains. This highlights the important link between the deep and observable atmosphere due to cloud mixing and non-equilibrium chemistry; while the deep atmosphere is not observationally probable, understanding the thermodynamical properties of this hidden region is imperative in retrieving the correct distribution of cloud and hence optical properties of the upper observable atmosphere.

5 Summary and further work

Clouds formed easily in the atmospheres of HD 209458 b and HD 189733 b due to large extents of their atmospheres having temperatures below the condensation temperature of the majority of our silicate-based growth species. The exception was the deepest regions of the hot HD 209458 b atmosphere, which unlike the standard profile of HD 209458 b, did not see a temperature inversion. Thus, in the deepest atmosphere, temperatures exceeded that which would allow for nucleation and particle surface growth. A banded cloud structure formed as a result of increased particle numbers and cloud particle volume at higher latitudes. This effect was a combined result of cloud particle and gas advection from the equator to higher latitudes, the latter leading to increased saturation and hence efficient cloud formation. This meridional variation in cloud coverage can only be seen in 3D (or potentially 2D) models, highlighting the importance of including the full horizontal advection that 1D models are unable to capture.

In cloud-free simulations, the atmosphere maintains its thermal balance by interior thermal emission and the absorption of stellar flux. The presence of our cloud had a significant effect on the thermal profile of the planet. Abundant sub-micron sized particles efficiently Mie scattered a large proportion of the incoming stellar short-wave flux. With agreatly reduced stellar absorption, the energy budget of the atmosphere decreased and subsequently temperatures decreased by up to 250 K. The reduced atmospheric temperature and opaque cloud blocking the emerging intrinsic planet flux lead to a significant decrease in the thermal emission. Correspondingly, the increased albedo in the visual, from scattering, boosted the returned stellar flux, making the planet much brighter in the Kepler bandpass.

Since there was no zonal asymmetry in the distribution of cloud particle number density or volume, with cloud extending in equal abundance across the day and night side, there was no detectable offset in the visual phase curve. The presence of meteorologyon the planet however was detectable by both the increased albedo in the visual and also the variable thermal flux originating from the advection of the disorganised cloud structure in the equatorial jet. In addition, the modification of the atmosphere’s thermal profile changed the wind speeds on the planet. Particularly, a widening of the equatorial jet was seen, leading to prograde flow at higher latitudes.

Sing et al. (2016) indicated a weak-cloudy atmosphere for HD 209458 b and cloudy and strongly hazy conditions for HD 189733 b. We found that both are cloudy at the observationally important millibar pressure level, with the higher gravity of HD 189733 b not settling cloud particles any more efficiently than in HD 209458 b. The moderately cooler temperatures in HD 189733 b did allow for cloud to form more abundantly however, with the particle number density increasing by an order of magnitude over HD 209458 b, making HD 189733 the “cloudier” planet in both particle number and volume. Future studies will include synthetic transmission spectra⁴ from the UM output which will assist in comparison with the observations of exoplanet atmospheres, such as those from Sing et al. (2016). This will, in addition to our already included emission spectra and phase curves, allow for a more complete understanding of cloud’s observational signatures. We also endeavour to improve the consistency between chemical phases which will allow for us to identify detectable species (such as the Sodium and Potassium) above the cloud photosphere, in the model spectra.

The success of coupling such a powerful cloud model to the UM naturally invites further investigations into exoplanetary nephology. Our chosen planets occupy a similar temperature range and therefore applying the model to a variety of atmospheric temperatures would be the natural progression; what clouds form and where, when the atmosphere is increasingly irradiated? Further studies will include the simulation of such super-heated atmospheres to explore the day-night condensation-evaporation cycle due to the large horizontal temperature contrasts. The application of the model to Brown Dwarf atmospheres would also be a useful test and may help us to determine why the temperature sequence in hot Jupiters does not correlate with cloud abundance and yet appears to do so for Brown Dwarfs.

This sophisticated coupled model may also be able to tackle further questions raised by observations. For example, WASP-121b has evidence of TiO and VO absorption in transmission (Evans et al. 2016), and water features that appear in emission on the day-side but absorption on the night-side (Evans et al. 2017) both pointing to the presence of a thermal inversion at low pressures. Such a different thermal profile is likely to have an effect on the presence of cloud. Could a forward 3D model with radiative clouds predict the presence of gas phase TiO and VO at the relevant pressures?

A number of improvements to our model are required in order to tackle, more accurately, these aforementioned questions. These additions are not an added layer of complexity but instead represent the minimum framework for a physically and chemically complete model; for example, no atmosphere is likely to contain only the limited condensible species in our current work, and some observations may be entirely explained by photochemistry. Such a complete model will be crucial in making the most meaningful comparisons with observations of exoplanetary atmospheres, and we list the most pressing improvements below.

• Introducing additional condensates is an essential step in improving the radiative effects from cloud coverage. While we have included a number of silicate species, in addition to titanium dioxide, there are many condensates missing from our model that are predicted to play a significant role in cloud formation.

• Cloud formation is dependent on, amongst other factors, the availability of gas phase species, which in turn depends on both the prescribed metallicity and the local thermochemical conditions. In future studies, we must explore both a changing metallicity and carbon-oxygen ratio to determine cloud sensitivities to this expanding parameter space; doing so will be useful to those in the planet formation community who are attempting to constrain their own models from these values.

• Our cloud model tracks the depletion of elements and solves for the key condensation reactants using a simple equilibrium scheme. These depleted or enriched element abundances are not coupled to our radiative transfer scheme which instead uses a parametrised calculation of the abundances of key gas phase absorbers (Burrows & Sharp 1999). We aim to improve chemical consistency in our model by introducing our own coupled flexible Gibbs energy minimisation scheme which will simultaneously provide gas phase abundances for both cloud formation and radiative transfer.

• Drummond et al. (2018b) also provides us with the option to couple a chemical relaxation scheme to our cloud model. This will provide more accurate determinations of the abundances of condensible gases, but also allow us to simultaneously assess the contributions of both clouds and the cloudless chemical instability theory (Tremblin et al. 2015). Such a process will be both challenging and time-consuming but a necessary step for full chemical consistency.

• Highly irradiated planets such as HD 209458 b and HD 189733 b may well have upper atmospheres with active photochemistry. With the importance of the Rayleigh scattering in observed spectra, it will be prudent to investigate the role of photochemically produced hydrocarbons. Adding a photochemical network will help to determine whether the Rayleigh scattering is set primarily by these complex species, tiny suspended condensates in our existing model, or a combination of both.

• The importance of particle sizes on the cloud’s radiative effects also invites us to consider implementing a more physical size distribution, rather than to assume a per-cell uniform particle size. In addition to the more accurate Mie calculations from such an applied size distribution, there will be improved accuracy to the drift velocities. While particle size is a crucial parameter in the cloud’s optical properties, the shape of the cloud particles also plays a role. Currently, particles are assumed to condense into perfect spheres. Realistically, dusty particles may have irregular shapes, much like ice crystal in Earth’s atmosphere. The assumption of spherical particulates may underestimate the scattering strength and must therefore be removed to understand the tendency of an atmosphere to net absorb or scatter.

Our results are an exciting introduction to modelling radiatively active clouds in 3D, but a number of model enhancements are in the pipeline. In addition, a much wider variety of planetary atmospheres should be explored to understand the change in cloud formation processes in super-hot gas giant to cool terrestrial atmospheres. While we present the most sophisticated 3D-coupled model of a cloudy atmosphere, it is complementary to other 1D, 2D and 3D schemes – both coupled and post-processed, which are required to tackle the complex nature of this field.

Acknowledgements

S.L. is funded by and thankful to support from the Leverhulme Trust. The calculations for this paper were performed on the University of Exeter Supercomputer, a DiRAC Facility jointly funded by STFC, the Large Facilities Capital Fund of BIS, and the University of Exeter. Material produced using Met Office Software. BD acknowledges funding from the European Research Council (ERC) under the European Unions Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement no. 336792. D.S.A. acknowledges support from the NASA Astrobiology Programme through the Nexus for Exoplanet System Science. EKHL acknowledges support from the Universities of Oxford and Bern through the Bernoulli Fellowship Programme. The authors would also like to thank the anonymous referee for their insightful report. The research data supporting this publication are openly available from the University of Exeter’s institutional repository at: https://doi.org/10.24378/exe.323.

References

All Tables

All Figures

	Fig. 1 The eight-stage DIHRT-UM cloud formation and evolution model flow. Heating rate solutions to the gas and cloud radiative transfer are applied to the UM and feedback onto the dynamics. Absorption and scattering coefficients are found independently for the cloud and gas phase but the radiative transfer equations are only solved once for the combined coefficients.
In the text

	Fig. 2 Horizontally averaged pressure–temperature profiles of our three atmosphere models, prior to modification by the cloud radiative feedback, at tcloud = 0 days (see Table 1 and Sect. 2.6). Dashed lines show condensation curve data for enstatite (MgSiO3) and forsterite (Mg2SiO4).
In the text

	Fig. 3 Zonally averaged temperature [K] (left) and zonal wind speed [m s−1] (right) of both the hot HD 209458 b (upper) and standard HD 209458 b (lower) atmospheres at tcloud = 0 days (see Table 1 and Sect. 2.6).
In the text

	Fig. 4 Horizontal temperature [K] distribution of both the hot HD 209458 b (upper), standard HD 209458 b (middle), and HD 189733 b (lower) atmospheres at tcloud = 0 days (see Table 1 and Sect. 2.6). Plots at constant pressures of 102 Pa and 2 ×102 Pa for HD 209458 b and HD 189733 b respectively.
In the text

	Fig. 5 Zonally averaged atmospheric profiles of HD 189733 b at tcloud = 0 days (see Table 1 and Sect. 2.6).
In the text

	Fig. 6 Time evolution of the horizontally averaged cloud particle number density, nd (left) and mean cloud particle volume, ⟨V ⟩ (right) during the transparent cloud stage of the hot HD 209458 b model at between tcloud = 0 and tcloud = 50 days (see Table 1 and Sect. 2.6).
In the text

	Fig. 7 Hot HD 209458 b slices at P = 102 Pa (upper), P = 103 Pa (middle), and P = 105 Pa (lower) for the cloud particle number density, nd [cm−3] and horizontal wind velocity [m s−1] (left) and temperature [K] (right). Sampled at tcloud = 50 days (see Table 1 and Sect. 2.6), prior to cloud radiation coupling.
In the text

	Fig. 8 Zonal mean of the cloud particle number density, log(nd) [cm−3] for hot HD 209458 b (upper), standard HD 209458 b (middle) and HD 189733 b (lower) atmospheres at tcloud = 50 days (see Table 1 and Sect. 2.6).
In the text

	Fig. 9 Chemistry in the hot HD 209458 b atmosphere: depletion of elemental abundances and volume contribution of dust species at tcloud = 50 days (see Table 1 and Sect. 2.6). Panel a: Zonal and meridional mean, density scaled, elemental abundances of Ti, Si, and Mg. The initial, undepleted values are the solar element abundances. Panel b: Zonal and meridional mean dust volume fractions.
In the text

	Fig. 10 Zonal average of the (left) mean cloud particle size, ⟨a⟩ [μm] and (right) mean cloud particle volume, ⟨V ⟩ [cm−3], for the hot HD 209458 b atmosphere at the end of the transparent cloud stage at tcloud = 50 days (see Table 1 and Sect. 2.6).
In the text

	Fig. 11 Thermal, dynamical and cloud properties of the hot HD 209458 b atmosphere after tcloud = 100 days (see Table 1 and Sect. 2.6) with 50 days of radiatively active clouds.
In the text

	Fig. 12 Longitudinally averaged gas temperature [K] (left) and particle number density (right) for the hot HD 209458 b model at the end of the transparent cloud stage at tcloud = 50 days (see Table 1 and Sect. 2.6) and the end of the radiatively active stage at tcloud = 100 days. For each time, three profiles are plotted that cover latitudes of ϕ = 0° (black), ϕ = 45° (blue) and ϕ = 60° (red).
In the text

	Fig. 13 Meridional mean at λ = 180° of the cloudunitless optical depth. The scattering component of the optical depth are shown by the solid lines and the total opticaldepth (scattering and absorption) from the cloud is shown by the dashed lines.
In the text

	Fig. 14 Upwards (returned) stellar flux (W m−2) for λ = 180°.
In the text

	Fig. 15 Spectral energy distributions of the active, cloudy, hot HD 209458 b (black line) and cloud free (blue line) simulations. The spectrum is taken at secondary eclipse (ϕ = 180°) and includes the reflected stellar and emitted thermal flux. Cloud free absorption is from Na and K.
In the text

	Fig. 16 Visual phase curve of the active, cloudy, hot HD 209458 b simulation for flux at 1 au received in the Kepler bandpass between 400 and 800 nm. Black points correspond to six complete orbits after the simulation end at tcloud = 100 days (see Table 1 and Sect. 2.6). Blue points follow the flux from a cloud free simulation.
In the text

	Fig. 17 Phase curves of the active, cloudy, hot HD 209458 b simulation for flux at 1 au received between 3.5 and 4.5 μm. Black dots correspond to six complete orbits after the simulation end at tcloud = 100 days. The red line is the flux from the thermal emission of the atmosphere only and is to be compared with the first orbit (uppermost black dots), and does not included reflected stellar flux.
In the text

	Fig. 18 Upwards thermal flux (upper) and cloud particle number density (lower) at P = 1 Pa in the hot HD 209458 b atmosphere after tcloud = 106 (left) and 106.5 (right) days. Latitude has been confined to the equatorial jet (− 25° < ϕ < 25°).
In the text

	Fig. 19 Horizontally averaged growth velocity, χnet for the active hot and standard HD 209458 b atmospheres at tcloud = 100 days (see Table 1 and Sect. 2.6).
In the text

Fig. 20 Vertical transport of cloud, considering the settling velocities (left) and the cloud mixing parameter (right). Left panel: Horizontally averaged drift velocities, vdr [m s−1], for both the active hot and standard HD 209458 b atmospheres at tcloud = 100 days (see Table 1 and Sect. 2.6). Right panel: Meridional mean of the cloud mixing parameter (ratio of vertical windspeed to drift velocity) as a function of pressure. The red line follows the night-side at λ = 0° and the blue line follows the day-side at λ = 180° with solid as HD 209458 b and dashed as HD 189733 b.

In the text

	Fig. 21 Longitudinally averaged gas temperature [K] (left) and particle number density (right) for the standard HD 209458 b model at the end of the transparent cloud stage at tcloud = 50 days (see Table 1 and Sect. 2.6) and the radiatively active stage at tcloud = 100 days. For each time, three profiles are plotted that cover latitudes of ϕ = 0° (black), ϕ = 45° (blue), and ϕ = 90° (red).
In the text

	Fig. 22 Thermal, dynamical and cloud properties of HD 189733 b after tcloud = 100 days (see Table 1 and Sect. 2.6) with 50 days of active clouds. Black regions correspond to areas with values below the colour-bar minimum.
In the text

	Fig. 23 Longitudinally averaged gas temperature [K] (left) and particle number density (right) for thestandard HD 189733 b model at the end of the transparent cloud stage at tcloud = 50 days (see Table 1 and Sect. 2.6) and the final active stage at tcloud = 100 days. For each time, three profiles are plotted that cover latitudes of ϕ = 0° (black), ϕ = 45° (blue) and ϕ = 90° (red).
In the text