# Convective self-aggregation as a two-regime damped gravity wave

Gilles Bellon, Aurélien Ribes, Benoit Meyssignac, Olivier Geoffroy, Grégory

Faye, Pascal Noble

► To cite this version:

Gilles Bellon, Aurélien Ribes, Benoit Meyssignac, Olivier Geoffroy, Grégory Faye, et al. Convective self-aggregation as a two-regime damped gravity wave. 2025. hal-05166602

HAL Id: hal-05166602

https://hal.science/hal-05166602v1

Preprint submitted on 17 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Convective self-aggregation as a two-regime damped gravity wave

Gilles Bellon<sup>a*</sup>, Aurélien Ribes<sup>a</sup>, Benoit Meyssignac<sup>b</sup>, Olivier Geoffroy<sup>a</sup>,

Grégory Faye<sup>c</sup>, Pascal Noble<sup>d</sup>

<sup>a</sup>Centre National de Recherches Météorologiques,
Centre National de la Recherche Scientifique – Météo France – Université de Toulouse, Toulouse, France.

<sup>b</sup>Laboratoire d'Etudes en Géophysique et Océanographie Spatiales,
Université de Toulouse, CNRS – CNES – UPS – IRD, Toulouse, France.

<sup>c</sup>Institut de Mathématiques de Toulouse; UMR5219, Université de Toulouse; CNRS, UPS, F-31077 Toulouse, France.

<sup>d</sup>Institut de Mathématiques de Toulouse; UMR5219, Université de Toulouse; CNRS, INSA, F-31077 Toulouse, France.

*Correspondence to: CNRM, 42 avenue Gaspard Coriolis, 31057 Toulouse, France. E-mail: gilles.bellon@meteo.fr

Convective self-aggregation is the spontaneous spatial organisation of deep-convective clouds into a limited region surrounded by drier, convectively-inhibited regions, which occurs in many models. We propose a simple, piecewise linear model for self-aggregation based on primitive equations. In this model, each atmospheric column is in one of two possible thermodynamic regimes, deep-convective or convectively inhibited, and the thermodynamics in each regime is linearised. The model simulates aggregated and non-aggregated steady states and reproduces many properties of self-aggregation as simulated by kilometre-resolution models. In particular, it exhibits an hysteresis with multiple, aggregated and non-aggregated equilibria and a similar sensitivity to convective inhibition, domain size, and boundary-layer radiative cooling in the dry region as in kilometre-resolution simulations. These results suggest that a self-aggregated state can be considered as a simple gravity wave with two phases: one convective and one convectively inhibited.

Key Words: Convective self-aggregation; gravity wave; convective inhibition; hysteresis; multiple equilibria.

Received ...

## 1. Introduction

Convective self-aggregation is a phenomenon well-documented in atmospheric model simulations, in which deep convective clouds regroup spontaneously in a fraction of the domain surrounded by

a drier region devoid of deep convection (see reviews by Wing et al. 2017; Muller et al. 2022). By “spontaneously”, we mean that this spatial organisation of convection occurs despite horizontally uniform initial and boundary conditions. The main example
of such configurations is the radiative-convective equilibrium (RCE), in which surface and top-of-the-atmosphere conditions are uniform and atmospheric exports from the domain are set to zero. Convective self-aggregation was first documented in limited-area, kilometre-resolution models* (Held et al. 1993) and later in general-circulation models (Held et al. 2007). It has now been widely documented via the Radiative-Convective Equilibrium model intercomparison project (Wing et al. 2018, 2020).

Analysis of RCE simulations has revealed a plurality of processes contributing to self-aggregation: solar and terrestrial radiation, wind-induced surface fluxes, convective downdraughts, shallow circulations. How important or crucial each of these processes is depends on models, model configurations, and/or boundary conditions (Wing et al. 2017, and references therein). Some factors are robust: for example, the boundary-layer radiative cooling in the dry region due to low clouds and/or to the clear-sky effect of a very dry free-troposphere is crucial, because it forces a shallow circulation that transport energy up-gradient from the low-energy dry regions to the high-energy convective regions (Bretherton et al. 2005; Muller and Held 2012; Muller and Bony 2015; Yang 2018). A moisture feedback also contributes to self-aggregation: environmental moisture in the convective region favours convection there through the entrainment into convective updraughts (Tompkins and Semie 2017).

A few robust sensitivities of self-aggregation have been identified: the larger the domain and the coarser the horizontal resolution, the more likely self-aggregation is to occur (Bretherton et al. 2005; Muller and Held 2012; Muller and Bony 2015; Wing et al. 2020). Interestingly, many sensitivities are non-linear and exhibit a hysteresis: the parameter interval for which a non-aggregated RCE can be obtained overlaps with the interval for which an aggregated RCE can be obtained (Bretherton et al. 2005; Muller and Held 2012; Muller and Bony 2015; Coppin and Bony 2015; Bellon and Coppin 2022). For a parameter in the overlapping interval, the stationary state reached by the model is a function of initial conditions only.

Many simple models attempt to explain self-aggregation. Some focus on the instability of the non-aggregated RCE to the

interaction of radiative-convective diabatic processes with the circulation in single-column models (Sobel et al. 2007; Sessions et al. 2010; Emanuel et al. 2014; Beucler and Cronin 2016; Beucler et al. 2018; Kuang 2018) or two-column models (Nilsson and Emanuel 1999; Raymond and Zeng 2000). Others focus on scale interaction such as between cold pools and convective clouds (Haerter 2019) or small-scale gravity waves and convective envelopes (Yang 2021). Models based on a reaction-diffusion equation applied to precipitable water presents self-aggregation as a coarsening process (Craig and Mack 2013; Windmiller and Craig 2019; Biagioli and Tompkins 2023). The last ones and some two-column models exhibit a hysteresis, but the others do not, or are not designed to, reproduce the non-linear behaviour documented in kilometre-resolution models.

An aggregated RCE can be considered as a system with one convective region and one convectively-inhibited, dry region. There is a significant tradition of modelling the tropics using two-column models, to study the tropical climate as a whole (Pierrehumbert 1995; Bellon et al. 2003), the Walker circulation (Kelly and Randall 2001; Bretherton and Sobel 2002), or, as pointed above, the instability of the non-aggregated RCE (Nilsson and Emanuel 1999; Raymond and Zeng 2000). Here we revisit this two-column framework to model the RCE and to attempt to reproduce the non-linear behaviour of kilometre-resolution models.

We base our model on the primitive equations and try to keep it conceptually as simple as possible. We use linearised dynamic and thermodynamic equations in the convective and dry regions and the only non-linearity in our model is the onset of convection at the boundary between the dry and the moist regions. Section 2 describes the model and the method of solution. Section 3 displays the solutions and their sensitivity to the parameters. Finally, Section 4 provides our conclusion.

## 2. Model formulation and solution

### 2.1. Basic equations

We will consider a two-dimensional horizontal-vertical plane $(x, z)$. The non-rotating primitive equations linearised about a thermodynamically stratified state of rest can be written in terms

<footnote>*also known as Cloud-Resolving Models (CRM) or sometimes as Storm-Resolving Models (SRM)</footnote>
of specific momentum (u, w) (density times velocity):

$$ \partial_x u + \partial_z w = 0, \quad (1) $$

$$ \partial_t u = -\partial_x p' - \epsilon u, \quad (2) $$

$$ 0 = -\partial_z p' + \rho_0 \frac{\theta'_\rho}{\theta_{\rho_0}} g, \quad (3) $$

$$ \partial_t \theta'_{\rho} + w \partial_z \theta_{\rho_0} = Q_{\theta} - \mu \theta'_{\rho}, \quad (4) $$

with the usual notations: u the horizontal momentum, w the vertical momentum, p the pressure, ρ the density, θρ = ρRd/p (pρ/p)Rd/Cρ the density (or virtual) potential temperature, with pρ a reference pressure, and Qθ the reference-state diabatic heating rate of density potential temperature. g is the acceleration of gravity, Rd is the gaz constant for dry air and Cρ is the heat capacity of air. ϵ and μ are the coefficients of, respectively, the linear Rayleigh damping and the linear Newtonian cooling, which account for the dissipation of, respectively, momentum and thermal energy. For simplicity, they are assumed to be independent of altitude. The subscript "0" indicates the reference-state values and the prime indicates the deviation from this reference state. Equation (1) is the continuity equation expressing the conservation of mass, (2) is the horizontal momentum equation, (3) is the hydrostatic approximation in which we use the additional approximation ρ'/ρ0 ≈ -θ'ρ/θρ0, and (4) is the thermodynamic equation. Without rotation and with their damping terms, these linear equations form the simplest model of atmospheric dynamics: they describe damped gravity waves, possibly forced by a non-zero diabatic heating Qθ.

Combining equations (1)-(3) yields the following fluid-dynamics equation:

$$ (\partial_t + \epsilon) \partial_z^2 w = \partial_x^2 b, \quad (5) $$

with b = ρ0(θ'ρ/θρ0)g the buoyancy.

The thermodynamic equation (4) can be rewritten in terms of buoyancy:

$$ \partial_t b + N^2 w = Q - \mu b \quad (6) $$

with Q = ρ0/θ0 g Qθ the reference-state diabatic heating in buoyancy units and N the Brunt-Väisälä frequency:

$$ N^2 = \frac{g}{\theta_{\rho_0}} \partial_z \theta_{\rho_0}. $$

In the atmosphere, N ≈ 10⁻² Hz, which corresponds to a period of less than 2 minutes. We set N = 10⁻² Hz uniform and non-dimensionalise the system with the following scales:

$$ [x] = 1200 \text{ km}, $$

$$ [z] = 12 \text{ km}, $$

$$ [t] = \frac{1}{N} \frac{[x]}{[z]} = 10^4 \text{ s} \approx 3 \text{ hr}, $$

$$ [u] = 10 \text{ kg m}^{-2} \text{ s}^{-1}, $$

$$ [w] = \frac{[z]}{[x]} [u] = 0.1 \text{ kg m}^{-2} \text{ s}^{-1}, $$

$$ [b] = \frac{[x]}{[z]} N[w] = 0.1 \text{ kg m}^{-2} \text{ s}^{-2}, $$

$$ [Q] = N^2[w] = 10^{-5} \text{ kg m}^{-2} \text{ s}^{-3}. $$

We used the large aspect ratio [x]/[z] to obtain a time scale [t] much larger than 1/N since we are interested in a phenomenon that occurs on time scales of days or weeks rather than minutes.

Using the same notations for the non-dimensional variables as previously for the dimensional ones, we obtain:

$$ (\partial_t + \epsilon) \partial_z^2 w = \partial_x^2 b, \quad (7) $$

$$ (\partial_t + \mu)b + w = Q, \quad (8) $$

which just eliminates the Brunt-Väisälä frequency from the equation set (5)-(6).

## 2.2. Two thermodynamic regimes

For each horizontal location, we consider that either deep convection occurs or that it is inhibited. These two convective regimes differ by their thermodynamic processes. At a deep-convective location, convective mixing, condensation, cloud and clear-sky radiative processes interact to adjust the atmosphere towards the thermodynamic equilibrium profiles in the absence of domain-scale atmospheric circulation; by definition, these are the non-aggregated RCE profiles. At a convectively inhibited location, the convective mixing is limited to the
boundary layer and radiative processes are dominant; the thermodynamic processes adjust the thermodynamic profiles to radiative equilibrium profiles controlled by the prescribed surface temperature. In the following, we will refer to the regime with deep convection as the convective regime and the convectively inhibited regime as the dry regime.

In our equations, these two thermodynamic regimes translates into different diabatic heating rates $Q$ and different Newtonian damping rates $\mu$. We will use the superscripts "c" to indicate the convective regime and "d" to indicate the dry regime. $Q^d$ and $Q^c$ depend on the choice of reference state used for linearisation. We assume that a non-aggregated RCE (i.e., with convection everywhere and no domain-scale circulation: $(u, w) = (0, 0)$) exists, for example for a small domain on the horizontal, and we use this non-aggregated RCE profile as our reference state.

In the convective regime, the balance between the diabatic heating due to radiative and convective processes is by definition the same as the one in the reference state so it dictates the total diabatic heating to be zero: $Q^c = 0$. To choose $\mu^c$, we apply a simple version of the relaxed quasi-equilibrium principle (Arakawa and Schubert 1974; Betts 1986; Moorthi and Suarez 1992; Emanuel et al. 1994; Yano and Plant 2012), which stipulates that convection eliminates convective instability very rapidly with an adjustment time orders of magnitude shorter than other processes. We therefore consider that convective processes dominate the adjustment to the reference profile where convection is active and that they damp buoyancy anomalies very efficiently, which corresponds to a large damping rate $\mu^c$. We choose $\mu^c = 1$ in the control case, which is equivalent to a relaxation time scale of about 3 hours, similar to the values in the literature (Arakawa and Schubert 1974; Betts and Miller 1986; Emanuel et al. 1994).

In the dry regime, the troposphere relaxes towards radiative equilibrium, which is cooler than the non-aggregated RCE. Since the non-aggregated RCE profile is our reference state, the diabatic processes in the dry regime create negative buoyancy perturbations with respect to the reference state. If we note $\Delta b^d(z)$ the (negative) difference between the radiative equilibrium profile and reference RCE profile, we get the reference-state diabatic heating in the dry regime $Q^d = \mu^d \Delta b^d < 0$. It corresponds to a cooling because, with respect to a convective reference

state, inhibiting convection removes a latent and free-tropospheric cloud-radiative heating and also allows more radiative cooling from the boundary layer because of the reduced greenhouse effect of the clear-sky and drier free troposphere above. We will represent the convective latent and cloud-radiative heating by a combination of a first-baroclinic ("convective") and second-baroclinic ("stratiform") heating profiles (Tao et al. 2006; Hagos et al. 2010; Bellon and Bony 2020) and the boundary-layer radiative cooling by a heat sink confined to the first two kilometres of the atmosphere (up to the inversion at $z_i = 1/6$):

$$Q^d = -Q_{cv} [\sin(\pi z) - \alpha_s \sin(2\pi z)] \\ -Q_{bl} \sin^2 \left[ \pi \min \left( \frac{z}{z_i}, 1 \right) \right]. \quad (9)$$

$Q_{cv}$ is the amplitude of the convective-cloud heating (latent and free-tropospheric cloud radiative heating); $\alpha_s$ is the ratio between the amplitudes of the stratiform heating and the convective heating, which controls how top-heavy the profile of convective-cloud heating is; $Q_{bl}$ is the amplitude of the boundary-layer radiative cooling. Figure 1 shows these profiles for the control values of parameters ($Q_{cv}$, $\alpha_s$, $Q_{bl}$) = (5 · 10⁻², 0.5, 2 · 10⁻²). This control value of $Q_{cv}$ corresponds to 60 W m⁻² over the depth of the troposphere or a precipitation of about 2 mm/day; and this control value of $Q_{bl}$ corresponds to a cooling of about 0.2 K/day in average over the atmospheric boundary layer. Since radiative adjustments are much slower than the convective adjustment, the thermodynamic damping coefficient in the dry regime is therefore much smaller than in the convective regime. In the control case, we choose $\mu^d = 0.02$, which corresponds to a time scale of about a week, an intermediate value between the time scales of Newtonian cooling in the free troposphere (a few weeks) and in the boundary layer (a few days) used in other models (e.g., Bellon and Sobel 2008).

Convection can be suppressed by a lower-tropospheric stable stratification. The canonical metric of such a stratification is convective inhibition (CIN), the kinetic energy necessary for a surface air parcel to reach, by adiabatic uplift, the level of free convection, where it becomes positively buoyant. Considering that our model does not represent the boundary-layer processes explicitly, CIN is too refined a metric to be used. Instead, we
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Profile $Q^d(z)$ of the reference-state, non-dimensional diabatic heating in the dry region for the control parameters (black line), and its decomposition into boundary-layer radiative contribution (dashed orange) and convective-cloud contribution (dashed blue). The latter is also decomposed into its convective and stratiform components (dotted blue).</figcaption></figure>

consider the buoyancy stratification between the surface and the top of the boundary layer at $z_i$ to be a rough equivalent of convective inhibition and use a threshold $\delta b_i$ on the buoyancy difference between $z = z_i$ and $z = 0$ to determine whether the atmospheric column is convective. Since the anomaly of $b$ is zero at the surface, the criterion results in:

Convection at location $x$ and instant $t$ is:
active if $b(x, z_i, t) \le \delta b_i$,
inhibited if $b(x, z_i, t) > \delta b_i$.

A small value of the threshold $\delta b_i$ corresponds to easy convective inhibition, i.e., difficult deep convection. Conversely, a large $\delta b_i$ corresponds to difficult inhibition or easy convection.

## 2.3. A two-column model

In designing our horizontal domain, we allow for one region to be convectively suppressed and one region to be convectively active. We consider the simplest geometrical configuration: a domain of size $2X$ ($x \in [-X, X]$) in which convection is possibly inhibited in a central region $[-\xi, \xi]$ and active outside of that region. The location $\xi$ of the transition between the dry and the convective region is an unknown, to be determined by the solution; if $\xi = 0$, convection occurs everywhere and is not aggregated; if $\xi > 0$, there is a region where convection is inhibited and the domain is in a self-aggregated state.

We impose periodic boundary conditions at $x = \pm X$, similarly to the boundary conditions in RCE simulations with limited-area, kilometre-resolution models. Therefore, our configuration can represent damped gravity waves with one dry and one deep-convective regions per wavelength.

As the geometry is symmetric with respect to $x = 0$, we consider flows with the same symmetry and we look for solutions over $[0, X]$. The symmetry and periodicity imply the following geometrical boundary conditions for $b$:

$$\partial_x b(0, z, t) = 0 \text{ and } \partial_x b(X, z, t) = 0.$$

Mass conservation imposes continuity of $b$ and $\partial_x b$ for $x = \xi$.

In the vertical, we will consider the troposphere as represented by $z \in [0, 1]$, with a flat surface and a rigid lid at the top. This constrains the vertical momentum:

$$w(x, 0, t) = 0 \text{ and } w(x, 1, t) = 0.$$

## 2.4. Method of solution for the equilibrium

If $\delta b_i \ge 0$, the model has a trivial solution, convective everywhere and motionless: $\xi_e = 0$, $b_e = 0$, $w_e = 0$, with the subscript "e" denoting the equilibrium solution. It verifies the condition for convection everywhere ($\forall x, b_e(x, z_i) = 0 \le \delta b_i$) and it is therefore a non-aggregated RCE.

We can also search for an equilibrium with $\xi_e > 0$. We first assume that such an equilibrium exists and derive the equilibrium field $b_e(x, z; \xi)$ that is a steady solution to equations (7)-(8) and all the geometrical and physical boundary conditions except the convection-inhibition criterion. A posteriori, we examine whether any such solution verifies the convection-inhibition criterion.

In (7) and (8), the only dependency in $z$ is the second derivative $\partial_z^2$, which has eigenmodes in sin form. Because of the surface and top conditions on $w$, we choose a base of vertical modes in $\sin(m\pi z)$ to decompose our variables and forcings:

$$w = \sum_{m=1}^{\infty} w_m(x, t) \sin(\tilde{m}z), \quad (10)$$

with $\tilde{m} = m\pi$. Buoyancy $b$ and heating $Q$ can be decomposed similarly. In practice, we will use a truncation at 100 vertical modes in our computations.
Equations (7) and (8) can be rewritten for each vertical mode, which yields for the steady state ($\partial_t = 0$):

$$w_{em}(x) = -\frac{1}{\epsilon\tilde{m}^2}\partial_x^2 b_{em}(x), \quad (11)$$

$$b_{em}(x) = \frac{Q_m - w_{em}(x)}{\mu}, \quad (12)$$

in which $\mu$ and $Q_m$ depend on the thermodynamic regime. Equation (12) shows the equilibrium $b_e$ as the sum of a direct response to the forcing $Q_\mu$ and a circulation feedback $-w_e/\mu$. Equation (11) shows that the circulation feedback acts as a horizontal diffusion on $b$, with an intensity that decreases with the vertical wavenumber $\tilde{m}$: the circulation response is strongest for a first-baroclinic buoyancy perturbation and decreases for higher baroclinic modes.

Equations (11) and (12) can be combined into one second-order linear equation for $b_{em}$ for each mode and each region:

$$\partial_x^2 b_{em} - \epsilon \mu^a \tilde{m}^2 b_{em} = -\epsilon \tilde{m}^2 Q_m^a, \quad (13)$$

with $a = d$ for $x \in [0, \xi]$ and $a = c$ for $x \in [\xi, X]$.

For each region, $Q_m^a/\mu^a$ is a particular solution of (13) and the solutions of the corresponding homogeneous equation which verify the geometrical boundary conditions are:

$$\cosh(\tilde{m}\sqrt{\epsilon\mu^d}x) \quad \text{for } x \le \xi, \quad (14)$$

$$\cosh(\tilde{m}\sqrt{\epsilon\mu^c}[x-X]) \quad \text{for } x \ge \xi. \quad (15)$$

Note that this particular solution corresponds to the thermodynamic equilibrium (RCE in the convective region, radiative equilibrium in the dry region) while the solution to the homogeneous equation encapsulates the circulation feedback.

The continuity conditions at $x = \xi$ determines the solution for each vertical mode. We choose to write it as a function of the normalised horizontal coordinate $\hat{x} = x/X$:

$$\begin{aligned} b_{em}(\hat{x}; \hat{\xi}) &= \Delta b_m^d \left( 1 + \frac{\nu^2}{D_m} \tanh(\tilde{m}\lambda\nu[\hat{\xi} - 1]) \frac{\cosh(\tilde{m}\lambda\nu^{-1}\hat{x})}{\cosh(\tilde{m}\lambda\nu^{-1}\hat{\xi})} \right) \quad \text{for } \hat{x} \le \hat{\xi}, \\ &= \frac{\Delta b_m^d}{D_m} \tanh(\tilde{m}\lambda\nu^{-1}\hat{\xi}) \frac{\cosh(\tilde{m}\lambda\nu[\hat{x} - 1])}{\cosh(\tilde{m}\lambda\nu[\hat{\xi} - 1])} \quad \text{for } \hat{x} \ge \hat{\xi}, \end{aligned} \quad (16)$$

with $\hat{\xi}$ the fraction of the domain covered by the dry region:

convective-inhibition criterion, i.e., if and only if:

$$\hat{\xi} = \frac{\xi}{X},$$

$$b_e(\hat{x}, z_i; \hat{\xi}_e) > \delta b_i \quad \text{for } \hat{x} < \hat{\xi}_e,$$

$$b_e(\hat{\xi}_e, z_i; \hat{\xi}_e) = \delta b_i, \quad (17)$$

and:

$$\text{and } b_e(\hat{x}, z_i; \hat{\xi}_e) < \delta b_i \quad \text{for } \hat{x} > \hat{\xi}_e.$$

$$\Delta b_m^d = \frac{Q_m^d}{\mu^d}, \quad \lambda = X \epsilon^{\frac{1}{2}} (\mu^c \mu^d)^{\frac{1}{4}}, \quad \nu = \left( \frac{\mu^c}{\mu^d} \right)^{\frac{1}{4}},$$

Therefore, we look for roots $\hat{\xi}_e$ of the function

$$D_m = \tanh(\tilde{m}\lambda\nu^{-1}\hat{\xi}) - \nu^2 \tanh(\tilde{m}\lambda\nu[\hat{\xi} - 1]).$$

$$f(\hat{\xi}) = b_e(\hat{\xi}, z_i; \hat{\xi}) - \delta b_i \quad (18)$$

Note that $Q_m^c = 0$ simplifies the expressions and that $\Delta b^d = \sum_{m=1}^\infty \Delta b_m^d \sin \tilde{m}z$ is the difference of buoyancy profiles between the radiative equilibrium and the reference profile (i.e., the non-aggregated RCE profile).

An equilibrium of our two-column model with $\hat{\xi} = \hat{\xi}_e > 0$ exists if $b_e(\hat{x}, z; \hat{\xi}_e) = \sum_{m=1}^\infty b_{em}(\hat{x}; \hat{\xi}_e) \sin(\tilde{m}z)$ verifies the

and check whether the corresponding functions $b_e(\hat{x}, z_i; \hat{\xi}_e)$ verify the inequalities (17).
## 2.5. Method of solution for equilibrium stability

The previous section provides a method to find the equilibria of our model but does not determine their stability, i.e., whether small perturbations from the equilibrium will lead to the system returning to or drifting away from equilibrium. For $\delta b_i > 0$, the trivial equilibrium with $\xi_e = 0$, $b_e = 0$, $w_e = 0$ is obviously stable: perturbations with $b$ smaller than $\delta b_i$ at $z = z_i$ are damped.

Assessing the stability of an equilibrium with $\hat{\xi}_e > 0$ requires a more complex derivation; we summarise it here and a detailed derivation is available in the Supporting Information. Again, we study the time evolution of a small perturbation, but the motion of the boundaries between dry and convective regions during this evolution makes it impractical to use equations (7) and (8). Instead, we perform a change of horizontal coordinate with the goal of setting the boundaries between the dry and convective regions to fixed locations in the new coordinate system and we study the time evolution of a small perturbation with the equations in this coordinate system.

For a perturbation of any kind from a given equilibrium with $\hat{\xi}_e > 0$, the symmetry of the equilibrium with respect to $x = 0$ is lost and we note $\xi_-$ the x-coordinate of the left boundary of the dry region and $\xi_+$ the x-coordinate of its right boundary. We define a new horizontal coordinate $\eta$ related to $x$ as follows:

$$x = XP(\eta, \sigma, \delta), \qquad (19)$$

$$\text{with } P = \frac{1}{1 - \hat{\xi}_e^2} \left[ (1 - \delta)\eta^3 - \hat{\xi}_e\sigma\eta^2 + (\delta - \hat{\xi}_e^2)\eta + \hat{\xi}_e\sigma \right],$$

with $\delta = (\xi_+ - \xi_-)/2\xi_e$ the normalized size of the dry region and $\sigma = (\xi_+ + \xi_-)/2\xi_e$ the normalized location of the centre of the dry region. With this change of coordinate, the dry region corresponds to $|\eta| < \hat{\xi}_e$ and the moist region to $\hat{\xi}_e < |\eta| \le 1$. We introduce the functions $B(\eta, z, t)$ and $W(\eta, z, t)$ that match the functions $b$ and $w$ in the new coordinate system:

$$B(\eta, z, t) = b(x, z, t) \quad \text{and} \quad W(\eta, z, t) = w(x, z, t).$$

We have the following relationships between the partial derivatives of $B(\eta, z, t)$ and the partial derivatives of $b$: $\partial_x b = \partial_x \eta \partial_\eta B$, $\partial_x^2 b = (\partial_x \eta)^2 \partial_\eta^2 B + \partial_x^2 \eta \partial_\eta B$, $\partial_z b = \partial_z^2 B$,

and $\partial_t b = \partial_t B + \partial_t \eta \partial_\eta B$. The same relationships between partial derivatives of $W$ and $w$ apply. We can now rewrite equations (7) and (8) in terms of $B$ and $W$:

$$(\partial_t + \epsilon + \partial_t \eta \partial_\eta) \partial_z^2 W = (\partial_x \eta)^2 \partial_\eta^2 B + \partial_x^2 \eta \partial_\eta B, \quad (20)$$

$$(\partial_t + \mu + \partial_t \eta \partial_\eta) B + W = Q. \qquad (21)$$

With $\eta$ defined by (19), we have:

$$\begin{aligned} \partial_t \eta &= -\frac{1}{P'} (P_\sigma \partial_t \sigma + P_\delta \partial_t \delta), \\ \partial_x \eta &= \frac{1}{X P'}, \quad \text{and} \quad \partial_x^2 \eta = -\frac{P''}{X^2 (P')^3}, \end{aligned}$$

in which, for brevity, the primes indicate the partial derivatives with respect to $\eta$, $\partial_\sigma P = P_\sigma$, and $\partial_\delta P = P_\delta$.

For the even steady state, we have $\xi_+ = -\xi_- = \xi_e$ and the series of vertical modes composing the buoyancy equilibrium field $B_e$ is given by replacing $\hat{x}$ with $P$ in (16). We now consider a small perturbation from this steady state. We note this perturbation with a breve such that: $B = B_e + \check{B}$, $W = W_e + \check{W}$, $\sigma = \check{\sigma}$, and $\delta = 1 + \check{\delta}$. Noting that $P'(\eta, 0, 1) = 1$ and $P''(\eta, 0, 1) = 0$, equations (20)-(21) are linearised as follows:

$$(\partial_t + \epsilon) \partial_z^2 \check{W} - (P_\sigma \partial_t \check{\sigma} + P_\delta \partial_t \check{\delta}) \partial_\eta \partial_z^2 W_e = \frac{1}{X^2} \left[ \partial_\eta^2 \check{B} - 2(P'_\sigma \check{\sigma} + P'_\delta \check{\delta}) \partial_\eta^2 B_e - (P''_\sigma \check{\sigma} + P''_\delta \check{\delta}) \partial_\eta B_e \right], \quad (22)$$

$$(\partial_t + \mu) \check{B} - (P_\sigma \partial_t \check{\sigma} + P_\delta \partial_t \check{\delta}) \partial_\eta B_e + \check{W} = 0. \qquad (23)$$

These linear equations have solutions in the exponential form $\check{B} = \check{B}_0 \exp(\Gamma t)$, with $\Gamma \in \mathbb{C}$, and equivalent expressions for $\check{W}$, $\check{\sigma}$, and $\check{\delta}$. We introduce these expressions in equations (22)-(23), use the decomposition in vertical modes (10) and eliminate the vertical-speed variables, which yields one equation for each vertical mode $\check{B}_m$ and each region as follows:

$$\begin{aligned} \partial_\eta^2 \check{B}_m - \tilde{m}^2 X^2 (\Gamma + \epsilon) (\Gamma + \mu^a) \check{B}_m = & \\ [P_\sigma'' \check{\sigma} + P_\delta'' \check{\delta} - \tilde{m}^2 X^2 \Gamma (\Gamma + \epsilon + \mu^a) (P_\sigma \check{\sigma} + P_\delta \check{\delta})] \partial_\eta \mathcal{B}_{em} & \\ + 2\epsilon \mu^a \tilde{m}^2 X^2 (P'_\sigma \check{\sigma} + P'_\delta \check{\delta}) \mathcal{B}_{em}, & \end{aligned} \quad (24)$$
with $a = d$ for $|\eta| \le \hat{\xi}_e$ and $a = c$ for $\hat{\xi}_e \le |\eta| \le 1$. All terms related to the steady state have been expressed in terms of $B_{em} = B_{em} - Q_m/\mu$. ($P_\sigma\hat{\sigma} + P_\delta\delta)\partial_\eta B_{em}$ is a particular solution of this equation in both regions and the solutions to the corresponding homogeneous equation are in the form $e^{\pm\Lambda_m^a\eta}$ with $\Lambda_m^a = \tilde{m}X\sqrt{(\Gamma + \epsilon)(\Gamma + \mu^a)}$ in each region.

We now apply the continuity conditions at $\eta = \pm\hat{\xi}_e$ and $\pm 1$ to the general solutions in each region in order to obtain the expression of $\tilde{B}_m$ for a given $\Gamma$. Using the symmetries of $B_{em}$ (even and periodic), $P_\sigma$ (even) and $P_\delta$ (odd), and the specific values $P_\sigma(1) = P_\delta(1) = 0$ and $P_\sigma(\hat{\xi}_e) = P_\delta(\hat{\xi}_e) = \hat{\xi}_e$ to simplify the equations, we obtain:

$$ 
\begin{aligned}
\tilde{B}_m &= (P_\sigma\hat{\sigma} + P_\delta\delta)\partial_\eta B_{em} + \hat{\xi}_e \left[ \partial_\eta^2 B_{em} \right] \times \left( \hat{\sigma} \frac{\sinh(\eta) \sinh \Lambda_m^c (|\eta| - 1)}{D_{1m} \sinh Y_m^c} + \delta \frac{\cosh \Lambda_m^c (|\eta| - 1)}{D_{2m} \cosh Y_m^c} \right) \quad \text{for } \hat{\xi}_e \le |\eta| \le 1, \\
&= (P_\sigma\hat{\sigma} + P_\delta\delta)\partial_\eta B_{em} + \hat{\xi}_e \left[ \partial_\eta^2 B_{em} \right] \times \left( \hat{\sigma} \frac{\sinh \Lambda_m^d \eta}{D_{1m} \sinh Y_m^d} + \delta \frac{\cosh \Lambda_m^d \eta}{D_{2m} \cosh Y_m^d} \right) \quad \text{for } |\eta| \le \hat{\xi}_e,
\end{aligned}
\quad (25) $$

with $[\partial_\eta^2 B_{em}]$ the jump in $\partial_\eta^2 \tilde{B}_{em}$ at $\eta = \hat{\xi}_e$, $Y_m^c = \Lambda_m^c(\hat{\xi}_e - 1)$, $Y_m^d = \Lambda_m^d \hat{\xi}_e$, $D_{1m} = \Lambda_m^d \coth Y_m^d - \Lambda_m^c \coth Y_m^c$, and $D_{2m} = \Lambda_m^d \tanh Y_m^d - \Lambda_m^c \tanh Y_m^c$.

Finally, the convective-inhibition criterion, which stipulates that $B(\pm\hat{\xi}_e, z_i, t) = \delta B_t$, or $\tilde{B}(\pm\hat{\xi}_e, z_i, t) = 0$ at all times, gives two constraints that determine the acceptable values of $\Gamma$. These constraints can be combined into:

$$ \hat{\sigma} f_1(\Gamma) = 0, \quad (26) $$

$$ \delta f_2(\Gamma) = 0, \quad (27) $$

with:

$$ f_1(\Gamma) = \partial_\eta B_e(\hat{\xi}_e, z_i) + \sum_{m=1}^{\infty} \frac{[\partial_\eta^2 B_{em}]}{D_{1m}} \sin \tilde{m} z_i = 0, \quad (28) $$

$$ f_2(\Gamma) = \partial_\eta B_e(\hat{\xi}_e, z_i) + \sum_{m=1}^{\infty} \frac{[\partial_\eta^2 B_{em}]}{D_{2m}} \sin \tilde{m} z_i = 0. \quad (29) $$

A perturbation with a given $\Gamma$ is possible if $f_1(\Gamma) = 0$ and/or $f_2(\Gamma) = 0$. If only $f_1(\Gamma) = 0$, the perturbation must verify $\delta = 0$, i.e. the size of the dry region is not changed by this perturbation, and (24) shows the perturbation is odd in $\eta$. If only $f_2(\Gamma) = 0$, the perturbation must verify $\hat{\sigma} = 0$, i.e., the centre of the dry region is not displaced by this perturbation, and the perturbation is even.

A perturbation verifying $f_1(\Gamma) = 0$ and/or $f_2(\Gamma) = 0$ will grow if the real part of $\Gamma$ is strictly positive. If $\Gamma = 0$, the perturbation will stay constant. If the real part of $\Gamma$ is strictly negative, the perturbation will decay.

If all possible perturbations are associated with $\Gamma$ with strictly negative real parts, all perturbations will decay and the equilibrium is stable. If there is one possible perturbation for $\Gamma$ with a strictly positive real part, it is enough to destabilize the equilibrium. If there is a possible perturbation with $\Gamma = 0$, the equilibrium is neutrally stable.

### 3. Results

#### 3.1. Multiple equilibria in the control case

For our control case, we take a domain of non-dimensional horizontal size $2X^o = 1$ (1200 km), with the subscript “o”. indicating the control value. For the mechanical damping, we choose $\epsilon^o = 0.12$ (corresponding to one day) similarly to previous studies considering damped gravity waves (Kuang 2008; Daleu et al. 2015, 2016). Table 1 summarizes the control values of the parameters and constants.

In the control case, the trivial solution $\hat{\xi}_e = 0$, $b_e = 0$, $w_e = 0$ is a stable equilibrium since $\delta b_i^o > 0$. As noted in the previous section, it is a non-aggregated equilibrium, with convection everywhere in the domain.

To look for other solutions with $\hat{\xi}_e > 0$, we compute function $f$ using equation (18) and plot it in Figure 2. It exhibits two roots, and the associated equilibrium buoyancy field decreases

<table><caption>Table 1. Parameter values in the control case.</caption><thead><tr><th>z<sub>i</sub></th><th>X<sup>o</sup></th><th>&epsilon;<sup>o</sup></th><th>&mu;<sup>co</sup></th><th>&mu;<sup>do</sup></th><th>&delta;b<sub>i</sub><sup>o</sup></th><th>Q<sub>cv</sub><sup>o</sup></th><th>&alpha;<sub>s</sub><sup>o</sup></th><th>Q<sub>bl</sub><sup>o</sup></th></tr></thead><tbody><tr><td>1/6</td><td>0.5</td><td>0.12</td><td>1</td><td>0.02</td><td>10<sup>-3</sup></td><td>0.05</td><td>0.5</td><td>0.02</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Function $f = b_e(\hat{\xi}, z_i; \hat{\xi}) - \delta b_i$ in the control case. $f(\hat{\xi}) = 0$ corresponds to the equilibria of the model with $\hat{\xi} > 0$. Equilibrium solutions $\hat{\xi}_e = \hat{\xi}_u$ and $\hat{\xi}_e = \hat{\xi}_s$ are indicated, respectively, by a red and a blue dots.</figcaption></figure>

monotonously with $\hat{x}$ at $z_i$ in both cases. There are therefore two aggregated equilibria which exhibit a dry, convectively suppressed region and a convective region additionally to the trivial, non-aggregated equilibrium. One of these has $\hat{\xi}_e = \hat{\xi}_u \approx 0.035$ (a few percents of the domain) and the other has $\hat{\xi}_e = \hat{\xi}_s \approx 0.71$ (almost three quarters of the domain). In order to investigate the stability of these equilibria, we study functions $f_1$ and $f_2$ given by (28) and (29). Figure 3 shows the real and imaginary parts of functions $f_1$ and $f_2$ for $\hat{\xi}_e = \hat{\xi}_u$. Both the real and imaginary parts of $f_2$ are zero for $\Gamma \approx 0.91 > 0$, which shows that there is a growing perturbation mode with an even structure in $B$ which destabilizes the equilibrium with $\hat{\xi}_e = \hat{\xi}_u$. There is no such a mode for $\hat{\xi}_e = \hat{\xi}_s$ (not shown), so the equilibrium with $\hat{\xi}_e = \hat{\xi}_s$ is stable. For both equilibria, there is a neutrally stable, odd perturbation mode for $\Gamma = 0$ ($f_1(0) = 0$ can be computed analytically), which expresses the fact that the equilibrium is invariant by translation.

The equilibrium fields of buoyancy and momentum in both aggregated equilibria are qualitatively similar, differing mostly in terms of amplitude and $\hat{\xi}_e$. In the following section, we analyse the stable aggregated equilibrium in more details, and the remarks apply to the unstable aggregated equilibrium as well.

## 3.2. Details of the stable aggregated equilibrium

Figure 4 shows the fields of diabatic heating $Q_e$, buoyancy $b_e$, horizontal momentum $u_e$, and vertical momentum $w_e$ for the stable aggregated equilibrium with $\hat{\xi}_e = \hat{\xi}_s$. The diabatic cooling $Q_e$ (fig. 4a) follows the profile shown in Figure 1 in

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Functions $f_1(\Gamma)$ and $f_2(\Gamma)$ for $\hat{\xi}_e = \hat{\xi}_u$ in the control case. $f_1(\Gamma) = 0$ and/or $f_2(\Gamma) = 0$ correspond to the eigenmodes of the model linearised about the equilibrium. Such modes are indicated with a dot.</figcaption></figure>

the dry, convectively inhibited region ($\hat{x} < \hat{\xi}_s$) and is zero in the convective region ($\hat{x} \ge \hat{\xi}_s$). This heating forces buoyancy anomalies $b_e$ (fig. 4b) and therefore a circulation (fig. 4c,d) which limits and smoothes the buoyancy gradients. The buoyancy anomalies are negative everywhere except near the surface and the top of the boundary layer in the dry region. These positive anomalies result from adiabatic subsidence warming $w_e$ being larger than the diabatic cooling $Q_e$ there. The positive buoyancy anomaly at the top of the boundary layer ensures convective inhibition in the dry region. The horizontal buoyancy gradients $\partial_x b_e$ are weak compared to the gradients of diabatic heating $Q_e$, as expected in the free troposphere (Bretherton and Smolarkiewicz 1989). These gradients are sufficient to support a domain-wide circulation with subsidence in the dry, convectively inhibited region ($\hat{x} < \hat{\xi}_s$) and ascending motion in the convective region ($\hat{x} \ge \hat{\xi}_s$).

As a result from the buoyancy pattern with sign changes in the boundary layer, this overturning circulation has two components: one shallow with horizontal flow converging in the convective region near the surface and return flow around the top of the boundary layer and the other deep, with additional horizontal convergence in the moist region in the middle troposphere and divergence in the upper troposphere. The shallow circulation is typical of self-aggregation and its role in transporting energy up-gradient has been identified as critical for self-aggregation (Bretherton et al. 2005; Muller and Held 2012; Muller and Bony 2015; Holloway and Woolnough 2016).
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Equilibrium fields of (a) $Q_e$, (b) $b_e$, (c) $u_e$, and (d) $w_e$ for $\hat{\xi}_e = \hat{\xi}_s$ in the control case. A line contour for $b_e = \delta b_i$ is also plotted in panel b. The altitude of the boundary-layer top $z_i$ and the location of the transition $\hat{\xi}_e = \hat{\xi}_s$ between the dry region $[0 \ \hat{\xi}_s]$ and the convective region $[\hat{\xi}_s \ 1]$ are indicated by dotted gray lines.</figcaption></figure>

The equilibrium buoyancy field $b_e(x, z; \hat{\xi}_s)$ in Figure 4b can be exactly decomposed into the sum of contributions from convective, stratiform, and boundary-layer radiative heating by computing the contributions of the different terms of diabatic heating in the dry region $Q^d$ given in equation (9) to $b_e(x, z; \hat{\xi}_s)$, keeping $\hat{\xi}_s$ fixed as determined by the sum of all three contributions. Each contribution includes the direct thermodynamic response and the dynamical feedback to that response. Figure 5 shows this decomposition, with the profiles of the three contributions to $Q^d$ (fig. 5a,c,e) and the buoyancy response to each contribution (fig. 5b,d,f). Buoyancy responses to the convective and stratiform coolings follow the same vertical structure as the corresponding diabatic heating, respectively first- and second-baroclinic (fig. 5a-d), and create buoyancy responses with small horizontal gradients. The response to the boundary-layer radiative cooling is more complex (fig. 5e,f); it creates the shallow circulation which is responsible for the positive anomalies

of $b_e$ near the surface and the top of the boundary layer in the dry region and for most of the total horizontal gradients of $b_e$.

These differences in dynamical responses result from dependence of the intensity of the dynamical feedback on the vertical wavenumber $\tilde{m}$. Equation (11) shows that, for a given buoyancy anomaly, the intensity of the vertical motion is inversely proportional to $\tilde{m}^2$: higher-order vertical modes have a smaller dynamical response than low-order modes and the dynamical response $w_e$ to a given perturbation is vertically smoother than the perturbation itself. As a result, the dynamical warming response to the boundary-layer radiative cooling extends significantly upwards and downwards from the radiative-cooling maximum and becomes larger than the radiative cooling near the surface and the top of the boundary layer, which explains the positive buoyancy anomalies there. Such smoothing effect does not appear in the response to convective and stratiform heating because each of them only involves one vertical mode. Equations (11) and (12) show that the dynamical feedback acts as a horizontal diffusion
<figure><img src="image_7.png" /><figcaption>**Figure 5.** (a,c,e) Profiles of the contributions to the dry-region anomaly $Q^d$ of diabatic heating respective to the reference state and (b,d,f) the corresponding buoyancy responses for (a,b) convective, (c,d) stratiform, and (e,f) boundary-layer radiative heating.</figcaption></figure>

on a buoyancy anomaly. The decrease of the intensity of the dynamical response $w_{em}$ with increasing vertical wavenumber $\tilde{m}$ also means that this diffusive effect is smaller for higher-order vertical modes than low-order vertical modes, and therefore the response to the convective and stratiform heating ($m=1$ and $2$) are horizontally smoother, with smaller horizontal gradients, than the response to higher-order vertical modes included in the boundary-layer radiative cooling.

Both the stratiform heating and the dynamical response to the boundary-layer radiative cooling contribute to the positive buoyancy anomaly at altitude $z_i$ in the dry region, and therefore

to convective inhibition there, but the latter contributes most to the buoyancy gradients between dry and convective regions at this altitude, and therefore to the transition between convectively inhibited region and the convective region. Our model therefore highlights the role of the shallow circulation in maintaining convective inhibition in parallel with its known role in up-gradient energy transport.

The equilibrium fields for $\hat{\xi}_e = \hat{\xi}_u$ are very similar to the fields shown in Figure 4 for $\hat{\xi}_e = \hat{\xi}_s$, but with much smaller amplitude and horizontal extents (not shown).

### 3.3. Hysteresis in the sensitivity to convective inhibition

The most interesting sensitivity of our model is that to the convective inhibition threshold $\delta b_i$. Physically, decreasing $\delta b_i$ favours convective inhibition and makes convection more difficult. Figure 6a shows the function $f$ for different values of $\delta b_i$, keeping all other parameters at their control values. Since $b_e$ does not depend on $\delta b_i$, varying this parameter shifts the function $f$ by a constant. For values of $\delta b_i$ above a threshold $\delta b_i^n \approx 7.8 \cdot 10^{-3}$, there is no solution to $f(\hat{\xi}) = 0$; only the trivial equilibrium with $\hat{\xi}_e = 0$ exists. For $0 < \delta b_i < \delta b_i^n$, $f$ has two positive roots $\hat{\xi}_u$ and $\hat{\xi}_s$ as in the control case; there are three equilibria: one stable, non-aggregated equilibrium with $\hat{\xi}_e = 0$, one unstable, aggregated equilibrium with $\hat{\xi}_e = \hat{\xi}_u$, and one stable, aggregated equilibrium with $\hat{\xi}_e = \hat{\xi}_s$. For $\delta b_i < 0$, $\hat{\xi}_u$ is negative and $f$ has one single positive root $\hat{\xi}_s$; also, the trivial equilibrium is impossible because the equilibrium buoyancy field $b_e = 0$ corresponds to inhibition everywhere rather than convection everywhere. There is only one realistic equilibrium with $\hat{\xi}_e = \hat{\xi}_s$. For $\delta b_i = 0$, $\hat{\xi}_u = 0$: the trivial and unstable equilibria coincide. Any small positive buoyancy perturbation from this equilibrium would cause localised convective inhibition and diabatic cooling, and subsequently a large circulation feedback because the amplitude of steady dynamical response is proportionately larger for localised heating than for large-scale heating (Reboredo and Bellon 2022), and this response will expand the convectively inhibited region.

Figure 6b shows the different equilibrium $\hat{\xi}_e$ as a function of $\delta b_i$: zero (black line), $\hat{\xi}_u$ (red line) and $\hat{\xi}_s$ (blue line). It exhibits
a hysteresis between a critical point at $\delta b_i = 0$ and a saddle node bifurcation at $\delta b_i = \delta b_i^n \approx 7.8 \cdot 10^{-3}$.

Our model's sensitivity to convective inhibition is qualitatively similar to the sensitivity of kilometre-resolution models: the easier the convective inhibition (or the more difficult deep convection is), the more self-aggregation-prone the model. Of course, in these models, there is no simple parameter controlling convective inhibition as in our simple model, but the sensitivity to rain re-evaporation in the boundary layer (Jeevanjee and Romps 2013; Muller and Bony 2015; Holloway and Woolnough 2016), to subsidence (Bellon and Coppin 2022), and to convective entrainment (Tompkins and Semie 2017) all follow this qualitative rule, and many of these sensitivities exhibit a hysteresis. Our simple model therefore provides a simple interpretation: if inhibition is easy, the circulation response to the diabatic forcing resulting from the inhibition of convection easily maintains a

buoyancy anomaly in the lower free troposphere of the dry region that is sufficient to maintain convective inhibition there, and therefore the existence of aggregated equilibria.

### 3.4. Other sensitivities

#### 3.4.1. Sensitivity to boundary-layer radiative cooling

In kilometre-resolution models, one important mechanism for aggregation involves the low-level radiative cooling in regions without deep convection, which causes a shallow circulation and up-gradient energy transport (Bretherton et al. 2005; Muller and Held 2012; Muller and Bony 2015; Holloway and Woolnough 2016). When this cooling is removed, self-aggregation rarely occurs. Figure 5c also shows that, in our model, the buoyancy response to the boundary-layer radiative cooling contributes significantly to convective inhibition in the dry region. It therefore seems interesting to investigate the sensitivity of our simple model to this cooling.

Figure 7 shows the sensitivity to its amplitude $Q_{bl}$, keeping all other parameters at their control values. The trivial non-aggregated equilibrium with $\hat{\xi}_e = 0$ exists for all values of $Q_{bl}$ because the boundary-layer radiative cooling is applied only in the dry region, which does not exist in non-aggregated states. By contrast, the two aggregated equilibria with $\hat{\xi}_e > 0$ only exist for $Q_{bl} > Q_{bl}^n \approx 1.7 \cdot 10^{-3}$. Self-aggregation is impossible without significant low-level cooling in the dry region, because this cooling induces a shallow circulation that warms the lower free-troposphere, as shown in Figure 5c, and contributes to maintaining buoyancy above the inhibition threshold $\delta b_i$ at $z_i$ in that region.

For large $Q_{bl}$, both $\hat{\xi}_u$ and $\hat{\xi}_s$ tend towards constant asymptotes, respectively zero and around 0.88. In the limit $Q_{bl} \to \infty$, both the convective inhibition threshold $\delta b_i$, and the convective and stratiform heating $Q_{cv}$ are negligible in front of $Q_{bl}$. The limit of $\hat{\xi}_u$ is its limit for $\delta b_i = 0$ as documented in previous section. As for $\hat{\xi}_s$, its limit is non-zero because the buoyancy equilibrium response to the boundary-layer radiative cooling includes a positive buoyancy anomaly at the top of the boundary layer (as seen in Figure 5c) and permits the existence of a convectively inhibited region†.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a) Functions *f* for different values of $\delta b_i$ and (b) equilibrium values of $\hat{\xi}_e$ as a function of $\delta b_i$: zero (black), $\hat{\xi}_u$ (red) and $\hat{\xi}_s$ (blue). The control parameter value is $\delta b_i^o = 10^{-3}$ and the saddle node bifurcation between the two aggregated equilibria occurs for $\delta b_i^n \approx 7.8 \cdot 10^{-3}$.</figcaption></figure>

<footnote>†unlike in the limit $Q_{cv} \to \infty$, which is equivalent to negligible inhibition threshold and boundary-layer radiative cooling, for which there is no solution</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Equilibrium values of $\hat{\xi}_e$ as a function of $Q_{bl}$: zero (black), $\hat{\xi}_u$ (red) and $\hat{\xi}_s$ (blue). The control parameter value is $Q_{bl}^o = 0.02$ and the saddle-node bifurcation occurs for $Q_{bl}^n \approx 1.7 \cdot 10^{-3}$.</figcaption></figure>

Again, our model's sensitivity is similar to the sensitivity of kilometre-resolution models: self-aggregation is impossible for small boundary-layer radiative cooling in the dry region. And our model provides an additional physical insight: the shallow circulation in response to the boundary-layer radiative cooling contributes significantly to the maintenance of convective inhibition in the dry region through subsidence warming. Our model's sensitivity to the boundary-layer radiative cooling does differ from the sensitivity of kilometre-resolution models in one aspect: while these models have only one stable aggregated stationary state for realistic, large boundary-layer radiative cooling, our model exhibits multiple equilibria even for $Q_{bl} \to \infty$. Nevertheless, $\hat{\xi}_u$ tends towards zero in this limit and the unstable equilibrium is close to the trivial equilibrium. A small perturbation from the trivial equilibrium could move the model state into the basin of attraction of the stable aggregated equilibrium and start a drift of the model state towards the stable aggregated equilibrium. Internal small-scale variability, which is not represented in our model, could create such perturbations, in which case the trivial equilibrium would only be metastable. Including internal variability, maybe as stochastic perturbations, might make our model's behaviour more similar to kilometre-resolution models' in this respect.

### 3.4.2. Sensitivity to stratiform heating

Figure 5b shows that the buoyancy response to the stratiform heating contributes significantly to convective inhibition in the dry region; the ratio between stratiform and convective heating $\alpha_s$ is therefore a parameter of interest as well. Figure 8 shows our model's sensitivity to this parameter, keeping all other parameters at their control value. The trivial non-aggregated equilibrium with $\hat{\xi}_e = 0$ exists for all values of $\alpha_s$ for the same reasons it exists for all values of $Q_{bl}$: these parameters impact the equilibrium only if convection is inhibited somewhere in the domain which is not the case in the trivial equilibrium. The two aggregated equilibria exist only for $\alpha_s$ above a threshold $\alpha_s^n \approx 0.19$. Below this threshold, the contribution of stratiform heating to the boundary-layer-top buoyancy anomaly is insufficient to inhibit convection. For $\alpha_s^n < \alpha_s < \alpha_s^X \approx 0.60$, both aggregated equilibria exist. And for $\alpha_s = \alpha_s^X$, the stratiform heating is large enough to inhibit convection everywhere and $\hat{\xi}_s = 1$: this is an unrealistic case corresponding to a convectively-stable radiative equilibrium. For $\alpha_s > \alpha_s^X$, only the stable non-aggregated equilibrium and the unstable aggregated equilibrium exist.

This sensitivity suggests a possible influence of the vertical profile of diabatic heating associated with convective clouds on convective self-aggregation. If the profile is not sufficiently top-heavy, the lower-free-tropospheric buoyancy response to

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Equilibrium values of $\hat{\xi}_e$ as a function of $\alpha_s$: zero (black), $\hat{\xi}_u$ (red) and $\hat{\xi}_s$ (blue). The control parameter value is $\alpha_s^o = 0.5$, and the thresholds for the saddle-node bifurcation and the critical point are $\alpha_s^n \approx 0.19$ and $\alpha_s^X \approx 0.60$.</figcaption></figure>

<footnote>$\hat{\xi}_e > 0$ because the buoyancy equilibrium response to the convective and stratiform heating is negative at the top of the boundary layer.</footnote>
convective inhibition might not be warm enough to inhibit convection anywhere in the domain.

### 3.4.3. Sensitivity to domain size

With kilometre-resolution models, domain size is a constraint on self-aggregation: for small domains, self-aggregation does not occur, while for large domains it almost always does (Muller and Held 2012; Muller and Bony 2015; Bellon and Coppin 2022; Patrizio and Randall 2019). Figure 9 shows our model's sensitivity to the horizontal domain size $X$. The trivial non-aggregated equilibrium with $\xi_e = 0$ exists for all values of $X$ for the same reasons it exists for all values of $Q_{bl}$ and $\alpha_s$. The self-aggregated equilibria exist only for $X$ above a threshold $X^n \approx 0.22$ (which corresponds to a 500-km periodic domain). The steady-state circulation response acts as horizontal diffusion on buoyancy perturbations, and this diffusion is more efficient at small scales (Reboredo and Bellon 2022; Bellon and Reboredo 2022). For small $x < X^n$, this response is too efficient to maintain significant buoyancy anomalies in the dry region, and in particular to maintain anomalies larger than $\delta b_i$ at boundary-layer top, and therefore to maintain convective inhibition.

For $X \to +\infty$ (or at planetary scale), $\xi_u$ and $\xi_s$ scales in $X^{-1}$, which corresponds to dry region of finite absolute sizes $\xi_u \to 1.3 \times 10^{-2}$ and $\xi_s \to 7.8$ within an infinite convective region. In this limit, the circulation feedback decreases with $x$ away from $\xi_e$ in both regions. In the convective region, the troposphere tends to a non-aggregated state of rest for $x \to \infty$. The size of the dry region is limited to maintain a circulation feedback sufficient to support convective inhibition all the way to the centre of the dry region (at $x=0$). As for the sensitivity to boundary-layer radiative cooling (see Section 3.4.1), our model's sensitivity to domain size is similar to the sensitivity of kilometre-resolution models', with no self-aggregation for small domains. Again, the main difference is that our model's non-aggregated equilibrium exists for all domain size while self-aggregation is systematic for very large domains with most kilometre-resolution models. But, again, since $\xi_u$ is very small for large domains, our model's trivial non-aggregated equilibrium is likely to be metastable if internal variability is included.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Equilibrium values of $\xi_e$ as a function of domain size $X$: zero (black), $\xi_u$ (red) and $\xi_s$ (blue). The control parameter value is $X^0 = 0.5$, and the thresholds for the saddle node bifurcation is $X^n \approx 0.22$.</figcaption></figure>

Note that here we are limiting this study to the case with one dry region and one convective region. Since our aggregated solutions are periodic, they provide multiple solutions for a domain of given size $2X$, one with $N$ dry regions and $N$ convective regions for each integer $N$ between 1 and $X/X^n$, by replicating $N$ times the solution with one dry region and one convective region for domain size $2X/N$ along the $x$-axis. These solutions have the same stability as the corresponding solution on $2X/N$: our model provides only a lower bound on the size of the combination of one dry and one convective region, or an upper bound on the number of dry and convective regions within a large domain, not a full scale selection, whereas kilometre-resolution models tend to select $N=1$ for square domains at least up to very large scales (Patrizio and Randall 2019) or a small $N$ for elongated domains (Wing and Cronin 2016; Wing et al. 2020). If, as we expect, the trivial equilibrium is metastable in the presence of internal variability, the stable aggregated equilibrium is probably also metastable for very large domains, since additional dry regions could appear in the convective region far away from the initial dry region, where the state of the atmosphere tends towards the metastable trivial equilibrium. This suggests that some mechanisms of scale selection might come from scale interaction. Another possibility is that non-linearities at domain scale contribute to scale selection.

### 3.4.4. Sensitivity to the damping rates

Rayleigh and Newton damping rates set the spatial scale of the steady-state circulation. In Gill's seminal study (Gill 1980),
these rates are equal and the sensitivity to spatial scales and the sensitivity to the damping rate are related (Reboredo and Bellon 2022). For the Gill circulation with differing Rayleigh and Newton damping rates, the geometric mean of the two is the main control on the scale of the circulation response (Wu et al. 2001). In our case, it is a little more complex as we have one Rayleigh damping rate $\epsilon$ and two Newtonian rates $\mu^c$ and $\mu^d$. Equation (16) shows that these damping rates appear in the solutions through:

• the radiative-equilibrium buoyancy anomaly $\Delta b_m^d = Q_m^d / \mu^d$;

• the ratio between the Newtonian rates via $\nu = \sqrt[4]{\mu^c / \mu^d}$;

• $\lambda = X \sqrt{\epsilon} \sqrt[4]{\mu^c \mu^d}$.

Consequently, our model's sensitivity to the Rayleigh damping rate $\epsilon$ is identical to the sensitivity to $X^2$ documented in the previous section, with no aggregated equilibrium for small $\epsilon$ and asymptotes $\hat{\xi}_u$ and $\hat{\xi}_s \sim \epsilon^{-1/2}$ for $\epsilon \to \infty$. Physically, since the amplitude of the circulation response is proportional to $\epsilon^{-1}$, the horizontal diffusion of buoyancy anomalies by the vertical transport is very efficient for small $\epsilon$ and decreasing $\epsilon$ reduces the positive buoyancy anomalies at boundary-layer top in the dry region. For small $\epsilon$, these anomalies are insufficient to maintain convective inhibition and no aggregated equilibrium can exist. For large $\epsilon$, the buoyancy anomalies decrease fast with $x$ away from the boundary between the dry and convective regions, which limits the size of the dry region and allows the part of the domain near $x = X$ to be essentially unaffected by the circulation and close to the reference state of rest.

Our model's sensitivities to the Newtonian rates $\mu^c$ and $\mu^d$ are more complex. For the radiative damping rate $\mu_d$, there is a contribution from the sensitivity of $\Delta b_m^d$, for which increasing $\mu^d$ is equivalent to decreasing all $Q_m^d$, which in turn is equivalent to increasing $\delta b_i$; this contribution to the sensitivity of the solutions to $\mu^d$ is therefore the same as the sensitivity to $\delta b_i$ documented in Section 3.3.

There is a contribution of both the convective and radiative damping rates $\mu^c$ and $\mu^d$ through their geometrical mean $\mu = \sqrt{\mu^c \mu^d}$, which is the same as the sensitivity to $\epsilon$ or $X^2$. There is no aggregated equilibrium for small $\mu$ and asymptotes $\hat{\xi}_u$ and $\hat{\xi}_s \sim \mu^{-1/2}$ for $\mu \to \infty$. The physical interpretation of this

sensitivity is similar to the interpretation of the sensitivity to $\epsilon$. The horizontally diffusive effect of the circulation on the buoyancy anomaly scales with $\mu^{-1}$. For small $\mu$, the buoyancy anomalies are smoothed efficiently and, in the dry region, insufficient to maintain convective inhibition. For large $\mu$, the diffusion is weak and the buoyancy anomalies do not extend very far into the convective region, the dry region interacts only with a limited part of the convective region.

Finally, there is a contribution of their ratio $\mu^c / \mu^d$, which we illustrate in Figure 10. While the trivial, non-aggregated equilibrium exists for all values of $\mu^c / \mu^d$, the aggregated equilibria exist only for a finite interval of this ratio, spanning 4 orders of magnitude. This interval includes 1 which shows that self-aggregation in this model occurs even if the characteristic time scales of adjustment of convective and radiative processes are the same. Still, for self-aggregation to be possible, the convective adjustment cannot be one order of magnitude slower than the radiative adjustment, while it can be three orders of magnitude faster. Physically, if the dry radiative adjustment is significantly faster than the convective adjustment ($\mu^d > \mu^c$), the fast adjustment of the buoyancy in the dry region towards the radiative-equilibrium profile $\Delta b_m^d$ requires intense subsidence to maintain a positive buoyancy anomaly at boundary-layer top and inhibit convection; but the ascent in the convective region is small because it only compensates the slow convective adjustment. To ensure mass conservation, the size of the dry region has to be small but this size compensation has a limit beyond which

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Equilibrium values of $\hat{\xi}_e$ as a function of the ratio of the Newtonian damping rates $\mu^c / \mu^d$: zero (black), $\hat{\xi}_u$ (red) and $\hat{\xi}_s$ (blue). The control parameter value is $\mu^{co} / \mu^{do} = 50$ and the saddle node bifurcations are for $\mu^c / \mu^d \approx 0.19$ and $2.3 \cdot 10^3$.</figcaption></figure>
aggregated equilibria cannot be sustained. On the contrary, if the convective adjustment is overwhelmingly faster than the radiative adjustment ($\mu^c \gg \mu^d$), the buoyancy anomalies are negligible in the convective region and the buoyancy gradients are too small to sustain a significant overturning circulation and maintain convective inhibition in the dry region, and aggregated equilibria cannot be sustained either.

## 4. Conclusion

Our model is conceptually simple: it is based on linear primitive equations and combining one of the simplest representation of dynamics, a damped gravity wave, with linear representations of moist and dry thermodynamics, together with a simple criterion for convective inhibition to determine where each thermodynamic regime applies. The model exhibits multiple equilibria, with one stable non-aggregated RCE and one stable aggregated RCE. It suggests that a self-aggregated RCE can be considered as a two-regime damped gravity wave, with one convective phase and one convectively inhibited phase. Each physical process (convection, radiation, and dynamics) by itself is dissipative. It is the difference between the adjustment thermodynamic profiles in the convective phase (non-aggregated RCE profile) and in the dry phase (radiative equilibrium profile) which permits the development of a steady-state circulation. In turn, subsidence warming maintains convective inhibition in the dry phase and cooling by ascent reinforces convection in the convective phase. Our model results are robust to the choice of parameters, with multiple equilibria for a significant interval of each parameter.

Despite its extreme physical simplicity, our model is able to reproduce and explain many of the sensitivities of convective self-aggregation previously described with kilometre-resolution models. It exhibits a non-linear sensitivity with a hysteresis qualitatively similar to kilometre-resolution models. The sensitivity to the parameter $\delta b_i$ modulating the criterion for convective inhibition exhibit a parameter interval with only one non-aggregated RCE and a parameter interval with only one aggregated RCE on opposite sides of the hysteresis. The general sensitivity to the ease of convection (or the stringency

of convective inhibition) is similar to the sensitivity of kilometre-resolution models: the easier the convection, the less favourable for self-aggregation.

Other sensitivities, such as the sensitivity to the boundary-layer radiative cooling in the dry region, or to the domain size, present some similarity to the sensitivity of kilometre-resolution models as well: self-aggregation does not occur for small boundary-layer radiative cooling or for small domains: only the non-aggregated equilibrium exists. It is also reasonable to presume that if some representation of internal variability was included in our model, its behaviour would be even more similar to that of kilometre-resolution models: for large domains or for large boundary-layer radiative cooling in the dry region, the non-aggregated equilibrium would be metastable rather than stable as in our results and self-aggregation would be systematic as in kilometre-resolution-model simulations.

Compared to previous models of self-aggregation, our model emphasizes the interplay between dynamics and convection (and its inhibition) rather than the intricate thermodynamic feedbacks as in column models of the instability of the non-aggregated RCE. For example, our model includes some of the radiative features identified as crucial in these studies (e.g., *Emanuel et al. 2014*) but as given heating profiles rather than interactively determined. Our model also bears some resemblance with the coarsening models (*Craig and Mack 2013; Biagioli and Tompkins 2023*), since vertical transport amounts to a horizontal diffusion in our model, so that it is mathematically equivalent to a reaction-diffusion system.

Our model, as others (e.g., *Yang 2021*), uses a simple, abrupt convective inhibition (or triggering) criterion and this is its only non-linearity. It is sufficient to grant our model a complex non-linear behaviour. Such a crucial role of triggering has been documented in more complex models (e.g., *Bellon and Sobel 2010*). Conceptually, this attributes a large role to activation control of convection (*Mapes 1997*) in the occurrence of self-aggregation.

There are significant limitations to our model in its current version. First, our method of solution can only compute steady states and cannot describe the transients, and therefore the
initiation of self-aggregation. We intend to develop a time-
evolving version of the model to investigate these transients in
the future. Second, our model includes only one thermodynamic
variable, buoyancy; while this is an advantage in terms of
presenting the simplest model possible, it excludes all possible
feedback differentiating between temperature effects and moisture
effects. Including a second thermodynamic variable in the model
is another promising extension of this work.

Data Availability Statement

Data sharing not applicable.

Funding Statement

This work is supported by the French national program LEFE
(Les Enveloppes Fluides et l'Environnement). G. Faye and P.
Noble acknowledge support from ANR project HEAD under grant
agreement ANR-24-CE40-3260.

Conflict of Interest Statement

None

Acknowledgements

The authors thank D. Saint-Martin for his methodological
insights and J.-F. Coulombel for facilitating our interdisciplinary
collaboration.

References

Arakawa A, Schubert WH. 1974. Interaction of a cumulus cloud ensemble with the large-scale environment, Part I. *Journal of the Atmospheric Sciences* **31**: 674–701.

Bellon C, Le Treut H, Ghil M. 2003. Large-scale and evaporation-wind feedbacks in a box model of the tropical climate. *Geophysical Research Letters* **30**(22): 2145, doi:10.1029/2003GL017895.

Bellon G, Bony S. 2020. Tropical and Subtropical Cloud Systems. *Clouds and Climate*: 251–278doi:10.1017/9781107447738.009.

Bellon G, Coppin D. 2022. Sensitivity of Convective Self-Aggregation to Subsidence. *Journal of Advances in Modeling Earth Systems* **14**(12): e2021MS002830, doi:10.1029/2021MS002830. eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2021MS002830.

Bellon G, Reboredo B. 2022. Scale Sensitivity of the Gill Circulation. Part II: Off-Equatorial Case. *Journal of the Atmospheric Sciences* **79**(1): 19–30, doi:10.1175/JAS-D-21-0068.1.

Bellon G, Sobel AH. 2008. Poleward-propagating intraseasonal monsoon disturbances in an intermediate-complexity axisymmetric model. *Journal of the Atmospheric Sciences* **65**: 470–489.

Bellon G, Sobel AH. 2010. Multiple Equilibria of the Hadley Circulation in an Intermediate-Complexity Axisymmetric Model. *Journal of Climate* **23**: 1760–1778.

Betts AK. 1986. A new convective adjustment scheme. Part I: Observational and theoretical basis. *Quarterly Journal of the Royal Meteorological Society* **112**(473): 677–691, doi:10.1002/qj.49711247307.

Betts AK, Miller MJ. 1986. A new convective adjustment scheme. Part II: Single column tests using GATE wave, BOMEX, ATEX and arctic air-mass data sets. *Quarterly Journal of the Royal Meteorological Society* **112**(473): 693–709, doi:10.1002/qj.49711247308.

Beucler T, Cronin T, Emanuel K. 2018. A Linear Response Framework for Radiative-Convective Instability. *Journal of Advances in Modeling Earth Systems* **10**(8): 1924–1951, doi:10.1029/2018MS001280. eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2018MS001280.

Beucler T, Cronin TW. 2016. Moisture-radiative cooling instability. *Journal of Advances in Modeling Earth Systems* **8**(4): 1620–1640, doi:10.1002/ 2016MS000763.

Biagioli G, Tompkins AM. 2023. A Dimensionless Parameter for Predicting Convective Self-Aggregation Onset in a Stochastic Reaction-Diffusion Model of Tropical Radiative-Convective Equilibrium. *Journal of Advances in Modeling Earth Systems* **15**(5): e2022MS003231, doi:10.1029/2022MS003231. eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2022MS003231.

Bretherton CS, Blossey PN, Khairoutdinov M. 2005. An energy-balance analysis of deep convective self-aggregation above uniform SST. *Journal of the Atmospheric Sciences* **62**(12): 4273–4292, doi:10.1175/JAS3614.1.

Bretherton CS, Smolarkiewicz PK. 1989. Gravity Waves, Compensating Subsidence and Detrainment around Cumulus Clouds. *Journal of Atmospheric Sciences* **46**(6): 740 – 759, doi:10.1175/1520-0469(1989)046(0740:GWCSAD)2.0.CO;2.

Bretherton CS, Sobel AH. 2002. A simple model of a convectively coupled walker circulation using the weak temperature gradient approximation. *Journal of Climate* **15**(20): 2907–2920, doi:10.1175/1520-0442(2002)015(2907:ASMOAC)2.0.CO;2.

Coppin D, Bony S. 2015. Physical mechanisms controlling the initiation of convective self-aggregation in a General Circulation Model. *Journal of Advances in Modeling Earth Systems* **7**(4): 2060–2078, doi:10.1002/ 2015MS000571.

Craig GC, Mack JM. 2013. A coarsening model for self-organization of tropical convection. *Journal of Geophysical Research: Atmospheres* **118**(16): 8761–8769, doi:10.1002/jgrd.50674. eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/jgrd.50674.

Daleu CL, Plant RS, Woolnough SJ, Sessions S, Herman MJ, Sobel A, Wang S, Kim D, Cheng A, Bellon G, others. 2015. Intercomparison of methods
of coupling between convection and large-scale circulation: 1. Comparison over uniform surface conditions. *Journal of Advances in Modeling Earth Systems* 7(4): 1576-1601.

Daleu CL, Plant RS, Woolnough SJ, Sessions S, Herman MJ, Sobel A, Wang S, Kim D, Cheng A, Bellon G, others. 2016. Intercomparison of methods of coupling between convection and large-scale circulation: 2. Comparison over nonuniform surface conditions. *Journal of advances in modeling earth systems* 8(1): 387-405.

Emanuel K, Wing AA, Vincent EM. 2014. Radiative-convective instability. *Journal of Advances in Modeling Earth Systems* 6(1): 75-90, doi:10.1002/2013MS000270.

Emanuel KA, David Neelin J, Bretherton CS. 1994. On large-scale circulations in convecting atmospheres. *Quarterly Journal of the Royal Meteorological Society* 120(519): 1111-1143, doi:10.1002/qj.49712051902.

Gill AE. 1980. Some simple solutions for heat-induced tropical circulation. *Quarterly Journal of the Royal Meteorological Society* 106(449): 447-462, doi:10.1002/qj.49710644905.

Haerter JO. 2019. Convective Self-Aggregation As a Cold Pool-Driven Critical Phenomenon. *Geophysical Research Letters* 46(7): 4017-4028, doi:10.1029/2018GL081817. ArXiv: 2002.03474 Publisher: Wiley Online Library.

Hagos S, Zhang C, Tao WK, Lang S, Takayabu YN, Shige S, Katsumata M, Olson B, L'Ecuyer T. 2010. Estimates of Tropical Diabatic Heating Profiles: Commonalities and Uncertainties. *Journal of Climate* 23(3): 542-558, doi:10.1175/2009JCLI3025.1. Publisher: American Meteorological Society Section: Journal of Climate.

Held IM, Hemler RS, Ramaswamy V. 1993. Radiative-convective equilibrium with explicit two-dimensional moist convection. *Journal of the Atmospheric Sciences* 50(23): 3909-3927, doi:10.1175/1520-0469(1993)050:3909:RCEWET)2.0.CO;2.

Held IM, Zhao M, Wyman B. 2007. Dynamic Radiative-Convective Equilibria Using GCM Column Physics. *Journal of the Atmospheric Sciences* 64(1): 228-238, doi:10.1175/JAS3825.11. Publisher: American Meteorological Society Section: Journal of the Atmospheric Sciences.

Holloway CE, Woolnough SJ. 2016. The sensitivity of convective aggregation to diabatic processes in idealized radiative-convective equilibrium simulations. *Journal of Advances in Modeling Earth Systems* 8(1): 166-195, doi:10.1002/2015MS000511.

Jeevanjee N, Romps DM. 2013. Convective self-aggregation, cold pools, and domain size. *Geophysical Research Letters* 40(5): 994-998, doi:10.1002/grl.50204.

Kelly MA, Randall DA. 2001. A two-box model of a zonal atmospheric circulation in the tropics. *Journal of Climate* 14(19): 3944-3964, doi: 10.1175/1520-0442(2001)014(3944:ATBMOA)2.0.CO;2.

Kuang Z. 2008. Modeling the interaction between cumulus convection and linear gravity waves using a limited-domain cloud system-resolving model. *Journal of the Atmospheric Sciences* 65(2): 576-591, doi:10.1175/2007JAS2399.1.

Kuang Z. 2018. Linear Stability of Moist Convecting Atmospheres. Part I: From Linear Response Functions to a Simple Model and Applications to Convectively Coupled Waves. *Journal of the Atmospheric Sciences* 75(9): 2889-2907, doi:10.1175/JAS-D-18-0092.1. Publisher: American Meteorological Society Section: Journal of the Atmospheric Sciences.

Mapes B. 1997. Equilibrium vs. Activation Control of Large-Scale Variations of Tropical Deep Convection. *The Physics and Parameterization of Moist Atmospheric Convection*: 321-358doi:10.1007/978-94-015-8828-7_13.

Moorthi S, Suarez MJ. 1992. Relaxed Arakawa-Schubert: a parameterization of moist convection for general circulation models. *Monthly Weather Review* 120(6): 978-1002, doi:10.1175/1520-0493(1992)120:0978:RASAPO)2.0.CO;2.

Muller C, Bony S. 2015. What favors convective aggregation and why? *Geophysical Research Letters* 42(13): 5626-5634, doi:10.1002/2015GL064260.

Muller C, Yang D, Craig G, Cronin T, Fildier B, Haerter JO, Hohenegger C, Mapes B, Randall D, Shamekh S, Sherwood SC. 2022. Spontaneous Aggregation of Convective Storms. *Annual Review of Fluid Mechanics* 54(Volume 54, 2022): 133-157, doi:https://doi.org/10.1146/annurev-fluid-022421-011319. Publisher: Annual Reviews Type: Journal Article.

Muller CJ, Held IM. 2012. Detailed investigation of the self-aggregation of convection in cloud-resolving simulations. *Journal of the Atmospheric Sciences* 69(8): 2551-2565, doi:10.1175/JAS-D-11-0257.1.

Nilsson J, Emanuel KA. 1999. Equilibrium atmospheres of a two-column radiative-convective model. *Quarterly Journal of the Royal Meteorological Society* 125(558): 2239-2264, doi:10.1256/smsqj.55813.

Patrizio CR, Randall DA. 2019. Sensitivity of Convective Self-Aggregation to Domain Size. *Journal of Advances in Modeling Earth Systems* 11(7): 1995-2019, doi:https://doi.org/10.1029/2019MS001672.

Pierrehumbert RT. 1995. Thermostats, radiator fins, and the local runaway greenhouse. *Journal of the Atmospheric Sciences* 52(10): 1784-1806, doi: 10.1175/1520-0469(1995)052(1784:TRFATL)2.0.CO;2.

Raymond DJ, Zeng X. 2000. Instability and large-scale circulations in a two-column model of the tropical troposphere. *Quarterly Journal of the Royal Meteorological Society* 126(570): 3117-3135, doi:10.1256/smsqj.57006.

Reboredo B, Bellon G. 2022. Scale sensitivity of the Gill circulation, Part I: equatorial case. *Journal of the Atmospheric Sciences* 79: 1-17.

Sessions SL, Sugaya S, Raymond DJ, Adam H. 2010. Multiple equilibria in a cloud resolving model. *Journal of Geophysical Research - Atmospheres*: 1-51.

Sobel AH, Bellon G, Bacmeister JT. 2007. Multiple equilibria in a Single-Column Model of the tropical atmosphere. *Geophysical Research Letters* 34: L22804, doi:10.1029/2007GL031320.

Tao WK, Smith EA, Adler RF, Haddad ZS, Hou AY, Iguchi T, Kakar R, Krishnamurti TN, Kummerow CD, Lang S, Meneghini R, Nakamura K, Nakazawa T, Okamoto K, Olson WS, Satoh S, Shige S, Simpson J,
Takayabu Y, Tripoli GJ, Yang S. 2006. Retrieval of Latent Heating from TRMM Measurements. *Bulletin of the American Meteorological Society* **87**(11): 1555-1572, doi:10.1175/BAMS-87-11-1555. Publisher: American Meteorological Society Section: Bulletin of the American Meteorological Society.

Tompkins AM, Semie AG. 2017. Organization of tropical convection in low vertical wind shears: Role of updraft entrainment. *Journal of Advances in Modeling Earth Systems* **9**(2): 1046-1068, doi:10.1002/2016MS000802.

Windmiller JM, Craig GC. 2019. Universality in the Spatial Evolution of Self-Aggregation of Tropical Convection. *Journal of the Atmospheric Sciences* **76**(6): 1677-1696, doi:10.1175/JAS-D-18-0129.1. Publisher: American Meteorological Society Section: Journal of the Atmospheric Sciences.

Wing AA, Cronin TW. 2016. Self-aggregation of convection in long channel geometry. *Quarterly Journal of the Royal Meteorological Society* **142**(694): 1-15, doi:10.1002/qj.2628.

Wing AA, Emanuel K, Holloway CE, Muller C. 2017. Convective Self-Aggregation in Numerical Simulations: A Review. In: *Shallow Clouds, Water Vapor, Circulation, and Climate Sensitivity*, Pincus R, Winker D, Bony S, Stevens B (eds). Springer International Publishing: Cham, pp. 1-25, doi:10.1007/978-3-319-77273-8_1.

Wing AA, Reed KA, Satoh M, Stevens B, Bony S, Ohno T. 2018. Radiative-convective equilibrium model intercomparison project. *Geoscientific Model Development* **11**(2): 793-813, doi:10.5194/gmd-11-793-2018. Publisher: Copernicus GmbH.

Wing AA, Stauffer CL, Becker T, Reed KA, Ahn MS, Arnold NP, Bony S, Branson M, Bryan GH, Chaboureau JP, De Roode SR, Gayatri K, Hohenegger C, Hu IK, Jansson F, Jones TR, Khairoutdinov M, Kim D, Martin ZK, Matsugishi S, Medeiros B, Miura H, Moon Y, Müller SK, Ohno T, Popp M, Prabhakaran T, Randall D, Rios-Berrios R, Rochetin N, Roehrig R, Romps DM, Ruppert Jr JH, Satoh M, Silvers LG, Singh MS, Stevens B, Tomassini L, van Heerwaarden CC, Wang S, Zhao M. 2020. Clouds and Convective Self-Aggregation in a Multimodel Ensemble of Radiative-Convective Equilibrium Simulations. *Journal of Advances in Modeling Earth Systems* **12**(9): e2020MS002138, doi:10.1029/2020MS002138. eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1029/2020MS002138.

Wu Z, Sarachik ES, Battisti DS. 2001. Thermally driven circulation under Rayleigh friction and Newtonian cooling: analytic solutions. *Journal of the Atmospheric Sciences* **58**: 724-741.

Yang D. 2018. Boundary Layer Diabatic Processes, the Virtual Effect, and Convective Self-Aggregation. *Journal of Advances in Modeling Earth Systems* **10**(9): 2163-2176, doi:10.1029/2017MS001261.

Yang D. 2021. A Shallow-Water Model for Convective Self-Aggregation. *Journal of the Atmospheric Sciences* **78**(2): 571-582, doi:10.1175/JAS-D-20-0031.1. Publisher: American Meteorological Society Section: Journal of the Atmospheric Sciences.

Yano JI, Plant RS. 2012. Convective quasi-equilibrium. *Reviews of Geophysics* **50**(4), doi:10.1029/2011RG000378.