New Critical Length for the Onset of Self‐Aggregation of Moist Convection

Convective self‐aggregation (CSA) in an idealized modeling framework is key to understanding the role of clouds. To investigate the existence of characteristic length of CSA onset, we conducted systematic cloud‐resolving simulations, with a scope covering the horizontal domain size and resolution. In the high‐resolution simulation, CSA can occur with a square domain larger than ~500 km. Based on the competition between two near‐surface horizontal divergent flows, we discuss the characteristic length existence. While the flow induced by radiative cooling in the subsidence region acts as positive feedback for moisture perturbation and scales with the domain size, the other flow induced by evaporative cooling of falling rain in the convective region acts as negative feedback and does not scale. The study suggests characteristic length existence for the organization of moist convection, even in real‐world conditions.


Introduction
Moist convection variability in the tropical atmosphere and its relationship with weather and climate are a subject of long-standing interest in atmospheric science research. Tropical moist convections often organize into a mesoscale convective system or a larger-scale cloud cluster. Since such organized cloud systems result in a large amount of precipitation and redistribute water and energy, the processes leading to the convective organization are important for understanding the weather-climate relationship. A paradigm of radiative-convective equilibrium (RCE), in which convective overturning develops in response to instability produced by radiative cooling, is a useful framework for investigating moist convective organization (Hartmann et al., 2019;Nakajima & Matsuno, 1988;Randall et al., 1994). In numerical simulations under the RCE framework, moist convection was organized into an aggregated cloud system even under constant sea surface conditions (Held et al., 1993;Tompkins & Craig, 1998a). This organization process is called convective self-aggregation (CSA) and has been intensively studied using a wide variety of models in recent years (Arnold & Randall, 2015;Bretherton et al., 2005, hereafter BBK05;Emanuel et al., 2014;Haerter, 2019;Hohenegger & Stevens, 2016;Khairoutdinov & Emanuel, 2013;Satoh & Matsuda, 2009;Tompkins, 2001b;Windmiller & Craig, 2019;Yanase & Takemi, 2018). It has recently been pointed out that CSA alters even mean fields (BBK05; Bony et al., 2015Bony et al., , 2016Cronin & Wing, 2017;Drotos et al., 2020;Holloway et al., 2017;Mauritsen & Stevens, 2015;Noda et al., 2019;Tobin et al., 2012;Wing, 2019). Thus, we now recognize that CSA is an essential player in controlling the Earth's climate.
Although many RCE studies have used cloud-resolving models, Muller and Held (2012, hereafter MH12) raised a serious concern related to the dependence of CSA on numerical configurations. MH12 investigated the sensitivity of CSA onset to horizontal domain size and resolution. They found that CSA is unlikely to occur from a homogeneous initial moisture field with a high resolution or small domain. Other studies also reported the existence of a critical domain size for the onset of CSA (e.g., BBK05; Jeevanjee & Romps, 2013, hereafter JR13), and the resolution dependence has also been reported again (e.g., Muller & Bony, 2015, hereafter MB15;Becker et al., 2018). JR13 discussed the domain size dependence in terms of the negative effect of a cold pool on CSA. Evaporative cooling below precipitating clouds induces cold pools, which spread horizontally and homogenize the boundary layer moisture. They also showed that CSA occurs at any domain size when the evaporative cooling effect is artificially removed, although it becomes weaker as the domain size decreases. Although other studies have confirmed the role of cold pools and evaporative cooling in the boundary layer (e.g., MB15;Holloway & Woolnough, 2016;Yang, 2018aYang, , 2019, it is still unclear how cold pools control CSA in full-physics simulations. In consideration of the radiation process, MB15 emphasized the importance of spatial variation of longwave radiative cooling on CSA. They implied that CSA does not occur in the simulation characterized by high-resolution and small domains due to insufficient cloud amounts and the low radiative cooling associated with it. By artificially imposing spatial variation of radiative cooling, MB15 showed that CSA occurs even at higher resolutions and smaller domains. Although their results suggested the importance of radiative cooling, the necessity of such an additional diabatic forcing for the onset of CSA in the high-resolution simulation is still debatable. The sensitivity of CSA to the model configuration itself also remains a problem in the RCE framework (Wing et al., , 2018. If we regard an atmospheric model as a mathematical model with uncertainties resulting from physics parameterization being neglected, it is reasonable to assume that the characteristics of the simulated atmosphere converge with the reduction in the grid spacing. Nonetheless, we should consider the domain size as a parameter that constrains the upper limit of the spatial scale for the variability of the simulated atmosphere. Hence, the question to be clarified is whether moist convection remains randomly distributed (scattered state) or is preferably organized into the larger scale (aggregated state) at a sufficiently high resolution as the domain size increases. By sweeping the domain size as a parameter while maintaining a high resolution, we can investigate (1) whether a critical length exists or not and (2) what the critical length is. The critical length is the minimum length for the onset of CSA. Several studies have investigated the way in which the statistical properties of the simulated atmosphere change and converge, as the resolution increases in a practical problem (e.g., Bryan et al., 2003;Kajikawa et al., 2016;Khairoutdinov et al., 2009;Miyamoto et al., 2013;Panosetti et al., 2018;Scheufele, 2014;Sueki et al., 2019;Yashiro et al., 2016). However, concerning the problem of CSA, the convergence characteristic of the critical domain size against the resolution has not yet been systematically examined.
According to MH12, CSA did not spontaneously occur within a domain range of less than 500 km at a resolution higher than 2 km (see Figure 6a of MH12). In this study, we extended the domain size to clarify whether CSA occurred in numerical simulations of moist convection with grid spacings less than 2 km. Considering the possibility of the existence of a critical domain size between scattered and aggregated states for a given grid spacing, we investigated the variation and convergence of the critical domain size with increasing resolution. For this purpose, we conducted a series of RCE simulations within a wide range of domain size and grid spacing parameters. Furthermore, we discussed a mechanism for the determination of the critical length for the onset of CSA, focusing on the negative feedback of the cold pool and the competing positive feedback.

Method
We conducted a series of numerical simulations using SCALE-RM Version 5.3.3, a regional atmospheric model constructed with Scalable Computing for Advanced Library and Environment Sato et al., 2015; SCALE) without the Coriolis force.
We imposed a doubly periodic lateral boundary condition on a square horizontal computational domain. In the vertical, we set the grid spacing at 75 m near the surface, 200 m at a height of about 4 km, 600 m at a height of about 12 km, and at 1 km at the top of the model at 24 km. The imposed conditions included a sponge layer condition above the 18 km height and a sea surface condition with a constant sea surface temperature of 300 K. subgrid-scale turbulent process, and the bulk method using a universal function (Beljaars & Holtslag, 1991;Wilson, 2001) as the surface flux scheme. For the radiative processes, we used the model simulation radiation transfer code (MSTRNX; Sekiguchi & Nakajima, 2008). The solar insolation parameters were the same as those used in previous studies (MH12 ;Tompkins & Craig, 1998a).
We performed a series of experiments on RCE sensitivity to horizontal grid spacing and the length of the square domain. Hereinafter we abbreviated the horizontal grid spacing and the length of the square domain as H (in meter) and L (in kilometer), respectively. We then referred to each experiment with its horizontal grid spacing and length; for example, an experiment with a grid spacing of 2,000 m and domain size of 192 km was referred to as H2000L192. We first ran a control experiment using H2000L96. Its vertical profile in the initial condition was based on the tropical squall line case in the Global Energy and Water Cycle Cloud System Study (GCSS) without horizontal wind. A potential temperature perturbation of 0.1 K was added for spinning up. We obtained a homogeneous distribution of precipitable water (PW) with mean values of approximately 35 kg m −2 in the statistical equilibrium state after 100 day integration in the control experiment. The vertical thermodynamic profiles averaged over the last day were used as the initial conditions for the other experiments. This procedure is similar to that used in previous studies (e.g., BBK05; Hohenegger & Stevens, 2016;Patrizio & Randall, 2019;.
The experiments are summarized in supporting information Table S1. The integration time for most experiments was 40 or 50 days; however, it was 13 days for H500L960 due to the computational resources constraint. Despite the relatively short time taken to fully reach an equilibrium, the CSA onset can be sufficiently distinguished, as shown in the next section.

Results
To verify the reproducibility of the results of MH12 with our model, we first checked four typical experiments: with two different horizontal resolutions (H1000 and H2000) and two different domain sizes (L96 and L384). Since CSA is characterized by the separation of moist and dry regions, the temporal evolution of the horizontal distribution of PW can be considered a primary indicator of CSA. As shown in Figure 1d, relatively dry contiguous areas, henceforth referred to as dry patches, emerged and expanded with time in H2000L384, while the PW field remained nearly homogeneous in H1000L384, H1000L96, and H2000L96 (Figures 1c, 1e, and 1f); convection was aggregated in H2000L384 and scattered in the others. This sensitivity of the CSA onset to the domain size and grid spacing is consistent with the results of MH12. Thus, the question is whether it remains scattered or transitions to the aggregated state as the domain size increases at a grid spacing of 1,000 m. Figure 1a indicates the transition to the aggregated state in H1000L960. CSA also occurs in H2000L960 (Figure 1b), as expected from the results of MH12.
Previous studies reported that the CSA onset is associated with the expansion of dry patches (e.g., BBK05; Hohenegger & Stevens, 2016;. To detect CSA in its early stage, we use a dry patch fraction where PW is lower than a threshold, as in Hohenegger and Stevens (2016). We employed a threshold of 30 kg m −2 . Our conclusion is unchanged with the threshold between 23 and 32 kg m −2 . Figure 2a shows the time evolution of the fractions of dry patches in all experiments. We define a case as aggregated if the fraction of the dry patch increases with time and as scattered if the fraction of the dry patch remains at approximately zero. The maximum fractions in all the scattered cases are less than 0.05, while all the aggregated cases have a fraction of greater than 0.25 at the end of the integration. The robustness of the regime classification is ensured by checking the horizontal distribution of PW for all simulations and the time evolution of other metrics, such as the interquartile range and the standard deviation of PW ( Figures S1-S3). Based on the above criterion, Figure 2b shows the RCE regime in the H-L parameter space. For any grid spacing, the RCE regime transitions from scattered to aggregated by increasing the domain size. In Figure 2b, the regime boundary is drawn and the line may be divided into three parts: I, II, and III. Lines II and III are concordant with the regime boundary in Figure 6a of MH12; that is, the tilt of the regime boundary drastically changes at a grid spacing of approximately 2 km. However, the linear extrapolation of Line II to smaller grid spacing cannot explain the aggregated regime simulated with the high-resolution cases in this study (e.g., H1000L560 and H500L560). In H1000 and H500, there should be a regime boundary somewhere between L384 and L560, and a new regime boundary line, namely, Line I, appears at a high resolution. As Line I indicates, the critical domain size converges around a value of 500 km as the resolution increases.

Geophysical Research Letters
Thus, Line II is not the lower limit of the grid spacing for the spontaneous occurrence of CSA (i.e., the perspective of MH12), but a sharp increase in the critical domain size connecting Lines I and III. Our results update the regime diagram of MH12 consistently.
We further explored the mechanism of the CSA onset. As is well known, the CSA onset is usually associated with the expansion of the dry patch. Although some positive feedback mechanisms that enhance the dry patch must surpass the negative feedback mechanisms, such as cold pools, for the CSA onset, it is still debatable why and under what circumstances a dry patch can expand. Because the effect of the cold pool is expected to be concentrated at the near surface, we note the change rate of the near-surface water vapor by transport. Figure 3 shows the horizontal flux convergence of water vapor and the effective stream function of BBK05, for the experiments the same as in Figure 1. All panels show a negative (drying) tendency from the surface to a height of several hundred meters at high PW ranks (blue ellipses in the leftmost panels). Phenomenologically, it corresponds to the occurrence of precipitation clouds in the relatively moist region; below the cloud, evaporatively driven cold pools spread the water vapor horizontally. In the scattered cases (e.g., H1000L96 and H1000L384), a bipolar structure is exhibited near the surface, that is, a pair of a negative tendency in the moister region and a positive tendency in the drier region. This means the transport of water vapor from the moist to the dry region (e.g., blue arrows in the leftmost panels), which acts as negative feedback, suppressing horizontal heterogeneity. This result supports the theory of the homogenizing effect of cold pools proposed by JR13. This can also be seen as the clockwise circulation near the surface. Contrarily, a tripolar structure appears in the aggregated cases, that is, in addition to the bipolar structure, another negative tendency emerges at low PW ranks (e.g., red ellipse in H1000L960). The transport of water vapor from the dry to the moist region (e.g., the red arrow in H1000L960) acts as positive feedback. Phenomenologically, it reflects horizontal divergent flow in the boundary layer of the radiatively driven compensating subsidence area, which is usually relatively dry. This corresponds to the counterclockwise circulation. The above results suggest that the horizontal flux of water vapor near the surface plays a crucial role in determining the regime. Previous studies have reported that, in the aggregated state, the large-scale circulation strengthens as the domain size increases (e.g., JR13; Arnold & Randall, 2015;Patrizio & Randall, 2019). Nevertheless, this is the first study to directly show that there is a qualitative difference, compared to the scattered state, in the competition of the opposite feedbacks in the early stages of CSA.

Discussion
Based on the above result, we explain why a positive feedback enhancing on the dry patch emerges with the increase in domain size and what factors set the critical length for the CSA onset. Wing and Cronin (2016), in their remoistening length theory, showed that the speed of wind during horizontal divergent flow in the boundary layer of the horizontally continuous subsidence area depends on the distance from the center point of the area. Let us consider the case in which the subsidence area takes a circular form. The downward vertical velocity in the subsidence region at the top of the boundary layer height z pbl is denoted as w sub . The vertical velocity at the ground surface is zero. The horizontal wind speed in the boundary layer u s linearly increases with the radius r. In this case, u s r ð Þ ¼ w sub 2z pbl r assuming u s (0) = 0. Intuitively, in a two-dimensional doubly periodic space that contains multiple points (convective elements), a continuous subsidence area can be imagined as a circle that does not contain any point. The wind speed u s can be regarded as one measure of positive feedback, whose negative feedback counterpart is considered the wind speed u c of a cold pool. If u s is negligibly small compared to u c , the cold pool can spread and homogenize the moisture field. Conversely, if u s becomes larger than or comparable to u c , the spread of the cold pool stops or vanishes. To roughly estimate their competition, if we use typical values of 5.0 × 10 −3 m s −1 for w sub and 5.0 × 10 2 m for z pbl (Tompkins & Craig, 1998b), then u s becomes 2.5 m s −1 at 5.0 × 10 5 m. Considering a typical value of 2-5 m s −1 for u c under RCE (Tompkins, 2001a), u s and u c become the same order when r is several hundreds of kilometers. This qualitatively explains that a characteristic length scale exists in hundreds of kilometers, where the two opposite feedbacks reverse. The maximum attainable size of the continuous subsidence area decreases as the domain size decreases. This is the reason CSA does not occur with a smaller domain size. Moreover, this idea is consistent with MB15, who argue that low-level diabatic processes and the associated moisture transports are essential for the CSA onset. Our study supports their view and adds a new insight that the competition between the opposite feedbacks sets the critical length. This interpretation can also explain previous sensitivity experiments, such as switching off the evaporation of rainfall and strengthening radiative cooling in the dry region by considering them as cases in which u c = 0 and strong u s , respectively.
Next, we explain why the critical length sharply changes at approximately 2 km resolution. This issue can be discussed by comparing H2000L384 (aggregated case) with H1000L384 (scattered case). We found that the low-cloud amount in the subsidence region during Days 11-20 was larger in H2000L384 than in H1000L384 ( Figures S4d and S4c), which is consistent with the results of MB15 and Yanase and Takemi (2018). This suggests that the low-resolution experiments represent a larger amount of low cloud, the enhanced radiative cooling strengthens the positive feedback and u s , and the critical length becomes smaller than that of high-resolution experiments.
Based on the above discussion, Figure 4 provides a schematic of the competition mechanism of two feedbacks from a phenomenological viewpoint. We initially suppose that a convective cloud occurs at some distance from the center of the subsidence area. Subsequently, a cold pool is generated in the boundary layer below the preceding convective cloud. In the scattered case (Figure 4a), the cold pool spreads over the preceding subsidence area and homogenizes the boundary layer moisture. This process corresponds to the blue arrow shown in Figure 3. In this situation, the next convective cloud can occur in the preceding subsidence area. Note that this cycle occurs everywhere in the whole domain, and a given region becomes a subsidence area or convective area by time. Thus, a state of statistical equilibrium is maintained where the temporal average equals the horizontal average. Conversely, in the aggregated case (Figure 4b), the spread of the cold pool is impeded at a certain large subsidence area. In this situation, the next convective cloud avoids the Figure 3. Horizontal flux convergence of water vapor in PW rank-height coordinate (shaded). The PW rank is the ascending PW average for each horizontal block with a size of (48 km) 2 . The experiments are the same as in Figure 1, and the averaged period is the previous 10 days of each corresponding panel in Figure 1. The effective stream function is shown by black dashed (solid) contours for counterclockwise (clockwise). The interval is set to 0.01 kg m −2 s −1 for L384 and is changed to be proportional to the number of PW rank for others. See the main text for the description of colored arrows and ellipses.

Geophysical Research Letters
preceding large subsidence area. Thus, the subsidence area can be maintained where the dry patch can grow without being moistened by cold pools. This process corresponds to the red arrow shown in Figure 3a. Note that for the onset of CSA, the existence of at least one such place is sufficient. Eventually, these opposite flows become balanced in the aggregated case, as described by Windmiller and Hohenegger (2019).

Conclusion and Future Perspectives
We conducted a series of cloud-resolving simulations under RCE to investigate the CSA onset dependence on domain size and resolution. We found that CSA occurs in high-resolution large-domain simulations. This result leads to an updated regime diagram. While the critical domain size, above which CSA occurs, depends on the resolution and sharply changes across the horizontal grid spacing of approximately 2 km, it appears to converge toward approximately 500 km as the grid spacing decreases. The regime change mechanism with the critical domain size can be interpreted based on the competition of opposite effects on the moisture variance in the boundary layer. As the domain size increases, the positive feedback from the radiatively driven near-surface horizontal divergent flow strengthens. When the domain size exceeds a critical length, the positive feedback overwhelms the negative feedback due to the evaporatively driven cold pools in the convective region.
Finally, the critical length of approximately 500 km, corresponding to the order of the meso-α scale, suggests the existence of a characteristic length in the interaction of cumulus clouds with the large-scale environment. To quantitatively understand the critical length from the perspective of moisture, it is necessary to include other effects such as surface flux and microphysical processes. Further, it should be clarified how the horizontal winds contributing to moisture transport are maintained from a dynamical perspective. The dynamical balance in the aggregated state has recently been revealed by considering the pressure gradient and frictional forces (e.g., Arnold & Putman, 2018;Patrizio & Randall, 2019;Yang, 2018b). However, it is less understood with respect to the CSA onset owing to the difficulty in describing the time-dependent transient process. Naumann et al. (2017) investigated the relationship between radiative cooling and low-level circulation. Their argument may be a clue to represent the radiative process in the competing feedbacks related to the CSA onset. Further refinement of the concept of critical length determination presented in this study will lead to a deeper understanding of the physical processes related to convective organization, even in real-world conditions. In addition, there may be other larger characteristic lengths of CSA, which should be explored in the several thousands of kilometers scale domain, as investigated by Patrizio and Randall (2019). Using such a larger domain will allow a clearer discussion on the physics of CSA.

Data Availability Statement
The following computational resources were used: the K-computer, the RIKEN R-CCS through the HPCI System Research project (Project ID: hp170323); supercomputer of ACCMS, Kyoto University; and Oakbridge-CX, The University of Tokyo. The authors are grateful to Team SCALE for providing SCALE Version 5.3.3, which is freely available at http://r-ccs-climate.riken.jp/scale/download/archives/under the 2-Clause BSD license. Similarly, the model configuration files, scripts, and data for developing the figures used for this study are publicly available in the Open Science Framework (https://osf.io/wdf3y/).