© P. Auclair-Desrotour and K. Heng 2020

1 Introduction

Both observations and theoretical studies achieved over the last decade suggest that rocky planets can be found around stars of different masses and they represent a large fraction of the population of exoplanets (e.g. Cumming et al. 2008; Mordasini et al. 2009; Howard 2013; Ronco et al. 2017). Many of these planets are located in the habitable zone of their host stars (e.g. Kopparapu et al. 2013), which is basically the region where the incident stellar flux and greenhouse effect enable them to sustain surface liquid water. This is particularly the case for planets detected around brown dwarfs and very-low-mass stars (e.g. Payne & Lodato 2007; Raymond et al. 2007; Kopparapu et al. 2017).

Three out of seven of the Earth-sized rocky planets hosted by the TRAPPIST-1 ultra-cool dwarf star (planets e, f, and g) that were discovered in 2017 are within the habitable zone (Gillon et al. 2017; Grimm et al. 2018). Recently, two temperate Earth-mass planet candidates were detected around the M-dwarf Teegarden’s Star by the CARMENES spectrograph using radial-velocity measurements (Zechmeister et al. 2019), while an Earth-sized planet orbiting in the habitable zone of the TOI-700 red dwarf was discovered by the TESS observatory (Gilbert et al. 2020; Rodriguez et al. 2020). This planet (TOI-700 d) is the first of its kind detected by TESS and it might harbor temperate surface conditions (Suissa et al. 2020). The number of discovered temperate rocky planets will keep growing in the future with the upcoming transit searches of the TESS (Barclay et al. 2018) and PLATO (Ragazzoni et al. 2016) observatories, which opens a new era of atmospheric characterisation and motivates a better understanding of the mechanisms that govern planetary climates and surface conditions.

Most planets located in the habitable zone of sub-Solar-mass stars, and particularly those orbiting very-low-mass dwarfs such as TRAPPIST-1, are expected to be circularised and tidally locked in spin-orbit synchronous rotation (e.g. Kasting et al. 1993) or in a spin-orbit resonance (e.g. Correia et al. 2014) over short timescales. As rocky planets locked into synchronous rotation exhibit permanent dayside and nightside, their climate is affected by a strong gradient of thermal forcing along the axis connecting the sub-stellar and anti-stellar points.

Since it is continually heated by the stellar incident flux, the dayside is hotter than the nightside, which is cooled by infrared radiation towards space. If the nightside temperature goes below the condensation temperature of greenhouse gases present in the atmosphere, then the nightside becomes a cold trap for these gases, which condense and form an ice sheet at the planet surface. This triggers atmospheric collapse (e.g. Joshi et al. 1997; Heng & Kopparla 2012; Wordsworth 2015): the greenhouse effect diminishes and the atmosphere cools down, thereby accelerating the condensation process in the cold trap and leading the atmospheric composition, climate, and surface conditions to change radically. It is thus crucial to characterise preliminarily the atmospheric stability of rocky planets against collapse to better constrain the surface conditions of the Earth-like exoplanets found in the habitable zone of their host star.

Kasting et al. (1993), who investigated the habitability conditions around main sequence stars, initially raised the question of atmosphericstability. This question was first addressed using 3D global climate models (GCMs), which are codes that solve the coupled angular momentum, mass conservation, radiative transfer, and energy equations in three dimensions on the surface of a sphere. In their pioneering works, Joshi et al. (1997) and Joshi (2003) examine the case of rocky planets orbiting M-dwarfs and hosting CO₂ /H₂O atmospheres, while Merlis & Schneider (2010) characterise the large-scale circulation patterns of Earth-like planets tidally locked in spin-orbitsynchronisation in various regimes.

A series of studies continued on this path, such as Heng et al. (2011a) and Heng & Vogt (2011), who characterise the atmosphericdynamics of a hypothetical tidally locked Earth and treat the case of the Gliese 581g super-Earth as a scaled-up version of Earth. By considering the effects of moist thermodynamics, Leconte et al. (2013) study the circulation patterns of close-in extrasolar terrestrial planets in the presence of clouds. They thus show evidence of the climate moist bistability, which is the fact that water can be either vaporised – which acts to increase the runaway greenhouse effect – or captured in permanent cold traps, leading the atmosphere to collapse. Carone et al. (2016) examine how surface friction affects the atmosphericstability, while Turbet et al. (2018) computed stability diagrams for TRAPPIST-1 planets assuming N₂ /CO₂ atmosphericmixtures. Similarly, Wordsworth (2015) and Koll & Abbot (2016) treat the case of CO, CO₂, and N₂ -dominated dry atmospheres, and Kang & Wordsworth (2019) investigate the effect of shortwave absorption on collapse.

Meanwhile, with the rise of next-generation space telescopes, such as the James Webb Space Telescope (JWST; Deming et al. 2009), the mechanism of global heat redistribution itself was examined in order to decipher the upcoming thermal phase curves of tidally locked extrasolar rocky planets and hot Jupiters (e.g. Seager & Deming 2009; Cowan & Agol 2011; Selsis et al. 2011; Perez-Becker & Showman 2013; Koll & Abbot 2015; Koll & Komacek 2018; Koll 2019). Although their goals are slightly different, these works are complementary with those dealing with atmospheric collapse.

While they are convenient to treat the complex non-linear physics governing planetary climates, GCMs are far too costly to explore a broad range of the parameter space, with any single climate simulation requiring days or weeks of CPU time. Furthermore, the asymptotic regimes of the steady states reached after convergence may actually be characterised using a relatively small set of nondimensional control parameters (e.g. Koll & Abbot 2015). This motivated the development of simplified analytic theories reducing both the physics complexity and the dimensionality of the problem.

First, 1D analytic models were introduced to characterise the thermal structure of planetary atmospheres. Guillot (2010) derived temperature-pressure (T−P) profiles in the limit of the “double-gray” approximation and pure absorption. Robinson & Catling (2012) considered T− P profiles in the presence of convection. Heng et al. (2012) generalised the work of Guillot (2010) to include isotropic scattering. Parmentier & Guillot (2014) used the “picket fence” model of Chandrasekhar (1935) and Mihalas (1978) to generalise the work of Guillot (2010) by including the simplified non-grey radiative transfer due to spectral lines versus continua. Heng et al. (2014) generalised the work of Heng et al. (2012) to include non-isotropic/anisotropic scattering, and also demonstrated that the governing equations of these T−P profiles and the two-stream solutions have a common origin. Mohandas et al. (2018) generalised the work of Parmentier & Guillot (2014) to include isotropic scattering.

Second, the analytic theory of global heat redistribution was developed incrementally through a diversity of approaches starting from analytical scalings of the atmospheric stability of super-Earths orbiting M-dwarf stars (e.g. Heng & Kopparla 2012). Laying the foundation of the 0-D theory, Wordsworth (2015) proposed a two-layer model that self-consistently includes both radiative transfer and sensible heat exchanges with planetary surface due to convection within the dayside planetary boundary layer (PBL). In this box model, the atmosphere is treated as a globally isothermal and optically thin layer, and the complex features of the three-dimensional general circulation are ignored. Despite these simplifications, the model approximately matches results obtained from GCM simulations for CO₂ -dominated planets. Particularly, it convincingly captures the behaviour of the collapse (or critical) pressure p_crit – that is the minimum surface pressure for atmospheric stability – at high stellar fluxes.

Comparable progress was achieved by Koll & Abbot (2016) regarding the 1D theory. Building on the early work by Robinson & Catling (2012), these authors developed a two-column radiative-convective model that takes the effect of subsidence (i.e. downwelling flows) on the nightside temperature profile into account. This model is thus referred to as a radiative-convective-subsiding (RCS) model in the following. By performing a series of simulations with a GCM, Koll & Abbot (2016) showed that the large-scale atmospheric circulation acts as a global heat engine (a similar work was done by Koll & Komacek 2018, in the case of hot Jupiters), and used this statement to derive scalings for the typical wind speeds parametrising heat fluxes. They thus recovered the stability diagrams of CO₂ -dominated atmospheres previously obtained by Wordsworth (2015) from GCM simulations.

In the present work, we aim to consolidate the analytic theory by building mainly on the developments made by Wordsworth (2015) and Koll & Abbot (2016). Particularly, our goal is to elucidate analytically how the collapse pressure is affected by the non-linear pressure dependence of the atmospheric optical thickness at low stellar fluxes, by dayside convection, and by large-scale advection. Therefore, we adopt the 0-D approach proposed by Wordsworth (2015) and introduce incrementally in the model the effects of longwave and shortwave absorption and scattering, turbulent heat exchanges with the planet surface – which we call “sensible heating” – and advective heat transport by a stellar and anti-stellar cell, which is the regime of slowly rotating planets. Following Koll & Abbot (2016), scalings of wind speeds are derived in the framework of heat engine theory. With every step, intermediate results are compared with those obtained in earlier studies. The full model – which encompasses radiative transfer, dayside convection, and large-scale circulation – is formally written as one single equation, given by Eq. (78), and controlled by a small set of dimensionless parameters. This model is used to compute stability diagrams, and finally compared to a 1D two-column radiative-convective-subsiding-upwelling (RCSU) model that we develop by introducing dayside upwelling flows in the RCS model of Koll & Abbot (2016).

In Sect. 2, we detail the physical setup of the 0-D model and discuss the main assumptions. In Sect. 3, we develop a two-layer radiative model of global heat redistribution including absorption and scattering in the two-stream, dual-band, and grey gas approximations. We then derive analytically the shortwave and longwave flux profiles, radiative transmission functions, and equilibrium temperatures. In Sect. 4, the dayside sensible heating is included in the model. In Sect. 5, we relax Wordsworth’s globally isothermal atmosphere assumption and make the day-night temperature gradient depend on the heat transport due to stellar and anti-stellar circulation. In Sect. 6, the model is used to compute stability diagrams for CO₂ -dominated atmospheresderived from the case treated by Wordsworth (2015). We thus show that the model captures the evolution of the collapse pressure at low stellar fluxes observed in GCM simulations, and the fact that the atmospheric stability increases with sensible heating and decreases as the planet’s size increases. In Sect. 7, we introduce our two-column RCSU model and discuss the role played by the atmospheric structure and, finally in Sect. 8 we summarise the conclusions and future works.

Fig. 1 Diagram of the two-layer radiative-convective-advective model detailed in the present study. The diagram shows the processes involved in dayside-nightside heat redistribution that are taken into account in the model. Radiative net fluxes (F− ≡ F↑− F↓) are upwards in the convention used, meaning that the shortwave net flux F−S is negative.

2 Physical setup and main assumptions

We introduce in this section the main features of the 0-D model developed in the study as well as frequently used notations. This provides a global overview of the mechanisms involved in day-night heat redistribution of terrestrial planets that we take into account (Fig. 1). These mechanisms are detailed in the next sections.

2.1 Tidal locking in 1:1 spin-orbit resonance

We consider the simplified case of a tidally-locked rocky planet, of mass M_p and radius R_p, synchronously moving on a coplanar and circular orbit of period P_⋆. This regime is mainly relevant for planets orbiting small main-sequence stars, such as K and M dwarfs, since the habitable zone of the host star is within the tidal lock radius in this case (Kasting et al. 1993; Edson et al. 2011). This result may be derived from a simple scale analysis.

First, the tidal lock radius is usually estimated for dry planets from the formula (Peale 1977; Kasting et al. 1993; Dobrovolskis 2009; Edson et al. 2011), (1)

Here, M_⋆ is the stellar mass, P_rot;0 the original rotation period of the planet, t_pf the time period from formation, and Q the tidal quality factor accounting for tidal dissipation in the planet’s interior (the smaller Q and the larger tidally dissipated energy; Goldreich & Soter 1966). All quantities are in CGS units (centimetres, grams, seconds) in the formula.

Second, the typical radius of the habitable zone r_HZ can be roughly defined as the star-planet distance at which the black body equilibrium temperature of the planet has a given value T_HZ. The equilibrium temperature is given as a function of the stellar flux received by the planet F_⋆, (2)

the parameter σ_SB = 5.670367 × 10⁻⁸ W m⁻² K⁻⁴ (Mohr et al. 2016) being the Stefan-Boltzmann constant, and the symbol ≡ meaning “defined by”. The above definition of r_HZ thus leads to . The stellar flux is also expressed as a function of the orbital radius a and the stellar bolometric luminosity L_⋆, (3)

By combining the two expressions of F_⋆, it follows that (4)

The dependence of the stellar bolometric luminosity of main-sequence stars on the stellar mass has been quantified by empirical formulae, such as (Barnes et al. 2008) (5)

with , the notations M_⊙ and L_⊙ referring to the Solar mass and luminosity, respectively. A linear regression of on the range − 1 ≤ μ ≤ 0 yields (6)

and thus . By comparing this scaling law to (Eq. (1)), we observe that the tidal lock radius (r_T) increases slower with the stellar mass than the characteristic radius of the habitable zone (r_HZ). This means that the probability to be tidally locked in synchronous rotation for a planet orbiting in the habitable zone decays as the stellar mass increases, which is illustrated for instance by Kasting’s plot (Kasting et al. 1993, Fig. 16).

In the circular coplanar regime, there is no seasonal variation of the local incident stellar flux because of the absence of obliquity and eccentricity. The stellar forcing at the top of the planet’s atmosphere is thus invariant with time. The spin rotation of the planet is defined by the spin angular velocity Ω = n_⋆, where n_⋆ ≡ 2π∕P_⋆ is the orbital frequency. Hence, Ω is determined using the third Keplerian law, (7)

where G designates the universal gravitational constant.

Besides, considering both the expression of the stellar flux received by the planet, given by Eq. (3), and the fact that the stellar bolometric luminosity L_⋆ of main-sequence stars may be written as a function of the stellar mass (see Eq. (5)), we remark that the planet’s angular velocity is actually fully constrained by the stellar mass and flux through the relationship (8)

and thus cannot be taken as a free parameter. Particularly, in simulations that we perform with general circulation models to benchmark the theory, the planet’s rotation rate is specified as a function of M_p, M_⋆, and F_⋆ using the preceding expression with the scaling law of L_⋆ given by Eq. (5).

2.2 Radiative transfer

The planet’s atmosphere is considered in the standard shallow atmosphere framework (Vallis 2006)¹, where the thickness of the fluid layer is assumed to be far smaller than horizontal scales, comparable to the planet’s radius in order of magnitude. In this case, the horizontal propagation of radiation may be neglected and radiation only propagates upwards and downwards. Upwelling fluxes are denoted by F_↑ and downwelling fluxes by F_↓. This sets the bases of the two-stream approximation (e.g. Heng 2017, Sect. 3.1), which is assumed to derive radiative net fluxes F₋ ≡ F_↑− F_↓ in the next section.

In addition to the two-stream approximation, we consider that the frequency spectra of the stellar and planetary radiative fluxes do not overlap, which is the so-called dual-band approximation (Heng 2017, Sect. 4.1). The radiative fluxes emitted by the star and the planet are thus split into two decoupled components, which are refereed to as the “shortwave” and “longwave” fluxes, and subscripted by S and L, respectively.While this simplification holds for rocky planets hosted by Sun-like stars, where the bodies’ surface temperatures are separated by more than one decade in orders of magnitude, it tends to be less adapted to those hosted by ultra-cool dwarf stars, such as TRAPPIST-1, although the stellar temperature (≈ 3000 K) is still far higher than temperate planetary atmospheres (≈300−500 K). In these cases, the stellar and planetary radiative fluxes partly overlap. Cool stars are also those that are more likely to host tidally-locked exoplanets in their habitable zone than Sun-like stars, as stated in Sect. 2.1. We should thus bear in mind the limitations of the dual-band approximation when applying the model to such systems.

Finally, grey gas opacities are assumed, meaning that atmospheric opacities in the longwave and shortwave are wavelength-independent in the model (Heng 2017, Sect. 4.1). In each band, absorption is described by a unique parameter that we call “effective opacity” to emphasise the fact that it does not really correspond to a mean opacity owing to the existing correlation between radiative spectra and opacity lines (see e.g. Wordsworth 2015, Fig. 10). Because of this correlation and of the strong wavelength dependence of opacity lines, this approximation appears as one of the strongest regarding radiative transfers. Particularly, as stated by Leconte et al. (2013) from the comparison of GCM simulations using different approaches, the cooling of the planet’s surface is considerably underestimated when grey opacities are used, meaning that the nightside temperature is overestimated.

This is what motivated the strong effort made to include refined treatments of radiative transfer (e.g. the k-correlated distribution method; Lacis & Oinas 1991) in GCMs used to study the heat redistribution on tidally-locked exoplanets (e.g. Leconte et al. 2013; Wordsworth 2015). In the present work however, we have to admit the grey gas assumption as the price to pay for the simplification of the theoretical analysis. For exhaustive discussions of its limitations, as well as those of the shallow-water, two-stream, and dual-band approximations, one may consult the reference books by Seager (2010), Pierrehumbert (2010), and Heng (2017). We note that, in analytic developments, we follow the conventions and notations employed by Heng (2017) in Chapters 3 and 4.

As shown by Koll & Abbot (2016) through the development of a two-column 1D model, the self-consistent treatment of the coupling between radiative transfer and the atmospheric structure is a major source of complexity, mainly because of the degeneracies affecting boundary conditions. As our goal is to introduce simplified atmospheric dynamics in the analytical theory, we make the choice to favour the atmospheric circulation over the atmospheric structure, noting that the atmospheric dynamics already induces non-negligible mathematical complications. Hence, following the early work by Wordsworth (2015), we opt for a zero-dimensional two-layer model where the atmosphere is considered as isothermal across the vertical direction. However, we relax the well-mixed atmosphere approximation made by Wordsworth (2015), and assume that dayside and nightside atmospheric temperatures, denoted by T_a;d and T_a;n, are not the same in the general case. Similarly, the planet surface is characterised by its dayside and nightside temperatures, T_s;d and T_s;n (Fig. 1).

2.3 Atmospheric heat transport

The four temperatures of the system are coupled together by the mechanisms of heat transport taken into account in the model, which are ofthree types: (i) radiative transfers along the vertical direction (shortwave and longwave absorption, radiation, and scattering); (ii) surface-atmosphere turbulent exchanges in the dayside convective boundary layer; and (iii) day-night horizontal heat transport. Each of these components may defined separately. Owing to the above simplifications, radiative transfers are described by simple analytic solutions, which are derived in Sect. 3. For the sensible heating due to turbulent heat transport, among the different prescriptions existing in literature, we choose to follow that proposed by Koll & Abbot (2016), which treats the atmosphere as a heat engine where the convective flow acts against friction in the surface boundary layer. This part of the theory is detailed in Sect. 4.

The large-scale day-night heat transport turns out to be the most complex effect to include since there is no general theory of atmosphericheat flux to our knowledge. This is due to the fact that the mechanisms contributing to the heat transport cannot be decoupled and are spanning over the three dimensions at planetary scales. Nevertheless, as demonstrated by Koll & Abbot (2015), the involved processes may be disentangled by making use of the Buckingham-Pi theorem (Buckingham 1914), which unravels the nondimensional parameters governing the system.

Among these parameters, we find the nondimensional Rossby deformation length, (9)

the length L_Ro being the Rossby radius of deformation (see e.g. Vallis 2006, Sect. 3.8.2), and c_wave the characteristic speed of gravity waves, defined hereafter. The nondimensional Rossby deformation length governs the circulation regime of the atmosphere. In the case of slow rotators, , meaning that waves can propagate planetwide. Because the effect of rotation is weak, the circulation that develops in this regime is approximately symmetric with respect to the axis connecting the substellar and anti-stellar points, the day-night flow being driven by the balance between advection and pressure-gradient acceleration (e.g. Leconte et al. 2013). In this stellar and anti-stellar circulation, high-altitude winds blow from the dayside to the nightside and low-altitude wind from the nightside to the dayside (e.g. Merlis & Schneider 2010; Leconte et al. 2013).

We illustrate the slow rotators regime by performing simulations with the Open-source General Circulation Model THOR (Mendonça et al. 2016; Deitrick et al. 2019), which solves the general non-hydrostatic equations within a spherical shell using the icosahedral grid. In these calculations, simple grey gas radiative transfers are used and two cases treated by early works are reproduced: a 0.1-bar CO₂ -dominated atmosphere at T_eq = 279 K (Wordsworth 2015, Fig. 2), and a 0.5-bar N₂-dominated atmosphere at T_eq = 400 K (Koll & Abbot 2016, Fig. 4). Figure 2 shows instantaneous snapshots of the zonally-averaged temperatureand vertical wind speeds in both cases. These quantities are plotted as functions of latitude and pressure in the reference frame where the north pole is located at the substellar point and the south pole at the anti-stellar point. We emphasise that this latitude coordinate is different from the usual latitude, which is defined from the planet axis of rotation. As we observe that a steady cycle has been reached in both simulations after 1000 days, snapshots are taken at this date. The plots particularly highlight the day-night temperature gradient (left panels) and the dayside convective cell (right panels), which is characterised by strong upward flows in the substellar region.

When the Rossby deformation radius becomes less than the planet radius (), the circulation regime changes and the atmospheric dynamics starts developing super-rotation, that is strong eastward equatorial jets (Showman & Guillot 2002; Showman & Polvani 2011; Thrastarson & Cho 2010; Heng & Vogt 2011; Heng et al. 2011b; Leconte et al. 2013; Mendonça 2020). As shown by Showman & Polvani (2011), the emergence of super-rotation is due to the formation of standing, planetary-scale equatorial Rossby and Kelvin waves (i.e. waves restored by the Coriolis acceleration; see e.g. Lee & Saio 1997), which feed zonal flows by continually pumping angular momentum from the mid-latitudes towards the equator. The latitude width of the equatorial band where super-rotating flows are accelerated through this mechanism scales as the equatorial Rossby radius (L_Ro), and thus decays as the planet angular velocity increases. As the Rossby and Kelvin waves responsible for the super-rotation are caused by the latitudinal variations in radiative heating, the strength of equatorial jets increases with the day-night temperature difference (Showman & Polvani 2011).

The speed of gravity waves that intervenes in the non-dimensional Rossby deformation length depends on the atmospheric structure. Denoting by g the surface gravity, R_GP the ideal gas constant, and the mean molecular weight of the atmosphere, we first introduce the specific gas constant and the pressure scale height, given as a function of temperature T by (10)

Gravity waves are restored by the Archimedean force associated with the fluid buoyancy (Gerkema & Zimmerman 2008). The strength of this force is quantified by the Brunt-Väisälä frequency N, which, in a dry stably stratified atmosphere, is given by (11)

where C_p designates the heat capacity per unit mass of the gas at constant pressure, and z the altitude. In terms of H and N, the typical speed of gravity waves is simply written as c_wave ~ HN (e.g. Koll & Abbot 2015; Leconte et al. 2013). Thus, substituting these parameters by the above expressions (Eqs. (10) and (11)) in Eq. (9) and assuming an isothermal atmosphere yields (12)

In the present work, we only consider the heat advected by a day-night stellar and anti-stellar circulation, which corresponds to the asymptotic regime of slowly rotating planets (). We do not include other mechanisms of day-night heat transport, such as that due to the propagation of gravity waves, which was already studied in early works (Koll & Abbot 2015, 2016). Thus, in Sect. 5, a rough approximation of the day-night heat flow is derived from a scale analysis in the heat engine framework (Koll & Abbot 2016) using for wind speeds the prescription given by Koll & Komacek (2018) in the case where circulation is balancing Rayleigh drag. As predicted by dimensional analyses (e.g. Leconte et al. 2013, Sect. 3.2.1), this heat flow is scaled by the ratio t_rad ∕t_adv between the radiative and advective timescales (e.g. Showman & Guillot 2002) (13)

where p_s designates the atmospheric surface pressure and V_adv the typical speed of stellar and anti-stellar advective flows. This accounts for the fact that the heat transported by weak flows (t_rad ∕t_adv ≪ 1) is radiated towards space before being advected to the nightside. Conversely, strong flows (t_rad ∕t_adv ≫ 1) efficiently transport heat since the advected fluid parcels reach the nightside before being radiatively cooled. In the present work, the scaling law chosen for the advected heat flux is different from that used to describe thermal exchanges with the surface on the dayside in order to emphasise the distinction between the two involved mechanisms in the general case. These mechanisms could nevertheless be linked to each other by assuming the same scalings for the associated heat fluxes.

Finally, we adopt dry thermodynamics as a first convenient step to set the basis of the model, ignoring thereby the role played by an ocean or moisture on the day-night heat redistribution. We note however that this role may be non-negligible. Edson et al. (2011) showed for instance that the presence of a slab ocean considerably affects the stability of the planet’s atmosphere with respect to collapse, by inducing an additional heat transport by latent heat flux in the atmosphere and by heat diffusion in the ocean itself. Water vapour also leads to the formation of clouds, which affects the albedo of the atmosphere as well as its thermal structure and general circulation (Leconte et al. 2013). While not taken into account in the present work, these aspects may be included in the model in future studies.

Fig. 2 Instantaneous snapshots of zonally-averaged temperature (left panels) and vertical wind speed distributions (right panels) as functions of latitude and pressure with north pole at substellar point. Top: 0.1-bar CO2 -dominated atmosphere at Teq = 279 K (Wordsworth 2015, Fig. 2). Bottom: 0.5 bar N2-dominated atmosphere at Teq = 400 K (Koll & Abbot 2016, Fig. 4). Latitude is given in degrees (horizontal axis). This coordinate does not correspond to the usual latitude but to a latitude coordinate in the system of spherical coordinates where the north pole is the substellar point and the south pole the anti-stellar point. The associated colatitude is the stellar zenithal angle in this system of coordinates. Pressure is given in bars (vertical axis), temperature in K, and vertical velocity in m s−1. The white area at the substellar point indicates the region where surface pressure is lower than the spherically averaged surface pressure due to the updraft. The dotted red line indicates the transition between rising and subsiding flows, where the vertical velocity is zero.

3 A two-layer grey radiative model

Radiative transfer is the first building block of the model. We establish analytically in this section the vertical profiles of net radiative fluxes in the short- and longwave for the isothermal atmosphere. These profiles provide the integrated transfer functions that are used to derive the atmospheric and surface power budget equations.

3.1 Two-stream analytic solution

In the two-stream, dual-band, and gray-gas approximations (see Sect. 2.2), the longwave and shortwave fluxes decouple and only propagate upwards or downwards. As it is governed by the same set of equations, the shortwave case is readily deduced from the longwave case that we treat hereafter. The longwave total and net fluxes are defined as F_+L ≡ F_↑L + F_↓L and F_−L ≡ F_↑L − F_↓L, respectively. They are governed by the two-stream Schwarzschild equations, which may be written as the system of first order partial differential equations (Heng 2017)

where ∂_X designates the partial derivative with respect to X, τ_L the optical depth of the atmosphere in the longwave, and β_L0 the scattering parameter. The optical depth characterises the optical thickness of the atmosphere and varies between 0 (top of the atmosphere) and τ_s;L (planet’s surface). The scattering parameter indicates the fraction of absorbed flux with respect to the scattered component and takes its values between 0 and 1 (β_{L 0} = 1 corresponds to pure absorption and β_L0 = 0 to pure scattering; see e.g. Heng 2017, Sect. 4.3). Eliminating the total flux in the preceding system yields the second order ordinary differential equation (16)

the notation designating the total derivative of some quantity X.

As the isothermal approximation is assumed, the atmospheric temperature profile is approximated by the constant T_a. In this particular case, the right-hand member of Eq. (16) vanishes, which leads to the analytic expressions of the total and net fluxes,

and of the upwelling and downwelling fluxes,

where and are the two integration constants of the solution, and and the couplingcoefficients defined by (21)

The parameters characterise the coupling of radiative fluxes due to scattering. In the case of pure absorption (β_{L 0} = 1), there is no coupling ( and ). Conversely, in the case of pure scattering (β_L0 = 0), the couplingis strong (). Similar expressions may be derived for the shortwave band by assuming T_a = 0 in Eqs. (17–20). The shortwave fluxes (F_+S, F_−S, F_↑S, and F_↓S) are parametrised by the shortwave optical depth τ_S (such that 0 ≤ τ_S ≤ τ_s;S, τ_s;S being the shortwave optical depth at the planet’s surface), scattering parameter β_{S 0}, coupling coefficients , and integration constants and .

The four integration constants of the solution are determined by assuming two boundary conditions for each wavelength domain. At the top of the atmosphere, the shortwave downwelling flux is the stellar incident flux F_i and the longwave downwelling flux is zero because of the absence of thermal source at the infinity, which is written as F_↓S= F_i and F_↓L = 0. At the planet surface, a fraction of the shortwave downwelling flux is reflected while the remaining part is absorbed by the ground and re-emitted in the longwave. Denoting by A_s the albedo of the planet’s surface in the shortwave, and F_s the flux emitted by the planet’s surface in the longwave, these conditions are expressed as F_↑S= A_sF_↓S and F_↑L = F_s. We thus end up with the net fluxes in the short- and longwave, (22)

and (23)

where we have introduced the atmospheric and surface fluxes, F_a and F_s, and the vertically-integrated transmission functions (24)

The vertical profiles of the short- and longwave net fluxes given by Eqs. (22) and (23) are plotted in Fig. 3 for various atmospheric optical depths. One may observe here that increasing the opacity tends to confine the absorption region in the upper atmosphere in the shortwave band, while it tends to neutralise the net flux within the atmosphere in the longwave.

Fig. 3 Net radiative fluxes in the short- and longwave (horizontal axis) as functions of optical depths (vertical axis) in case of isothermal atmosphere. Fluxes profiles are computed from the analytic formulae given by Eqs. (22) and (23) and normalised by the black body equilibrium flux (the equilibrium temperature Teq being givenby Eq. (2)). Optical depths are normalised by their value at the planet’s surface, namely τs;S for the shortwave and τs;L for the longwave. Shortwave (blue lines) and longwave (red lines) net fluxes are plotted for various optical thicknesses, the optical depths at the planet’s surface taking the values τs;S = 0.1, 1, 10 (from light to dark blue lines) and τs;L = 0.1, 1, 10 (from light to dark red lines). Pure absorption is assumed (βL0 = βS0 = 1). The incident, atmospheric, and surface fluxes are arbitrarily set to Fi = 4Feq, Fa = 0.7Feq, and Fs = 2Feq, respectively.

3.2 Surface and atmospheric energy budgets

The derivedanalytic solution is now used to establish the local energy budgets of the planet’s surface and atmosphere at thermodynamical equilibrium. By considering the upper (τ_L = τ_S = 0) and lower (τ_L = τ_s;L and τ_S = τ_s;S) boundaries, we obtain the fluxes radiated towards space at the top of the atmosphere and towards the atmosphere at the planet’s surface, respectively. These fluxes are expressed as (25) (26)

and (29)

which brings out the coefficient parameterising the surface and atmospheric energy balances,

with (similarly, in the short wavelength domain, we may define ). These coefficients are expressed as (34) (35) (36) (37) (38)

and (39)

Thus, introducing the heat power per unit area due to surface turbulent exchanges on the dayside (F_sen) and day-night advection (F_adv) – specified in Sects. 4 and 5, respectively – the local energy balances of the planet’s surface and atmosphere are formulated as

We note that Eqs. (40) and (41) are the counterparts of Eqs. (14) and (15) in Wordsworth (2015), respectively. Following Pierrehumbert (2011) and Wordsworth (2015), the notation F_sen refers to “sensible heat flux”, which designates an energy flux corresponding to a change of temperature of the fluid².

As illustrated by Fig. 1, the dayside and nightside hemispheres of the planet are characterised by different surface and atmospheric temperatures (T_s;d, T_a;d, T_s;n, and T_a;n). This leads to average the preceding energy budgets over each of the two hemispheres. Hence, we introduce the dayside and nightside hemisphere-averaged surface fluxes, denoted by B_s;d and B_s;n, the blackbody fluxes radiated by the atmosphere, B_a;d and B_a;n, and the averaged heat flows, these quantities being defined by (see Wordsworth 2015, Eqs. (16) and (17)) (44)

Besides, the incident stellar flux is given as a function of the stellar zenithal angle ϕ by (45)

Thus, averaged over the dayside and nightside hemispheres, the system of Eqs. (40)–(43) becomes

This set of equations is the counterpart of Eqs. (21)–(23) in Wordsworth (2015) extended outside of the optically thin limit and with an additional energy flux between the dayside and nightside hemispheres ().

3.3 Purely radiative case

Before investigating the role played by sensible heating and advection, it is worth examining the purely radiative case. For simplification, we place ourselves in the regime treated by Wordsworth (2015), where day-night heat transfers are strong enough to homogenise the atmospheric temperature (T_a;d = T_a;n = T_a and B_a;d = B_a;n = B_a). In this framework, the sensible heating is neglected (), and the advected heat flow does not intervene in the equilibrium, meaning that the dayside and nightside atmospheric budgets may be combined into a globally averaged atmospheric budget, (50)

Solving this equation together with Eqs. (46) and (48) yields the surface and atmospheric equilibrium temperatures (51) (52) (53)

which capture the non-linear dependence of radiative transfer on the atmospheric optical thickness in the general case, and highlight asymptotic regimes.

Assuming pure absorption in the longwave (β_L0 = 1) and transparency in the shortwave (τ_s;S = 0), we recover the behaviour derived by Wordsworth (2015) in the optically thin limit (). Conversely, as the optical depth increases, both the atmosphere and surface temperatures converge towards the black body temperature, as showed by Table 1. We note that the surface temperature of rocky planets in the optically thick limit are expected to be higher than that predicted by the model, because the atmospheric structure cannot be approximated by an isothermal temperature profile in this limit. For instance, Venus’ troposphere is characterised by a strong temperature gradient (e.g. Seiff et al. 1985), which is far closer to the adiabatic profile (N = 0) than to the isothermal profile. Consequently, Venus’s mean surface temperature is T_s ≈ 740 K whereas its blackbody equilibrium temperature in the absence of albedo is T_eq = 327.3 K (we have taken the value F_⋆ = 2601.3 W m⁻² for the Solar bolometric flux).

The effect that absorption and scattering have on the steady state is explored in Fig. 4, where the atmospheric and surface temperatures are plotted against the longwave opacity for various scattering parameters and shortwave opacities. We first consider the quasi-transparent limit in the shortwave (τ_s;S = 10⁻⁴), which corresponds to τ_s;S ≪ τ_s;L. As the longwave optical thickness increases, the equilibrium state switches from the optically thin regime, where , to the optically thick one, where the three temperatures reach a plateau. The interplay between short- and longwave optical depths defines another dimension in the parameter space. When τ_s;L ≲ τ_s;S, the atmosphere has to increase its temperature to evacuate the power absorbed in the shortwave. As a consequence, the atmospheric temperatureincreases as the optical thickness in the longwave decays, which is the behaviour observed in Fig. 4 (right panels).

Scattering in the shortwave and in the longwave have opposite effects on temperatures. When scattering in the shortwave is present (β_{S 0} = 0.5, middle panels), it induces an anti-greenhouse cooling. A fraction of the incident stellar flux is scattered back to space, meaning that the planet absorbs less energy than in the case of pure absorption. As a consequence, the atmosphere and surface temperatures are smaller. Conversely, the presence of scattering in the longwave (β_{L 0} = 0.5, bottom panels) generates an effect known as scattering greenhouse effect (e.g. Heng 2017, Sect. 4.7, p. 71). Thermal emissions of the surface and atmosphere towards space are scattered back to the planet, leading to the observed increase of the atmospheric, dayside and nightside surface temperatures with respect to the case of pure absorption.

4 Inclusion of sensible heating

As a second step, the dayside sensible heating is included in the model. While we still assume a uniform atmospheric temperature (T_a;d = T_a;n = T_a), the surface turbulent heat flux is now expressed as a function of the dayside surface and atmospheric temperatures. We introduce first the expression of turbulent heat exchanges between the surface and the atmosphere within the planetary boundary layer (PBL), and second the scaling law of the horizontal velocity parametrising this flux, following Koll & Abbot (2016), where the atmospheric circulation is treated as a heat engine (Fig. 5).

4.1 Scaling of turbulent heat exchanges

The turbulent heat exchanges due to convection in the dayside PBL may be quantified in the framework of the mixing length theory (Holton 1973, 5th edition, Sect. 8.3.5, p. 268), as detailed in Appendix A. The hemisphere-averaged sensible heat flux is written as (e.g. Pierrehumbert 2011, Eq. (6.11), p. 396) (54)

where is the surface density, C_D the bulk drag coefficient accounting for the strength of friction, and V_sen the typical horizontal wind speed quantifying the strength of the convective cell that develops in the vicinity of the substellar point. If the ground is warmer than air, heat is carried away from the ground () at a rate proportional to the temperature difference. Conversely, if the ground is cooler than the air, the atmosphere warms the ground (). We note that the convective cell vanishes in case of stable stratification, which implies V_sen = 0 and . Thus, we ignore the sensible flux on the nightside. On the dayside however, although radiative transfer is described in the isothermal atmosphere approximation for simplification, we assume that the planetary boundary layer is convective, so that in the general case.

The bulk drag coefficient weakly depends on the surface and boundary layer properties. In the case of a neutrally buoyant boundary layer, following the empirical scaling law proposed by von Kármán (1930) for the mixing length (Appendix A), it is expressed as a function of the surface roughness length z₀, the thickness of the turbulent boundary layer z_h, and the von Kármán (1930) constant , as (e.g. Esau 2004, Eq. (7)) (55)

The bulk drag coefficient is of order C_D ~ 10⁻³ (a typical value of C_D over oceans is C_D ≈ 1.5 × 10⁻³; e.g. Holton 1973, Sect. 8.3.1), although it may be larger over rough ground. In the present work, this parameter is set to C_D = 3.4 × 10⁻³, which is the value used by Wordsworth (2015) and corresponding to z₀ ~ 10⁻² m and z_h ~ 10 m.

The characteristic horizontal wind speed varies with the surface and atmospheric temperatures. To take this dependence into account, Wordsworth (2015) derived a relationship between V_sen, T_a, and T_s;d from the thermodynamic equation by assuming the Weak Temperature Gradient (WTG) approximation, that is considering that winds are driven by small horizontal temperature gradients. This equation (Eq. (43) in the article) was combined with Eq. (54) to close the system of equations derived from hemisphere-averaged energy budgets. Later, by performing a series of numerical simulations using a GCM, Koll & Abbot (2016) showed that the atmospheric circulations of rocky planets resemble heat engines. Particularly, they established that the heat engine theory better fit horizontal wind speeds than the relationship derived by Wordsworth (2015) in the explored region of the parameter space. We thus choose to follow here the prescription proposed by Koll & Abbot (2016).

In this framework, the dayside convection is idealised as a single overturning cell between the substellar point and cooler regions. This cell is driven by the surface-atmosphere temperature difference, and the heat engine is fed by the absorption of radiative fluxes. The atmosphere absorbs heat near the dayside surface at the temperature T_s;d (hot reservoir) and emits it to space at the temperature T_a;d (cold reservoir) in the upper regions of the atmosphere. During the cycle, a fluid parcel works against friction in the boundary layer, which is the place where the energy is dissipated. This work is formulated as (Bister & Emanuel 1998; Koll & Abbot 2016). In the case of an isentropic cycle (i.e. composed of adiabatic reversible processes solely), W is related to the amount of power per unit area available to drive atmospheric motion Q_in by Carnot’s theorem (e.g. Schroeder 2000) (56)

the factor being the atmosphere’s thermodynamic efficiency. The typical horizontal wind speed of the convective cell is thus expressed as (57)

where we have introduced the efficiency coefficient ε_sen accounting for the non-isentropic nature of the thermodynamic cycle. For an isentropic cycle, ε_sen = 1, which corresponds to the idealised Carnot’s heat engine. In reality, the thermodynamic cycle is not isentropic because additional dissipative processes, such as diffusion, induce an irreversible production of entropy, which decreases the efficiency of the heat engine. As a consequence, the effective amount of energy available to drive atmospheric motion is smaller than its theoretical estimate, and ε_sen < 1. Koll & Abbot (2016) hence noted that the typical wind speeds given by their GCM simulations could be twice smaller than those predicted by the theory, that is ε_sen = 1∕2.

The power per unit area available to drive atmospheric motion (Q_in) has to be specified from an ad hoc prescription. As it was benchmarked against GCM simulations, we follow that proposed by Koll & Abbot (2016), which is formulated as (58)

with , the black body equilibrium temperature T_eq being defined in Eq. (2).

The first factor of the expression given by Eq. (58) is just the hemisphere-averaged power received by the planet on the dayside. The second factor is the vertically-integrated atmospheric transmission function in the shortwave. The third factor is a function of optical depth scaling the fraction of flux emitted by the surface in the longwave that is absorbed by the atmosphere. Assuming τ_s;L ≪ 1 leads to 1 − e^−τ_s;L ≈ τ_s;L and we recover the usual transfer function in the optically thin limit (Pierrehumbert 2011).

Fig. 4 Atmospheric, dayside and nightside surface temperatures in pure radiative regime. Temperatures are plotted as functions of the logarithm of τs;L in the case of pure absorption (top panels), scattering in the shortwave (middle panels) and scattering in the longwave (bottom panels). Left: quasi-transparence in the shortwave (τs;S =10−4). Middle: small shortwave opacity (τs;S = 10−2). Large shortwave opacity (τs;S = 100). The atmospheric temperature (solid black line), dayside surface temperature (dashed red line), nightside surface temperature (dashed blue line), and black-body equilibrium temperature (dotted grey line) are plotted using Eqs. (51)–(53) and (2), respectively. Parameters values: As = 0.2, and F⋆ = 1366 W m−2.

Fig. 5 Diagram of the dayside atmospheric heat engine. The heat engine is driven by dayside heating and cooling to space. A fluid parcel works against friction in the boundary layer, which limits the speed of the flow.

4.2 Role played by sensible heating

The formulation of the sensible heat flux given by Eqs. (54) and (57) closes the system formed by energy budget equations and allows us to determine the state of equilibrium of the system. It is convenient here to adopt the formalism of Wordsworth (2015) since we follow the same approach as this early work. Thus, introducing the normalised temperature, radiative flux, and sensible heat flux, (59)

we transform Eqs. (46), (50), and (54) into

where the dimensionless parameter L_sen is defined as (63)

The parameter L_sen controls the intensity of sensible heating with respect to radiative heating on the dayside. The case L_sen = 0 correspondsto pure radiative equilibrium (no turbulent heat exchanges), while L_sen ≫ 1 indicates a very strong turbulent mixing within the boundary layer, which leads to T_a = T_s;d. The system of Eqs. (60)–(62) is finally reduced to one single equation. By combining Eqs. (60) and (61) we first express the normalised stellar flux as a function of the normalised temperature, (64)

Substituted in Eqs. (61) and (62), this expression of yields (65)

which is the counterpart of Eq. (45) in Wordsworth (2015) with the heat engine approach proposed by Koll & Abbot (2016). This equation can be solved by using a combination of the dichotomy and secant methods (e.g. Press et al. 2007, p. 449). The obtained is then substituted in Eq. (64) and Eqs. (60)–(62) to compute successively , , the atmospheric temperature, and the dayside and nightside surface temperatures. Numerical calculations are performed using TRIP (Gastineau & Laskar 2011).

Figure 6 shows the evolution of the atmospheric temperature (T_a), and the dayside (T_s;d) and nightside (T_s;n) surface temperatures with the atmospheric opacity in the longwave – quantified by τ_s;L – and the control parameter of sensible heating L_sen. In this example, pure absorption is assumed (β_L0 = β_S0 = 1), the atmosphere is quasi-transparent in the shortwave (τ_s;S = 10⁻⁴), and parameters are set to A_s = 0.2 and F_⋆ = 1366 W m⁻², which is the stellar flux received by the Earth. The value of L_sen in the typical case of a 1-bar CO₂-dominated atmosphere (green triangle) is estimated by taking p_s = 1 bar, R_s = 189 J kg⁻¹ K⁻¹, C_p = 650 J kg⁻¹ K⁻¹, C_D = 3.4 × 10⁻³, and ε_sen = 1∕2.

The figure highlights the role played by sensible heating. As L_sen increases, the thermodynamic state of equilibrium switches from the purely radiative regime (L_sen ≪ 1) to a strongly convective regime (L_sen ≫ 1), where atmospheric and dayside surface temperatures are homogenised by sensible exchanges. Both in the optically thin (τ_s;L ≪ 1) and thick (τ_s;L ≫ 1) limits, we retrieve for a similar evolution with L_sen as that shown by Fig. 7 of Wordsworth (2015), which indicates that the two different prescriptions used to estimate the typical horizontal wind speed (heat engine and scaling analysis using the thermodynamic equation) give comparable results. As expected, increasing the atmospheric opacity tends to accentuate the greenhouse effect, and thus to increase the nightside surface temperature.

Fig. 6 Atmospheric, dayside and nightside surface temperatures as functions of logarithm of control parameter Lsen. Left: τs;L = 0.01 (quasi-transparent). Middle: τs;L = 0.1 (small optical thickness). Right: τs;L = 1 (Earth-like optical thickness). In all cases, the atmosphere is quasi-transparent in the shortwave (τs;S =10−4), pure absorption is assumed (βL0 = βS0 = 1), As = 0.2, and F⋆ = 1366 W m−2. The atmospheric temperature is designated by the solid black line, the dayside surface temperature by the dashed red line, the nightside surface temperature by the dotted blue line, and the black body equilibrium temperature (Eq. (2)) by the dotted grey line. The value of Lsen for a 1-bar CO2-dominated atmosphere is indicated by the green triangle.

5 Inclusion of large-scale advection

The last physical ingredient that has to be introduced in the model is a mechanism responsible for heat transport from the dayside to the nightside. As discussed in Sect. 2.3, heat can be advected by large-scale atmospheric flows or transported by gravity waves (e.g. Koll & Abbot 2015, 2016), which propagate in stably stratified fluid layers (Gerkema & Zimmerman 2008). In the present work, we assume that heat transport is due to atmospheric circulation solely and we focus on the regime of slow rotators (), where mean flows are symmetric with respect to the axis connecting the stellar and anti-stellar points. From now on, we relax the uniform temperature approximation and consider that T_a;d≠T_a;n in the general case. First, we derive a scaling of the heat flux associated with advection F_adv as a function of the dayside and nightside atmospheric temperatures by using the heat engine theory (Fig. 7). Second, we use this expression to derive and solve the equilibrium equation in the general case.

Fig. 7 Diagram of large-scale stellar and anti-stellar atmospheric circulation as planetary heat engine. The heat engine is driven by dayside heating and nightside cooling to space. In this study, the atmospheric circulation is balancing Rayleigh drag, which causes frictional energy dissipation.

5.1 Scaling of the advected heat flux

To quantify the atmospheric heat transport, we proceed to a standard scale analysis (e.g. Pierrehumbert 2011, Sect. 9.2.3, p. 610). The heat per unit mass contained in a fluid parcel starting from the dayside is C_p T_a;d. When the parcel moves to the nightside, it decays to C_pT_a;n. The energy released by the parcel thus corresponds to . Moreover, the mass flow that goes across the annulus separating the day- and nightside is proportional to the air column mass p_s ∕g, to the size of the annulus 2πR_p, and to a typical windspeed V_adv. According to these considerations, the hemisphere-averaged heat flux due to atmospheric circulation may be written as (66)

where ε_circ is an efficiency coefficient related to the three-dimensional geometry of the flow. The wind speed V_adv has now to be specified as a function of the atmospheric temperatures and of the system parameters. In order to clearly separate the mechanisms responsible for day-night heat transport on the one hand, and for sensible heat exchanges with the surface on the other hand, we choose to adopt a scaling law for the advection velocity different from that used to parametrise the convective cell (Eq. (57)). We note however that the two mechanisms could be linked to each other by adopting the same scaling for both velocities.

Koll & Komacek (2018) showed that the large-scale atmospheric circulation of hot Jupiters can be modelled as planetary heat engines and derived the typical wind speeds of mean flows in this framework. The reasoning applied to hot Jupiters in their work holds for rocky planets, and is very similar to that followed in Sect. 4.1 to derive the expression of the sensible heating flux used in the model (Eq. (62)). Here, the hot and cold reservoirs of the heat engine are the dayside and nightside parts of the atmosphere, and the temperature gradient driving the flow is thus horizontal instead of being vertical. In steady state, the work generated by differential heating and cooling is balanced by dissipation of kinetic energy, which determines the typical wind speed of the flow. In the present study, we assume that the atmospheric circulation is balancing Rayleigh drag in the terminator region. As a consequence, wind speeds scale as (e.g. Koll & Komacek 2018, Eq. (11)) (67)

where designates Carnot’s thermodynamic efficiency for the large-scale circulation, ε_adv the efficiency coefficient associated with the non-isentropic nature of the cycle, Q_in the theoretical amount of power available to drive atmospheric motion, given by Eq. (58), and t_drag the drag timescale, that is the timescale over which winds are linearly damped by Rayleigh drag. Substituting V_adv by Eq. (67) in Eq. (66), introducing the total efficiency factor ε_eff ≡ ε_circε_adv, and extracting the dependence on the dayside and nightside atmospheric temperatures yields (68)

where we retrieve the scaling given by earlier studies (e.g. Leconte et al. 2013). We remind ourselves here that t_rad and t_adv are the radiative and advective timescales defined by Eq. (13). The ratio t_rad ∕t_adv quantifies the ability of the stellar and anti-stellar circulation to transport heat from the dayside to the nightside. The case t_rad ∕t_adv = +∞ corresponds to the regime treated by Wordsworth (2015) and Sect. 4 of the present work, where advection is so strong that the atmosphere is horizontally well mixed and its temperature homogenised. Conversely, if t_rad ∕t_adv ≪ 1, advection is not efficient in warming the nightside, which leads to strong day-night temperature gradients. We note that the drag timescale (t_drag) and the efficiency factor (ε_eff) appearing in Eq. (68) are unknown a priori. Both these parameters and the scaling law used for should eventually be constrained with the help of more sophisticated models solving the coupled momentum and thermodynamic equations, such as GCMs typically.

5.2 Steady state equation in the general case

We can now write the system of equations describing the steady state in the general case, where radiative transfer, sensible heating and large-scale circulation are taken into account (Fig. 1). In accordance with the formalism employed in Sect. 4.2, we introduce the normalised temperature and heat flux, (69)

Equations (46), (50), and (49) thus become

where the normalised sensible heating and advection fluxes are expressed as

In the preceding equation, we have introduced the non-dimensional parameter (75)

which controls the intensity of heating due to atmospheric circulation with respect to radiative cooling. This parameter is complementary with L_sen, defined by Eq. (63), which controls the intensity of sensible heating with respect to radiative cooling on dayside.

The system of Eqs. (70)–(74) may be reduced to one single equation by proceeding similarly as in Sect. 4.2. First, is expressed as a function of and by substituting by Eq. (74) in Eq. (72), (76)

Second,Eqs. (70) and (71) are combined to eliminate , and the resulting equation is combined with Eq. (74) to eliminate . This yields the expression of as a function of , (77)

Third, is substituted by Eq. (77) in Eq. (76), and both and are substituted by their expressions as functions of in the normalised fluxes (Eqs. (73) and (74)), which allows us to reduce the system of Eqs. (70)–(74) to a single equation of . By using Eq. (70), we get (78)

This equation can be solved using a combination of the dichotomy and secant methods (Press et al. 2007), as previously done for Eq. (65). However, we note that the denominator of Eq. (77) is negative if is less than a minimum value , meaning that is not defined in this range. This minimum value simply results from the fact that the advected heat flux cannot be greater than the heat absorbed on the dayside hemisphere. As the heat transport by large-scale advection increases as decays (see Eq. (74)), a maximum of corresponds to a minimum of , which is . We note that if the efficiency of advection is low (L_adv ≪ 1), while if the efficiency is high (L_adv ≫ 1). As a consequence, one has to preliminarily calculate by finding the zero of the denominator of in Eq. (77). As a second step, the solution of Eq. (78) is sought within the interval considering the fact that the nightside atmospheric temperature cannot be greater than the dayside atmospheric temperature in this model. Although there is no evidence for the existence and unicity of solutions from an analytical point of view, we verify them numerically in the studied range of the parameter space (see Appendix B).

However, as L_adv → +∞, the dayside and nightside atmospheric temperatures tend to join together, meaning that very small temperature variations lead to huge variations of the flux associated with atmospheric circulation (). In other words, the day-night atmospheric temperature difference do not contain information any more in the strong circulation limit, and one should drop Eqs. (72) and (74), which corresponds to the case treated in Sect. 4.2 and by Wordsworth (2015). This degeneracy has repercussions on the numerical solution of Eq. (78), which cannot be performed beyond a certain upper bound of L_adv ≫ 1 (typically L_adv ~ 10²−10³). We also note that the logarithm of the distance to is a better appropriate coordinate than itself in the root-finding procedure because of the strong variation of the function given by Eq. (78) in the vicinity of (see Fig. B.1).

Figure 8 shows the evolution of the normalised temperature (Eq. (59)) as a function of the parameters controlling sensible heating (L_sen) and large-scale heat transport (L_adv) in the optically thin case of Fig. 6 (τ_s;L = 0.1, middle panels). We observe here the transition between the asymptotic limits of weak (L_adv ≪ 1) and strong (L_adv ≫ 1) atmosphericcirculations. In the regime of strong advection, the solution resembles that obtained by Wordsworth (2015) (see Fig. 7, top panel). Conversely, while the circulation weakens, a larger amount of heat remains on the planet’s dayside. As a consequence, the temperature difference between the surface and the atmosphere decays.

In Fig. 9, the four temperatures of the system in steady state are plotted as functions of L_sen for L_adv = 0.01, 1, 100 in the optically thin (τ_s;L = 0.1) and thick (τ_s;L = 1.0) regimes. This figure shows that the nightside surface temperature strongly depends on L_adv, while it is not very sensitive to L_sen. As discussed above, increasing the efficiency of heat transport by the atmospheric circulation leads the nightside atmospheric and surface temperature to increase. In the regime of strong circulation (L_adv ≫ 1), the atmosphericnightside and dayside temperatures join together and we recover the results obtained in Fig. 6, where L_adv = +∞. The dependence of the nightside temperature on the atmospheric circulation is increased by greenhouse effect. In the optically thick case the temperature evolution between the cases L_adv = 0.01 and L_adv = 100 is larger than in the optically thick one.

Fig. 8 Normalised temperature as function of logarithm of parameter Lsen for various values of Ladv. The parameter Ladv, which controls advective heat transport by the stellar and anti-stellar circulation, varies in the range (from dark blue to green). The atmosphere is quasi-transparent in the shortwave (τs;S =10−4), pure absorption is assumed (βL0 = βS0 = 1), As = 0.2, and F⋆ = 1366 W m−2, similarly as in Fig. 6. Here, the atmosphere is optically thin in the longwave: τs;L = 0.1, which corresponds to the middle panels of Fig. 6.

6 Atmospheric stability

The analytical developments made in the preceding sections allow us to characterise the stability of the atmosphere, that is its ability to conserve its greenhouse gases and, more generally, to remain unchanged in composition. In the following, we use the model (Eq. (78)) in parallel with that proposed by Wordsworth (2015, Eq. (45) of the article) to quantify the effects of the different processes participating to global heat redistribution on atmospheric stability.

The collapse is triggered by the condensation of one of the gases present in the atmospheric mixture, which occurs from the moment that the temperature of the air goes below the condensation temperature of the gas. The gas then condensates and forms an ice sheet at the planet’s surface in the cold trap. This modifies the radiative properties of the atmosphere, as well as the atmospheric circulation, which are both coupled to the atmospheric gas mixture. In the standard case of greenhouse gases, collapse acts as a positive feedback by leading the atmosphere to cool down, and favouring thereby the condensation of other gases.

In the simplified approach of the present work, the complex dynamics of the collapse is not treated and the cold trap is assumed to be the whole nightside surface. We thus consider that, because of atmospheric circulation, a fluid parcel will end up cooling down at the nightside surface temperature, and we do not investigate the timescale of the process consequently. Under these hypotheses, the atmospheric stability is simply determined by the ratio between the nightside surface temperature and the lowest condensation temperature among the greenhouse gases present in the mixture, following early analytical studies (Wordsworth 2015; Koll & Abbot 2016).

In order to better understand the role played by the different involved processes, we use as a reference case the case studied by Wordsworth (2015), which is an Earth-sized planet hosting a CO₂-dominated atmosphere. In this case, the atmosphere is assumed to be stable with respect to atmospheric collapse if , the condensation temperature of CO₂ in K being given, below the triple point (p < 5.18 × 10⁵ Pa), by (Fanale et al. 1982; Wordsworth et al. 2010; Wordsworth 2015) (79)

and, beyond the triple point (p ≥ 5.18 × 10⁵ Pa), by (80)

To specify the relation between the short- and longwave optical thicknesses and surface pressure, we use a linear law of the form (Wordsworth 2015, Eq. (12)) (81)

where κ_S and κ_L are the supposed constant shortwave and longwave opacities, respectively. The longwave opacity is set to κ_L = 5 × 10⁻⁵ m² kg⁻¹, which is the value used by Wordsworth (2015) in GCM simulations. This value corresponds to τ_s;L ≈ 1 for a 1-bar atmosphere. The shortwave opacity is set to κ_S = 10⁻⁹ m² kg⁻¹, that is a value corresponding to the quasi-transparent regime (κ_S ≪ κ_L). We note that pure absorption is assumed in the reference case (β_S0 = β_L0 = 1). Finally, the parameters characterising heat engines are specified from the observations made by Koll & Abbot (2016) and Koll & Komacek (2018) using GCM simulations. As mentioned above, the efficiency of the heat engine responsible for sensible heating was noticed to be twice smaller than the theoretical value (Koll & Abbot 2016). Thus, it is set to ε_sen = 0.5.

As the timescale of Rayleigh drag and the efficiency parameter of advection both scale the inter-hemispheric heat transport due to atmospheric circulation (see Eq. (75)), only one of these two parameters is needed to quantify the flux. Thus we arbitrarily fix t_drag = 10 days, and choose the efficiency parameter ε_adv so that L_adv ≫ 1, which corresponds to the asymptotic regime of the thermally homogenised atmosphere studied by Wordsworth (2015). The parameters values used in the reference case are gathered in Table 2.

Fig. 9 Atmospheric, dayside and nightside surface temperatures as functions of control parameters Lsen and Ladv for optically thin (τs;L = 0.1, top panels) and thick (τs;L = 1.0, bottom panels) atmospheres in longwave. Left: Ladv = 0.01 (weak circulation). Middle: Ladv = 1.0 (medium-strength circulation). Right: Ladv = 100 (strong circulation). In all cases, the atmosphere is quasi-transparent in the shortwave (τs;S =10−4), pure absorption is assumed (βL0 = βS0 = 1), As = 0.2, and F⋆ = 1366 W m−2, similarly as in Fig. 6. The dayside atmospheric (solid pink line) and surface (dashed red line) temperatures, as well as the nightside atmospheric (dotted cyan line) and surface (dashed blue line) temperatures are plotted as a function of the logarithm of Lsen. The notation log10 designates the decimal logarithm.

6.1 Optical opacity favours collapse at low stellar fluxes

Our investigations start with the role played by longwave opacity on atmospheric collapse. By performing simulations with a GCM with correlated-k radiative transfer, Wordsworth (2015) computed a diagram of the atmospheric stability as a function of stellar flux and surface pressure in the reference case of the CO₂-dominated atmosphere for an Earth-sized planet and a 10 M_⊕-super-Earth (Wordsworth 2015, Fig. 12). This diagram shows that the thresholds between stable and unstable regions of the parameter space predicted by GCM simulations and the analytic theory in the optically thin limit diverge from each other at low stellar fluxes, while they match well at high stellar fluxes. The divergence is patent around F_⋆ ≈ 500 W m⁻², where the critical pressure increases drastically as the stellar flux decays.

Following Wordsworth (2015), we interpret this feature as a limitation of the optically thin approximation assumed in this early work. In order to better understand it, we compute the nightside temperature both with Wordsworth’s model (designated by the acronym “W2015”), and our model (“This work”), which captures in a simplified way the non-linear dependence of the atmospheric thickness onthe surface pressure. Figure 10 shows the results of these calculations in the purely radiative case (L_sen = 0, L_adv = +∞) and in the case where sensible heating is introduced (L_sen≠0), refereed to as “General case”. As mentioned above, this later case corresponds to the asymptotic regime treated by Wordsworth (2015), where the day-night heat transport is efficient (L_adv ≫ 1), leading to T_a;n ≈ T_a;d. In this regime, the values chosen for the timescale of Rayleigh drag (t_drag) and the efficiency parameter of advection (ε_adv) do not affect the planet’s thermal state of equilibrium and the associated atmospheric stability from the moment that L_adv ≫ 1. Besides, the computed nightside temperature and atmospheric stability correspond to upper bounds since they can only decay if the efficiency of the day-night heat transport is modified³.

The values of parameters used in the two cases are given by Table 2. Nightside temperatures as well as the corresponding stability diagrams are plotted as function of the stellar flux and surface pressure in logarithmic scales. For comparison, the incident stellar fluxes of TRAPPIST-1 f (Gillon et al. 2017) and a hypothetical tidally locked Venus are indicated in the bottom right panel of the figure (green triangles).

We first focus on stability diagrams (Fig. 10, bottom panels). Both in the purely radiative and in the general cases, the two models diverge at low stellar fluxes. While the evolution of the critical pressure with the stellar flux is well described by a power law in the optically thin limit (W2015), it undergoes a radical change of behaviour around F_⋆ ≈ F_⊕ when the non-linear dependence of the atmospheric optical thickness on surface pressure is taken into account (this work). This change may be explained by comparing the scalings of the nightside temperature with surface pressure given by the two models (Fig. 10, top and middle panels). Below p_s ≈ 1 bar, the atmospheric opacity is linear with κ_L, as predicted by Eq. (53) in the purely radiative case (see Table 1). Combining the fact that in the optically thin limit with the expression of the condensation temperature of CO₂ given by Eq. (79) yields the scaling (82)

where K is the condensation temperatureof carbon dioxide at 1 mbar. We remark that this scaling law depends on the relation between the longwave optical thickness and pressure. In the case where this relation is a power law of the form with n >0, the preceding scaling becomes . For instance, if the opacity increases due to pressure broadening, then n = 2 (Robinson & Catling 2012; Pierrehumbert 2011), meaning that the dependence of the critical pressure on the stellar flux is weaker in this case than in the studied case.

Beyond p_s ≈ 1 bar, the atmospheric opacity stops growing linearly with κ_L. As a consequence, as surface pressure increases, which leads to the observed dependence inversion. We note that the change of behaviour predicted by the model exceeds that derived from GCM simulations since the critical pressure tends to decay as surface pressure increases for p_s ≳ 3 bar. This may be related to the limitations of the isothermal approximation made to derive radiative transfer functions in Sect. 3, which is not appropriate to model the structure of thick atmospheric layers as discussed in the case of Venus. In spite of these limitations, the isothermal approximation turns out to be sufficient to capture the effect of large longwave opacities on the atmospheric stability at low stellar fluxes, and recover the associated behaviour of the critical pressure highlighted by GCM simulations.

Besides,we observe that the prediction of the model for atmospheric stability in Fig. 10 corresponds to an upper limit since stability diagrams were computed in the asymptotic regime previously treated by Wordsworth (2015), which maximises day-night heat transport (L_adv ≫ 1). Decreasing the efficiency parameter of advection would decrease the atmospheric stability and widen the region of the parameter space where collapse may occur as detailed in Sect. 6.3. In the light of these considerations, the obtained stability diagrams (Fig. 10, bottom panels) predict atmospheric collapse on TRAPPIST-1 f, which agrees with the results obtained by Wolf (2017) using 3D GCM simulations.

Fig. 10 Nightside temperature (K) predicted by analytical theory in this work (top panels) and in Wordsworth (2015) (middle panel), and associated comparative stability diagrams (bottom panels) as functions of logarithms of stellar flux (horizontal axis) and surface pressure (vertical axis). Left: purely radiative case as defined in Sect. 3.3 (Lsen = 0, Ladv = +∞). Right: general case (i.e., with sensible heating and heat transport by stellar and anti-stellar atmospheric circulation). The stellar flux is normalised by the modern Earth’s value F⊕ = 1366 W m−2, and the surface pressure is given in bar. In stability diagrams, the acronyms “TW” and “W2015” are used for “This work” and “Wordsworth (2015)”, respectively. Colours indicate the stability of the steady state predicted by the two models: both models predict stability (orange-red areas); the present model predicts stability while W2015 predicts collapse (orange areas); the present model predicts collapse while W2015 predicts stability (light blue areas); both models predict collapse (blue areas). The incident stellar fluxes of TRAPPIST-1 f and a hypothetical tidally locked Venus are designated by green triangles. These plots are obtained by solving Eqs. (78) and (45) of Wordsworth (2015) with parameters values given by Table 2.

6.2 Dayside sensible heating increases stability

As a secondstep, we investigate how the dayside sensible heating affects the atmospheric stability. To do so, we place ourselves in the asymptotic regime studied by Wordsworth (2015), where the atmosphere has the same temperatureon the dayside and nightside hemispheres. This allows us to ignore the effects of advection for the moment (L_adv = +∞) and to use the simplified equation given by Eq. (65) in the root-finding procedure instead of the general one (Eq. (78)). We solve this equation for various values of the efficiency parameter controlling the intensity of sensible heating (ε_sen) and proceedsimilarly with Wordsworth’s model (Wordsworth 2015, Eq. (45)).

Figure 11 shows the stability diagrams derived from these calculations. The general case (i.e., with sensible heating) is compared with the purely radiative case treated in Sect. 3.3 (Fig. 11, left column), and with the general case of W2015 (Fig. 11, right column). In addition, outcomes computed from 3D GCM simulations for the 1 M_⊕-planet in this early study (Wordsworth 2015, Fig. 12, top panel) are included (Fig. 11, top right panel). Violet dots indicate simulations where the atmosphere remained stable, and blue dots simulations where collapse occurred.

As discussed in Sect. 4, increasing the intensity of sensible heating acts to warm the atmosphere up. Consequently, the critical pressure decays as the efficiency of the dayside convective heat engine increases. However, this evolution is bounded. From the moment that L_sen ≫ 1, T_a;d ≈ T_s;d and p_crit reaches a minimum that is not sensitive to L_sen, as observed for the nightside temperature in Figs. 6 and 9. One may show analytically that the maximum amplitude of critical pressure variations is of 0.6 decades in the studied case.

First, assuming that T_a;d = T_s;d in the strongly convective asymptotic regime (L_sen ≫ 1), we obtain from Eq. (64) that (83)

which yields the nightside temperature in this regime, (84)

If the atmosphere is transparent in the shortwave and optically thin in the longwave, then this expression simplifies to (85)

The critical optical depth below which collapse occurs is determined by the equality . This critical optical depth is denoted by in the purely radiative regime and in the regime dominated by sensible heating, and the corresponding critical pressures are denoted by and , respectively. Thus, assuming that , and using both Eqs. (53) and (85), we get the ratio between these two boundaries (86)

that is roughly 0.6 decades. This corresponds to what is observed for ε_sen = 1 in Fig. 11 (top left panel).

In W2015, no efficiency parameter such as ε_sen intervene to modulate the sensible heat flux except a geometrical coefficient (χ), which is implicitly set to 1 here. This is why the critical pressure predicted by W2015 is the same in all panels of Fig. 11. Interestingly, setting ε_sen = 1.0 reproduces the asymptotic behaviour of the critical pressure derived from W2015 at high stellar fluxes, although the two models use different approaches to scale the horizontal wind speed at planet’s surface. Moreover, we note that the behaviour of the critical pressure predicted by Eq. (78) at low fluxes matches GCM simulations fairly well in this case.

Fig. 11 Comparative stability diagrams as functions of logarithms of stellar flux (horizontal axis) and surface pressure (vertical axis) for various intensities of sensible heating, εsen = 10−4, 10−2, 100 (from bottom to top). Left: comparison between the general case (Lsen > 0) and the purely radiative case characterised in Sect. 3.3 (Lsen = 0). Right: comparison between this work and the model developed by Wordsworth (2015) in thegeneral case. The stellar flux is normalised by the modern Earth’s value F⊕ = 1366 W m−2, and the surface pressure is given in bar. The acronyms “TW”, “RC”, and “W2015” are used for “This work” (with sensible heating), “Radiative case” (without sensible heating), and “Wordsworth (2015)”, respectively. Colours indicate the stability of the steady state predicted by the two models, similarly as in Fig. 10. These plots are obtained by solving Eq. (65) (Ladv = +∞) and Eq. (45) of Wordsworth (2015) with parameters values given by Table 2. Violet dots (top right panel) indicate the region of atmospheric stability computed by Wordsworth (2015) from GCM simulations (Fig. 12), while small blue dots indicate the region of collapse.

6.3 Increasing the planet’s size decreases stability

In the text above, the atmosphere is assumed to be homogeneous in temperature (T_a;n = T_a;d). We now assume that the day-night atmospheric temperature distribution is controlled by the mechanism of heat transport by atmospheric circulation, which depends on dayside and nightside temperatures in return. This allows us to study the link between the intensity of advective heat transport – scaled by the efficiency parameter ε_eff – and the atmospheric stability. Particularly, we aim to show how the 0-D theory, in spite of its limitations, may give us some insight to better understand the decrease of atmospheric stability that is observed as the planet size increases (e.g. Wordsworth 2015, Fig. 12).

In order to benchmark the predictions of the model against outcomes computed from GCM simulations, we reproduce the case treated by Wordsworth (2015) and thus consider two planets: the Earth-sized planet characterised by Table 2 and a super-Earth of 10 M_⊕ and 1.88 R_⊕. As GCM simulations were performed using correlated-k distributions for radiative transfers in this early work, there is no evident value for the equivalent gray-gas longwave opacity. We thus set the absorptivity to κ_L = 8.5 × 10⁻⁵ m² kg⁻¹ for a pedagogical purpose, this value leading the analytic model to match GCM simulations for a simple value of ε_eff, as shown inthe following. The equation of the steady state in the general case (Eq. (78)) is then solved for various values of the efficiency parameters controlling sensible heating (ε_sen) and heat transport by atmospheric circulation (ε_eff).

Diagrams comparing the atmospheric stability of the two planets are plotted in Fig. 12 in the absence of sensible heating (ε_sen = 0) and in the regime of strong dayside convection (ε_sen = 1.0), for (this later case is treated by solving Eq. (65) instead of Eq. (78)). Besides, we include the data computed from 3D GCM simulations (Wordsworth 2015, Fig. 12) in the figure (middle right panel). Violet dots designate simulations where the atmosphere remained stable in both cases, while red dots indicate those where the atmosphere of the 1 M_⊕-planet only remained stable. The figure shows that, in all cases, the atmospheric stability decays as the planet size increases.

We first consider the purely radiative regime (Fig. 12, top left panel). In this case, the observed difference between the Earth-sized planet and the super-Earth comes from the dependence of the longwave atmospheric optical depth on the planet mass and radius that is hidden in the surface gravity (Eq. (81)). As the planet size increases, the atmospheric optical depth decays, and so does the greenhouse effect. Consequently, the nightside surface temperature decays. The difference between the two planets does not change when sensible heating is included (top right panel) because this flux does not depend on the planet size at first order in the scaling given by Eq. (63) (there is a dependence hidden in Q_in but it is weakened by the 1∕3 exponent in the scaling of the horizontal wind speed). However, we observe that the gap of collapse pressure between the two planets predicted by GCM simulations is larger than that associated with the aforementioned dependence. Among the multiple reasons that may be invoked to explain this difference, we examine the hypothesis where it is related to a decrease of the efficiency of heat redistribution mechanisms through the example of large-scale advection.

Except in the asymptotic limit of strong horizontal mixing (L_adv ≫ 1), the intensity of the heat flux due to large-scale advection is modulated by the efficiency parameter ε_eff, leading the atmospheric stability to decay with ε_eff (Fig. 12, from top to bottom). Furthermore, in the derived scaling of the heat transport by stellar and anti-stellar atmospheric circulation, the heat flux is quantified by (Eq. (75)) and thus decays as the planet size increases. As a consequence, the intensity of the heat flux from the dayside to the nightside decays as one switches from the Earth-sized planet to the super-Earth, which tends to widen the gap between the two planets (Fig. 12, bottom panels). This effect is stronger in the transition regime (L_adv ~ 1) than in asymptotic ones, as shown by the evolution of the nightside and dayside temperatures ratio plotted in Fig. 8.

The case (i.e. ε_eff ≈ 3 × 10⁻³, middle panels) highlights the divergence between the analytic model and GCM simulations (Wordsworth 2015, Fig. 12). First, we note that the model tends to underestimate atmospheric stability at low stellar fluxes and high pressures, which is partly due to the simplified temperature profile used to derive radiative fluxes, as discussed in Sect. 6.1. Second, the critical pressure predicted by the analytic model scales similarly with the stellar flux for both planets, while it is not the case in GCM results. This is obviously due to the numerous simplifications made in the present work, where the processes responsible for heat redistribution are described using rough scalings. Despite these limitations, the model approximately matches GCM simulations and captures the decrease of atmospheric stability caused by the increase of the planet size. It should be possible to improve the analytical theory with better scalings of the fluxes in future studies.

Fig. 12 Comparative stability diagrams for 1 M⊕ and 10 M⊕ planets as functions of logarithms of stellar flux (horizontal axis) and surface pressure (vertical axis). The efficiency parameter controlling advective heat transport εeff take the values (from bottom to top). Left column: case without sensible heating (εsen = 0). Right column: case with efficient sensible heating (εsen = 1.0). The stellar flux is normalised by the modern Earth’s value F⊕ = 1366 W m−2, and the surface pressure is given in bar. The acronyms “E” and “SE” are employed for “Earth” and “Super-Earth”, and designate the 1 M⊕ and 10 M⊕ planets, respectively. Colours indicate the stability of the steady state obtained in the two cases, similarly as in Fig. 10. These plots are obtained by solving Eq. (65) (Ladv = +∞, top panels) and Eq. (78) (other panels) with κL = 8.5 × 10−5 m2 kg−1 and parameters values given by Table 2 for the Earth-sized planet. In the case of the 10 M⊕ planet, the planet’s radius and mass are set to Mp = 10 M⊕ and Rp= 1.88 R⊕, respectively. Data computed by Wordsworth (2015) from GCM simulations (Fig. 12) are included in the middle right panel for comparison and indicate the region of atmospheric stability for the 1 M⊕ planet (red dots) and both planets (violet dots), and the region of atmospheric collapse (small blue dots).

7 Role played by atmospheric structure

In our zero-dimensional approach, the implications of the atmospheric structure on the nightside temperature are ignored, given that atmospheric temperature profiles are reduced to dayside and nightside bulk atmospheric temperatures. This approximation enabled us to derive radiative fluxes analytically, and to capture in a simplified way the non-linear dependence of the atmospheric opacity on the shortwave and longwave optical depths. In reality, the atmospheric structure is not the same on the dayside and on the nightside, and differs from the idealised isothermal temperature profile, as highlighted by analytical models (e.g. Robinson & Catling 2012; Koll & Abbot 2016) and GCM simulations (e.g. Leconte et al. 2013; Wordsworth 2015). Particularly, the evolution of the temperature profile with the stellar zenithal angle appears clearly in simulations performed using THOR (Fig. 2, left panels). The substellar region is characterised by a quasi-adiabatic temperature gradient, which is due to the strong convection generated by the absorption of stellar heating near the ground. On the nightside, the structure inversion that may be observed indicates that the atmosphere is stable with respect to convection.

As mentioned above, the role played by the atmospheric structure was already investigated in the framework of the analytical theory. For instance Robinson & Catling (2012) proposed a one-dimensional radiative convective model to predict the atmospheric structure of a wide range of planets. Koll & Abbot (2016) refined this approach in the case of tidally locked rocky planets by including in a simplified way the effect of subsidence on the nightside temperature profile and surface temperature, which led them to develop a two-column radiative-convective-subsiding (RCS) model. In this early work, the dayside temperature profile was supposed to be shaped by convection, and was thus specified in calculations. This assumption holds in the case of convective tropospheres, which corresponds to the regime of strongly irradiated atmospheres transparent in the shortwave. However, it may break down if the absorption near the ground is not large enough to generate strong convection. In view of these limitations, the prospect of extending the theory to optically thicker atmospheres in the shortwave motivates us to go past this assumption.

Therefore, although the core of the present study is the zero-dimensional model, we make a first attempt to relax the specification of the dayside atmospheric structure in one-dimensional models for the purpose of this section. We aim to give here some insight about the transition from zero-dimensional to one-dimensional models where both the dayside and nightside temperature profiles are computed from atmospheric opacities. This leads us to develop a simplified two-column radiative-convective-subsiding-upwelling (RCSU) model by introducing dayside upwelling winds in the RCS model of Koll & Abbot (2016). The RCSU model is detailed as a first step, and used as a second step to compute dayside and nightside atmospheric temperature profiles and study atmospheric stability against collapse.

7.1 A two-column RCSU model

To treat the coupling between the atmospheric structure and mean flows, one should integrate self-consistently the momentum, mass conservation, thermodynamic and radiative transfer equations, which is achieved by three-dimensional GCM but cannot be envisaged in the simplified framework of one-dimensional models. We thus introduce vertical winds by using the thermodynamic equation solely, following the method by Koll & Abbot (2016). Written in pressure coordinates, this equation reads (Vallis 2006) (87)

the notation∇_p designating the horizontal gradient over an isobar, V the horizontal velocity vector, the pressure velocity (i.e. the vertical velocity in pressure coordinates; stands for the material derivative), and the vertical diffusive power flux. To relate pressure to optical depth consistently with the early work by Koll & Abbot (2016), we use the standard power law (88)

where the exponent n specifies how optical depths increases with pressure. Introducing the dimensionless parameter β ≡ and substituting the pressure coordinate by the optical depth in the longwave, Eq. (87) thus becomes (89)

In steady state, ∂_tT = 0. Besides, we ignore the diffusion term assuming that radiative transfers and advective heat transport predominates. In the general case, the horizontal advection term (V ⋅∇_pT) is comparable to the vertical advection term (), and should be taken into account. However, we ignore it for simplicity following Koll & Abbot (2016), since including horizontal advection would require to specify ad hoc horizontal wind speeds and temperature distributions, which goes beyond the scope of a simplified two-column approach. Nevertheless, one may think about clever ways to introduce this component in the modelling in future studies. The above approximations yield (90)

where the complex atmospheric circulation is now reduced to one single quantity, ω, scaling the strength of vertical flows. We notice that ω > 0 corresponds to subsidence, and ω < 0 to updraft, since pressure and optical depth decay with altitude.

In this one-dimensional model, we consider a regime similar to that treated in Sect. 4 in the framework of the zero-dimensional RC model. The heat transport from the dayside to the nightside is supposed to be efficient, and the stratosphere is consequently globally isothermal. However, the two cells introduced in the zero-dimensional model to characterise sensible heating and large-scale heat transport by atmospheric circulation separately are now encompassed into one single stellar and anti-stellar cell (see Fig. 5), so that updraft takes place within the dayside substellar region and subsidence on the nightside. Moreover, we assume uniform vertical profiles for pressure velocity both in the updraft and subsidence regions. The corresponding amplitudes are denoted by ω_up and ω_sub, respectively.The conservation of mass thus yields the relationship (91)

the notations A_up and A_sub designating the updraft and subsidence areas, respectively. As shown by snapshots of vertical wind speeds obtained from GCM simulations in typical cases (Fig. 2, right panels), there is an asymmetry between rising and sinking air motions. Particularly, air rises rapidly near the substellar point where heat absorption is the strongest, while it sinks slowly over a large area in cooler atmospheric regions. This asymmetry is quantified by the ratio A_up∕A_sub, which is a fixed parameter of the model. As it characterises the strength of convection, the pressure velocity of rising air in the substellar region is related to the typical horizontal wind speed V_sen introduced previously in the expression of the sensible heating flux (Eq. (54)) and quantified using the heat engine theory (Eq. (57)). A simple scale analysis based on the conservation of mass yields (92)

With a similar analysis, we derive the typical time for a fluid parcel to rise up and to subside, (93)

Because of the asymmetry between rising and sinking motions, t_sub ≫ t_up. For simplification, we assume in the following that and A_sub∕A_up = t_sub∕t_up ~ 10 and ignore the subsidence happening on the dayside (see Fig. 2, right panels). The steady state thermodynamic equation thus reads, for the dayside updraft region, (94)

and, for the dayside subsidence region, (95)

where t_rad is the timescale of radiative cooling defined by Eq. (13). We note that the dayside net flux in the right-hand member of Eq. (94) is the sum of the longwave and shortwave components. For the later, we use the analytical solution given by Eq. (22), which is not limited to the case of isothermal profiles since it does not depend on temperature.

Combined with the day- and nightside net flux equations (Eq. (16)), Eqs. (94) and (95) form the system of equations governing the global heat redistribution and dayside and nightside temperature profiles. For both the dayside and nightside, two boundary conditions are required to integrate the radiative transfer equation, and one to integrate the thermodynamic equation, meaning that six boundary conditions have to be specified. As may be noticed, the simplified dynamics adopted to derive Eqs. (94) and (95) generates a singularity at τ_L = 0 in their left-hand members. For this reason, upper boundary conditions cannot be applied at τ_L = 0, and we have to set arbitrarily the optical depth of the top of the atmosphere τ_{top; L} > 0. We verify a posteriori that results do not vary with τ_{top; L} from the moment that τ_{top; L} ≪ τ_s;L. We remark that Koll & Abbot (2016) could preliminarily calculate the level of the tropopause to which they applied nightside boundary conditions by assuming a convective dayside temperature profile. As the dayside temperature profile is not fixed in the present work, the pressure level of the nightside upper boundary cannot be determined in a similar way. Thus, upper boundary conditions have to be applied to the fixed level τ_{top; L}.

On the dayside substellar area of rising air, we assume that there is no radiative source at infinity in the longwave, and subtract the amount of power available to drive atmospheric motion from the outgoing longwave radiation, assuming that dissipation occurs outside of the area. In the scaling of the typical wind speed V_sen given by Eq. (57), the power available to drive atmospheric motion is quantified by , where Q_in is the theoretical power defined by Eq. (58), and ε_sen = 1∕2 the efficiency coefficient evaluated by Koll & Abbot (2016) for the wind speed using GCM simulations⁴. We remind here that the day-night heat transport by atmospheric circulation is not taken into account in this section. The stratosphereis considered as globally isothermal, meaning that dayside and nightside atmospheric temperatures are equal at the upper boundary. The dayside surface temperature is defined by the blackbody radiation of the surface, meaning that T_a;d = T_s;d at surface. On the nightside, we assume the absence of downwelling longwave flux at infinity, as on the dayside. Besides, following Koll & Abbot (2016), we consider that the net flux in the longwave vanishes at the planet surface on the nightside – which is a consequenceof stable stratification (no sensible heat flux) – and that the stratosphere is horizontally isothermal. Mathematically, these conditions read, for the dayside,

and, for the nightside,

The boundary value problem is solved over a fixed range of optical depths by means of a relaxation method. Starting from initial temperatureprofiles, the flux and temperature equations are integrated iteratively until convergence, an adaptive relaxation coefficient depending on the history of residuals being applied every iteration to stabilise the convergence process. At a given step, the dayside and nightside radiative transfer equations are integrated first with a finite difference scheme, by means of the shooting method (Press et al. 2007, Sect. 18.1) and Thomas’s algorithm (Press et al. 2007, Sect. 2.4), respectively. Then, quantities parametrising the dynamics are calculated using auxiliary equations. Finally, the dayside and nightside thermodynamic equations are integrated, which yields new temperature profiles. The old temperature profiles are then incremented by the temperature difference between the old and new solutions weighted by the relaxation coefficient. This defines the temperature profiles for the next iteration, and so on.

Fig. 13 Temperature structure and fluxes profiles calculated with radiative-convective-subsiding-upwelling (RCSU) model in reference case for various surface pressures as functions of pressure levels in logarithmic scale. Left: ps = 0.01 bar and τs;L = 0.01. Middle: ps = 0.1 bar and τs;L = 0.1. Right: ps = 1 bar and τs;L = 1. Dayside atmospheric (solid pink line), nightside atmospheric (dashed sky blue line), dayside surface (red square), nightside surface (blue square), and radiative equilibrium (dotted grey line) temperatures (K) are plotted in top panels. Dayside (shades of red) and nightside (shades of blue) radiative net, total, upwelling, and downwelling fluxes in the longwave, and net flux in the shortwave (dashed cyan line) are plotted in bottom panels. Fluxes are normalised by the black body equilibrium flux , where Teq is the black body equilibrium temperature defined by Eq. (2). Calculations are performed using parameters values given by Table 2 with n = 1, which corresponds to Fig. 2 (right panel) of Wordsworth (2015). Complementary parameters of the RCSU model are set to and Asub∕Aup = 10.

7.2 Temperatures and fluxes profiles

We benchmark the two-column RCSU model against simulations performed by Wordsworth (2015), who used the three-dimensional LMD GCM (Hourdin et al. 2006) with grey gas opacities. These simulations correspond to the reference case described by parameters values of Table 2 with n = 1 and p_s = 0.01, 0.1, 1 bar (see Wordsworth 2015, Fig. 2, right panel). We remind ourselves that the associated optical depths in the longwave are τ_s;L = 0.01, 0.1, 1, respectively, owing to the chosen effective longwave opacity (κ_L = 5 × 10⁻⁵ m² kg⁻¹). As may be noted, the substellar region is far smaller than the whole dayside hemisphere. Thus, the dayside atmospheric structure in this region is calculated by taking the value of the incident stellar flux at the substellar point, which is F_⋆. This differs both from Wordsworth (2015) and Koll & Abbot (2016), who calculated dayside averaged profiles (with an averaged incident flux of F_⋆∕2), and leads consequently to higher dayside temperatures and stronger fluxes.

Temperatures and fluxes profiles computed using our RCSU model are plotted in Fig. 13, which is the counterpart of Fig. 2 (right panel) in the study by Wordsworth (2015). For comparison, we indicate in plots the radiative equilibrium temperature T_rad, which corresponds to F_−L = 0, and is expressed in the absence of scattering as (e.g. Pierrehumbert 2011) (102)

In this equation, T_skin designates the skin temperature, which is temperature the outer regions of the atmosphere would have in the absence of in situ heating by stellar absorption (Pierrehumbert 2011, Sect. 3.6, p.169). In the case where the atmosphere is transparent to incident stellar radiation, , which is typically the temperature of the stratosphere.

We recover in Fig. 13 the behaviour described by the zero-dimensional RC model for surface temperatures, which are both increasing with the atmospheric optical thickness in the longwave because of greenhouse effect. Quantitatively, the predictions of the RCSU model for the nightside surface temperature match rather well the results obtained from GCM simulations, although the RCSU model tends to overestimate the nightside temperature with respect to the GCM in the case p_s = 1 bar. The RCSU model also captures the variation of atmospheric structure between the dayside and the nightside. While convection on the dayside leads to a quasi-adiabatic temperature profile, the atmosphere is cooler and stably stratified on the nightside, which induces a smaller temperature gradient. Particularly, we retrieve here the nightside structure inversion observed in GCM simulations (e.g. Leconte et al. 2013; Koll & Abbot 2016) and due to subsidence, as discussed by Koll & Abbot (2016).

Because of the assumed vertically-uniform pressure velocity, both dayside and nightside atmospheric temperatures decay with pressure in order to compensate the increase of the singular term of the simplified thermodynamic equations (Eqs. (94) and (95)). In reality, vertical winds vanish above the troposphere, and there is no singularity. The tropopause approximately corresponds to the pressure level where T_a;d and T_a;n join together. Above this pressure level, the temperature profile should approach the radiative equilibrium temperature profile T_rad (Eq. (102)) instead of decaying, as shown by GCM simulations (e.g. Wordsworth 2015, Fig. 2, right panel). In spite of the unrealistic tendency observed in the stratosphere, the nightside temperature profile derived from the RCSU model in the troposphere is qualitatively similar to that obtained by Koll & Abbot (2016) with their radiative-convective-subsiding model (Fig. 4 of the article). Thus we may now compare the predictions of the two models regarding atmospheric stability and test thereby the RCSU model.

7.3 Atmospheric stability

Considering the reference case of the study, we run grid calculations using the RCSU model to compute the nightside temperature of the planet, and determine how the stability of the steady state against collapse evolves withthe incident stellar flux and surface pressure. Consistent with previous sections, the atmospheric stability is simply obtained by comparing the nightside surface temperature to the condensation temperature of CO₂, given by Eqs. (79) and (80). We note that the boundary value problem becomes difficult to solve accurately in the optically thick limit because tiny changes of the upper boundary condition (τ_L = τ_{top; L}) in Schwarzschild equation (Eq. (16)) substantially affect radiative fluxes at large τ_L. As we start encountering convergence issues for τ_s;L ≳ 10, we set the upper bound of the surface pressure range to (bar).

Figure 14 shows the obtained stability diagram, on which data computed by Wordsworth (2015) from GCM simulations are also plotted. This diagram may be compared to those calculated using the zero-dimensional RC model (Fig. 11) and to the results obtained by Koll & Abbot (2016) with their RCS model (Koll & Abbot 2016, Fig. 12, left panel). Although it includes the day-night difference of atmospheric structure, the RCSU model does not match GCM simulations, which were well approached by the RCS model of Koll & Abbot (2016). Particularly, the model does not capture the behaviour caused by the non-linear dependence of temperature on the atmospheric optical thickness, whereas this behaviour was captured by the 0-D approach, as shown by Fig. 11 (top right panel). The scaling of the collapse pressure given by the RCSU model actually corresponds to the purely radiative case of the 0-D model, that is (Eq. (82)).

These limitations of the theory are apparently due to the oversimplified dynamics assumed to derive the differential equations of dayside and nightside temperatures (Eqs. (94) and (95)). Particularly, the pressure velocity (i.e. the velocity in pressure coordinates) is supposed uniform over the whole air column while it should tend to zero above the tropopause. As a consequence, the model does not describe properly the stratosphere. The inclusion of shortwave opacities in the two-column RCS approach thus still remain an open question for future studies. Moreover, temperature and winds profiles strongly vary with the stellar zenithal angle on the dayside, which cannot be properly modelled in the framework of a 1D theory. This encourages the development of 2D intermediate models treating the coupling between the large-scale atmospheric circulation and radiative transfers in a self-consistent way in the regime of slowly rotating planets.

Fig. 14 Stability diagram computed with one-dimensional radiative-convective-subsiding-upwelling (RCSU) model in reference case. The orange area indicate a stable atmosphere, and the blue area collapse in the RCSU model. Data computed by Wordsworth (2015) from GCM simulations (Fig. 12) are included for comparison. Violet dots indicate simulations where the atmosphere remained stable, while small blue dots indicate atmospheric collapse. Calculations are performed using parameters values given by Table 2 with n = 1, which corresponds to Fig. 2 (right panel) of Wordsworth (2015). Complementary parameters of the RCSU model are set to and Asub∕Aup = 10.

8 Discussion and conclusions

In order to better understand the mechanism of atmospheric collapse, which is of fundamental importance to characterise surface conditions of tidally locked rocky planets, we developed a hierarchy of models building on pioneering works of analytical theory in the spirit of Held (2005). We aimed thereby to enrich and to consolidate the theory by means of cross-comparisons between multiple approaches ranging from 0-D scalings to 3D GCM simulations. Models of increasing complexity were thus developed following a step-by-step method to treat the case of dry atmospheres hosted by slowly rotating rocky planets by including incrementally the couplings between radiative transfer, convection and large-scale circulation. In these models, radiative transfer is described using the two-stream, dual-band and grey gas approximations, while the heat transport by atmospheric circulation is treated in the framework of heat engine theory following Koll & Abbot (2016).

We first revisited the early study by Wordsworth (2015), which is based on a 0-D approach. Assuming an isothermal atmosphere,we derived analytically the short- and longwave net flux profiles in presence of absorption and scattering, as well as the atmospheric, dayside and nightside surface temperatures characterising the steady state in the purely radiative case. This enabled us to quantify at first order how scattering and shortwave absorption modify the equilibrium temperatures and atmosphericstability against collapse by acting on greenhouse effect. Typically, longwave and shortwave absorption both increase the greenhouse effect, and thus the atmospheric stability, which is not the case of scattering. While scattering in the shortwave induces anti-greenhouse cooling and favours collapse, scattering in the longwave tends to stabilise the atmosphere by generating scattering greenhouse effect. The model captures these dependences in a simplified way.

It also captures the particular behaviour of the collapse pressure (or critical pressure) at low stellar fluxes due to the non-linear dependence of the atmospheric thickness on the longwave opacity. Although highlighted by Wordsworth (2015) using GCM simulations, this behaviour was absent from the theory proposed by the author since the regime of optically thin atmospheres was solely considered in this early study. It results from the fact that the absorbed power ceases to grow linearly with the atmospheric optical depth at the planet’s surface as one leaves the optically thin regime. The present study shows that, in spite of being a rough approximation of the atmospheric structure, the isothermal atmosphere assumption is sufficient to capture the induced scaling of collapse pressure with stellar flux that may be observed in GCM simulations.

As a second step, we included in the model the turbulent exchanges due to convection within the dayside boundary layer in the asymptotic regime where the atmosphere is globally isothermal (efficient large-scale heat transport). The sensible heat flux between the atmosphere and surface was scaled considering the atmosphere as a heat engine controlled by the atmospheric and dayside surface temperatures. In this framework, we derived the steady-state equation that governs the global heat redistribution, and solved it numerically to recover the early results obtained by Wordsworth (2015) and Koll & Abbot (2016) in the case of CO₂-dominated atmospheres. Besides, we derived analytically a lower bound for the collapse pressure. This lower bound corresponds to the regime of strong convection, which maximises greenhouse effect.

As a third step, we ceased to treat the atmosphere as a globally isothermal layer and examined the case where the dominating heat transport mechanism is coupled with dayside and nightside temperatures. By way of an example, we supposed the heat to be advected from the dayside to the nightside by a stellar and anti-stellar atmospheric circulation, which is the typical dynamical regime of slow rotators. The associated heat flux was thus quantified from a scale analysis by treating the stellar and anti-stellar cell as a heat engine, as done for dayside convection. This simplified modelling accounts for the fact that the efficiency of heat transport by large-scale advection decays as the planet size increases. It allowed us to derive a generalised version of the steady state equation including the coupled effects of radiative transfer, sensible heating, and large-scale atmospheric circulation.

In GCM simulations, the atmospheric stability decays as the planet-size increases. Our generalised 0-D model approximately captures this behaviour and relate it to two conjugated effects: (i) for given gas opacity and surface pressure, the optical depth at the planet’s surface decays as the surface gravity increases, which weakens greenhouse effect; and (ii) the ratio of the advective and radiative timescales decay, which decreases the hemisphere-averaged amount of heat transported from the dayside to the nightside.

The 0-D approach may certainly be improved in future studies by introducing refined scalings for heat fluxes and mean flows, which would better capture the dependence of the collapse pressure on the physical properties of the star-planet system. However, the model already allows us to characterise asymptotic regimes using a small number of dimensionless control parameters, such as – for instance – L_sen and L_adv, which compare dayside sensible heating and large-scale heat transport with radiative cooling, respectively. Furthermore, the fact that the root-finding procedure solving the steady state equation can be massively parallelised makes it possible to widely explore the parameter space. This clearly encourages the development of analytic 0-D models as a complement of studies based on GCM simulations.

As a last step, we investigated the transition between 0-D and 1D models by studying how the atmospheric circulation affects the dayside and nightside atmospheric structure. Starting from the two-column radiative-convective-subsiding (RCS) model proposed byKoll & Abbot (2016), we developed a radiative-convective-subsiding-upwelling (RCSU) model. In this approach, the dayside and nightside temperature profiles are modified by vertical winds. While upwelling winds due to convection in the substellar region drive the temperature gradient towards the adiabatic profile, subsidence on the nightside causes a structure inversion leading to stable stratification. In order to extend the theory to the case of non-transparent atmospheres in the shortwave, we let the dayside temperature profile be integrated simultaneously with the nightside temperature profile, instead of specifying its shape, which was done by Koll & Abbot (2016).

From a quantitative point of view, the RCSU model does not match well the behaviour of the collapse pressure obtained from GCM simulations, and rather predict a scaling similar to that given by the 0-D model in the purely radiative case. This apparently results from the simplifications made in the dynamics to include the effect of vertical winds. Particularly, in both subsidence and updraft areas, the vertical wind speed was supposed to be uniform over the whole air column, although the nature of the dynamics is not the same in tropospheric and stratospheric layers. Moreover, horizontal advection was ignored notwithstanding the fact that it is comparable to vertical advection in the general case. Future studies will have to remedy to both of these rough approximations to derive self-consistently realistic temperature profiles and collapse pressure at the same time with the 1D two-column approach, which still remains an open question. Nevertheless, in spite of its limitations, the RCSU model is a promising first attempt that allowed us to recover the features identified by early works,and particularly the nightside structure inversion previously captured by the RCS model (Koll & Abbot 2016).

For future works, it is crucial to improve the treatment of the dynamics in the analytic theory since it is a major aspect of the problem. This suggests to develop an intermediate class of 2D atmospheric models describing the complexcoupling between mean flows and the thermodynamics in the asymptotic regime of slowly rotating tidally locked planets, the circulation being symmetric with respect to the star-planet axis in this regime. These simplified 2D GCMs would also be useful to quantify timescales associated with the atmospheric collapses, which cannot be achieved with steady states models, and requires a substantial computational effort when 3D GCMs are used.

In addition to the dynamics, one should evolve towards improved modelling of radiative transfers including wave- length-dependent opacities and accounting for the relationship between atmospheric absorption and stellar spectra, which is a blind spot of the dual-band approximation. This aspect is of fundamental importance when the stellar and planetary radiative fluxes overlap, which is typically the case for temperate exoplanets orbiting cool dwarf stars. TRAPPIST-1 planets are a representative example of such a configuration. This question could be addressed for instance by introducing simplified dynamics in 1D radiative transfer codes.

Similarly as their predecessors, the simplified 0-D and 1D analytical models detailed in the present work are limited to dry thermodynamics and slow rotation, which corresponds to the simplest physics. The presence of water is however a feature of utmost interest in the study of Earth-sized exoplanets located in the habitable zone of their host star because of its implications on climate and prebiotic chemistry. Thus, this feature should be introduced in the analytic theory in order to consolidate the conclusions of works essentially based on GCM simulations (e.g. Merlis & Schneider 2010; Leconte et al. 2013).

Particularly, the presence of a slab ocean induces an additional heat transport both by latent heat flux in the atmosphere and by heat diffusion in the ocean itself (Edson et al. 2011), and tends thereby to stabilise the atmosphere against collapse. Additionally, water vapour leads to the formation of clouds, which affects the atmospheric albedo, thermal structure, and general circulation in a complex way. To complete the theory, similar inclusions may be done for the heat transport by propagating gravity waves, and super-rotation in the regime of rapid rotators. Most of these aspects were discussed in early works (see e.g. Leconte et al. 2013).

Finally, we shall emphasise that steady states do not necessarily exist, and if so, are not necessarily stable. For example, Edson et al. (2011) showed evidence of the bistability characterising the dependence of the wind speeds on the planet spin by highlighting an abrupt transition between two regimes. This transition leads to a significant change of velocity. Regarding this type of questions, the analytic theory has a strong asset with respect to GCMs since it offers the opportunity to characterise formally the existence and stability of steady states.

Acknowledgements

The authors thank the reviewer, D. D. B. Koll, for helpful comments that contributed to improve the manuscript. They acknowledge financial support from the European Research Council via the Consolidator grant EXOKLEIN (grant number 771620). K.H. also acknowledges partial financial support from the National Swiss Foundation, the PlanetS National Center of Competence in Research, the Center for Space & Habitability and the MERAC Foundation. This research has made use of NASA’s Astrophysics Data System.

Appendix A Turbulent heat exchanges

At the planet surface, the heat flux resulting from turbulent friction is written as (A.1)

where is the density of the atmosphere at the planet’s surface, T′ and are fluctuations of temperature and radial velocity, respectively.

Fluctuations of the radial flow are induced by turbulence due to the interaction between the dominating horizontal flow and the planet surface. Thus, we can write the energy transport in Eq. (54) as , where V_* and T_* correspond to frictional velocity and temperature, respectively. Using the flux gradient theory with the mixing length hypothesis, these quantities are defined as (A.2)

where z is the altitude with respect to the surface, ℓ the turbulence mixing length, and and the average horizontal velocity and temperature, respectively.

The turbulence mixing length is calculated from the empirical scaling law of Blackadar (1962), (A.3)

the parameter being the von Kármán constant and ℓ₀ the maximum attainable mixing length in the boundary layer. By substituting Eq. (A.3) into Eq. (A.2), and integrating the equations, we obtain

where z₀ is the roughness height. Thus, considering that z is the thickness of the boundary layer and introducing the bulk drag coefficient (A.6)

we can write the vertical turbulent heat flux (from the surface to the atmosphere) as (A.7)

By assuming ℓ₀ = +∞ (no fixed maximum for the turbulence mixing length), we recover the scaling law of the mixing length initially proposed by von Kármán on the basis of laboratory experiments, (von Kármán 1930), and the classical expression for the bulk drag coefficient given by Eq. (55). The parameters z and z₀ are free parameters in the general case.They are set to z = 10 m and z₀ = 1 cm in Wordsworth (2015).

Appendix B Solving the steady state equation in the general case

Equation (78) defines a root-finding problem. This problem is formulated as (B.1)

where f is the function defined by the left-hand member of Eq. (78), that is

To solve Eq. (B.1), we first have to determine the range of where is defined, which corresponds to the range where the denominator of the expression given by Eq. (77) is strictly positive. This is a preliminary zero-finding problem, which is treated using the standard secant method (e.g. Press et al. 2007). We thus obtain the lower bound of the range, namely . Equation (B.1) is then solved by using a combination of the dichotomy and secant method within the interval .

We remark that f is subject to very sharp variations inthe vicinity of while it smoothly varies with in the rest of the interval. As a consequence, the coordinate is not well appropriate to the root-finding procedure. We rather use the logarithm of the normalised distance, (B.2)

the notation log₁₀ designating the decimal logarithm. The problem thus becomes , where . Finally, to deal with the possible issues resulting from the variations of f over several orders of magnitude, we apply the root finding method to the problem (B.3)

where ○ is the function composition operator and h the smoothing function defined as a function of the variable X by (B.4)

the notation sign referring to the sign function ( if X < 0 or X > 0, respectively, and ).

For an illustrative purpose, f is plotted in Fig. B.1 as a function of for differentvalues of the nondimensional parameters controlling atmospheric circulation, L_sen and L_adv. In these calculations, the atmosphere is optically thick in the longwave (τ_s;L = 10) and thin in the shortwave (τ_s;S = 0.1). The figure shows that there is one unique solution to the steady state equation in each of the treated cases. Besides, roots are found to be very close to with values ranging within the interval . This interval depends on the atmospheric opacity, and increases when the optical depth in the longwave decays.

Fig. B.1 Left-hand member of Eq. (78) as function of normalised temperature difference in logarithm scale for Lsen, Ladv = 0.01, 1, 100. The functionis denoted by f. The parameters Lsen and Ladv are thenon-dimensional control parameters of sensible heating and advected heat transport defined by Eqs. (63) and (75), respectively. Parameters values: βL 0 = βS0 = 1, τs;L = 10; τs;S = 0.1; F⋆ = 1366 W m−2; As = 0.2.

References

All Tables

All Figures

	Fig. 1 Diagram of the two-layer radiative-convective-advective model detailed in the present study. The diagram shows the processes involved in dayside-nightside heat redistribution that are taken into account in the model. Radiative net fluxes (F− ≡ F↑− F↓) are upwards in the convention used, meaning that the shortwave net flux F−S is negative.
In the text

Fig. 2 Instantaneous snapshots of zonally-averaged temperature (left panels) and vertical wind speed distributions (right panels) as functions of latitude and pressure with north pole at substellar point. Top: 0.1-bar CO2 -dominated atmosphere at Teq = 279 K (Wordsworth 2015, Fig. 2). Bottom: 0.5 bar N2-dominated atmosphere at Teq = 400 K (Koll & Abbot 2016, Fig. 4). Latitude is given in degrees (horizontal axis). This coordinate does not correspond to the usual latitude but to a latitude coordinate in the system of spherical coordinates where the north pole is the substellar point and the south pole the anti-stellar point. The associated colatitude is the stellar zenithal angle in this system of coordinates. Pressure is given in bars (vertical axis), temperature in K, and vertical velocity in m s−1. The white area at the substellar point indicates the region where surface pressure is lower than the spherically averaged surface pressure due to the updraft. The dotted red line indicates the transition between rising and subsiding flows, where the vertical velocity is zero.

In the text

Fig. 3 Net radiative fluxes in the short- and longwave (horizontal axis) as functions of optical depths (vertical axis) in case of isothermal atmosphere. Fluxes profiles are computed from the analytic formulae given by Eqs. (22) and (23) and normalised by the black body equilibrium flux (the equilibrium temperature Teq being givenby Eq. (2)). Optical depths are normalised by their value at the planet’s surface, namely τs;S for the shortwave and τs;L for the longwave. Shortwave (blue lines) and longwave (red lines) net fluxes are plotted for various optical thicknesses, the optical depths at the planet’s surface taking the values τs;S = 0.1, 1, 10 (from light to dark blue lines) and τs;L = 0.1, 1, 10 (from light to dark red lines). Pure absorption is assumed (βL0 = βS0 = 1). The incident, atmospheric, and surface fluxes are arbitrarily set to Fi = 4Feq, Fa = 0.7Feq, and Fs = 2Feq, respectively.

In the text

Fig. 4 Atmospheric, dayside and nightside surface temperatures in pure radiative regime. Temperatures are plotted as functions of the logarithm of τs;L in the case of pure absorption (top panels), scattering in the shortwave (middle panels) and scattering in the longwave (bottom panels). Left: quasi-transparence in the shortwave (τs;S =10−4). Middle: small shortwave opacity (τs;S = 10−2). Large shortwave opacity (τs;S = 100). The atmospheric temperature (solid black line), dayside surface temperature (dashed red line), nightside surface temperature (dashed blue line), and black-body equilibrium temperature (dotted grey line) are plotted using Eqs. (51)–(53) and (2), respectively. Parameters values: As = 0.2, and F⋆ = 1366 W m−2.

In the text

	Fig. 5 Diagram of the dayside atmospheric heat engine. The heat engine is driven by dayside heating and cooling to space. A fluid parcel works against friction in the boundary layer, which limits the speed of the flow.
In the text

Fig. 6 Atmospheric, dayside and nightside surface temperatures as functions of logarithm of control parameter Lsen. Left: τs;L = 0.01 (quasi-transparent). Middle: τs;L = 0.1 (small optical thickness). Right: τs;L = 1 (Earth-like optical thickness). In all cases, the atmosphere is quasi-transparent in the shortwave (τs;S =10−4), pure absorption is assumed (βL0 = βS0 = 1), As = 0.2, and F⋆ = 1366 W m−2. The atmospheric temperature is designated by the solid black line, the dayside surface temperature by the dashed red line, the nightside surface temperature by the dotted blue line, and the black body equilibrium temperature (Eq. (2)) by the dotted grey line. The value of Lsen for a 1-bar CO2-dominated atmosphere is indicated by the green triangle.

In the text

	Fig. 7 Diagram of large-scale stellar and anti-stellar atmospheric circulation as planetary heat engine. The heat engine is driven by dayside heating and nightside cooling to space. In this study, the atmospheric circulation is balancing Rayleigh drag, which causes frictional energy dissipation.
In the text

Fig. 8 Normalised temperature as function of logarithm of parameter Lsen for various values of Ladv. The parameter Ladv, which controls advective heat transport by the stellar and anti-stellar circulation, varies in the range (from dark blue to green). The atmosphere is quasi-transparent in the shortwave (τs;S =10−4), pure absorption is assumed (βL0 = βS0 = 1), As = 0.2, and F⋆ = 1366 W m−2, similarly as in Fig. 6. Here, the atmosphere is optically thin in the longwave: τs;L = 0.1, which corresponds to the middle panels of Fig. 6.

In the text

Fig. 9 Atmospheric, dayside and nightside surface temperatures as functions of control parameters Lsen and Ladv for optically thin (τs;L = 0.1, top panels) and thick (τs;L = 1.0, bottom panels) atmospheres in longwave. Left: Ladv = 0.01 (weak circulation). Middle: Ladv = 1.0 (medium-strength circulation). Right: Ladv = 100 (strong circulation). In all cases, the atmosphere is quasi-transparent in the shortwave (τs;S =10−4), pure absorption is assumed (βL0 = βS0 = 1), As = 0.2, and F⋆ = 1366 W m−2, similarly as in Fig. 6. The dayside atmospheric (solid pink line) and surface (dashed red line) temperatures, as well as the nightside atmospheric (dotted cyan line) and surface (dashed blue line) temperatures are plotted as a function of the logarithm of Lsen. The notation log10 designates the decimal logarithm.

In the text

Fig. 10 Nightside temperature (K) predicted by analytical theory in this work (top panels) and in Wordsworth (2015) (middle panel), and associated comparative stability diagrams (bottom panels) as functions of logarithms of stellar flux (horizontal axis) and surface pressure (vertical axis). Left: purely radiative case as defined in Sect. 3.3 (Lsen = 0, Ladv = +∞). Right: general case (i.e., with sensible heating and heat transport by stellar and anti-stellar atmospheric circulation). The stellar flux is normalised by the modern Earth’s value F⊕ = 1366 W m−2, and the surface pressure is given in bar. In stability diagrams, the acronyms “TW” and “W2015” are used for “This work” and “Wordsworth (2015)”, respectively. Colours indicate the stability of the steady state predicted by the two models: both models predict stability (orange-red areas); the present model predicts stability while W2015 predicts collapse (orange areas); the present model predicts collapse while W2015 predicts stability (light blue areas); both models predict collapse (blue areas). The incident stellar fluxes of TRAPPIST-1 f and a hypothetical tidally locked Venus are designated by green triangles. These plots are obtained by solving Eqs. (78) and (45) of Wordsworth (2015) with parameters values given by Table 2.

In the text

Fig. 11 Comparative stability diagrams as functions of logarithms of stellar flux (horizontal axis) and surface pressure (vertical axis) for various intensities of sensible heating, εsen = 10−4, 10−2, 100 (from bottom to top). Left: comparison between the general case (Lsen > 0) and the purely radiative case characterised in Sect. 3.3 (Lsen = 0). Right: comparison between this work and the model developed by Wordsworth (2015) in thegeneral case. The stellar flux is normalised by the modern Earth’s value F⊕ = 1366 W m−2, and the surface pressure is given in bar. The acronyms “TW”, “RC”, and “W2015” are used for “This work” (with sensible heating), “Radiative case” (without sensible heating), and “Wordsworth (2015)”, respectively. Colours indicate the stability of the steady state predicted by the two models, similarly as in Fig. 10. These plots are obtained by solving Eq. (65) (Ladv = +∞) and Eq. (45) of Wordsworth (2015) with parameters values given by Table 2. Violet dots (top right panel) indicate the region of atmospheric stability computed by Wordsworth (2015) from GCM simulations (Fig. 12), while small blue dots indicate the region of collapse.

In the text

Fig. 12 Comparative stability diagrams for 1 M⊕ and 10 M⊕ planets as functions of logarithms of stellar flux (horizontal axis) and surface pressure (vertical axis). The efficiency parameter controlling advective heat transport εeff take the values (from bottom to top). Left column: case without sensible heating (εsen = 0). Right column: case with efficient sensible heating (εsen = 1.0). The stellar flux is normalised by the modern Earth’s value F⊕ = 1366 W m−2, and the surface pressure is given in bar. The acronyms “E” and “SE” are employed for “Earth” and “Super-Earth”, and designate the 1 M⊕ and 10 M⊕ planets, respectively. Colours indicate the stability of the steady state obtained in the two cases, similarly as in Fig. 10. These plots are obtained by solving Eq. (65) (Ladv = +∞, top panels) and Eq. (78) (other panels) with κL = 8.5 × 10−5 m2 kg−1 and parameters values given by Table 2 for the Earth-sized planet. In the case of the 10 M⊕ planet, the planet’s radius and mass are set to Mp = 10 M⊕ and Rp= 1.88 R⊕, respectively. Data computed by Wordsworth (2015) from GCM simulations (Fig. 12) are included in the middle right panel for comparison and indicate the region of atmospheric stability for the 1 M⊕ planet (red dots) and both planets (violet dots), and the region of atmospheric collapse (small blue dots).

In the text

Fig. 13 Temperature structure and fluxes profiles calculated with radiative-convective-subsiding-upwelling (RCSU) model in reference case for various surface pressures as functions of pressure levels in logarithmic scale. Left: ps = 0.01 bar and τs;L = 0.01. Middle: ps = 0.1 bar and τs;L = 0.1. Right: ps = 1 bar and τs;L = 1. Dayside atmospheric (solid pink line), nightside atmospheric (dashed sky blue line), dayside surface (red square), nightside surface (blue square), and radiative equilibrium (dotted grey line) temperatures (K) are plotted in top panels. Dayside (shades of red) and nightside (shades of blue) radiative net, total, upwelling, and downwelling fluxes in the longwave, and net flux in the shortwave (dashed cyan line) are plotted in bottom panels. Fluxes are normalised by the black body equilibrium flux , where Teq is the black body equilibrium temperature defined by Eq. (2). Calculations are performed using parameters values given by Table 2 with n = 1, which corresponds to Fig. 2 (right panel) of Wordsworth (2015). Complementary parameters of the RCSU model are set to and Asub∕Aup = 10.

In the text

Fig. 14 Stability diagram computed with one-dimensional radiative-convective-subsiding-upwelling (RCSU) model in reference case. The orange area indicate a stable atmosphere, and the blue area collapse in the RCSU model. Data computed by Wordsworth (2015) from GCM simulations (Fig. 12) are included for comparison. Violet dots indicate simulations where the atmosphere remained stable, while small blue dots indicate atmospheric collapse. Calculations are performed using parameters values given by Table 2 with n = 1, which corresponds to Fig. 2 (right panel) of Wordsworth (2015). Complementary parameters of the RCSU model are set to and Asub∕Aup = 10.

In the text

Fig. B.1 Left-hand member of Eq. (78) as function of normalised temperature difference in logarithm scale for Lsen, Ladv = 0.01, 1, 100. The functionis denoted by f. The parameters Lsen and Ladv are thenon-dimensional control parameters of sensible heating and advected heat transport defined by Eqs. (63) and (75), respectively. Parameters values: βL 0 = βS0 = 1, τs;L = 10; τs;S = 0.1; F⋆ = 1366 W m−2; As = 0.2.

In the text