© P. Auclair-Desrotour et al. 2022

1 Introduction

Launched recently from Kourou’s spaceport in French Guiana, the James Webb Space Telescope (Deming et al. 2009) is on the point of unravelling the features of exoplanetary atmospheres at resolutions never before reached. With the current or upcoming transit searches of the TESS (Barclay et al. 2018) and PLATO (Ragazzoni et al. 2016) observatories, this telescope will accelerate the dynamics initiated by previous space missions by populating the continuum of extrasolar planets and constraining the properties of the detected objects. Many of these planets are rocky planets in close-in star–planet systems, notably planets that orbit brown dwarfs and very low-mass stars (e.g. Payne & Lodato 2007; Raymond et al. 2007; Kopparapu et al. 2017) such as the seven Earth-sized planets hosted by the TRAPPIST-1 ultra-cool dwarf star (Gillon et al. 2017; Grimm et al. 2018). Therefore, it is crucial to better understand the mechanisms governing their climate, atmospheric circulation, and surface conditions.

Tidal locking in 1:1 spin–orbit resonance is the most probable final spin state of planets in close-in star-planet systems (Goldreich 1966). This evolution results from the action of the gravitational tides raised by the perturbing tidal potential of the star. Because of dissipative mechanisms, the tidal response of the planet is delayed with respect to the perturber. As a consequence, the resulting tidal torque tends to drive the spin towards the configuration where the star is motionless in the frame of reference rotating with the planet. This spin state corresponds to spin-orbit synchronisation and is reached when the spin angular frequency of the planet, Ω, equalises its orbital frequency, n_⋆.

Additionally, gravitational tides act to decrease both the obliquity and the eccentricity of the planet, which is thereby driven towards the equilibrium configurations of coplanarity and circularity (Hut 1980, Hut 1981) unless it spirals towards the star until being engulfed by it if the system is very close (Hut 1981; Levrard et al. 2009). Asynchronous final spin states may also exist. For instance, eccentric orbits maintained by orbital resonances in a multiple-planet system lead to spin-orbit resonances of higher degrees where the planet can be trapped (Correia et al. 2014; Auclair-Desrotour et al. 2019a). Similarly, it has been shown that significant thermal tides generated by stellar irradiation are able to prevent Venus-like planets from reaching spin–orbit synchronisation by inducing a torque opposed to the solid tidal torque (Gold & Soter 1969); Ingersoll & Dobrovolskis 1978; Leconte et al. 2015; Auclair-Desrotour et al. 2019b).

The probability of a planet being tidally locked in 1:1 spin– orbit resonance with temperate surface conditions is determined from the interplay between two radii: the tidal-lock radius, r_T, and the radius of the habitable zone, r_HZ. While the tidal-lock radius indicates the size of the region where planets are likely to be tide-locked in spin–orbit synchronisation, the radius of the habitable zone corresponds to the typical star-planet distance at which a planet can sustain liquid water at its surface (Kopparapu et al. 2013). By assuming that the planet behaves as a black body, and by writing the stellar luminosity as a function of the stellar mass with the empirical formula given by Barnes et al. (2008), it can be shown that $r_{\text{HZ}}\propto M_{*}^{1.32}\text{for\hspace{0.17em}}{M}_{*}\lesssim 1$ (Auclair-Desrotour & Heng 2020), whereas the tidal-lock radius scales as $r_{\text{T}}\propto M_{*}^{1/3}$ (Peale 1977; Kasting et al. 1993; Dobrovolskis 2009; Edson et al. 2011). Thus, the size of the habitable zone radius decays faster than the tidal-lock radius with decreasing stellar mass, which means that planets located in the habitable zone have a greater chance of being tide-locked if they orbit low-mass stars than if they orbit Sun-like stars (Kasting et al. 1993).

For planets orbiting low-mass M stars, tide-locking times are actually very short, and even extremely short in the case of lava planets (e.g. 55 Cancri e, Kepler 10b), with maximum values reaching just a few million years. For instance, the time required for the cool planet LHS1140 b to become tide-locked is about 14 million years (Pierrehumbert & Hammond 2019), which is low compared with the typical ages of planetary systems. This strongly suggests that planets orbiting in the habitable zone of very low-mass stars such as TRAPPIST 1d are tide-locked in the 1:1 spin–orbit resonance. Therefore, the rotation rate of these planets is well constrained and is given by Ω = n_⋆. Additionally, the strength of tidal forces makes the existence of orbital configurations with significant eccentricities or obliquities unlikely. Such planets can thus be reasonably supposed to be close to the stable equilibrium configurations of coplanarity and circularity.

In this configuration, the planet exhibits a permanent day-side and nightside centred on the star–planet axis. The dayside is irradiated by the incident stellar flux while the nightside is radiatively cooled, and the energy absorbed on the dayside is transported towards the nightside by mean flows, which act to decrease the horizontal temperature gradient (e.g. Showman & Guillot 2002; Leconte et al. 2013; Pierrehumbert & Hammond 2019). The differential thermal forcing induced by spin–orbit synchronisation plays a crucial role in the evolution of the planet’s climate. Typically, as shown by the pioneering work of Joshi et al. (1997), the thermal state of synchronously rotating rocky planets is determined by the interplay between the night-side surface temperature and the condensation temperature of greenhouse gases. From the moment that the condensation temperature of the gas exceeds the surface temperature, the nightside acts as a cold trap. The greenhouse gas initially present in the atmosphere condenses and forms an ice sheet on the surface, which induces a temperature decrease in return. This triggers a positive feedback that cools down the atmosphere until the gas has been fully condensed, which is called atmospheric collapse (e.g. Joshi et al. 1997; Heng & Kopparla 2012; Wordsworth 2015). In the opposite case, the atmosphere is said to be stable against collapse, its composition remaining unchanged.

Over the past decade, a substantial effort has been made to characterise the climate of tide-locked rocky planets using a broad range of modelling approaches, from simplified analytical greenhouse models (e.g. Heng & Kopparla 2012; Wordsworth 2015; Auclair-Desrotour & Heng 2020) to 3D general circulation model (GCM) simulations (e.g. Merlis & Schneider 2010; Leconte et al. 2013; Carone et al. 2014; Haqq-Misra et al. 2018; Ding & Pierrehumbert 2020; Sergeev et al. 2020; Turbet et al. 2018), including intermediate semi-analytical or numerical approaches (e.g. Yang & Abbot 2014; Koll & Abbot 2016; Song et al. 2021) that cannot be listed here in an exhaustive way. Although based on robust methodologies, most of these approaches cannot be related self-consistently to one anotherdue to major differences in modelling choices. These discrepancies raise the two questions of how to disentangle the possible causes of different predictions between two models and how to assess the epistemic value of a given model. They can be reformulated into the more concrete question of how to consistently characterise the climate of tide-locked planets from multiple modelling approaches. This major concern was explicitly formulated by Held (2005), who argued for the need of model hierarchies on which to base one’s understanding in climate modelling. Such hierarchies appear as the only way to close the gap between idealised modelling and high-end simulations as they allow the essence of each particular source of complexity to be captured.

The aim of the present work is to tackle these questions from the angle of atmospheric stability against collapse. In the continuity of a former study on the atmospheric stability of tide-locked rocky planets (Auclair-Desrotour & Heng 2020), we have developed a multi-dimensional model hierarchy that we call a general circulation meta-model (GCMM) in order to bridge the gap between the analytic theory of planetary climates and simulations performed with 3D GCMs. This model hierarchy is based on a systematic bottom-up approach in the spirit of Held (2005).

We need to specify the sense given here to meta-modelling. By meta-model, we mean that the model ought to be able to exactly reproduce the setups of both simplified greenhouse models and GCMs – as well as the configurations in between – with the same intrinsic theoretical background. In that sense, such models are possible instances of the meta-model, which can generate any of them. Hence, the so-defined GCMM allows the effects of mechanisms that are either strongly coupled in standard GCMs or ignored in simplified analytic models to be disentangled. These effects are added or subtracted as a function of the number of degrees of freedom of the model. Increasing the number of degrees of freedom amounts to adding key sources of complexity.

Typically, radiative 0D models are essentially based on radiative exchanges between the planet’s surface and the atmosphere (e.g. Wordsworth 2015; Auclair-Desrotour & Heng 2020). At the next level of complexity, 1D models take the coupling between radiative transfer and the atmospheric structure into account (e.g. Robinson & Catling 2012). Two-column – or 1.5D – models are the minimum setup to self-consistently couple the large-scale day-night overturning circulation with radiative transfer and the atmospheric structure (e.g. Yang & Abbot 2014; Koll & Abbot 2016). This coupling is refined at the level of 2D GCMs, which allow the interaction between physical mechanisms – clouds, turbulent diffusion in the planetary boundary layer (PBL), convection – and mean flows to be calculated self-consistently in the slow rotation regime (e.g. Song et al. 2021). Finally, 3D GCMs complete the picture by introducing Coriolis effects and non-axisymmetric flows where super-rotation can develop (e.g. Leconte et al. 2013; Carone et al. 2014; Haqq-Misra et al. 2018; Ding & Pierrehumbert 2020; Sergeev et al. 2020; Turbet et al. 2021).

Thus, the essential function of a GCMM is to model all these levels of complexity at the same time so that the roles played by the different mechanisms involved in the planet’s climate can be clearly separated. As, to our knowledge, such a model has not yet been developed, the present work is a first attempt to design a GCMM dedicated to the study of tide-locked rocky planets. For simplicity, we confine ourselves to the slow rotation regime and ignore Coriolis effects in the dynamics. This allows us to avoid the mathematical complications related to the 3D geometry and to speed up calculations. Our GCMM is thus designed to generate models ranging from 0D configurations to 2D configurations. Additionally, we opt for a finite-volume method to solve the hydrostatical primitive equations (HPEs), following the approaches detailed by Yao & Stone (1987) and implemented in standard finite-volume GCMs such as the LMDZ (Hourdin et al. 2006) or THOR (Mendonça et al. 2016) GCMs. As the finite-volume method divides the atmosphere into cells, it is appropriate to describe the radiative energy balance models on which the analytic theory is built. For instance, the one-cell configuration (1 × 1 grid) corresponds to the single-layer isothermal atmosphere of Wordsworth’s model (Wordsworth 2015). Similarly, increasing numbers of horizontal and vertical grid intervals generate the 1D (1 × 50 grid), 1.5D (2 × 50 grid), and 2D (32 × 50 grid) model setups.

In order to minimise the size of the parameter space, we confine ourselves to the dry case in this study. The effects of moisture (clouds, latent heat transport, sedimentation, and surface condensation or evaporation) are ignored. Radiative transfer is described in the double-grey approximation, meaning that the fluxes are divided into two bands – shortwave and longwave – characterised by effective absorption parameters (e.g. Sect. 4.1 of Heng 2017). We also make the two-stream approximation and consider that radiative fluxes only travel upwards and downwards (e.g. Heng 2017, Sect. 3.1). In addition to radiative transfer, the vertical turbulent diffusivity induced by the interactions between mean flows and the planet surface in the PBL is taken into account by making use of a model based on the mixing length theory (Holtslag & Boville 1993). Finally, the thermal diffusion in the soil is included in the diffusion scheme of the GCMM with a 1D finite-difference model following the method described by Wang et al. (2016). As shown by earlier studies (Wordsworth 2015; Koll & Abbot 2016; Auclair-Desrotour & Heng 2020), the three abovementioned physical ingredients (circulation, radiative transfer, and turbulent diffusion) predominantly determine the nightside surface temperature and, thereby, the atmospheric stability against collapse. Table 1 summarises the physics described by the studied instances of the meta-model and the THOR 3D GCM, with the latter used to benchmark the former. Physical mechanisms are gradually captured by the grid formats that characterise the models as the number of degrees of freedom increases.

In Sect. 2, we introduce some of the control parameters and analytical scaling laws that characterise the climate and circulation regime of tide-locked planets. In Sect. 3 we detail the main features of the GCMM and the physical setup of the studied Earth-like and pure CO₂ atmospheres. Section 4 introduces the four instances of the meta-model used in this work: 0D, 1D, 1.5D, and 2D. In Sect. 5, we run grid simulations for these instances to characterise the atmospheric stability of Earth-sized synchronous planets against collapse as a function of the stellar flux and surface pressure. Particularly, this vertical inter-comparison highlights the influence of the PBL on climate, day-night advection, and surface conditions. In Sect. 6, we investigate the limitations of the zero-spin rate approximation assumed in this approach by running simulations with the THOR 3D GCM. We show that the approximation holds from the moment that the dimensionless equatorial Rossby deformation length is greater than 2. Finally, in Sect. 7 we summarise the conclusions of the study.

2 Preliminary scalings

The circulation regime and thermal state of equilibrium of tide-locked planets is controlled by a few parameters and scaling laws derived either from the shallow water approximation (e.g. Vallis 2006; Pierrehumbert & Hammond 2019) or from the weak temperature gradient (WTG) approximation (e.g. Pierrehumbert 1995; Sobel et al. 2001). One ought to recall these scalings before introducing the physical setup of the numerical approach.

2.1 Circulation regimes of synchronous planets

If we assume that the planet surface is isotropic, all the physics and dynamics that govern the atmospheric circulation are symmetric about the star–planet axis except Coriolis terms. As a consequence, the circulation regime is essentially characterised by one control parameter depending on the planet’s spin angular velocity, which determines whether mean flows are bidimensional and symmetric about the star–planet axis or if they are sufficiently deviated by the planet’s rotation as to become 3D.

Such a parameter appears naturally in analyses making use of the Buckingham-Pi theorem (Buckingham 1914) in the primitive equations of fluid dynamics and thermodynamics (e.g. Koll & Abbot 2015). In the present study, following Leconte et al. (2013) and Auclair-Desrotour & Heng (2020), we define the dimensionless equatorial Rossby deformation length ${\tilde{L}}_{\text{Ro}}$ from the equatorial Rossby deformation radius L_Ro as (e.g. Menou et al. 2003) $${\tilde{L}}_{\text{Ro}}\equiv \frac{L_{\text{Ro}}}{R_{\text{p}}}=\sqrt{\frac{c_{\text{wave}}}{2\text{Ω}{R}_{\text{p}}}},$$(1)

where R_p designates the planet radius, and c_wave the speed of horizontally propagating gravity waves. The dimensionless equatorial Rossby deformation length is essentially the square root of the distance – in radius unit – that fast gravity waves can travel before they feel the Coriolis effects and geostrophically adjust (Vallis 2006).

If ${\tilde{L}}_{\text{Ro}}>1$, the Coriolis effects are small and the two-way force balance between advection and pressure-gradient accelerations leads to a day-night overturning circulation symmetric about the star–planet axis (Leconte et al. 2013; Pierrehumbert & Hammond 2019; Hammond & Lewis 2021). In this regime, the dynamics of mean flows is the same in all planes containing the star–planet axis, with high-altitude winds blowing from the day-side to the nightside and near-surface winds blowing from the nightside to dayside. This essentially corresponds to the steady state expected in the WTG theory, where small Coriolis forces, friction, and non-linearities make the heat advection unable to annihilate completely the day-night temperature gradient (Sobel et al. 2001; Hammond & Pierrehumbert 2017; Pierrehumbert & Hammond 2019).

Conversely, for ${\tilde{L}}_{\text{Ro}}\lesssim 1$, Showman & Polvani (2011) demonstrated that 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) favours the emergence of super-rotation by pumping angular momentum from the mid-latitudes towards the equator. In this regime, the equatorial Rossby deformation radius (L_Ro) essentially corresponds to the latitudinal width of the produced eastward equatorial jet, and mean flows take the form of the Matsuno-Gill standing wave pattern (Matsuno 1966b; Gill 1980; Showman & Polvani 2011; Tsai et al. 2014).

The dimensionless equatorial deformation length introduced in Eq. (1) can be related to the atmospheric parameters by considering the properties of gravity waves. Gravity waves are restored by the Archimedean force associated with the fluid buoyancy and their typical speed is given by c_wave = HN, where H designates the characteristic vertical scale length of the atmosphere and N the Brunt-Vaisala frequency, which scales the strength of the atmospheric stratification against convection (Gerkema & Zimmerman 2008). In a dry, stably stratified atmosphere, this frequency is expressed as (Gerkema & Zimmerman 2008) $$N^2=\frac{g}{T}\left(\frac{g}{C_p}+\frac{\partial T}{\partial z}\right),$$(2)

where we have introduced the gravity g, the heat capacity per unit mass of the gas C_p, the temperature T, and the partial derivative operator over the altitude $\frac{\partial }{\partial z}$. In the idealised case of the vertically isothermal atmosphere (constant temperature), $N=g/\sqrt{C_{p}T}$, and the vertical scale is the pressure height (Vallis 2006, Sect. 1.4), $$H\equiv \frac{R_{\text{d}}T}{g},$$(3)

where R_d ≡ R_GP/M_a designates the specific gas constant for dry air, R_GP being the ideal gas constant and M_a the mean molecular weight of the atmosphere. Thus, in this configuration (e.g. Leconte et al. 2013), $${\tilde{L}}_{\text{Ro}}=\sqrt{\frac{R_{\text{d}}{T}^{1/2}}{2\text{Ω}{R}_{\text{p}}{C}_p^{1/2}}},$$(4)

which highlights the fact that the circulation regime depends on the planet’s spin rotation, thermal state, and atmospheric composition.

The dimensionless equatorial Rossby deformation length conveys exactly the same information as the WTG parameter Λ introduced in the WTG theory (see e.g. Pierrehumbert & Hammond 2019), which is defined as the Rossby radius of deformation – distinct from the equatorial deformation radius – normalised by the planet radius (Vallis 2006, Sect. 3.8.2). The two parameters are linked together by the relationship ${\tilde{L}}_{\text{Ro}}=\sqrt{\Lambda /2}$ (Pierrehumbert & Hammond 2019), meaning that either the former or the latter can be chosen to characterise the circulation regime. In the present study, we confine ourselves to the configuration of the WTG theory $\left({\tilde{L}}_{\text{Ro}}>1\right)$ and consider thereby that mean flows are symmetric about the star–planet axis.

In addition with the slow and fast rotators regimes, there exists a third dynamical state that is proper to intermediate stellar cases in the range of 3000–3300 K and that is described as the Rhines rotation regime (Haqq-Misra et al. 2018; Sergeev et al. 2020). This regime is related to the Rhines length, which determines the maximum extent of zonally elongated turbulent structures (Rhines 1975). It occurs when the non-dimensional Rossby deformation radius is greater than one but the non-dimensional Rhines length is less than one (Haqq-Misra et al. 2018). The slow rotation and fast rotation regimes occur when both the non-dimensional Rhines length and Rossby deformation radius are greater than or less than one, respectively. The Rhines rotation regime is not considered here, meaning that we focus on the configuration where both the non-dimensional Rhines length and Rossby deformation radius are greater than one.

2.2 Thermal states predicted by radiative box models

Over the past decade, analytic solutions and scalings characterising the thermal state of equilibrium of tide-locked planets have been obtained both for hot Jupiters (Komacek & Showman 2016); Zhang & Showman 2017; Koll & Komacek 2018), lava planets (Hammond & Pierrehumbert 2017), and cooler rocky planets orbiting in the habitable zone of their host star (Showman et al. 2013; Wordsworth 2015; Koll & Abbot 2016; Pierrehumbert & Ding 2016; Koll 2022; Auclair-Desrotour & Heng 2020). Most of them were derived in the framework of the WTG theory and involve simplified atmospheric physics and structure. The present study builds on these works, and particularly those based on box model approaches, where the atmosphere and surface are reduced to large-scale energy reservoirs exchanging heat with each other (Wordsworth 2015; Koll & Abbot 2016; Auclair-Desrotour & Heng 2020). Although they are based on strong simplifications (isothermal atmosphere, large-scale averages, no self-consistent coupling between the dynamics and the thermodynamics), these models provide scalings that capture the behaviour of the thermal state of tide-locked rocky planets with a minimum set of physical parameters. Particularly, they lead to closed-form solutions for the nightside surface temperature T_n, which determines the whole atmospheric stability against collapse.

By considering a globally isothermal atmosphere interacting with dayside and nightside surface hemispheres, Wordsworth (2015) shows that the pure radiative equilibrium of the surface-atmosphere system corresponds, in the optically thin layer approximation (i.e. transparent in the visible and optically thin in the infrared), to the nightside temperature scaling (Wordsworth 2015, Eq. (29)) $$T_{\text{n;low}}\equiv T_{\text{eq}}{\left[\frac{1}{2}\left(1-A_{\text{s}}\right){\tau }_{\text{s;L}}\right]}^{\frac{1}{4}},$$(5)

which can be generalised to optically thick atmospheres with scattering (Auclair-Desrotour & Heng 2020). In the above expression, A_s designates the surface albedo in the shortwave, τ_s;L the longwave optical depth at planet’s surface, and T_eq the black body equilibrium temperature, which is defined as $$T_{\text{eq}}\equiv {\left(\frac{F_{\text{*}}}{4{\sigma }_{\text{SB}}}\right)}^{\frac{1}{4}},$$(6)

with F_⋆ the incident stellar flux and σ_SB = 5.670367 × 10⁻⁸ W m⁻² K⁻⁴ the Stefan-Boltzmann constant (Mohr et al. 2016).

Since it ignores all types of energy exchanges except radiative transfer, the estimate given by Eq. (5) can be interpreted as a lower bound for the nightside surface temperature of a rocky tide-locked planet. In reality, the strong convection generated by the thermal forcing of the atmosphere in the dayside PBL increases surface-atmosphere heat fluxes, which significantly affects the nightside temperature (Sergeev et al. 2020). The friction of the flow against the surface generates sensible heat exchanges in dry thermodynamics. Additionally, in moist atmospheres, surface evaporation generates latent heat exchanges, which results from the energy taken from or released in the fluid during the changes of phases of the component (Pierrehumbert 2010).

We consider here that the atmosphere is dry and thereby ignore the contribution of latent heat exchanges. Sensible heat exchanges can be introduced in the radiative box model by including, in the energy balance equations, the hemisphere-averaged sensible heat flux given by (e.g. Pierrehumbert 2010, Eq. (6.11), p. 396) $${\overline{F}}_{\text{sen}}=C_{\text{D}}{C}_p{\rho }_{\text{a}}{\upsilon }_{\text{sen}}\left(T_{\text{d}}-T_{\text{a}}\right),$$(7)

where ρ_a is the atmospheric density at planet’s surface, C_D the bulk drag coefficient characterising the strength of friction in the surface layer, and υ_sen the typical horizontal wind speed of the flow. Among these parameters, υ_sen accounts for the circulation, meaning that it cannot be self-consistently related to the thermal state of the system in this simplified approach. Nevertheless, it can be scaled from a dimensional analysis of the thermodynamic equation (e.g. Wordsworth 2015), or by making use of the heat engine theory (e.g. Koll & Abbot 2016; Koll & Komacek 2018; Auclair-Desrotour & Heng 2020). For instance, by modelling the overturning circulation as an ideal heat engine and using Carnot’s theorem (Bruhat 1968), Koll & Abbot (2016) found an upper bound of the dayside average surface wind speed (Koll & Abbot 2016, Eq. (12)), $${\upsilon }_{\text{sen}}={\left\{\left[T_{\text{d}}-{\left(1-A_{\text{s}}\right)}^{\frac{1}{4}}{T}_{\text{eq}}\right]\left(1-e^{-{\tau }_{\text{s;L}}}\right)\left(1-A_{\text{s}}\right)\frac{R_{\text{d}}{F}_{\text{*}}}{2C_{\text{D}}{p}_{\text{s}}}\right\}}^{\frac{1}{3}},$$(8)

which agrees well with the values obtained numerically from 3D GCM simulations (Koll & Abbot 2016); Koll & Komacek 2018).

For an isentropic cycle (i.e. an idealised Carnot’s heat engine), the weight of dayside sensible heating¹ relative to radiative heating is controlled by the dimensionless parameter (Auclair-Desrotour & Heng 2020, Eq. (63)) $$L_{\text{sen}}\equiv \frac{2C_p{C}_{\text{D}}{p}_{\text{s}}}{{\tau }_{\text{s;L}}{R}_{\text{d}}{F}_{\text{*}}}{\left[\frac{Q_{\text{in}}{R}_{\text{d}}}{C_{\text{D}}{p}_{\text{s}}}{\left(\frac{F_{\text{*}}}{2{\sigma }_{\text{SB}}}\right)}^{\frac{1}{4}}\right]}^{\frac{1}{3}},$$(9)

which is written here for an atmosphere optically transparent in the shortwave and thin in the longwave. The notation $Q_{\text{in}}\propto {\upsilon }_{\text{sen}}^3$ (e.g. Koll & Abbot 2016) designates the amount of power per unit area available to drive atmospheric motion. Looking at the zero-convection limit (L_sen = 0) we recover the purely radiative regime, while the opposite asymptotic regime (L_sen = +∞) implies that T_d = T_a and provides an upper bound for the night-side surface temperature of a tide-locked rocky planet (e.g. Auclair-Desrotour & Heng 2020, Eq. (85)), $$T_{\text{n;up}}\equiv T_{\text{eq}}{\left[2\left(1-A_{\text{s}}\right){\tau }_{\text{s;L}}\right]}^{\frac{1}{4}}=\sqrt{2}{T}_{\text{n;low}},$$(10)

which is valid in the optically thin layer approximation as well as Eq. (5). Therefore, for a globally isothermal and optically thin atmosphere, the nightside surface temperature falls into the interval $$T_{\text{n;low}}\le T_{\text{n}}\le \sqrt{2}{T}_{\text{n;low}}.$$(11)

However, we remark that T_d = T_a actually corresponds to an extreme regime that is never reached in standard configurations, and we thus consider T_n;up as a theoretical upper limit.

Similarly as the dayside convective planetary layer, the nightside atmospheric structure alters the nightside equilibrium temperature. Its effect can be quantified by relaxing the isothermal approximation and by dividing the atmosphere into dayside and nightside air columns, which is the essence of two-column models (Yang & Abbot 2014); Koll & Abbot 2016; Auclair-Desrotour & Heng 2020). As shown by Koll & Abbot (2016), the nightside subsidence induced by the day-night overturning circulation generates a temperature inversion in the lowest layers of the atmosphere if the subsidence timescale is slightly greater than the radiative cooling timescale. The resulting atmospheric structure leads to large day-night differences.

3 A general circulation meta-model (GCMM)

We introduce in this section the main features of the meta-model and the used physical setup.

3.1 Primitive equations

At a given time, t, the dynamical core of the GCMM solves the HPEs over the Cartesian rectangular domain defined by the colatitude θ ∈ [0°, 180°] of the tidally locked coordinates (TLCs; see Koll & Abbot 2015, Appendix B) and the mass-based vertical coordinate defined, in the absence of the topography, as (e.g. Yao & Stone 1987; Carone et al. 2014) $$\sigma \equiv \frac{p-p_{\text{t}}}{p_{\text{s}}-p_{\text{t}}}\in \left[0,1\right],$$(12)

where we have introduced the pressure p, the surface pressure p_s, and the pressure at the top of the atmosphere p_t. In these coordinates, θ = 0° and θ = 180° correspond to the sub-stellar and anti-stellar points, respectively, while σ = 0 and σ = 1 correspond to the top and the bottom of the atmosphere, respectively. While p_s and p vary over time and spatial coordinates, p_t is set to p_t = 0 in the model, which corresponds to the usual sigma coordinate σ = p/p_s. The vertical coordinate given by Eq. (12) is well suited to the study of the tide-locked planets since it follows the distortion of the atmosphere induced by the differential day-night thermal forcing: the pressure of an altitude level may differ by orders of magnitude between the dayside and nightside, which would possibly generate numerical issues with the altitude coordinate.

The relationship between the altitude (z) and the generalised vertical coordinate (σ) is contained in the so-called pseudodensity (e.g. Kasahara 1974), $$\mathfrak{p}\equiv -g\rho \frac{\partial z}{\partial \sigma },$$(13)

where ρ denotes the density. The pseudo-density is proportional to the mass contained in a generalised volume where the vertical dimension is not a length but an interval of the generalised coordinate σ. With the chosen mass-based coordinate (Eq. (12)) and the assumed hydrostatic balance, it is simply expressed as (e.g. Yao & Stone 1987) $$\mathfrak{p}=p_{\text{S}}-p_{\text{t}}.$$(14)

We note that p would be the density to a constant factor if the vertical coordinate were the altitude. Using the pseudo-density instead of the density allows the HPEs given further in the same form to be written for any chosen vertical coordinate.

The HPEs are the mass continuity equation (e.g. Kasahara 1974), $$\frac{\partial \mathfrak{p}}{\partial t}+{\nabla }_{\sigma }\cdot \left(\text{p}{\upsilon }_{\sigma }\right)+\frac{\partial }{\partial \sigma }\left(\mathfrak{p}\dot{\sigma }\right)=0,$$(15)

the horizontal momentum equation, $$\begin{array}{l}\frac{\partial }{\partial t}\left(\mathfrak{p}{\upsilon }_{\theta }\sin \theta \right)+\frac{1}{R_{\text{p}}}\frac{\partial }{\partial \theta }\left(\mathfrak{p}{\upsilon }_{\theta }^2\sin \theta \right)\hfill \\ +\frac{\partial }{\partial \sigma }\left(\mathfrak{p}{\upsilon }_{\theta }\dot{\sigma }\sin \theta \right)+\frac{1}{R_{\text{p}}}\mathfrak{p}\sin \theta \left[\frac{\partial \varphi }{\partial \theta }+\text{Θ}\frac{\partial E}{\partial \theta }\right]=\mathfrak{p}\sin \theta F_{\theta }.\hfill \end{array}$$(16)

the potential temperature equation, $$\frac{\partial \mathfrak{p}\text{Θ}}{\partial t}+{\nabla }_{\sigma }\cdot \left(\mathfrak{p}\text{Θ}{\upsilon }_{\sigma }\right)+\frac{\partial }{\partial \sigma }\left(\mathfrak{p}\text{Θ}\dot{\sigma }\right)=\mathfrak{p}\frac{Q}{E},$$(17)

and the hydrostatic equation combined with the ideal gas law, $$\frac{\partial \varphi }{\partial \sigma }+\text{Θ}\frac{\partial E}{\partial \sigma }=0,$$(18)

where we have introduced the horizontal velocity vector υ_σ ≡ υ_θ e_θ, the sigma-velocity $\dot{\sigma }\equiv \frac{\text{d}\sigma }{\text{d}t}$ (with $\frac{\text{d}}{\text{d}t}$ the material time derivative), the geopotential ϕ ≡ gz, the potential temperature Θ, the Exner function E (e.g. Vallis 2006, Sect. 3.9), and the horizontal divergence operator at constant σ, $$\nabla \cdot \equiv \frac{1}{R_{\text{p}}}{e}_{\theta }\cdot \frac{\partial }{\partial \theta }.$$(19)

The potential temperature and Exner function are defined as $$\begin{array}{ll}\text{Θ}\equiv T{\left(\frac{p}{p_{\text{ref}}}\right)}^{-K},\hfill & E\hfill \end{array}\equiv C_p\frac{T}{\text{Θ}}=C_p{\left(\frac{p}{p_{\text{ref}}}\right)}^K,$$(20)

where T is the temperature, p_ref a constant reference pressure set to p_ref = 1 bar, and κ ≡ R_d/C_p. We note that the Exner function is a proxy for pressure. It is used for convenience, as it facilitates the integration of the hydrostatic equation. In right-hand members of Eqs. (15)–(18), the source-sink terms are the force per unit mass and the heat power per unit mass Q. We note that the primitive equations are written in their conservative forms, which involve mass flows and mass-integrated quantities rather than the original quantities themselves. Besides, these equations are given here for any vertical coordinate varying monotonically with altitude for the sake of generality.

The non-dimensional HPEs derived from Eqs. (15)–(18) (see Appendix A) are solved for {p, υ_θ, Θ, ϕ} on a staggered Arakawa C grid (Arakawa & Lamb 1977) with uniformly spaced horizontal intervals, and σ-dependent vertical intervals refined near the model bottom and top (see Appendix B). Following the method implemented in the LMDZ GCM (Hourdin et al. 2006), the integration is done using a leapfrog scheme with a periodic predictor/corrector timestep. The source-sink terms associated with the physics, {F_θ, Q}, as well as other physical variables, are updated periodically using implicit schemes except radiative transfer equations, which are solved directly from the current thermodynamical state. In the 2D GCM-like configurations, integrating the HPEs on a discrete domain generates numerical instabilities that develop at grid scale and may strongly disrupt calculations. In GCMs, this concern is usually handled by introducing horizontal hyper-diffusion, which is ideally designed to efficiently damp the numerical instabilities at grid scale while leaving the mean flows unchanged (e.g. Thrastarson & Cho 2011). Therefore, we include in the model a fourth-order hyper-diffusion (or bi-harmonic diffusion) using an anisotropic super-diffusivity that vanishes at the sub-stellar and anti-stellar points (see Appendix C.1) in order to avoid the stability concerns associated with isotropic diffusion near the poles (Sect. 13.3 of Lauritzen et al. 2011). The corresponding hyper-diffusion terms for the momentum and temperature equations are given by $$F_{\text{diff};\upsilon }=-K_4{\sin }^{2\alpha }\left(\theta \right){\nabla }_{\sigma }^4{\upsilon }_{\theta },$$(21) $$F_{\text{diff};T}=-K_4{\sin }^{2\alpha }\left(\theta \right){\nabla }_{\sigma }^{4}T,$$(22)

where ${\nabla }_{\sigma }^4\equiv {\nabla }_{\sigma }^2{\nabla }_{\sigma }^2$ designates the second-order horizontal hyper-Laplacian operator, α the anisotropy exponent (α = 1 in the model), and K₄ the super-diffusivity defined by Eq. (C.5). Validation test simulations were run to verify that the mean flow and temperature distribution are insensitive to the hyper-diffusion scheme (see Fig. C.1).

In very hot cases, exponentially growing gravity waves propagating upwards and reflected by the upper boundary of the model can lead to extreme fluctuations of the dynamical quantities in the upper atmosphere. To address these instabilities, it can be necessary to use a sponge layer in addition with horizontal hyper-diffusion. In the present work, we introduce a linear Rayleigh friction sponge (e.g. Lauritzen et al. 2011, Sect. 13.4.5) following the formulation proposed by Polvani & Kushner (2002) for the vertical profile of the Rayleigh coefficient (see Appendix C.2). The sponge layer is thus modelled by a Rayleigh damping increasing with the altitude between a critical level σ = σ_SL and the top of the atmosphere, σ = 0, which tends to annihilate horizontal winds in the vicinity of the upper boundary. Finally, the strong convection generated by the thermal forcing on the dayside can induce super-adiabatic vertical temperature gradients (i.e. ∂Θ/∂z < 0), which is another source of numerical errors. This behaviour can be prevented by introducing a convective adjustment scheme in the model. The convective adjustment scheme regularises the atmospheric structure by dynamically correcting the tendencies to adjust super-adiabatic regions to adiabatic profiles. We implemented a standard scheme (e.g. Hourdin et al. 1993; Mendonça & Buchhave 2020) that can be activated when necessary (see Appendix D). However, we did not have to use the convective adjustment scheme nor the sponge layer in the simulations performed for this work.

3.2 Radiative transfer

In the model, radiative transfer is described through the double-grey approximation, which consists in (i) decoupling the stellar radiation (shortwave flux) and the planet radiation (longwave flux) – each band being characterised by an effective optical depth – and (ii) assuming that radiative fluxes only propagate upwards and downwards, which is the essence of the two-stream approximation (e.g. Heng 2017, Sects. 3.1 and 4.1). This allows fluxes in the shortwave and in the longwave to be calculated separately (see Appendix E). We denote upward and downward fluxes by F_↑ and F_↓ respectively. The equations governing the propagation of the wavelength-integrated total flux F₋ ≡ F_↑ + F_↓ and net flux F₊ ≡ F_↑ + F_↓ have the same formulation in both bands. They are written as (Heng 2017) $$\frac{\text{d}{F}_{+}}{\text{d}\tau }=\frac{1}{{\beta }_0}{F}_{-},$$(23) $$\frac{\text{d}{F}_{-}}{\text{d}\tau }={\beta }_0\left(F_{+}-2B\right),$$(24)

where B ≡ σ_SBT⁴ designates the black body radiation of the gas, which is zero in the shortwave since the atmosphere is assumed to radiate in the infrared only.

In these equations, τ designates the optical depth of the associated band, and β₀ the scattering parameter (β₀ = 1 for pure absorption; 0 < β₀ < 1 in the presence of scattering). The fluxes equations given by Eqs. (23) and (24) are solved numerically in the code by computing the transmission functions of each atmospheric layer as a first step, and by solving the boundary condition problem as a second step. We note that the solution obtained numerically for the single-layer atmosphere this way exactly corresponds to that derived in radiative box models based on the isothermal atmosphere approximation (see e.g. Wordsworth 2015; Auclair-Desrotour & Heng 2020). The optical depths in the shortwave τ_S and longwave τ_L are both assumed to scale linearly with pressure, and are defined as $$\begin{array}{cc}{\tau }_{\text{S}}\equiv \frac{k_{\text{S}}p}{g},& {\tau }_{\text{L}}\equiv \frac{k_{\text{L}}p}{g}\end{array},$$(25)

where we have introduced the effective absorption coefficients of the gas in the short- and longwave, κ_S and κ_L, respectively.

We remark that the optical depths defined by Eq. (25) depict horizontally averaged vertical profiles rather than local profiles varying as functions of the propagation angles of radiative fluxes, α_S and α_L. Therefore, the effective absorption coefficients κ_S and κ_L introduced in Eq. (25) are related both to the absorption properties of the gas and to the mean cosine of the propagation angles in the visible and in the infrared, $\overline{\cos {\alpha }_{\text{S}}}$ and $\overline{\cos {\alpha }_{\text{L}}}$. Typically, these coefficients are related to Wordsworth’s absorption coefficients (Wordsworth 2015, Eq. (12)) – denoted by ${\kappa }_{\text{S}}^{\text{W}}$ and ${\kappa }_{\text{L}}^{\text{W}}$ – by the equations² ${\kappa }_{\text{S}}={\kappa }_{\text{L}}^{\text{W}}/\overline{\cos {\alpha }_{\text{S}}}$ and ${\kappa }_{\text{L}}={\kappa }_{\text{L}}^{\text{W}}/\overline{\cos {\alpha }_{\text{L}}}$. Therefore, changing the value of the mean cosine of the propagation angle in this definition amounts to changing the value of the absorption coefficient in Eq. (25). One shall also bear in mind that the absorption coefficients defined in the double-grey approach are not fundamental parameters of the gas but parameters that mimic the averaged effect of highly frequency-dependent atmospheric opacities, as shown by Wordsworth et al. (2010a) for CO₂-dominated atmospheres.

Although it does not fully account for the complex physics of radiative transfer, the adopted double-grey approximation with average optical depths captures the dependence of optical depths on pressure, and it allows for fast numerical computation of radiative fluxes. The radiative transfer scheme may be refined in future works by using more sophisticated approaches such as the correlated-k distribution method (e.g. Lacis & Oinas 1991), but this goes beyond the scope of the present study.

3.3 Planetary boundary layer

In the PBL the shear instability generates turbulence, which acts to mix the flow. This turbulent mixing induces a vertical diffusion of momentum, heat, and potentially moisture in moist cases, near the planet’s surface. The associated eddy diffusivities are controlled by the gradient Richardson number Ri defined as (Vallis 2006) $$\text{Ri}\equiv \frac{g\frac{\partial \text{Θ}}{\partial z}}{\text{Θ}{\left(\frac{\partial {\upsilon }_{\theta }}{\partial z}\right)}^2},$$(26)

which characterises fluid stratification. The Richardson number is essentially the ratio of the production of turbulent energy due to the shear instability over the restoring force induced by buoyancy. For a given quantity X = υ_θ, C_pΘ (or q in moist cases, q being the specific humidity), the vertical diffusion equation is written, in the gradient-flux theory (e.g. Garratt 1994), as $$\frac{\text{d}X}{\text{d}t}=\frac{1}{\rho }\frac{\partial }{\partial z}\left[\rho K_X\frac{\partial X}{\partial z}\right],$$(27)

the upward diffusive flux being given by $$F_{\text{diff};X}=-\rho K_X\frac{\partial X}{\partial z}.$$(28)

In these equations, K_X is the eddy diffusivity associated with turbulent mixing for X. This parameter is a function of the mean fields. The lower boundary condition is a continuity condition determined by the exchanges with the planet’s surface, while the upper condition is a zero-flux condition. In a dry atmosphere, the tendencies for the momentum and potential temperature equations in the dry case are expressed from the diffusive terms given by Eq. (27) as $$\begin{array}{cc}{F}_{\theta }=\frac{\text{d}{\upsilon }_{\theta }}{\text{d}t},& Q={\left(\frac{p}{p_{\text{ref}}}\right)}^k\frac{\text{d}}{\text{d}t}\left(C_p\text{Θ}\right).\end{array}$$(29)

To calculate the eddy diffusivities, we make use of the formulation given by Holtslag & Boville (1993). In that model, the eddy diffusivities of Eq. (27) are expressed as functions of a mixing length scale ℓ, the local shear $\left|\frac{\partial {\upsilon }_{\theta }}{\partial z}\right|$, and the gradient Richardson number Ri defined by Eq. (26). They read (e.g. Holtslag & Boville 1993, Eq. (3.2)) $$K_X={\ell }^2\left|\frac{\partial {\upsilon }_{\theta }}{\partial z}\right|{ℱ}_X\left(\text{Ri}\right),$$(30)

where F_X (Ri) describes the functional dependence of K_X on the gradient Richardson number. The form of F_X is determined by the turbulent regime, which can be either stable (Ri ≥ 0) or unstable (Ri < 0). The mixing length is expressed as (Blackadar 1962) $$\ell =\frac{{\ell }_0\mathcal{K}z}{\mathcal{K}z+{\ell }_0},$$(31)

the parameter K ≈ 0.4 being the von Karman constant (e.g. Garratt 1994), and ℓ₀ the asymptotic length scale (ℓ ≈ ℓ₀ for Kz ≫ ℓ₀), which varies as a function of z (see Appendix F.1).

The turbulent friction of mean flows against the soil in the surface layer leads to sensible momentum and heat exchanges that are described in the form of parametrised surface fluxes. Denoting by M and H the subscripts for the momentum and heat components, respectively, we write the upward momentum and heat surface fluxes as $$F_{\text{M}}=-C_{\text{M}}{\rho }_{\text{SL}}\left|{\upsilon }_{\sigma ;\text{SL}}\right|{\upsilon }_{\theta ;\text{SL}},$$(32) $$F_{\text{H}}=-C_{\text{H}}{\rho }_{\text{SL}}{C}_p\left|{\upsilon }_{\sigma ;\text{SL}}\right|\left({\text{Θ}}_{\text{SL}}-{\text{Θ}}_{\text{s}}\right),$$(33)

where the subscripts s and SL refer to values at the surface and at the top of the surface layer³, respectively. The surface-layer exchange coefficients C_M and C_H are defined as (Holtslag & Boville 1993) $$C_{\text{M}}\equiv C_{\text{N}}{f}_{\text{M}}\left({\text{Ri}}_0\right),$$(34) $$C_{\text{H}}\equiv C_{\text{N}}{f}_{\text{H}}\left({\text{Ri}}_0\right).$$(35)

Here, C_N designates the neutral exchange coefficient (e.g. Holtslag & Boville 1993), $$C_{\text{N}}\equiv {\left[\frac{\mathcal{K}}{\text{1n}\left(1+z_{\text{SL}}/z_{\text{r}}\right)}\right]}^2,$$(36)

where z_r denotes the roughness height, while f_M and f_H are two functions of the bulk Richardson number, $${\text{Ri}}_0\equiv \frac{g{z}_{\text{SL}}\left({\text{Θ}}_{\text{SL}}-{\text{Θ}}_{\text{s}}\right)}{{\text{Θ}}_{\text{SL}}{\left|{\upsilon }_{\sigma ;\text{SL}}\right|}^2},$$(37)

which controls the stability of the surface layer. We note that the bulk Richardson number Ri₀ corresponds to the local gradient Richardson number Ri (Eq. (26)) characterising the surface layer. The functions f_M and f_H, as well as the functions F_X introduced in Eq. (30), are detailed in Appendix F.1. The physical tendencies resulting from turbulent diffusion are evaluated every physical time step using an implicit scheme (see Appendix F.2).

As it accounts for the dependence of eddy diffusivities on the gradient Richardson number, the above turbulent diffusion scheme describes both the regime of strong convection occurring on the dayside and the regime of stable stratification associated with the nightside temperature inversion (see Figs. 1 and 2 in the following). It thus captures the evolution of turbulent diffusivities in the PBL between the dayside, where they are high, and the nightside, where they are low. However, we note that the used standard formulation of turbulent diffusion is not sophisticated enough to account properly for the heat and momentum exchanges due to turbulent flows in the case of strong stratification (Ri ≫ 1). In this regime, the vertical momentum mixing continues even at relatively high Ri due to the momentum transport of vertically propagating internal gravity waves, which may increase the surface-atmosphere sensible heat exchanges and warm up the nightside surface (e.g. Sukoriansky et al. 2005; Joshi et al. 2020). Considering this effect, the used turbulent scheme might tend to underestimate the atmospheric stability against collapse. Nevertheless, we adopt it as a convenient compromise between simplicity and realism, which can be refined in the future by using more advanced turbulent diffusion models (e.g. Sukoriansky et al. 2005).

3.4 Soil heat transfer

The soil heat transfer is solved by using a classical 1D soil heat conduction approach (e.g. Hourdin et al. 1993; Wang et al. 2016). The evolution of the ground temperature, T, due to vertical diffusion is governed by the heat equation $$C_{\text{gr}}\frac{\partial T}{\partial t}=-\frac{\partial F_{\text{c}}}{\partial u},$$(38)

where u = −z is the depth from surface (u ≥ 0), C_gr the heat capacity of the ground per unit volume, and F_c the conductive flux propagating downwards. The conductive flux is expressed as $$F_{\text{c}}=-{\lambda }_{\text{gr}}\frac{\partial T}{\partial u},$$(39)

the parameter λ_gr being the thermal conductivity of the material. Both C_gr and λ_gr are assumed to be constants in the model.

We introduce the normalised pseudo-depth $$\tilde{u}\equiv u\sqrt{\frac{C_{\text{gr}}}{{\lambda }_{\text{gr}}}},$$(40)

which has dimensions of s^1/2. Expressed in terms of ũ, the vertical diffusion is controlled by one parameter solely, the thermal inertia $$I_{\text{gr}}=\sqrt{{\lambda }_{\text{gr}}{C}_{\text{gr}}}.$$(41)

The downward conductive flux given by Eq. (39) then becomes $$F_{\text{c}}=-I_{\text{gr}}\frac{\partial T}{\partial \tilde{u}},$$(42)

and the vertical diffusion equation simply reads $$\frac{\partial T}{\partial t}=\frac{{\partial }^{2}T}{\partial {\tilde{u}}^2}.$$(43)

The above equation is solved by means of the finite difference method between the surface (ũ = 0) and an inner boundary (ũ = ũ_bot) using an implicit time scheme. At the surface, a continuity condition is applied: the incoming heat fluxes equalise the outcoming fluxes. At the inner boundary, a zero-flux condition is applied, which enforces the assumption that the planet interior is in thermal equilibrium. These two conditions are respectively formulated as $$\begin{array}{ll}-F_{\text{c}}+F_{\downarrow }-F_{\uparrow \text{S}}-F_{\text{H}}-{\epsilon }_{\text{s}}{\sigma }_{\text{SB}}{T}_{\text{s}}^4=0\hfill & \text{at\hspace{0.17em}}\tilde{u}=0,\hfill \\ \hfill F_{\text{c}}=0& \text{at\hspace{0.17em}}\tilde{u}={\tilde{u}}_{\text{bot}},\hfill \end{array}$$(44)

where F_↓ designates the downward radiative flux (i.e. the sum of the shortwave and longwave contributions), F_↑S the reflected shortwave flux, T_s the surface temperature, and ε_s the surface emissivity in the longwave. The scheme used to solve the soil heat transfer is detailed in Appendix G. Similarly as the radiative transfer and turbulent diffusion equations, the soil heat transfer equation is discretised and integrated as a boundary condition problem by means of the tri-diagonal matrix algorithm (see Appendix H).

Several simulations were run in order to assess the sensitivity of mean fields to the heat transfer scheme. The obtained results are discussed in Appendix G. They show that mean flows and temperatures do not depend much on the ground thermal inertia that parametrises the soil heat transfer scheme (Eq. (41)). Particularly, the nightside surface temperature is almost insensitive to the value chosen for I_gr. Nevertheless, we recall that horizontal diffusion is ignored in the scheme since we focus of dry rocky planets, while it might play an important role for an ocean planet due to heat advection by oceanic flows, as discussed by Wordsworth (2015).

Fig. 1 Two-day averaged snapshots of pressure (left), temperature (middle), and vertical wind speed (right) distributions for the 1 bar-atmosphere of an Earth-sized tidally locked planet (Earth-like case of Table 2) with a stellar irradiation of 1366 W m−2 after convergence (t = 400 days). From bottom to top: OD (1 × 1 grid), 1D (1 × 50 grid), 1.5D (2 × 50 grid), and 2D (32 × 50 grid) instances of the meta-model. The sub-stellar point corresponds to θ = 0° and the anti-stellar point to θ = 180°.

Fig. 2 Two-day averaged potential temperature field in the Earth-like case of Table 2 with ps = 1 bar and F⋆ = 1366 W m−2. Plotted values belong to the interval [248 K, 376 K]. Outside of this interval, the colour is set to that of the closest bound. The 2D instance of the meta-model (32 × 50 grid) is used. The latitudes 90° and −90° correspond to the sub-and anti-stellar points, respectively.

3.5 Physical setup

In the following, we perform simulations for the two cases defined in Table 2: (i) and Earth-like atmosphere, and (ii) a pure CO₂ atmosphere. For the Earth-like case, we use the values given by Deitrick et al. (2020) in the synchronous Earth case (see Deitrick et al. 2020, Table 2). These values correspond to a tide-locked Earth-sized planet with an atmosphere having the thermodynamical properties of the Earth’s atmosphere. Following Wordsworth (2015), we assume that τ_L = 1 at p = 1 bar. Similarly, we fix τ_S = 0.01 at p = 1 bar to enforce the assumption that the atmosphere is optically thin in the visible. The effective absorption coefficients introduced in Eq. (25) are set accordingly. Besides, we consider the case of pure absorption (no scattering). The scattering parameter is therefore set to β₀ = 1 both for the shortwave and longwave. Assuming that the planet’s surface is made of bare rocks, we use the typical values commonly assumed for Venus-like soils (e.g. Lebonnois et al. 2010) to set the planet’s surface properties. Following Frierson et al. (2006), the roughness height is set to z_r = 3.21 × 10⁻⁵ m so that C_N = 10⁻³ for z_SL = 10 m.

The pure CO₂ case is similar to the Earth-like case except for the specific gas constant, R_d = 188.9 J kg⁻¹ K⁻¹ (calculated from Meija et al. 2016), the heat capacity per unit mass, C_p = 909.3 J kg⁻¹ K⁻¹ (evaluated for Τ = 350 Κ from Yaws 1996, Appendix E), and the absorption coefficient in the longwave, κ_L = 2.5 × 10⁻⁴ m² kg⁻¹. The latter is an effective value of κ_L for which the 2D instance of the GCMM presented in Sect. 4 approximately reproduces the stability diagram obtained by Wordsworth (2015) from 3D GCM simulations with correlated-k radiative transfer (Wordsworth 2015, Fig. 12), as shown in Sect. 5. We note that this value is less than that obtained by Wordsworth (2015) by adjusting the analytic solution given by Eq. (5) to GCM simulations (κ_L = 3.2 × 10⁻⁴ m² kg⁻¹ by taking into account the factor 2 difference between the definition of κ_L given by Eq. (12) of the article and ours, given by Eq. (25)). As highlighted by Sect. 5, this discrepancy is due to the fact that the PBL acts to warm up the nightside surface in the general case, which consequently requires a weaker greenhouse effect to reach the same thermal state.

4 Climate Regime: From 0D to 2D Models

The GCMM described in the preceding section can be used at the same time for various grid configurations, each of them being a possible instance of the meta-model. This allows the results obtained from a bench of models of various complexities, though with the same intrinsic theoretical background and physical setup, to be compared. Four instances of the meta-model are examined in the present study (Table 1). There are introduced below in ascending order of complexity, while the resulting pressure, temperature, and vertical wind speed snapshots are plotted in Fig. 1.

We note that, for legibility, the arrows representing wind speeds are uniformly spaced both horizontally and vertically in the figure, rather than being centred at the points where winds speeds are really evaluated. For instance, in the 1.5D model, the arrows are determined from a linear interpolation and located at θ = 45°, 135°, while horizontal wind speeds are evaluated at θ = 90°. In the range 180° < θ < 360°, the distributions are plotted by applying a symmetric transformation to the distributions calculated in the range 0° < θ < 180^°. The additional inner ring in temperature snapshots (middle column of Fig. 1) corresponds to the soil temperature.

The 0D model (1 × 1 grid) is the simplest configuration. Here the atmosphere is vertically and horizontally isothermal, and it exchanges heat with dayside and nightside isothermal surface hemispheres through radiative transfer only. This configuration corresponds exactly to that of the two-layer grey radiative model described in Sect. 3 of Auclair-Desrotour & Heng (2020), which provides closed-form solutions. In the optically thin limit, it reproduces the analytic model introduced by Wordsworth (2015).

The next level of complexity is the 1D model (1 × 50 grid). In this configuration, the vertical structure of the atmosphere is allowed to adjust with radiative transfer although it is still horizontally isothermal. Similarly as in the 0D configuration, the planet’s surface is divided into the dayside and nightside hemispheres, and the surface-atmosphere heat exchanges are purely radiative. The 1D model thus provides a more realistic vertical temperature profile than that assumed in the 0D model. This profile can be interpreted as an approximation of the globally averaged atmospheric structure. This instance of the GCMM appears as a simplified version of more sophisticated 1D models (e.g. Robinson & Catling 2012).

Then we introduce the 1.5D model (2 × 50 grid). This level corresponds to the two-column approach, where the atmosphere is modelled by dayside and nightside hemispherical air columns, similarly as the planet’s surface (e.g. Yang & Abbot 2014; Koll & Abbot 2016). Therefore, the vertical structure ceases to be horizontally uniform and differences appear between the day-side and the nightside. However, these differences are mitigated by the fact that the dayside and nightside PBLs are both controlled by the horizontal wind speed at the planet’s terminator (θ = 90°) here, whereas horizontal wind speeds strongly differ between the dayside and nightside in reality. Besides, this configuration allows the coupling between the thermodynamics and the day-night overturning circulation to be taken into account, particularly the contribution of heat advection to the planet’s thermal state of equilibrium.

Finally, the 2D model (32 × 50 grid) is the more complex instance of the meta-model. In this configuration, the 2D model behaves similarly as a 2D GCM (e.g. Song et al. 2021) or a 3D GCM in the slow rotation regime (e.g. Leconte et al. 2013; Haqq-Misra et al. 2018; Pierrehumbert & Hammond 2019; Turbet et al. 2021). Both the atmospheric structure and mean flows are fully resolved, which describes the global heat engine circulation. The PBL strongly differs between the dayside, where it is unstable (Ri < 0) due to convection (see e.g. Koll & Abbot 2016), and the nightside, where it is stable (Ri ≥ 0). Also, the circulation is highly asymmetric with respect to the terminator: as shown by Fig. 1 (top panels), it exhibits strong upwelling flows in a small region around the sub-stellar point, and weak subsidence over a large area that spreads from θ ≈ 60° to θ = 180° (anti-stellar point).

We recover the day-night evolution of the atmospheric structure in the potential temperature distribution shown by Fig. 2, which was obtained using the 2D model. In this figure, the latitude is measured with respect to the terminator, located at 0°. The sub- and anti-stellar point thus correspond to the latitudes 90° and −90°, respectively. The nightside stably stratified region is characterised by a positive vertical potential temperature gradient. On the dayside, the convective mixing induced by the thermal forcing makes the temperature gradient converge towards the adiabat. As a consequence, convective regions are indicated by vertically uniform profiles of potential temperature. We note that the thickness of the PBL grows as the latitude increases, and that it is maximal at the sub-stellar point.

Owing to spatial bi-dimensionality, the solver is relatively fast even for the 2D instance, where one day of simulation with a two-minute dynamical timestep is equivalent to 0.61 seconds of CPU time on one CPU. This allows grid simulations to be run, which is the purpose of the next section.

5 Stability diagrams

In this section we examine both the role played by several physical features in the atmospheric stability, and the sensitivity of model predictions to the simplifications made in the different approaches. To do so, we proceed to a vertical inter-comparison between the four models introduced in Sect. 4 for the two Earth-sized synchronous planets defined in Table 2. For each instance of the meta-model, simulations were performed on a 15 × 13 grid in the space of stellar flux and initial surface pressure, with 0.2F_⊕ < F_⋆ < 3F_⊕ and 0.01 bar ≤ p_s ≤ 10 bar, F_⊕ = 1366 W m⁻² being the Earth’s incident stellar flux. Starting from isothermal and zero-velocity initial conditions, simulations were run for a period t_run = min {max[t_rad, t_min], t_max} ranging between t_min = 300 and t_max = 30000 Earth days, t_rad being an empirical estimate of the timescale necessary to reach radiative equilibrium that includes the dependence of the radiative cooling timescale on surface pressure and stellar flux (e.g. Showman & Guillot 2002), $$t_{\text{rad}}=900\text{\hspace{0.17em}days}\times \left(\frac{p_{\text{S}}}{1\text{\hspace{0.17em}bar}}\right){\left(\frac{F_{★}}{1366\text{\hspace{0.17em}W\hspace{0.17em}}{\text{m}}^{-2}}\right)}^{-3/4}.$$(45)

At the end of simulations, two-day averaged distributions of mean fields were computed, as well as the resulting minimum (or nightside) surface temperature (T_n).

As highlighted by Wang & Wordsworth (2020) in the case of sub-Neptunes, the convergence time of 3D numerical simulations can be extremely long (t_run ~ 250 000 Earth days, typically) for massive atmospheres (p_s ≳ 80 bar) owing to the long radiative timescale of deep atmospheric layers. In the present study, the maximum surface pressure (10 bar) is smaller than that usually assumed for sub-Neptunes, the circulation described by our 2D model is simpler than that described by 3D GCMs – where complex structures and super-rotating zonal jets can emerge –, and the major part of the incident stellar flux reaches the planet’s surface, which tends to facilitate vertical energy transport. Therefore the resulting convergence times are expected to be smaller than those reached in the case of sub-Neptunes by one order of magnitude approximately. However, we verified a posteriori that both the energy balance between the outgoing longwave radiation and the absorbed stellar radiation on the one hand, and the circulation on the other hand, had reached a steady state at the end of several test simulations. Besides, by running grid simulations with larger values of t_run, we noticed that increasing this parameter did not affect the considered mean fields.

Following earlier studies (Wordsworth 2015; Koll & Abbot 2016; Auclair-Desrotour & Heng 2020), we assume that the greenhouse effect is mainly due to the presence of CO₂ in the atmosphere. The condensation temperature of CO₂ is given, in K, by (Fanale et al. 1982; Wordsworth et al. 2010b; Wordsworth 2015) $$T_{\text{cond},{\text{CO}}_2\left(p\right)}=\{\begin{array}{ll}\frac{3167.8}{23.23-1\text{n}\left(0.01p\right)}\hfill & \text{if\hspace{0.17em}}p<p_{\text{tr}},\hfill \\ 684.2-92.31\text{n}\left(p\right)+4.321{\text{n}}^2\left(p\right)\hfill & ifp\ge p_{\text{tr}},\hfill \end{array}$$(46)

where the partial pressure of the gas p is given in Pa, and p_tr = 5.18 × 10⁵ Pa designates the triple point pressure. The stability diagrams are therefore obtained by comparing the minimum surface temperature calculated from simulations with the condensation temperature of CO₂ at the planet’s surface, $T_{{\text{cond,CO}}_2}\left(\chi p_s\right)$, where χ designates the volume mixing ratio of CO₂. In the Earth-like case, the volume mixing ratio of CO₂ is set to the value of Earth at the beginning of the 21st century, namely χ = 370 ppm (e.g. Etheridge et al. 1996), while χ = 1.0 in the pure CO₂ case. The atmosphere is considered to be stable if $T_{\text{n}}\left(F_{\ast },p_{\text{s}}\right)>T_{\text{cond,}}{}_{{\text{CO}}_2}\left(\chi p_{\text{s}}\right)$ and unstable else. We remark that the collapse itself, when it occurs, is not described by the model since the changes of phases of CO₂ are not taken into account.

Figure 3 shows the simulation results. The minimum surface temperature is plotted in both cases as a function of the normalised stellar flux F_⋆/F_⊕ and the surface pressure p_s in logarithmic scale. The stability diagrams are plotted too, with large red dots indicating stability and small blue dots collapse. The collapse pressures associated with the lower and upper bounds of the nightside temperature given by Eqs. (5) and (10), are obtained by solving for ps the equations $$T_{\text{n;low}}\left(F_{★},p_{\text{S}}\right)=T_{{\text{cond,CO}}_2}\left(\chi p_{\text{S}}\right),$$(47) $$T_{\text{n;up}}\left(F_{★},p_{\text{S}}\right)=T_{{\text{cond,CO}}_2}\left(\chi p_{\text{S}}\right),$$(48)

and are denoted by p_C;low (orange dashed line) and p_C;up (pink dotted line), respectively.

With the 0D model, we recover the behaviour of the nightside temperature predicted by the closed-form solutions of Auclair-Desrotour & Heng (2020) in the purely radiative regime. This is due to the fact that the two models are actually the same in this configuration, the temperatures being computed numerically here instead of analytically. Compared with the other models, the 0D model tends to underestimate the nightside temperature in the high surface pressure regime. As the surface pressure increases, T_n reaches a plateau that corresponds to the planet’s equilibrium temperature, T_eq (F_⋆), given by Eq. (6), and no longer evolves with p_s. This unrealistic behaviour is a consequence of the isothermal approximation, which does not account for the strong vertical temperature gradient characterising thick atmospheres, especially their convective regions. In the low stellar flux regime, the 0D model captures the stability decrease observed for pure CO₂ atmospheres in the 3D GCM simulations performed by Wordsworth (2015). However, this feature is due to the isothermal temperature profile too since it vanishes from the moment that the atmospheric structure is allowed to adjust with radiative transfer, in models of higher dimensions. This effect is a caveat of the limitations of idealised models in explaining predictions of much more sophisticated 3D GCMs.

The 1D model exhibits the same behaviour as the 0D model for surface pressures of less than 1 bar, which corresponds to the regime where the vertically isothermal approximation holds. Beyond p_s ≈ 1 bar, the nightside surface temperature increases as a function of both the stellar flux and surface pressure, which makes the collapse pressure associated with T_n;low capture the threshold of the stability region for the whole stellar flux interval. The two-column configuration (1.5D model) relaxes the horizontally isothermal atmosphere approximation. As a consequence, the nightside temperature becomes dependent on the efficiency of the interhemispheric heat redistribution, and can thereby be less than the lower bound obtained in the horizontally isothermal atmosphere approximation, T_n;low. The coarse spatial resolution of the 1.5D model for the horizontal direction does not account for the strong convection generated in the sub-stellar region. The wind speed is therefore underestimated, and so are the strength of the overturning circulation and the heat advected from the dayside to the nightside. This leads to the observed stability decrease: the collapse pressure is approximately increased by ~25% with respect to p_C;low in the Earth-like case, and by ~80% in the pure CO₂ case.

Conversely, the 2D model predicts a wider stability region. The collapse pressure is lowered by ~ 10–40% with respect to p_C;low for F_⋆ ≳ 0.7F_⊕, which is significant, albeit less than the 75% maximum decrease predicted by the analytic theory (p_C;up/p_C;low = 1/4; see Auclair-Desrotour & Heng 2020, Eq. (86)) in the case of intense sensible heating (L_sen → + ∞; see Eq. (9)). This stability increase results from the effect of the PBL. The vertical turbulent diffusion generated by the friction of mean flows against the planet’s surface in the planetary layer acts both (i) to increase the thermal forcing of the atmosphere by intensifying sensible exchanges, and (ii) to enhance the day-night heat advection by strengthening the overturning circulation. As a consequence, the nightside surface is warmer by ~4–14 K in the vicinity of the threshold between the stability and collapse regions.

The nightside temperature increase induced by the PBL was quantified by running simulations in the 2D configuration without turbulent diffusion both for Earth-like and pure CO₂ atmospheres. In these simulations, the surface-atmosphere heat exchanges are induced by radiative transfer only, and there is no friction of mean flows against the surface. Figure 4 shows the resulting stability diagrams, as well as the corresponding night-side surface temperature difference between the cases with and without turbulent diffusion, denoted by the superscripts TD (turbulent diffusion) and NTD (no turbulent diffusion), respectively (top panels). In addition with the nightside surface temperature, we consider the day-night advection timescale t_adv, which is defined here as the mean period necessary for a fluid parcel to accomplish one full cycle of the day-night overturning circulation (see Appendix I), $$t_{\text{adv}}\equiv \frac{4R_{\text{p}}}{{\left[{\displaystyle {\int }_0^1\left|{\upsilon }_{\theta }\right|\text{d}\sigma }\right]}_{90°}},$$(49)

where the subscript 90° indicates that the integral of the mass flow rate is performed over the terminator annulus (θ = 90°). The day-night advection timescale also corresponds to the mean renewal time of the air contained in one atmospheric hemisphere (dayside or nightside). This timescale can be compared to the advection timescales introduced in earlier studies, such as the analytic expression of the lower bound obtained by Koll & Abbot (2016) from the heat engine theory (Koll & Abbot 2016, Eq. (12)), $$t_{\text{adv}}^{\text{KA}}\equiv \frac{R_{\text{p}}}{{\upsilon }_{\text{sen}}},$$(50)

where the typical speed υ_sen, given by Eq. (8), is a function of the stellar flux, the surface pressure, the surface albedo, the day-side surface temperature, the optical depth in the longwave at surface, the specific gas constant, and the bulk drag coefficient of the surface layer. The timescale $t_{\text{adv}}^{\text{KA}}$ is estimated by setting the bulk drag coefficient to the typical value C_D = 10⁻³ for convenience and by taking the maximum surface temperature for the dayside temperature. Thus, in addition with the nightside temperature difference mentioned above, we plot in Fig. 4 the day-night advection timescale in the case with turbulent diffusion (bottom panels), the ratio of this timescale over $t_{\text{adv}}^{\text{KA}}$, and the ratio between advection timescales in the cases with and without turbulent diffusion (middle panels).

We first consider the stability diagrams obtained in the absence of turbulent diffusion (Fig. 4, top and middle panels). Similarly as in the ID configuration, the threshold of the stability region coincides with the lower bound of T_n associated with the purely radiative regime in the radiative box model, namely p_C;low. This indicates that the bulk atmosphere is horizontally isothermal for p_s ≳ 0.1 bar, which corresponds to an efficient interhemispheric heat redistribution. As highlighted by temperature differences (Fig. 4, top panels), turbulent diffusion tends to warm up the nightside surface in the general case with, for instance a temperature increase of 6–28 K in the Earth-like case. However, this temperature increase does not vary monotonically with the stellar flux and surface pressure, but instead it exhibits a bi-modal behaviour with maxima and minima depending on surface pressure.

Particularly, turbulent diffusion in the PBL somehow counter-intuitively acts to decrease the nightside surface temperature instead of increasing it in a region centred on p_s ~ 1 bar for pure CO₂ atmospheres, with a minimum of −8 K for p_s = 1 bar and F_⋆ = 1.2F_⊕. A similar – although slightly smaller – negative difference (–5 K) was obtained by running grey gas simulations with the fully global 3D GCM THOR (Mendonça et al. 2016; Deitrick et al. 2020) in this configuration using the values given by Table 2 and assuming a zero-spin angular velocity, which tends to corroborate the prediction of the 2D model. This effect of the PBL can be analysed through the interplay between the day-night advection timescale given by Eq. (49) and the dayside radiative timescale, t_rad, which is the typical timescale needed for a warm fluid parcel located in the upper atmosphere to cool down radiatively. Both parameters are altered by the turbulent diffusion taking place within the PBL.

Notwithstanding the high pressure – and optically thick – regime, where strong convection develops, the effect of turbulent diffusion reaches a maximum around p_s ~ 0.1 bar for Earth-like atmospheres and around p_s ~ 0.03 bar for pure CO₂ atmospheres (Fig. 4, top panels). This maximum is consistent with the fact that the additional thermal forcing due to turbulent diffusion is all the more significant as the radiative absorption is weak, which tends to maximise the impact of turbulent diffusion for small optical thicknesses. However, the radiative timescale of the atmosphere scales as $t_{\text{rad}}\propto p_s/T_{\text{a}}^3$ (e.g. Showman & Guillot 2002, Eq. (10)) in the optically thin regime, meaning that the heat surplus provided by sensible exchanges is radiated towards space over timescales that become extremely short as the surface pressure tends to zero. Thus, in spite of the decay of t_adv induced by turbulent diffusion, a fluid parcel is radiatively cooled before being advected to the nightside by mean flows in the optically thin limit (t_rad ≪ t_adv), which mitigates the impact of the PBL in this regime and makes the temperature difference decay as p_s → 0.

Similarly, as the surface pressure increases, the atmosphere switches from the optically thin regime to the optically thick regime, with a transition occurring at lower pressures in the pure CO₂ case than in the Earth-like case due to the difference between optical depths in the infrared. In this transition regime, the PBL strongly affects the advection timescale, which is increased up to three times (Fig. 4, middle panels) and becomes thereby greater than the radiative timescale. As a consequence, less heat is advected towards the nightside than in the absence of PBL – although the PBL generates a heat surplus on the dayside – and the nightside temperature difference falls to the observed minimum valley, where it reaches negative values in the pure CO₂ case.

We note that this behaviour could be significantly altered by the presence of gas, dust, and aerosols inducing the so-called anti-greenhouse effect by increasing the shortwave scattering and absorption, as observed on Titan (McKay et al. 1991). The anti-greenhouse effect refers to the cooling of the planet surface resulting from the fact that a substantial part of the incident stellar flux is absorbed and re-radiated towards space in the infrared by the upper layers of the atmosphere, and that only a fraction of it reaches the surface (e.g. Pierrehumbert 2010). This effect is likely to play a major role on rocky planets with thick atmospheres similar to Venus, where only ~ 0.1–1% of the incident solar flux reaches the surface (Lacis 1975).

We now consider the evolution of the advection timescale itself, which is plotted in the case including turbulent diffusion in the PBL (Fig. 4, bottom panels). Similarly as the nightside surface temperature (Fig. 3), t_adv exhibits a relatively smooth dependence on the incident stellar flux and surface pressure, with values spanning over two orders of magnitude from ~8 days in the low pressure-high stellar flux regime to ~500 days in the high pressure-low stellar flux regime. These values are larger than those given by the analytic scaling law of Koll & Abbot (2016) by one order of magnitude owing to the difference in the used definitions: since it is computed from the horizontal wind speed in the sub-stellar region, the analytic advection timescale is necessarily smaller than the advection timescale defined by Eq. (49), which is computed from mass flows going through the terminator annulus. Notwithstanding this scaling factor, matches relatively well for both Earth-like and pure CO₂ atmospheres in the optically thin regime, where the ratio varies by a factor of two. The two quantities diverge from each other as the stellar flux decreases and the surface pressure increases, the model of the present study predicting larger advection timescales in this regime. However, we remark that this discrepancy is less significant for pure CO₂ atmospheres than for Earth-like atmospheres, which suggests that the thermodynamic and absorption properties of the gas affect the large-scale overturning circulation in a non-negligible way when the atmosphere is optically thick.

The strength of the overturning circulation can be characterised by examining the behaviour of the Eulerian mean stream function, which is defined in TLCs as (e.g. Pauluis et al. 2008) $$\text{Ψ}\equiv \frac{2\pi R_{\text{p}}}{g}{\displaystyle {\int }_p^{p_{\text{s}}}{\upsilon }_{\theta }\sin \theta \text{d}p.}$$(51)

The Eulerian mean stream function measures here the strength of longitude-averaged cells in the TLC. It accounts for the vor-ticity of the flow in a plane containing the planet-star axis. In the slow rotation regime, the large-scale cell of the predominant overturning circulation corresponds to a large region centred on a maximum of Ψ in absolute value, Ψ taking negative values owing to the flow direction (see e.g. Fig. 6 in the next section). The maximum value of Ψ over the atmospheric domain, defined as $${\text{Ψ}}_{\max }\equiv \max \left\{-\text{Ψ}\right\},$$(52)

can be scaled analytically, as shown by Innes & Pierrehumbert (2022) who establishedforsub-Neptunes the scaling law ${\Psi }_{\max }\propto F_{*}^{3/4}$. In the case of rocky planets, we need to account for the presence of the planet’s surface, which sizes the thickness of the atmosphere. The mass flow thus depends on the atmospheric mass in addition with the stellar flux. As the atmospheric mass is directly proportional to the surface pressure, we introduce for rocky planets a scaling law of the form $${\text{Ψ}}_{\max }^{\text{SL}}=\frac{R_{\text{p}}^2{F}_{\oplus }}{C_p{T}_{\oplus }}{\left(\frac{F_{★}}{F_{\oplus }}\right)}^{\alpha }{\left(\frac{p_{\text{s}}}{p_{\text{ref}}}\right)}^{\beta },$$(53)

where T_⊕ is the equilibrium temperature of Earth given by Eq. (6), p_ref = 1 bar a reference pressure, and α and β two real exponents to be defined. We remark that (α,·β) = (3/4, 0) corresponds to the scaling law obtained by Innes & Pierrehumbert (2022) for sub-Neptunes. However, for the Earth-sized rocky planets of the present study, we adopt the empirical values (α,β) = (1/2, 1), which seem to be more appropriate as discussed further.

Figure 5 shows the evolution Ψ_max and of the ratio ${\Psi }_{\max }/{\Psi }_{\max }^{\text{SL}}$ as functions of instellation and surface pressure for the Earth-like and pure CO₂ atmospheres defined in Table 2. In both cases the maximum of the Eulerian mean stream function varies over more than three orders of magnitude (Ψ_max ~ 10¹⁰−10¹³ kg s⁻¹), by increasing with F_⋆ and p_s monotonically. The scaling law given by Eq. (53) (with α = 1/2 and β = 1) approximately captures this behaviour although it can only be considered as a rough estimate of Ψ_max. We observe that the dependence of Ψ_max on the instellation and surface pressure evolves between the optically thin and optically thick regimes with transition zones located around p_s ~ 0.3 bar for Earth-like atmospheres and 0.1 bar for pure CO₂ atmospheres. The scaling law ${\Psi }_{\max }\propto F_{*}^{\text{1/2}}$ is appropriate to describe the dependence of Ψ_max on F_⋆ in the regime of optically thick atmospheres in the infrared, but not in the regime of optically thin atmospheres where increases faster with F_⋆. Conversely, the slow evolution of ${\Psi }_{\max }/{\Psi }_{\max }^{\text{SL}}$ with p_s at low surface pressures suggests that the scaling law Ψ_max ∝ p_s captures well the dependence of Ψ_max on the surface pressure in the regime of optically thin atmospheres, while this dependence is weaker in the regime of thick atmospheres.

Fig. 3 Stability diagrams of Earth-sized tidally locked planets with Earth-like (left panels) and pure CO2 (right panels) atmospheres, and the associated minimum surface temperature. Quantities are plotted as functions of the stellar flux normalised by the Earth’s stellar flux, F⋆/F⊕ (horizontal axis), and the surface pressure, ps, in logarithmic scale (vertical axis). From bottom to top: 0D (1 × 1 grid), 1D (1 × 50 grid), 1.5D (2 × 50 grid), and 2D (32 × 50 grid) instances of the meta-model. Large red dots indicate simulations where the atmosphere remained stable, while small blue dots indicate atmospheric collapse. The dashed orange (or dotted pink) line indicates the collapse pressure pC;low (or pC;up) corresponding to the lower (or upper) bound of the nightside surface temperature predicted by Wordsworth’s analytic model and given by Eq. (11).

Fig. 4 Stability diagrams of Earth-sized tidally locked planets hosting Earth-like (left panels) and pure CO2 (right panels) atmospheres, with turbulent diffusion (TD, four bottom panels) and without turbulent diffusion (NTD, four top panels), superimposed on the maps of various quantities. From bottom to top: day-night advection timescale with turbulent diffusion, day-night advection timescale with turbulent diffusion over the lower bound estimated analytically by Koll & Abbot (2016) using the heat engine theory, ratio of day-night advection timescales with and without turbulent diffusion, and nightside surface temperature difference between the cases with and without turbulent diffusion. The curves and symbols are the same as in Fig. 3.

Fig. 5 Maximum of Eulerian mean stream function (top) and its ratio over the scaling law (bottom) for Earth-sized planets hosting Earth-like (left panels) and pure CO2 (right panels) atmospheres. The two estimates of the Eulerian mean stream function are given by Eqs. (52) and (53), respectively. The curves and symbols are the same as in Fig. 3.

6 From slow to fast rotation

Since the GCMM is designed to study the asymptotic regime of slowly rotating planets, where mean flows are predominantly driven by the day-night temperature gradient, we have ignored the effects of rotation on the general circulation until now. As discussed in Sect. 2, the impact of these effects is mainly controlled by the non-dimensional equatorial Rossby deformation radius ${\tilde{L}}_{\text{Ro}}$ introduced in Eq. (1). The zero-spin rate limit treated by the model corresponds to ${\tilde{L}}_{\text{Ro}}=+\infty$. In reality Coriolis acceleration tends to deviate the divergent winds blowing between the sub-stellar and the anti-stellar points. This alters the large-scale circulation regime, which switches from the bi-dimensional Hadley cell described by the 2D model to the 3D structure characterising fast rotators, where super-rotation develops, as ${\tilde{L}}_{\text{Ro}}$ decays. In this section we investigate the limitations of the bidimensional approach by benchmarking the results obtained with the 2D instance of the GCMM against the simulations performed with the THOR GCM (Mendonça et al. 2016; Deitrick et al. 2020), which solves the 3D non-hydrostatic Euler equations on an icosahedral grid. Particularly, THOR fully accounts for the horizontal vortical component of mean flows that is ignored in the 2D model.

Four simulations were performed with THOR for the pure CO₂ atmospheres defined by Table 2. They each correspond to a given spin period Ρ taken in a range spanning from fast to slow rotators all things being equal (P = 10⁻³, 10⁻², 10⁻¹, 10⁰ P₀ with P₀ = 365 days). For comparison, one additional simulation was run with the 3D GCM in the slow rotator case (P = P₀) by assuming no turbulent diffusion or sensible surface-atmosphere heat exchanges. In these simulations, the values of the surface pressure (p_s = 0.18 bar) and stellar flux (F_⋆ = F_⊕) were chosen so that the studied cases are in the vicinity of the threshold between the stability and collapse regions predicted by the 2D instance of the meta-model (see Fig. 3, top right panel). The horizontal resolution of the icosahedral grid used in the model was set to ~4 degrees, and the atmosphere was divided into 40 vertical intervals logarithmically refined in the vicinity of the surface with a top altitude of 37000 m and a 2 m-thick lowest layer. For the sake of consistency, the simulations were run using the double-grey approximation for radiative transfer and a physical setup as close as possible as that implemented in the GCMM. After the circulation and radiative transfer have reached a state of equilibrium, the physical quantities were averaged over the longitude of the TLCs to transform the 3D fields into 2D fields similar to those displayed in Fig. 1. The bulk dimensionless equatorial Rossby deformation length (${\tilde{L}}_{\text{Ro}}$) was computed from the resulting mass-averaged temperature using the expression given by Eq. (4). In parallel of the simulations performed with THOR, a simulation corresponding to the zero-spin rate limit was run with the 2D instance of the meta-model.

Notwithstanding the vertical coordinates used in dynamical cores – altitude-based in THOR, and mass-based in the GCMM – the main differences between the two models essentially lie in the description of the surface thermal response, PBL, and numerical energy dissipation. In THOR, the surface temperature evolution is integrated using a 0D thermodynamic equation parametrised by an effective surface heat capacity ($C_{\text{s}}={10}^7{\text{JK}}^{-1}{\text{m}}^2$), while a ID soil heat transfer scheme is used in the GCMM (see Appendix G). In the treatment of the PBL, the two models closely follow the method proposed by Holtslag & Boville (1993), except for some parameters, which are set to different values. For instance, the asymptotic length scale characterising the evolution of the mixing length with altitude (Eq. (31)) for the heat equation is set to three times the function used for the momentum equation (Eq. (F.6)) in THOR, while the same function is used in both cases in the GCMM. As regards numerical energy dissipation, all THOR simulations utilised fourth-order horizontal hyper-diffusion and 3D divergence damping with a non-dimensional coefficient D_hyp = D_div = 0.002 (Mendonça et al. 2016; Deitrick et al. 2020, for descriptions of the hyper-diffusion scheme). We additionally used sixth-order vertical hyper-diffusion with a non-dimensional coefficient of D_ver = 5 × 10⁻⁴. Finally, a sponge layer was applied to the upper 25% of the model domain to damp spuriously reflected waves off the top boundary, with a damping timescale of 5000 s (see Mendonça et al. 2018; Deitrick et al. 2020, for the sponge layer description).

Figure 6 shows the simulation results. The dayside and nightside surface temperatures, as well as the day night surface temperature difference, are plotted as a function of the bulk equatorial Rossby deformation radius (top panels). Solid lines designate the temperatures obtained in TH0R 3D simulations, and the horizontal dashed lines the temperatures obtained in the 2D simulation performed with the meta-model. The dotted lines indicate the corresponding averaged dayside and nightside surface temperatures predicted by Wordsworth’s purely radiative greenhouse model (Wordsworth 2015). The temperatures obtained from the 3D simulation performed in the absence of turbulent diffusion are designated by the star symbol ⋆ In addition, the Eulerian mean stream function defined in Eq. (51) is plotted for the 2D model and the two extrema of the 3D model (P = 365 days and P = 0.365 days) as functions of pressure (or sigma coordinate) and the latitude of the TLCs, where the north and south poles correspond to the sub-stellar and anti-stellar points, respectively (bottom panels).

We first consider the evolution of surface temperatures with ${\tilde{L}}_{\text{Ro}}$ (Fig. 6, top panels). Whereas the dayside surface temperature predicted by 3D GCM simulations is hardly affected by the planet’s rotation, the nightside surface temperature increases monotonically with the equatorial Rossby deformation radius until it converges towards the asymptotic limit of slow rotators described by the 2D model. The effect of rotation is particularly significant for the minimum surface temperature, which varies by ~35 K between the two extremal cases. This difference results from the decay of day night advection provoked by the breaking of the day night overturning circulation in the fast rotator regime (${\tilde{L}}_{\text{Ro}}\lesssim 1$). The super-rotating equatorial jets induced by Coriolis effects do not compensate the overturning circulation in terms of mass flux crossing the terminator annulus, which makes the nightside surface temperature decrease. However, we remark that the nightside temperature stays close from the asymptotic limit of slow rotators from the moment that ${\tilde{L}}_{\text{Ro}}>2$, the difference of the cases P = 36.5 days and P = 365 days to this limit being less than 3 K. This suggests that the 2D model is relevant to describe the climate and large-scale circulation regime of the planet for equatorial Rossby deformation radii exceeding the critical value ${\tilde{L}}_{\text{Ro}}\approx 2$. Below this value, the circulation and heat redistribution are strongly affected by the planet’s rotation.

Additionally, we remark that the 2D model and the THOR GCM are in agreement with each other for the nightside surface temperature, which is greater by ~2 K only in 3D simulations. The discrepancy between the two models is more significant for the dayside average temperature, where the effects of the PBL in surface-atmosphere heat exchanges are predominant. As observed in the 2D model, the frictional interaction of mean flows with the planet’s surface taking place in the PBL tends to increase the atmospheric stability against collapse by warming up the nightside surface hemisphere. To understand whether the differences in PBL parameters between the 2D and 3D models could be the source of the discrepant dayside temperatures, one additional THOR simulation was run at the slowest rotation rate. In this simulation, the asymptotic scale lengths (for heat and momentum mixing in the boundary layer) were set equal to those used in the 2D model. The resulting values are shown as the square points in top panels of Fig. 6. Interestingly, these results are nearly indistinguishable from the simulation with different scale lengths used for the heat and momentum equations, indicating that the scenario is insensitive to the asymptotic scale lengths used in Eq. (31).

We now consider the snapshots of the Eulerian mean stream function (Fig. 6, bottom panels). We remark that the vertical coordinates used for the plots in the 2D and 3D cases differ from each other, which leads to slight distortions of the distributions. The sigma coordinate is such that σ = 1 everywhere at planet’s surface, while the surface pressure is less than the globally averaged surface pressure on the dayside, and greater on the nightside due to the day night thermal forcing gradient. As a consequence, isobars in sigma coordinates are slightly shifted downwards on the dayside and upwards on the nightside. The 3D simulation performed with THOR GCM in the slow rotation regime is in good agreement with the 2D simulation performed with the GCMM in the zero-spin rate limit (Fig. 6, bottom left and middle panels). The two snapshots both exhibit large day-night cells of comparable strengths, albeit slightly weaker in the 2D model. The snapshot corresponding to the fast rotation regime (bottom right panel) exhibits more complex features due to Coriolis effects although a weak day night cell still remains visible.

Fig. 6 Dayside and nightside surface temperatures vs. normalised equatorial Rossby deformation radius ${\tilde{L}}_{\text{Ro}}$ in the grey 3D GCM simulations performed with THOR. Top left: maximum (solid red line) and hemisphere-averaged (solid purple line) dayside surface temperatures. Top middle: minimum (solid blue line) and hemisphere-averaged (solid green line) nightside surface temperature. Top right: day night averaged (solid grey line) and extremal (solid black line) surface temperature differences. Dashed lines, dotted lines, star symbols, and square points indicate (i) the asymptotic surface temperatures computed using the 2D model, (ii) the hemisphere-averaged surface temperatures predicted by Wordsworth’s greenhouse model (W15; Wordsworth 2015), (iii) the surface temperatures obtained in the 3D simulation without turbulent diffusion in the slow rotator case, and (iv) those obtained in the simulation performed with the same asymptotic scale lengths in the PBL as in the GCMM, respectively. Bottom, from left to right: averaged snapshots of the Eulerian mean stream function (e.g. Pauluis et al. 2008) obtained from the 2D simulation (zero-spin rate limit), the 3D simulation for Ρ = 365 days (slow rotator), and the 3D simulation for Ρ = 0.365 days (fast rotator).

7 Concluding remarks

In this work we have developed a GCMM to bridge the gap between the analytic solutions provided by simplified greenhouse models for synchronous planets and the numerical simulations obtained from 3D GCMs in the asymptotic regime of slow rotators. This model hierarchy is based on a systematic bottom up approach in the spirit of Held (2005), wherein the number of degrees of freedom determines the key sources of complexity that are added or subtracted. The solver of the GCMM integrates the HPEs using the finite-volume method for arbitrary numbers of horizontal and vertical intervals, each configuration being an instance of the meta-model. Consistent with a previous analytical study (Auclair-Desrotour & Heng 2020), the physics implemented in the meta-model includes double-grey radiative transfer, turbulent diffusion in the PBL, and soil heat conduction. Particularly, the meta-model was designed so that the solutions obtained with the 0D instance exactly correspond to the analytic solutions of the purely radiative box model detailed in Sect. 3.3 of Auclair-Desrotour & Heng (2020).

As a first step, we proceeded to a vertical model inter-comparison by running grid simulations for four instances of the meta-model (0D, 1D, 1.5D, and 2D) in the cases of dry Earthlike and pure CO₂ atmospheres. In each case, we computed from simulations the nightside surface temperature of the planet as a function of the stellar flux and surface pressure as well as the resulting stability diagrams of the atmosphere against collapse. These diagrams were compared to the scaling laws predicted by the analytic theory. With the 0D and 1D instances of the meta-model, we recovered the stability diagrams predicted by simplified radiative models in the optically thin regime, which shows that the globally isothermal approximation used in these models is relevant in this regime. The 1.5D instance tends to underestimate the atmospheric stability, by predicting a collapse pressure 25% to 80% larger than that given by radiative box models. Conversely, the collapse pressure computed using the 2D instance of the meta-model is 10% to 40% smaller than the analytic estimate owing to the warming effect of the PBL, which is less than the theoretical 75% maximum decrease predicted by radiative models, albeit still significant.

As a second step, we investigated the role played by the PBL in the thermal state of equilibrium and atmospheric circulation of the planet by examining with the 2D instance of the GCMM how the turbulent diffusion taking place in the PBL alters the nightside temperature, the day night advection timescale, and the collapse pressure. We compared the advection timescale obtained from simulations with the analytic scaling law proposed by Koll & Abbot (2016). We observed that the turbulent diffusion taking place in the PBL increases the nightside surface temperature by 4–14 K around the threshold of the stability region. However, we found that the PBL can also contribute to cooling the nightside surface of the planet by acting on the day night advection timescale in the transition zone between optically thin and optically thick atmospheres. This result was corroborated a posteriori by 3D GCM simulations. The effect of the PBL on the large-scale circulation is complex and depends on the interplay between the advection and radiative timescales. The day night advection timescale estimated with the 2D model varies over two orders of magnitude in the studied domain of stellar fluxes and surface pressures, with values ranging between 8 and 500 days. In the optically thin regime its evolution matches the scaling law derived by Koll & Abbot (2016) relatively well, but it diverges from it at high pressures and low stellar fluxes. We noticed that this behaviour also depends on the thermodynamic and absorption properties of the atmosphere in the optically thick regime. In addition we empirically obtained, for the circulation of slowly rotating rocky planets, a scaling law analogous to that established analytically by Innes & Pierrehumbert (2022) for sub Neptunes.

As a third and final step, we characterised the limitations of the slow rotator approximation that forms the foundations of the meta-model by performing simulations with the THOR 3D GCM, which fully accounts for the effects of rotation on mean flows. We computed from these simulations the evolution of the planet’s dayside and nightside surface temperatures as functions of the dimensionless equatorial Rossby deformation radius, which controls the large-scale circulation regime of the planet. The results obtained with the 3D GCM were compared with the outcomes of the 2D instance of the meta-model. The 3D GCM simulations highlight the transition between the slow and fast rotation regimes. We found that the 2D model properly accounts for the climate and large scale atmospheric circulation from the moment that the normalised equatorial Rossby deformation radius is greater than the critical value ${\tilde{L}}_{\text{Ro}}\approx 2$, which corresponds to a broad region of the parameter space. In the slow rotation regime, the circulation and surface temperature predicted by the THOR 3D GCM and the 2D instance of the meta-model are similar.

This study is a first attempt to fill the continuum between analytic greenhouse models and 3D GCMs in a self-consistent way. The obtained results show that the meta-modelling approach is efficient at disentangling the mechanisms that determine the climatic state of the planet, which are narrowly coupled together in GCMs. This approach allows models of various complexities, albeit with the very same physical setup, to be run in parallel series so that each simulation can be interpreted using the others. As a consequence, the final outcome of the meta-model conveys information not only on the climate itself but also on the separated contributions of key physical ingredients such as radiative transfer, atmospheric structure, dynamics, and turbulent friction in the PBL. In addition to these diagnostic aspects, the meta-modelling approach appears as a robust method for refining the analytic theory of planetary climates given that it allows the relevance of solutions obtained at low spatial dimensionality to be assessed with consistently generated GCM numerical solutions. Finally, we note that the meta-modelling approach can be extended to moist atmospheres and rapidly rotating planets, but we leave that to the content of a future study.

Acknowledgements

The authors thank the anonymous referee for helpful comments that improved the manuscript, and Vera Matarese for inspiring discussions about the representational capacity of numerical simulations. They acknowledge financial support from the European Research Council via the Consolidator grant EXOKLEIN (grant number 771620). Kevin Heng 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. The authors gratefully acknowledge the open source softwares that made this work possible: numpy (Harris et al. 2020), matplotlib (Hunter 2007), scipy (Virtanen et al. 2020). THOR calculations were performed on UBELIX (http://www.id.unibe.ch/hpc), the HPC cluster at the University of Bern.

Appendix A Non-dimensional primitive equations

The model solves the HPEs given by Eqs. (15-18) in their non-dimensional form. Although it is not used in the present study, the moisture conservation equation is included in the solver and we give it here for the sake of generality. In its conservative form, this equation reads (e.g. Yao & Stone 1987) $$\frac{\partial }{\partial t}\left(\mathfrak{p}q\right)+{\nabla }_{\sigma }\cdot \left(\mathfrak{p}q{\upsilon }_{\sigma }\right)+\frac{\partial }{\partial \sigma }\left(\mathfrak{p}q\dot{\sigma }\right)=\mathfrak{p}\dot{q},$$(A.1)

where q designates the specific humidity of any tracer, and $\dot{q}$ the corresponding net evaporation or condensation rate per unit mass. Since the elemental unit of time of the discretised equations is the dynamical time step Δt used in the time-differencing scheme, it is convenient to normalise the time by Δt so that the time lapse between two dynamical time steps is always unity.

We then introduce the reference pressure p₀, temperature T₀, and specific humidity q_o, from which we can define the reference velocity υ₀, energy per unit mass e₀, density p₀, force per unit mass F₀, heat power per unit mass Q₀, net evaporation per unit mass ${\dot{q}}_0$, Exner function E₀, and potential temperature Θ₀, $$\begin{array}{lll}{\upsilon }_0\equiv \sqrt{C_p{T}_0},\hfill & e_0\equiv C_p{T}_0,\hfill & {\rho }_0\equiv \frac{p_0}{C_p{T}_0},\hfill \\ F_0\equiv \frac{\sqrt{C_p{T}_0}}{\text{Δt}},\hfill & Q_0\equiv \frac{C_p{T}_0}{\text{Δt}},\hfill & {\dot{q}}_0\equiv \frac{q_0}{\text{Δ}t},\hfill \\ E_0\equiv C_p{\left(\frac{p_0}{p_{\text{ref}}}\right)}^K,\hfill & {\Theta }_0\equiv \frac{C_p{T}_0}{E_0}.\hfill & \hfill \end{array}$$(A.2)

Besides, we make the horizontal coordinate vary within the range [0,1] like the vertical coordinate by re-normalising the colatitude (θ). The normalised variables and colatitude are therefore defined as $$\begin{array}{llll}{\tilde{p}}_{\text{s}}\equiv \frac{p_{\text{s}}}{p_0},\hfill & \tilde{p}\equiv \frac{p}{p_0},\hfill & \tilde{\mathfrak{p}}\equiv \frac{\mathfrak{p}}{p_0},\hfill & {\tilde{\upsilon }}_{\theta }\equiv \frac{{\upsilon }_{\theta }}{{\upsilon }_0},\hfill \\ \tilde{\dot{\sigma }}\equiv \text{Δ}t\dot{\sigma },\hfill & \tilde{T}\equiv \frac{T}{T_0},\hfill & \tilde{\Theta }\equiv \frac{\Theta }{{\Theta }_0},\hfill & \tilde{\rho }\equiv \frac{\rho }{{\rho }_0},\hfill \\ \tilde{\varphi }\equiv \frac{\varphi }{e_0},\hfill & {\tilde{F}}_{\theta }\equiv \frac{F_{\theta }}{F_0},\hfill & \tilde{Q}\equiv \frac{Q}{Q_0},\hfill & \tilde{q}\equiv \frac{q}{q_0},\hfill \\ \tilde{\dot{q}}\equiv \frac{\dot{q}}{{\dot{q}}_0},\hfill & \tilde{\theta }\equiv \frac{\theta }{\pi }.\hfill & \hfill & \hfill \end{array}$$(A.3)

As a second step, the spatial coordinates ($\tilde{\theta },\sigma$) are converted to grid coordinates⁴, (Y, Z). In this transformation, the normalised colatitude $\tilde{\theta }=\tilde{\theta }\left(Y\right)$ and vertical coordinate σ = σ (Ζ) are assumed to be monotonie functions of Y and Z, respectively. These functions are defined further so that the interval between the two boundaries of a finite volume is equal to 1 both in the horizontal and in the vertical directions, as showed by Fig. A.1. The transformation ($\left(\tilde{\theta },\sigma \right)\to \left(Y,Z\right)$) is defined at any point by the horizontal and vertical metric coefficients, C_Y and C_z, defined as $$\begin{array}{llllll}{c}_{\text{Y}}\equiv \frac{\partial \tilde{\theta }}{\partial Y},\hfill & \hfill & \hfill & \hfill & \hfill & c_{\text{Z}}\equiv \frac{\partial \sigma }{\partial Z},\hfill \end{array}$$(A.4)

which yields the change relations of partial derivatives $$\begin{array}{llllll}\frac{\partial }{\partial \tilde{\theta }}=\frac{1}{c_Y}\frac{\partial }{\partial Y},\hfill & \hfill & \hfill & \hfill & \hfill & \frac{\partial }{\partial \sigma }=\frac{1}{c_Z}\frac{\partial }{\partial Z},\hfill \end{array}$$(A.5) $$\begin{array}{llllll}\frac{\partial }{\partial Y}=c_Y\frac{\partial }{\partial \tilde{\theta }},\hfill & \hfill & \hfill & \hfill & \hfill & \frac{\partial }{\partial Z}=c_{\text{Z}}\frac{\partial }{\partial \sigma }.\hfill \end{array}$$(A.6)

Fig. A.1 Change of coordinate, ($\left(\tilde{\theta },\sigma \right)\to \left(Y,Z\right)$), from the physical spatial coordinates to the grid spatial coordinates. The grid coordinates vary in the range 0 ≤ Y ≤ M and 0 ≤ Ζ ≤ Ν, where M and Ν are two integers that do not need to be defined at this stage but will correspond to the numbers of horizontal and vertical grid levels in the finite-volume discretisation of the primitive equations. The metric is illustrated by the spacing variation between two isolines. The interval between two isolines in grid coordinates is of size 1 to mark the cells of the finite-volume method employed in the following.

We remark that the change of coordinate is just a dilation if the adopted horizontal and vertical spacings are uniform since c_Y and c_Z are constants in this case. We also notice that c_z < 0 with the chosen sigma coordinate given that σ decreases monotonically as the altitude increases.

To include the metric in the primitive equations, we introduce the respective covariant and contravariant horizontal velocities, $$\begin{array}{llllll}\widehat{\upsilon }\equiv c_Y{\tilde{\upsilon }}_{\theta },\hfill & \hfill & \hfill & \hfill & \hfill & \stackrel{⌣}{\upsilon }\equiv \frac{{\tilde{\upsilon }}_{\theta }}{c_Y},\hfill \end{array}$$(A.7)

the normalised area density (area per unit length), $$\tilde{\mathcal{A}}\equiv c_Y\sin \left(\theta \right)$$(A.8)

mass density⁵, $$\tilde{m}\equiv -\tilde{\mathfrak{p}}\tilde{\mathcal{A}}{c}_Z$$(A.9)

and horizontal and vertical mass flux densities, $$\begin{array}{lllll}V\equiv \tilde{m}\stackrel{⌣}{\upsilon },\hfill & \hfill & \hfill & \hfill & W\hfill \end{array}\equiv \tilde{m}\frac{\tilde{\dot{\sigma }}}{c_Z}.$$(A.10)

After the above manipulations, the system of HPEs given by Eqs. (15-18) and Eq. (A.1) become the normalised primitive equations $$\frac{\partial \tilde{m}}{\partial \tilde{t}}+b\frac{\partial V}{\partial Y}+\frac{\partial W}{\partial Z}=0,$$(A.11) $$\begin{array}{c}\frac{\partial }{\partial \tilde{t}}\left(\tilde{m}{\tilde{\upsilon }}_{\theta }\right)+b\frac{\partial }{\partial Y}\left(V{\tilde{\upsilon }}_{\theta }\right)+\frac{\partial }{\partial Z}\left(W{\tilde{\upsilon }}_{\theta }\right)\\ +\frac{b}{c_Y}\tilde{m}\left(\frac{\partial \tilde{\varphi }}{\partial Y}+\widetilde{\text{Θ}}\frac{\partial \tilde{E}}{\partial Y}\right)=\tilde{m}{\tilde{F}}_{\theta },\end{array}$$(A.12) $$\frac{\partial }{\partial \tilde{t}}\left(\tilde{m}\widetilde{\text{Θ}}\right)+b\frac{\partial }{\partial Y}\left(V\widetilde{\text{Θ}}\right)+\frac{\partial }{\partial Z}\left(W\widetilde{\text{Θ}}\right)=\frac{\tilde{m}\tilde{Q}}{\tilde{E}},$$(A.13) $$\frac{\partial \tilde{\varphi }}{\partial Z}+k\widetilde{\text{Θ}}+\frac{\partial \tilde{E}}{\partial Z}=0,$$(A.14) $$\frac{\partial }{\partial \tilde{t}}\left(\tilde{m}\tilde{q}\right)+b\frac{\partial }{\partial Y}\left(V\tilde{q}\right)+\frac{\partial }{\partial Z}\left(W\tilde{q}\right)=\tilde{m}\tilde{\dot{q}},$$(A.15)

where the normalised Exner function, derived from Eq. (20), is expressed as a function of the normalised pressure $\tilde{p}$, $$\tilde{E}\equiv \frac{E}{E_0}={\tilde{p}}^K.$$(A.16)

As may be noticed, with the normalisation adopted in Eq. (A.3), the dynamics of the non-dimensional HPEs is controlled by one unique dimensionless parameter, $$b\equiv \frac{{\upsilon }_0\text{Δ}t}{\pi R_{\text{p}}}=\frac{\sqrt{C_p{T}_0\text{Δ}t}}{\pi R_{\text{p}}},$$(A.17)

which is a global Courant number weighting the contribution of horizontal advection⁶. The density is not a variable of the primitive equations but can be evaluated using the perfect gas equation, $$\tilde{p}=K\tilde{\rho }\tilde{T},$$(A.18)

while the vertical velocity can be calculated from the conversion formula $${\upsilon }_z=\frac{{\upsilon }_{z}0}{\tilde{m}}\left[\tilde{m}\frac{\partial \tilde{\varphi }}{\partial \tilde{t}}+bV\frac{\partial \tilde{\varphi }}{\partial Y}+W\frac{\partial \tilde{\varphi }}{\partial Z}\right],$$(A.19)

where we have introduced the characteristic vertical velocity $${\upsilon }_{z}0\equiv \frac{e_0}{g\text{Δ}t}.$$(A.20)

Appendix B Dynamical core

The system of prognostic equations given by Eqs. (A.11-A.15) is integrated numerically using a standard finite-volume method. This appendix details the main features of the dynamical core of our model.

Appendix B.1 Grid

The primitive equations are discretised and solved on a staggered Arakawa C grid (Arakawa & Lamb 1977). The atmosphere is divided into elemental cells. Mass fluxes and velocities are evaluated at cell interfaces, and volume quantities ($\left(\tilde{m},\tilde{\Theta },\tilde{q},\tilde{\varphi },\tilde{E}\right)$) at cells centres, except the pressure, which is evaluated at horizontal cells interfaces. A diagram representing the grid is shown by Fig. B.1.

As a first step, the spatial domain defined by ($\left(\tilde{\theta },\sigma \right)\in {\left[0,1\right]}^2$) is divided into non-uniform intervals both along the horizontal and vertical directions. There are M intervals for the colatitude, and N intervals for the vertical coordinate, which represents M × N cells. The vertical spacing can be specified arbitrarily. It determines the coordinates of cell boundaries. We note that the σ levels of cell centres can be explicitly computed from the σ levels of cell boundaries if p_t = 0 (standard sigma coordinate) as discussed in Sect. B.2.

As as second step, the change of coordinates ($\left(\tilde{\theta },\sigma \right)\to \left(Y,Z\right)\in \left[0,M\right]\times \left[0,N\right]$) is applied, and the corresponding metric coefficients c_Y and c_Z are calculated. Typically, for a uniform grid, $\tilde{\theta }=Y/M$, and thus c_Y = 1/M everywhere. The case of non-uniform horizontal intervals is more complicated, since this implies $\tilde{\theta }$-dependent coefficients. To treat the general case, we adapt the method employed for the horizontal grid in the LMDZ GCM to the 1D case: each surface area is divided into two subsurface areas, as shown by Fig. B.2, and the C_Y are evaluated at subarea centres and boundaries (Sadourny 1975a,b). The c_Z are just the difference between two vertical levels. We introduce the difference and average operators, $${\delta }_X\psi \left(X\right)=\psi \left(X+\frac{1}{2}\right)-\psi \left(X-\frac{1}{2}\right),$$(B.1) $${\overline{\psi }}^X\left(X\right)=\frac{1}{2}\left[\psi \left(X+\frac{1}{2}\right)+\psi \left(X-\frac{1}{2}\right)\right],$$(B.2)

where X is a placeholder for Y or Z and ψ any variable. In this formalism, the c_Z are simply expressed as $c_Z={\delta }_Z\sigma$.

In the present work, we use a vertical grid refined near the surface and the upper boundary. The sigma coordinates of vertical levels are given by the function $$\sigma \left(x\right)=\frac{1}{2}\left[1+\cos \left(\pi x^a\right)\right],$$(B.3)

where x ≡ Z/N takes its values between 0 and 1. Here, the exponent a is a dilation coefficient that controls the growth or decay rate of vertical intervals in the vicinity of the lower and upper bounds. Introducing the vertical coordinate difference between the lowest model level and the surface, ∆σ_SL:, this parameter is defined as $$a\equiv \frac{\text{1n}[{\pi }^{-1}\arccos \left(1-2\text{Δ}{\sigma }_{\text{SL}}\right)}{\text{1n}\left(1/N\right)}.$$(B.4)

For an isothermal temperature profile near planet’s surface, $\text{Δ}{\sigma }_{\text{SL}}\approx z_{\text{SL}}/H$, where H is the local pressure height. In practice, we set $\text{Δ}{\sigma }_{\text{SL}}=3.0\times {10}^{-3}$.

Fig. B.1 Staggered grid. In the adopted finite-volume method, the horizontal axis is divided into M intervals and the vertical axis into N intervals, which represents M × N cells (here, M = 5 and N = 4). Horizontal levels are indexed by j and vertical levels by k. The left and right boundaries of the physical domain correspond to $\tilde{\theta }=0$ (sub-stellar point) and $\tilde{\theta }=1$ (anti-stellar point), respectively. The top and bottom boundaries correspond to σ = 0 (space) and σ = 1 (ground), respectively. The horizontal and vertical mass flows, V and W, are evaluated at vertical and horizontal cell interfaces, respectively, while the mass, temperature, geopotential, specific humidity, and Exner function (blue) are evaluated at cell centres. The physical domain is surrounded by ghost cells (grey), which are employed to facilitate vectorisation.

Appendix B.2 Exner Function

The calculation of the Exner function follows the method described by Arakawa & Lamb (1977) and Hourdin (1994), which is summarised here. Instead of computing E at the middle of the layer by extrapolating the pressure at the middle of the layer from the level pressures at the interfaces, which is computationally expensive, one rather uses a supplementary relationship between the interface levels and the Exner function. This relationship is directly derived from considerations about energy conservation principles within the atmospheric air column. Particularly, in the hydrostatic approximation, the internal and potential energy are proportional, which implies (Arakawa & Lamb 1977); Hourdin 1994) $${\displaystyle {\int }_0^{m_{\text{col}}}\varphi \text{d}m}={\displaystyle {\int }_0^{m_{\text{col}}}{R}_{\text{d}}T\text{d}m},$$(B.5)

where m = pdz is the infinitesimal parcel of mass per unit surface, and m_col the mass of the air column per unit surface. The expression of ф as a function of Θ and E is simply obtained by making use of the hydrostatic balance equation given, in grid vertical coordinates, by $$\frac{\partial \varphi }{\partial Z}+\text{Θ}\frac{\partial E}{\partial Z}=0,$$(B.6)

and reads $$\varphi ={\displaystyle {\int }_0^Z\frac{\partial \varphi }{\partial Z^{\prime }}\text{d}{Z}^{\prime }}=-{\displaystyle {\int }_0^{Z}K\text{Θ}\frac{\partial E}{\partial Z^{\prime }}\text{d}{Z}^{\prime }},$$(B.7)

which allows us, noticing that R_dT = κΘΕ, to rewrite Eq. (B.5) as $${\displaystyle {\int }_0^{m_{\text{col}}}\left[{\displaystyle {\int }_0^Z\text{Θ}\frac{\partial E}{\partial Z^{\prime }}\text{d}{Z}^{\prime }}+k\text{Θ}E\right]}\text{\hspace{0.17em}\hspace{0.17em}d}m=0.$$(B.8)

Fig. B.2 Subdivision of a surface area into two subsurface areas, A1 and A2 (light grey and grey regions). The big blue dots indicate the centres of horizontal intervals. The areas are calculated analytically by taking their dependence on colatitude into account. The configuration of the diagram corresponds to small horizontal intervals, where this dependence is approximately linear.

The discretised form of the first term is given by $${\displaystyle {\int }_0^{Z_l}\text{Θ}\frac{\partial E}{\partial Z}\text{d}{Z}^{\prime }}={\displaystyle \sum _{k=1}^k\left(\frac{{\text{Θ}}_k+{\text{Θ}}_{k-1}}{2}\right)}\left(E_k-E_{k-1}\right)+{\text{Θ}}_0\left(E_0-E_{\text{s}}\right),$$(B.9)

where Θ_l, E_l are the potential temperature and Exner function evaluated at the middle of the layer, indexed by l = 0, . . . , N − 1 (for N vertical intervals), and E_s = C_p (p/p_ref)^κ is the Exner function evaluated at the planet’s surface. By making use of the difference and average operators introduced in Eqs. (B.1) and (B.2), the preceding equation can be rewritten in the compact form $${\displaystyle {\int }_0^{Z_l}\text{Θ}\frac{\partial E}{\partial Z}\text{d}{Z}^{\prime }={{\displaystyle \sum _{l=0}^{N-1}\left[{\overline{\text{Θ}}}^Z{\delta }_{Z}E\right]}}_k,}$$(B.10)

with $$\begin{array}{lll}{\left[{\overline{\text{Θ}}}^Z{\delta }_{Z}E\right]}_k\hfill & =\left(\frac{{\text{Θ}}_k+{\text{Θ}}_{k-1}}{2}\right)\left(E_k-E_{k-1}\right),\hfill & k=1,\dots ,N-1,\hfill \\ {\left[{\overline{\text{Θ}}}^Z{\delta }_{Z}E\right]}_0\hfill & ={\text{Θ}}_0\left(E_0-E_s\right).\hfill & \hfill \end{array}$$(B.11)

Then, introducing the mass m_l per unit area of layer l and interchanging the sums in the double integral, we obtain $${\displaystyle \sum _{l=0}^{N-1}{m}_l{\displaystyle {\int }_0^{Z_l}\text{Θ}\frac{\partial E}{\partial Z}\text{d}{Z}^{\prime }={\displaystyle \sum _{l=0}^{N-1}{m}_l{{\displaystyle \sum _{k=0}^l\left[{\overline{\text{Θ}}}^Z{\delta }_{Z}E\right]}}_k,}}}$$(B.12) $$={\displaystyle \sum _{k=0}^{N-1}{\left[{\overline{\text{Θ}}}^Z{\delta }_{Z}E\right]}_k{\displaystyle \sum _{l=k}^{N-1}{m}_l.}}$$(B.13)

In the above equation, we recognise the atmospheric mass per unit surface above the k interfaces (k = 0‚…, N), which is expressed in the framework of the hydrostatic approximation as $${\displaystyle \sum _{l=k}^N{m}_l=\frac{p_k}{g}.}$$(B.14)

It follows $${\displaystyle \sum _{l=0}^{N-1}{m}_l}{\displaystyle {\int }_0^{Z_l}\text{Θ}\frac{\partial E}{\partial Z}d{Z}^{\prime }=\frac{1}{g}{\displaystyle \sum _{k=0}^{N-1}{\left[{\overline{\text{Θ}}}^Z{\delta }_{Z}E\right]}_k{p}_k.}}$$(B.15)

Finally, expanding the coefficients of the sum, we remark that this later may be rewritten $${{\displaystyle \sum _{k=0}^{N-1}\left[{\overline{\text{Θ}}}^Z{\delta }_{Z}E\right]}}_k{p}_k={\displaystyle \sum _{k=0}^{N-1}{\text{Θ}}_k{\left[{\overline{p{\delta }_{Z}E}}^Z\right]}_k,}$$(B.16)

with the conventions $${\left[{\overline{p{\delta }_{Z}E}}^Z\right]}_0=\frac{1}{2}\left(E_1-E_0\right)p_1+\left(E_0-E_{\text{s}}\right)p_0$$(B.17) $${\left[{\overline{p{\delta }_{Z}E}}^Z\right]}_{N-1}=\frac{1}{2}\left(E_{N-1}+E_{N-2}\right)p_{N-1}+\left(E_t-E_{N-1}\right)p_t.$$(B.18)

Thus, in the general case, the Exner function can be evaluated at cell centres by applying, at all vertical levels, the relation $${\overline{p{\delta }_{Z}E}}^Z=-kE{\delta }_{Z}p.$$(B.19)

The use of standard sigma coordinates appreciably reduces the complexity of the problem. Assuming pt = 0 and introducing the notation s = σ^κ, the preceding equation simplifies to $${\overline{\sigma {\delta }_{Z}E}}^Z=kS{\delta }_Z\sigma ,$$(B.20)

which allows the sigma coordinates of cell centre levels to be pre-calculated. Indexing these levels by k = o, …,N – 1 and the intermediate (or interface) σ levels by k = o,…, N , we have $$\begin{array}{rr}\hfill {\left[{\overline{\sigma {\delta }_{Z}s}}^Z\right]}_k& =\kappa s_k{\left[{\delta }_Z\sigma \right]}_{k}fork=1,\dots ,N-2,\hfill \\ \hfill {\left[{\overline{\sigma {\delta }_{Z}s}}^Z\right]}_0& =\frac{1}{2}\left(s_1-s_0\right){\sigma }_1+\left(s_0-s_{\text{s}}\right){\sigma }_0=\kappa s_0\left({\sigma }_1-{\sigma }_0\right),\hfill \\ \hfill {\left[{\overline{\sigma {\delta }_{Z}s}}^Z\right]}_{N-1}& =\frac{1}{2}\left(s_{N-1}-s_{N-2}\right){\sigma }_{N-1}=\kappa s_{N-1}\left({\sigma }_N-{\sigma }_{N-1}\right)\hfill \end{array}$$(B.21)

where s_s = 1 and s_t = 0 are the value of s at planet’s surface and atmospheric upper boundary, respectively.

In practice, the Exner function is evaluated at grid centre levels by solving an algebraic equation of the form AX = B. In the case of standard sigma coordinates (p_t = 0), this calculation is performed only once, when the grid is constructed. The algebraic system to solve then writes $$\left[\begin{array}{ccccc}{A}_{0,0}& A_{0,1}& & & \\ A_{1,0}& A_{1,1}& A_{1,2}& & \\ & A_{k,k-1}& A_{k,k}& A_{k,k+1}& \\ & & A_{N-2,N-3}& A_{N-2,N-2}& A_{N-2,N-1}\\ & & & A_{N-1,N-2}& A_{N-1,N-1}\end{array}\right]\left[\begin{array}{c}{s}_0\\ s_1\\ s_k\\ s_{N-2}\\ s_{N-1}\end{array}\right]=\left[\begin{array}{c}{B}_0\\ B_1\\ B_k\\ B_{N-2}\\ B_{N-1}\end{array}\right],$$(B.22)

where the s_k is the coordinate s = σ^κ evaluated at the mid-level of the k-layer, A the tri-diagonal matrix of coefficients $$A_{k,k}=\left(\frac{1}{2}+\kappa \right)\left({\sigma }_k-{\sigma }_{k+1}\right)\text{\hspace{0.17em}\hspace{0.17em}for\hspace{0.17em}}k=1,\dots ,N-2;$$(B.23) $$A_{0,0}=\left(1+\kappa \right){\sigma }_0-\left(\frac{1}{2}+\kappa \right){\sigma }_1,$$(B.24) $$A_{N-1,N-1}=\left(\frac{1}{2}+K\right){\sigma }_{N-1}-\left(1+K\right){\sigma }_N,$$(B.25) $$A_{k,k+1}=-A_{k+1,k}=\frac{1}{2}{\sigma }_{k+1},$$(B.26)

and B the vector of coefficients $$B_k=0\text{\hspace{0.17em}\hspace{0.17em}for\hspace{0.17em}\hspace{0.17em}}k=1,\dots ,N-2;$$(B.27) $$B_0={\sigma }_0{s}_{\text{s}},$$(B.28) $$B_{N-1}=-{\sigma }_N{s}_{\text{t}}.$$(B.29)