How weakened cold pools open for convective self-aggregation

In radiative-convective equilibrium (RCE) simulations, convective self-aggregation (CSA) is the spontaneous organization into segregated cloudy and cloud-free regions. Evidence exists for how CSA is stabilized, but how it arises favorably on large domains is not settled. Using large-eddy simulations (LES), we link the spatial organization emerging from the interaction of cold pools (CPs) to CSA. We systematically weaken simulated rain evaporation to reduce maximal CP radii, $R_{\text{max}}$, and find reducing $R_{\text{max}}$ causes CSA to occur earlier. We further identify a typical rain cell generation time and a minimum radius, $R_{\text{min}}$, around a given rain cell, within which the formation of subsequent rain cells is suppressed. Incorporating $R_{\text{min}}$ and $R_{\text{max}}$, we propose a toy model that captures how CSA arises earlier on large domains: when two CPs of radii $r_{i,j}\in[R_{\text{min}},R_{\text{max}}]$ collide, they form a new convective event. These findings imply that interactions between CPs may explain the initial stages of CSA.

unstable atmosphere, establishing connections between the loci where new convective cells emerge and loci at which the previous cells dissipated. In particular, new cells were suggested to be spawned by the CP gust front alone or by collisions between mobile gust fronts [18][19][20][21]. Inspired by the notion of CP interactions, conceptual work has formulated CPs as cellular automata [22][23][24][25], and CP representations have been incorporated into large-scale models [26].
Studies on convective self-aggregation often argue that sufficiently large domain size is required for the phenomenon to emerge [7,13], especially when the horizontal resolution is coarse, e.g., at least 500 km for resolutions finer than 2 km [10]. To examine this claim more closely, for deliberately small domain sizes and fine horizontal resolution, we here show that CSA sets in earlier when CPs are weakened through reductions in rain evaporation -that is, when the CP maximal radius, termed R max , is reduced. We track the CP gust fronts to motivate that loci of gust front collisions are preferable for subsequent convective rain cells -hence, that CP interactions are essential in organizing the precipitation field. Dependent on rain evaporation, we further detect a minimal distance R min , within which subsequent rain cells are unlikely to form, as well as a typical generation time. Using these findings, we build, simulate, and analyze a simple mathematical model, which helps understand cloud-free regions' formation. We explore the phase diagram of this model and find that, when initialized at a high density of rain cells, the population of subsequent rain cells will either remain high or transition to a low state with only a subregion covered by rain cells.
The change from high to low density occurs later for large R max or small domain sizes L. Large-eddy simulations. We carried out (96 km) 2 simulations of the convective atmosphere for up to five simulation days using the University of California, Los Angeles (UCLA) Large Eddy Simulator with sub-grid scale turbulence parametrized after Smagorinsky [27]. We combine with a delta four-stream radiation [28] and a two-moment cloud microphysics scheme [29]. Rain evaporation is accounted for by Seifert and Beheng (2006). All five simulations are identical, except that rain evaporation is varied within fractions {0, .1, .2, .6, 1.0} of its default. In the following, the corresponding simulations are labeled "Evap=0", "Evap=0.1", etc. We reduce computed rain evaporation by the corresponding factor -hence, less hydrometeor mass is converted to water vapor. Surface temperatures are set constant to 300 K, and insolation is fixed using a constant equatorial zenith angle of 50 • (constant 655 W m −2 ) [13].
Surface heat fluxes are computed interactively and depend on the vertical temperature and humidity gradients and horizontal wind speed, which is approximated using the Monin-Obukhov similarity theory. Temperature and humidity are initialized using a prior approximately 3-day spin-up using 400 m horizontal resolution (see Fig. S1). The horizontal model grid is regular, and periodic boundary conditions are applied in both lateral dimensions. Vertically, the model resolution varies from 100 m below 1 km, stretching to 200 m near 6 km and finally 400 m in the upper layers. There are 75 vertical levels. The Coriolis force and the mean wind were set to zero with weak random initial perturbations added as noise to break complete spatial symmetry. At each output time step of 10 min, instantaneous surface precipitation intensity, specific humidity, temperature, liquid water mixing ratio, outgoing long-wave radiation, and 3D velocities are output at various model levels. To explore resolution effects, we supplemented the simulations described above using Evap=1 by otherwise similar simulations at 1km, 2km, and 4km horizontal resolution, each using 480 × 480 horizontal grid boxes.
Tracking of cold pools. Cold pool (CP) gust fronts are tracked using tracer particles [24,31]. In any given time step, a sufficient number of tracers, n tr , are placed at the edges (8-neighborhood) of an existing surface precipitation patch, which is a spatially contiguous area of rainy pixels (intensity I > I 0 , with I 0 ≡ .5 mm h −1 ). The first tracers are placed when the patch is first detected, and this time defines the time origin for each CP (t = 0 in Fig. 2B). If the number of tracers placed is less than a maximal number n tr,max ≡ 100, i.e., n tr < n tr,max , then further tracers are placed at the patch edges during the subsequent timesteps until n tr,max is reached. If the rain patch disappears before n tr,max is reached, fewer tracers are used for that CP. This procedure was found sufficient in reliably tracking the gust fronts. Using a simple Euler method, the tracer particles are advected along with the radial velocity component in the lowest model level (z = 50 m) and reliably settle in the gust fronts surrounding each CP. As this tracking method is implemented to run "offline," it uses only the recorded discrete 10 min output time steps. In the initial stages of tracking each CP, this sometimes leads to particles having to "catch up" with the gust front (see Fig. 2A, for examples of such CPs). However, it is found that tracers consistently settle in the gust front after few time steps due to the strong radial velocity gradient. To transparently compare CP radii between the different simulations, which differ in time of rainfall onset (Fig. 2B), we first determine the time of onset for each simulation, that is, the time when the first surface rain patch with I > I 0 occurred. We then track all CP gust fronts during 1100 min. This interval was found sufficient to yield significant statistics on the spreading of each CP, but short enough so that not many CP collisions were encountered. Conversely, to study collision effects (Fig. 3), we used a late-stage (∼ 4 days after initialization) of the realistic simulation (Evap=1). As CPs are space-filling, any new CP inevitably collides with recent CPs in its surroundings.
Mathematical model. We initialize N 1 = L 2 /(10R 2 min ) points -the subscript indicates "generation one" -at random locations selected uniformly on a 2D domain of size L × L with double-periodic boundary conditions. All points grow into circles representing CPs with equal and constant radial speed. The expanding circles have centers at c i = [x i , y i ] and increasing radii r i . At the collision point, [x, y], of two circles belonging to the same generation g and having radii r i and r j , a new expanding circle belonging to the subsequent generation g + 1 emerges only if R min < r i , r j < R max . Hence, for two CPs to form a new circle, their centers' separation distance needs to be greater than 2R min and less than 2R max . Colliding circles continue expanding, and new circles start growing instantaneously; hence, all circles, at all times, expand with equal and constant speed.
We find the collision point by solving ( where dr is the distance from the two circles' rims to the collision point. Only collisions that fall onto the straight line between the two circle centers are allowed since that is the collision point with the highest momentum, yielding These three quadratic equations have three unknowns (x, y, dr), and two solutions. One of these two solutions can be ruled out because only the positive real solution is relevant to this model. Similar to the analytical approach in [24], we do not run the model strictly chronologically. To reduce simulation time, we take advantage of the fact that circles belonging to different generations cannot interact since they grow with equal and constant speed. Thereby, we calculate all collision points for each generation before proceeding to the next generation. Generally, the last circle in any generation will initiate later than the next generation's first circle.
Therefore, when moving to the next generation, we return to the time when the first circle of that generation was seeded.
The list of potential collision points is calculated for each generation and sorted incrementally by dr. We update the system by inserting circles at the collision points if R min < r i , r j < R max , and if the subsequent generation does not occupy the collision point. This process is simulated until circle generation number 500 is reached. and overall higher near-surface temperature (Fig. S2C), along with a systematically earlier onset of persistent dry patches, e.g., near day 2 for Evap=0.2 (Fig. 1C). This comparison underlines findings from Jeevanjee and Romps [6] and Muller and Bony [7], who reported that CPs hamper self-aggregation. The five experiments highlight that reducing rain evaporation weakens subsidence drying in the center of CPs (compare dark spots in Fig. 1A-B vs. C-D) and visibly reduces CP radii. We also note that intermediate values of evaporation appear to allow for a bandlike aggregation state, where rain cells form a quasi-one-dimensional chain around one of the horizontal dimensions ( Fig. 1C on day 4). When rain evaporation is entirely removed (Fig. 1E), any organizing effect through CPs is absent: one is left with a coarsening process akin to reaction-diffusion dynamics, [32] small impurities gradually merging into larger structures. The dominant feedback can be ascribed to the net warming in the troposphere induced by deep convective clouds, thus promoting further convective activity in existing precipitation cells' surroundings. The average CP radii as the CPs evolve after their emergence (lines) and the 90th radius percentile (dots). Note that CPs initially grow quickly but monotonically slow and that maximal CP radii increase with rain evaporation rate.
Measuring cold pool radius. Using a rain cell [33] and CP [24,31] tracking method, we seed tracer particles at the boundary of surface rain patches (see cartoon in Fig. 2A). We advect these tracers using the radial velocity field, forcing them to gather in pronounced convergence areas caused by the CP gust fronts (Details: Methods).
Superimposing the resulting pattern of tracers onto the near-surface vertical velocity field ( Fig. 2A) shows that the tracers indeed gather along the edge of each CP (subsident or featureless vertical wind field). It is also visually apparent that radii systematically increase with the evaporation rate. The radius increase is confirmed by plotting the time evolution of the average CP radii in each simulation (Fig. 2B). For Evap=0, radii simply show the radii of the corresponding surface rain cells, as, without CPs, there is no pronounced wind field to advect the tracers. For Evap=1, CPs typically expand to ≈ 11 km a few hours after initiation, a value comparable to previous simulation results found on various domain sizes [3,16] and observational findings [34][35][36].
New convective events are initiated in the vicinity of cold pool collisions. What is then the specific role of CPs in maintaining domain-wide convection? To explore this, first consider locations of rainfall at a particular time step of Evap=1 (Fig. 3A), the associated cloud-base vertical velocity (Fig. 3B), and specific humidity (Fig. 3C).
Updrafts form shortly before the onset of rainfall, whereas specific humidity becomes elevated earlier -in line with the analysis of RCE simulations, where a considerable buoyancy build-up before any subsequent convective event was reported. [21] Second, we determine gust front loci using CP tracer particles [24,31] (Fig. 3D). Many of these tracers are located at the intersection between two CPs. It is visually apparent that such loci coincide with enhanced humidity (Fig. 3D). To quantify that tracers lie in regions of pronounced updrafts, we verified that the q v (50 m)-histogram for all loci with tracers is shifted to markedly positive values (Fig. 3E). The histogram of peak humidity during rain event buildup (peak highlighted in Fig. 3C) shows a shift towards elevated values. Collecting, as a comparison, the specific humidity at CP gust fronts (Fig. 3E), it is found that this histogram is similarly shifted to moister values. In summary, loci of CP collisions do provide the positive humidity anomalies typical of subsequent convective events. New deep convective events are initiated at a specific distance away from earlier events. To quantify a possible suppression effect caused by a present rain cell's CP on subsequent cells forming within the surroundings, we examine whether rain events after the initial rain onset are spaced randomly or not. A non-random spacing would imply either suppression (larger distance) or activation (smaller distance), whereas a random spacing would speak against a direct spatial influence on subsequent rain cell formation. We thus identify all rain events within the first 12h after rain onset [37] allowing us to compare non-aggregating simulations with aggregating simulations (day 1 in Fig. 1). We measure each rain cell's distance to it's nearest rain event occurring within a time window ∆t.
For the control simulation (Fig. 4A), we find an inhibitory effect causing the nearest neighbor distance to be > 5 km for up to 8h. In contrast, it would linear decay on a log-log plot if events were distributed randomly in space (compare orange and blue curves in Fig. 4). We refer to this distance as R min ≈ 6 km and explain it by CPs being too negatively buoyant to initialize new convective cells within this distance [38,39]. When increasing the time window of included events from 8-12 hours, we find that this suppression effect disappears; that is, the distribution function matches the random one. On this time scale, the CPs associated with two rain events have time to grow larger than R min , collide, and trigger the formation of a new, closer, rain event belonging to the subsequent generation. We, therefore, interpret this time scale as the generation time of one CP. We find this time scale unchanged across simulations, while the spatial scale increases for coarser horizontal resolutions (Fig. 4).
A simple mathematical model captures the onset of self-aggregation. To understand the role of CP collisions, we introduce a model consisting of growing and colliding circles. The circles represent the gust fronts of CPs and they are initialized uniformly at random inside the model domain. The circle radius, r, initially set to zero, increases linearly until r = R max , hence the expansion velocity, v(r), is a step function v(r) = v 0 θ(R max − r) where v 0 is a constant and θ(R max − r) = 1 for r < R max , otherwise θ = 0. After a given circle reaches R max , it has no further effect and is removed. When two circles meet, both having their radii lie between R min and R max , they instantly produce a new circle at the first point of intersection. The reasoning is that in RCE, most new rain cells result from thermodynamic pre-conditioning near the gust front collision lines (Fig. 3 and Fuglestvedt and Haerter, 2020), and the delay between the collision time and the initiation of the resultant rain cell is so large (typically several hours) that direct forced lifting can be ruled out. In line with the findings in Fig. 4, CPs with r < R min are considered too negatively buoyant to initialize new CPs [38,39]. However, circles may collide with multiple other circles until they reach R max beyond which they have no role.
Mathematically, the collision dynamics allow us to categorize circles into generations: the initially seeded circles  constitute generation one (panels 1-2 in Fig. 5A,i where all circles are smaller than R min ). Each collision is between circles belonging to the same generation g, and the resultant new circles belong to the subsequent generation g + 1.
R min > 0 prevents singularities, that is, infinitely rapid replication. Given R min = 0, an initial random generation-one population N 1 of circles would yield N 2 = 2 N 1 (Details: Supplement), and subsequent growth of N g vs. g would be nearly exponential. For R min > 0 circles must grow beyond R min and only then replicate to yield generation-two events (note the different color shadings in Fig. 5A,i, third panel). In Fig. 5A,ii new circles of various generations are initiated throughout the domain with no obvious patterning. However, during Fig. 5A,iii a separation into a circlefilled (convecting) and a circle-free (non-convecting) sub-region occurs. Note the visual similarity with the numerical experiment in Fig. 1C-D. The number of circles N g in all simulations eventually drops from an initial high-N state to one with low N (Fig. 5B). The histogram of N , which is bimodal, confirms the notion of two distinct meta-stable states (Fig. 5C).
Given this finding, we explore how the two states depend on the independent model parameters L/R min and L/R max .
We find that the low-N state scales as N L = L/(2R min ), whereas the high-N state scales as N H = L 2 /(10R 2 min ), both independent of R max (Fig. 5D). The transition point occurs at N T ≈ 1.5L/R min The linear scaling N L ∼ L is commensurate with band-like, one-dimensional, structures (compare Fig. 5A,iii and Fig. 1C-D on day 2-4), whereas N H ∼ L 2 is in line with two-dimensional organization. By fitting the fraction of simulations in the high-N state to an exponential function (Fig. 5E), we show that a characteristic time exists when the simulations decay to the low state. Thereby, we find that the circle model predicts decreasing R max , increasing L, or increasing R min speed up the characteristic time when the transition occurs (Fig. 5F). Decreasing R max is in correspondence with the results presented in Figs. 1-2. Increasing L has previously been reported to facilitate self-aggregation [7,13]. The literature also matches that coarser horizontal resolution favors self-aggregation given our result in Fig. 4E-H where we show that coarser horizontal resolution results in increased R min [10,40].

VI. DISCUSSION AND CONCLUSION
Convincing theories have been proposed for Turing-like coarsening of the RCE atmosphere into moist and dry sub-regions [13,41,42]. Such reaction-diffusion systems describe local positive feedback, e.g., of moisture on itself.
In contrast, we here explicitly model the non-local, two-particle interaction resulting from interacting cold pool (CP) gust fronts -competing with the well-studied local moisture feedback. The notion that CP collisions form new convective events is well documented [43][44][45] and addressed in toy models [22,46]. Our study investigates how the CP radius influences CSA and builds on the finding that the total rain cell number is approximately constant across simulations (Fig. S3). Indeed, relatively constant rain cell numbers (Fig. S3) and rain intensities (Fig. S2B) are supported by radiation constraints on precipitation in RCE [47]. Our circle model implicates that large CPs, as formed by pronounced rain evaporation, become space-filling where the whole domain is filled by CPs and there is a connected patch through the domain among CPs from the same generation. From hexagonal close-packing, a lower radius bound for space-filling would be R max > L(3 √ 3N ) −1/2 ≈ 5.7 km, where L = 92 km is the side length of the simulation domain and N ≈ 50 is the number of CPs per generation (Fig. S3). For radii smaller than this, areas emerge that cannot be reached by newly initialized circles -a gap results, and the transition to CSA sets in.
When (realistically) departing from perfect close-packing, the required value of R max may lie somewhat highercommensurate with our findings (Fig. 2) and the transition to CSA between Evap=.6, where R max ≈ 12 km, and Evap=.2, where R max ≈ 8 km.
Our model simplifies CP expansion by assuming a step function v(r) = v 0 θ(R max − r). In reality, CPs initially grow quickly and their speed of spreading decreases gradually over the course of a few hours (Fig. 2B) [48,49]. Introducing a smoothly varying gust front speed into our model would require a time-dependent expansion speed factor and a numerically-approximate approach might then be more practicable [24]. The presented model further does not reach a final, fully-aggregated state, where a small fraction of the domain intensely convects indefinitely. This sustained activity might be obtained by adding spatial noise (displacing new circles slightly away from the exact geometric collision point) and systematically increased triggering probabilities for decreased overall rain area [46]. Extensions could include explicit incorporation of the "super-CP" [32] and radiatively driven CP [10,50], constituting the two components of the final large-scale circulation. This model extension would allow triggering events at the edges of the intensely-convecting supercell due to convective CPs colliding with the opposing radiatively driven CP. Such circulation feedbacks may well be essential in stabilizing the final steady state but may not be required to develop the first dry patches and their initial growth, which we have focused on in this work.

COMPETING INTERESTS
The authors declare no competing interests. [

Supplementary Information
Here, we analytically find the average number of circle collisions for neighbors in a system with randomly positioned cells, and we provide supplementary figures, including the initial conditions for the large-eddy simulations (Fig. S1), the time development of the domain-mean specific humidity variation, rain intensity, and temperature (Fig. S2), as well as the number of rainfall tracks (Fig. S3).
Direct circle collisions for random seeding. We seek to compute the number of direct circle collisions for an initial random population of circle centers, with circles emerging synchronously and at equal speed from all circle centers. By direct we mean, that the collision takes place along the line connecting the two circle centers involved.
Finding the number of such collisions is equivalent to finding the number of line-of-sight connections (also known as a Gabriel connection) among these points (Fig. S4).
Consider N points randomly seeded in a total area A = πR 2 . A line-of-sight connection between any two points r i and r j at distance l ≡ |r i − r j | exists if there are no points located inside the circle of radius l/2 centered at (r i + r j )/2 ( Fig. S4). This is the type of connection that we require for any two colliding circles to initialize the growth of a new circle in Fig. 5. One must further consider the probability of finding two points at distance l. For this purpose, define the density of points as ρ ≡ N/A = N/πR 2 . Now consider the infinitesimal area a(l) ≡ 2π l dl between two circles of radii l and l + dl. The number of points contained in this area is n(l)dl = ρ a(l) = 2N l R 2 dl. (S1) For any two points at a given distance l, we now consider the probability p(l), that none of the remaining N − 2 points lie within the circle of radius l/2: where is the probability that any single point is not inside the area enclosed by a circle of radius l/2. Now the total number of expected line-of-sight (LOS) connections for a fixed given point to any of the other points can be computed:    S4. A line-of-sight connection. Schematic illustrating points in a 2D domain. Two points, ri and rj, separated by a distance l have a line-of-sight connection given that no points are located inside the circle of radius l/2 with the two points on its rim. In mathematics, this concept is also known as Gabriel neighbors.